Reference Manual

Go to the documentation of this file.
00001 #include "RecoTracker/DebugTools/interface/TestSmoothHits.h"
00002 
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00005 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00006 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
00007 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00008 #include "RecoTracker/TrackProducer/interface/TrackingRecHitLessFromGlobalPosition.h"
00009 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
00010 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
00011 #include <TDirectory.h>
00012 #include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
00013 
00014 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
00015 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00016 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00017 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00018 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00019 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00020 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00021 
00022 typedef TrajectoryStateOnSurface TSOS;
00023 typedef TransientTrackingRecHit::ConstRecHitPointer CTTRHp;
00024 using namespace std;
00025 using namespace edm;
00026 
00027 TestSmoothHits::TestSmoothHits(const edm::ParameterSet& iConfig):
00028   conf_(iConfig){
00029   LogTrace("TestSmoothHits") << conf_<< std::endl;
00030   propagatorName = conf_.getParameter<std::string>("Propagator");   
00031   builderName = conf_.getParameter<std::string>("TTRHBuilder");   
00032   srcName = conf_.getParameter<std::string>("src");   
00033   fname = conf_.getParameter<std::string>("Fitter");
00034   sname = conf_.getParameter<std::string>("Smoother");
00035   mineta = conf_.getParameter<double>("mineta");
00036   maxeta = conf_.getParameter<double>("maxeta");
00037 }
00038 
00039 TestSmoothHits::~TestSmoothHits(){}
00040 
00041 void TestSmoothHits::beginRun(edm::Run & run, const edm::EventSetup& iSetup)
00042 {
00043  
00044   iSetup.get<TrackerDigiGeometryRecord>().get(theG);
00045   iSetup.get<IdealMagneticFieldRecord>().get(theMF);  
00046   iSetup.get<TrackingComponentsRecord>().get(propagatorName,thePropagator);
00047   iSetup.get<TransientRecHitRecord>().get(builderName,theBuilder);
00048   iSetup.get<TrajectoryFitter::Record>().get(fname, fit);
00049   iSetup.get<TrajectoryFitter::Record>().get(sname, smooth);
00050 
00051   file = new TFile("testSmoothHits.root","recreate");
00052   for (int i=0; i!=6; i++)
00053     for (int j=0; j!=9; j++){
00054       if (i==0 && j>2) break;
00055       if (i==1 && j>1) break;
00056       if (i==2 && j>3) break;
00057       if (i==3 && j>2) break;
00058       if (i==4 && j>5) break;
00059       if (i==5 && j>8) break;
00060       title.str("");
00061       title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
00062       hPullGP_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00063       title.str("");
00064       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
00065       hPullGP_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00066       title.str("");
00067       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
00068       hPullGP_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00069       title.str("");
00070       title << "Chi2Increment_" << i+1 << "-" << j+1;
00071       hChi2Increment[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00072 
00073       title.str("");
00074       title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
00075       hPullGM_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00076       title.str("");
00077       title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
00078       hPullGM_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00079       title.str("");
00080       title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
00081       hPullGM_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00082 
00083       title.str("");
00084       title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
00085       hPullGP_X_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00086       title.str("");
00087       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
00088       hPullGP_Y_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00089       title.str("");
00090       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
00091       hPullGP_Z_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00092 
00093       title.str("");
00094       title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
00095       hPullGP_X_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00096       title.str("");
00097       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
00098       hPullGP_Y_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00099       title.str("");
00100       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
00101       hPullGP_Z_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00102 
00103       if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
00104         //mono
00105         title.str("");
00106         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
00107         hPullGP_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00108         title.str("");
00109         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00110         hPullGP_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00111         title.str("");
00112         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00113         hPullGP_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00114 
00115         title.str("");
00116         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
00117         hPullGM_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00118         title.str("");
00119         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00120         hPullGM_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00121         title.str("");
00122         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00123         hPullGM_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00124 
00125         title.str("");
00126         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
00127         hPullGP_X_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00128         title.str("");
00129         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
00130         hPullGP_Y_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00131         title.str("");
00132         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
00133         hPullGP_Z_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00134 
00135         title.str("");
00136         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
00137         hPullGP_X_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00138         title.str("");
00139         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
00140         hPullGP_Y_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00141         title.str("");
00142         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
00143         hPullGP_Z_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00144 
00145         //stereo
00146         title.str("");
00147         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00148         hPullGP_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00149         title.str("");
00150         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00151         hPullGP_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00152         title.str("");
00153         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00154         hPullGP_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00155 
00156         title.str("");
00157         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00158         hPullGM_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00159         title.str("");
00160         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00161         hPullGM_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00162         title.str("");
00163         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00164         hPullGM_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00165 
00166         title.str("");
00167         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
00168         hPullGP_X_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00169         title.str("");
00170         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
00171         hPullGP_Y_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00172         title.str("");
00173         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
00174         hPullGP_Z_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00175 
00176         title.str("");
00177         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
00178         hPullGP_X_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00179         title.str("");
00180         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
00181         hPullGP_Y_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00182         title.str("");
00183         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
00184         hPullGP_Z_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00185       }
00186     }
00187   hTotChi2Increment = new TH1F("TotChi2Increment","TotChi2Increment",1000,0,100);
00188   hChi2_vs_Process  = new TH2F("Chi2_vs_Process","Chi2_vs_Process",1000,0,100,17,-0.5,16.5);  
00189   hChi2_vs_clsize  = new TH2F("Chi2_vs_clsize","Chi2_vs_clsize",1000,0,100,17,-0.5,16.5);
00190 }
00191 
00192 
00193 void TestSmoothHits::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00194 {
00195   LogTrace("TestSmoothHits") << "new event" << std::endl;
00196 
00197   iEvent.getByLabel(srcName,theTCCollection ); 
00198   hitAssociator = new TrackerHitAssociator(iEvent);
00199 
00200   TrajectoryStateCombiner combiner;
00201 
00202   for (TrackCandidateCollection::const_iterator i=theTCCollection->begin(); i!=theTCCollection->end();i++){
00203 
00204     LogTrace("TestSmoothHits") << "new candidate" << std::endl;
00205       
00206     const TrackCandidate * theTC = &(*i);
00207     PTrajectoryStateOnDet state = theTC->trajectoryStateOnDet();
00208     const TrackCandidate::range& recHitVec=theTC->recHits();
00209 
00210     //convert PTrajectoryStateOnDet to TrajectoryStateOnSurface
00211     
00212     
00213     DetId  detId(state.detId());
00214     TrajectoryStateOnSurface theTSOS=
00215       trajectoryStateTransform::transientState(state, &(theG->idToDet(detId)->surface()),theMF.product());
00216 
00217     if (theTSOS.globalMomentum().eta()>maxeta || theTSOS.globalMomentum().eta()<mineta) continue;
00218     
00219     //convert the TrackingRecHit vector to a TransientTrackingRecHit vector
00220     TransientTrackingRecHit::RecHitContainer hits;
00221     
00222     for (edm::OwnVector<TrackingRecHit>::const_iterator i=recHitVec.first;
00223          i!=recHitVec.second; i++){
00224       hits.push_back(theBuilder->build(&(*i) ));
00225     }
00226 
00227     //call the fitter
00228     std::vector<Trajectory> fitted = fit->fit(theTC->seed(), hits, theTSOS);
00229     //call the smoother
00230     std::vector<Trajectory> result; 
00231     for(std::vector<Trajectory>::iterator it = fitted.begin(); it != fitted.end(); it++) {
00232       std::vector<Trajectory> smoothed = smooth->trajectories(*it);
00233       result.insert(result.end(), smoothed.begin(), smoothed.end());
00234     }
00235     if (result.size()==0) continue;
00236     std::vector<TrajectoryMeasurement> vtm = result[0].measurements();
00237 
00238     TSOS lastState = theTSOS;
00239     for (std::vector<TrajectoryMeasurement>::iterator tm=vtm.begin(); tm!=vtm.end();tm++){
00240 
00241       TransientTrackingRecHit::ConstRecHitPointer rhit = tm->recHit();
00242       if ((rhit)->isValid()==0&&rhit->det()!=0) continue;
00243       LogTrace("TestSmoothHits") << "new hit" ;
00244 
00245       int subdetId = rhit->det()->geographicalId().subdetId();
00246       int layerId  = 0;
00247       DetId id = rhit->det()->geographicalId();
00248       if (id.subdetId()==3) layerId = ((TIBDetId)(id)).layer();
00249       if (id.subdetId()==5) layerId = ((TOBDetId)(id)).layer();
00250       if (id.subdetId()==1) layerId = ((PXBDetId)(id)).layer();
00251       if (id.subdetId()==4) layerId = ((TIDDetId)(id)).wheel();
00252       if (id.subdetId()==6) layerId = ((TECDetId)(id)).wheel();
00253       if (id.subdetId()==2) layerId = ((PXFDetId)(id)).disk();
00254       LogTrace("TestSmoothHits") << "subdetId=" << subdetId << " layerId=" << layerId ;
00255 
00256       double delta = 99999;
00257       LocalPoint rhitLPv = rhit->localPosition();
00258 
00259       std::vector<PSimHit> assSimHits = hitAssociator->associateHit(*(rhit->hit()));
00260       if (assSimHits.size()==0) continue;
00261       PSimHit shit;
00262       for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
00263         if ((m->localPosition()-rhitLPv).mag()<delta) {
00264           shit=*m;
00265           delta = (m->localPosition()-rhitLPv).mag();
00266         }
00267       }
00268 
00269       TSOS currentState = combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
00270       TSOS updatedState = tm->updatedState();
00271  
00272       //plot chi2 increment
00273       double chi2increment = tm->estimate();
00274       LogTrace("TestSmoothHits") << "tm->estimate()=" << tm->estimate();
00275       title.str("");
00276       title << "Chi2Increment_" << subdetId << "-" << layerId;
00277       hChi2Increment[title.str()]->Fill( chi2increment );
00278       hTotChi2Increment->Fill( chi2increment );
00279       hChi2_vs_Process->Fill( chi2increment, shit.processType() );
00280       if (dynamic_cast<const SiPixelRecHit*>(rhit->hit()))      
00281         hChi2_vs_clsize->Fill( chi2increment, ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size() );
00282       if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit()))    
00283         hChi2_vs_clsize->Fill( chi2increment, ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() );
00284      
00285       //test hits
00286       const Surface * surf = &( (rhit)->det()->surface() );
00287       LocalVector shitLMom;
00288       LocalPoint shitLPos;
00289        if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
00290         double rechitmatchedx = rhit->localPosition().x();
00291         double rechitmatchedy = rhit->localPosition().y();
00292         double mindist = 999999;
00293         float distx, disty;
00294         std::pair<LocalPoint,LocalVector> closestPair;
00295         const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) ((const GluedGeomDet *)(rhit)->det())->stereoDet();
00296         const BoundPlane& plane = (rhit)->det()->surface();
00297         for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++) {
00298           //project simhit;
00299           std::pair<LocalPoint,LocalVector> hitPair = projectHit((*m),stripDet,plane);
00300           distx = fabs(rechitmatchedx - hitPair.first.x());
00301           disty = fabs(rechitmatchedy - hitPair.first.y());
00302           double dist = distx*distx+disty*disty;
00303           if(sqrt(dist)<mindist){
00304             mindist = dist;
00305             closestPair = hitPair;
00306           }
00307         }
00308         shitLPos = closestPair.first;
00309         shitLMom = closestPair.second;
00310       } else {
00311         shitLPos = shit.localPosition();        
00312         shitLMom = shit.momentumAtEntry();
00313       }
00314       GlobalVector shitGMom = surf->toGlobal(shitLMom);
00315       GlobalPoint shitGPos = surf->toGlobal(shitLPos);
00316 //      if (dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())) {
00317 //      double rechitmatchedx = rhit->localPosition().x();
00318 //      double rechitmatchedy = rhit->localPosition().y();
00319 //      double mindist = 999999;
00320 //      double distx, disty;
00321 //      std::pair<LocalPoint,LocalVector> closestPair;
00322 //      const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) ((const GluedGeomDet *)(rhit)->det())->stereoDet();
00323 //      const BoundPlane& plane = (rhit)->det()->surface();
00324 //      for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
00325 //        const PSimHit& hit = *m;
00326 //        const StripTopology& topol = stripDet->specificTopology();
00327 //        GlobalPoint globalpos = stripDet->surface().toGlobal(hit.localPosition());
00328 //        LocalPoint localHit = plane.toLocal(globalpos);
00329 //        //track direction
00330 //        LocalVector locdir=hit.localDirection();
00331 //        //rotate track in new frame     
00332 //        GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
00333 //        LocalVector dir=plane.toLocal(globaldir);
00334 //        float scale = -localHit.z() / dir.z();          
00335 //        LocalPoint projectedPos = localHit + scale*dir;                 
00336 //        float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));         
00337 //        LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame    
00338 //        LocalVector localStripDir = LocalVector( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
00339 //        std::pair<LocalPoint,LocalVector> hitPair( projectedPos, localStripDir);
00340 //        distx = fabs(rechitmatchedx - hitPair.first.x());
00341 //        disty = fabs(rechitmatchedy - hitPair.first.y());
00342 //        double dist = distx*distx+disty*disty;
00343 //        if(sqrt(dist)<mindist){
00344 //          mindist = dist;
00345 //          closestPair = hitPair;
00346 //        }
00347 //      }
00348 //      shitLPos = closestPair.first;
00349 //      shitLMom = closestPair.second;
00350 //       } else {
00351 //      shitLPos = shit.localPosition();        
00352 //      shitLMom = shit.momentumAtEntry();
00353 //       }
00354 //       GlobalVector shitGMom = surf->toGlobal(shitLMom);
00355 //       GlobalPoint shitGPos = surf->toGlobal(shitLPos);
00356 
00357       GlobalVector tsosGMom = currentState.globalMomentum();
00358       GlobalError  tsosGMEr(currentState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00359       GlobalPoint  tsosGPos = currentState.globalPosition();
00360       GlobalError  tsosGPEr = currentState.cartesianError().position();
00361 
00362       GlobalPoint rhitGPos = (rhit)->globalPosition();
00363       GlobalError rhitGPEr = (rhit)->globalPositionError();
00364 
00365       double pullGPX_rs = (rhitGPos.x()-shitGPos.x())/sqrt(rhitGPEr.cxx());
00366       double pullGPY_rs = (rhitGPos.y()-shitGPos.y())/sqrt(rhitGPEr.cyy());
00367       double pullGPZ_rs = (rhitGPos.z()-shitGPos.z())/sqrt(rhitGPEr.czz());
00368       //double pullGPX_rs = (rhitGPos.x()-shitGPos.x());
00369       //double pullGPY_rs = (rhitGPos.y()-shitGPos.y());
00370       //double pullGPZ_rs = (rhitGPos.z()-shitGPos.z());
00371       
00372       LogTrace("TestSmoothHits") << "rs" << std::endl;
00373 
00374       title.str("");
00375       title << "PullGP_X_" << subdetId << "-" << layerId << "_rs";
00376       hPullGP_X_rs[title.str()]->Fill( pullGPX_rs );
00377       title.str("");
00378       title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs";
00379       hPullGP_Y_rs[title.str()]->Fill( pullGPY_rs );
00380       title.str("");
00381       title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs";
00382       hPullGP_Z_rs[title.str()]->Fill( pullGPZ_rs );
00383 
00384       double pullGPX_tr = (tsosGPos.x()-rhitGPos.x())/sqrt(tsosGPEr.cxx()+rhitGPEr.cxx());
00385       double pullGPY_tr = (tsosGPos.y()-rhitGPos.y())/sqrt(tsosGPEr.cyy()+rhitGPEr.cyy());
00386       double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z())/sqrt(tsosGPEr.czz()+rhitGPEr.czz());
00387       //double pullGPX_tr = (tsosGPos.x()-rhitGPos.x());
00388       //double pullGPY_tr = (tsosGPos.y()-rhitGPos.y());
00389       //double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z());
00390 
00391       LogTrace("TestSmoothHits") << "tr" << std::endl;
00392 
00393       title.str("");
00394       title << "PullGP_X_" << subdetId << "-" << layerId << "_tr";
00395       hPullGP_X_tr[title.str()]->Fill( pullGPX_tr );
00396       title.str("");
00397       title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr";
00398       hPullGP_Y_tr[title.str()]->Fill( pullGPY_tr );
00399       title.str("");
00400       title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr";
00401       hPullGP_Z_tr[title.str()]->Fill( pullGPZ_tr );
00402 
00403       double pullGPX_ts = (tsosGPos.x()-shitGPos.x())/sqrt(tsosGPEr.cxx());
00404       double pullGPY_ts = (tsosGPos.y()-shitGPos.y())/sqrt(tsosGPEr.cyy());
00405       double pullGPZ_ts = (tsosGPos.z()-shitGPos.z())/sqrt(tsosGPEr.czz());
00406       //double pullGPX_ts = (tsosGPos.x()-shitGPos.x());
00407       //double pullGPY_ts = (tsosGPos.y()-shitGPos.y());
00408       //double pullGPZ_ts = (tsosGPos.z()-shitGPos.z());
00409 
00410       LogTrace("TestSmoothHits") << "ts1" << std::endl;
00411 
00412       title.str("");
00413       title << "PullGP_X_" << subdetId << "-" << layerId << "_ts";
00414       hPullGP_X_ts[title.str()]->Fill( pullGPX_ts );
00415       title.str("");
00416       title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts";
00417       hPullGP_Y_ts[title.str()]->Fill( pullGPY_ts );
00418       title.str("");
00419       title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts";
00420       hPullGP_Z_ts[title.str()]->Fill( pullGPZ_ts );
00421 
00422       double pullGMX_ts = (tsosGMom.x()-shitGMom.x())/sqrt(tsosGMEr.cxx());
00423       double pullGMY_ts = (tsosGMom.y()-shitGMom.y())/sqrt(tsosGMEr.cyy());
00424       double pullGMZ_ts = (tsosGMom.z()-shitGMom.z())/sqrt(tsosGMEr.czz());
00425       //double pullGMX_ts = (tsosGMom.x()-shitGMom.x());
00426       //double pullGMY_ts = (tsosGMom.y()-shitGMom.y());
00427       //double pullGMZ_ts = (tsosGMom.z()-shitGMom.z());
00428 
00429       LogTrace("TestSmoothHits") << "ts2" << std::endl;
00430 
00431       title.str("");
00432       title << "PullGM_X_" << subdetId << "-" << layerId << "_ts";
00433       hPullGM_X_ts[title.str()]->Fill( pullGMX_ts );
00434       title.str("");
00435       title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts";
00436       hPullGM_Y_ts[title.str()]->Fill( pullGMY_ts );
00437       title.str("");
00438       title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts";
00439       hPullGM_Z_ts[title.str()]->Fill( pullGMZ_ts );
00440 
00441       if (dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())) {
00442         //mono
00443         LogTrace("TestSmoothHits") << "MONO HIT" << std::endl;
00444         auto m = dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->monoHit();
00445         CTTRHp tMonoHit = 
00446           theBuilder->build(&m);
00447         if (tMonoHit==0) continue;
00448         vector<PSimHit> assMonoSimHits = hitAssociator->associateHit(*tMonoHit->hit());
00449         if (assMonoSimHits.size()==0) continue;
00450         const PSimHit sMonoHit = *(assSimHits.begin());
00451         const Surface * monoSurf = &( tMonoHit->det()->surface() );
00452         if (monoSurf==0) continue;
00453         TSOS monoState = thePropagator->propagate(lastState,*monoSurf);
00454         if (monoState.isValid()==0) continue;
00455 
00456         LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
00457         GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
00458         LocalPoint monoShitLPos = sMonoHit.localPosition();
00459         GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
00460 
00461         GlobalVector monoTsosGMom = monoState.globalMomentum();
00462         GlobalError  monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00463         GlobalPoint  monoTsosGPos = monoState.globalPosition();
00464         GlobalError  monoTsosGPEr = monoState.cartesianError().position();
00465 
00466         GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
00467         GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
00468 
00469         double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x())/sqrt(monoRhitGPEr.cxx());
00470         double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y())/sqrt(monoRhitGPEr.cyy());
00471         double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z())/sqrt(monoRhitGPEr.czz());
00472         //double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x());
00473         //double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y());
00474         //double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z());
00475 
00476         title.str("");
00477         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
00478         hPullGP_X_rs_mono[title.str()]->Fill( pullGPX_rs_mono );
00479         title.str("");
00480         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
00481         hPullGP_Y_rs_mono[title.str()]->Fill( pullGPY_rs_mono );
00482         title.str("");
00483         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
00484         hPullGP_Z_rs_mono[title.str()]->Fill( pullGPZ_rs_mono );
00485 
00486         double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x())/sqrt(monoTsosGPEr.cxx()+monoRhitGPEr.cxx());
00487         double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y())/sqrt(monoTsosGPEr.cyy()+monoRhitGPEr.cyy());
00488         double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z())/sqrt(monoTsosGPEr.czz()+monoRhitGPEr.czz());
00489         //double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x());
00490         //double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y());
00491         //double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z());
00492 
00493         title.str("");
00494         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
00495         hPullGP_X_tr_mono[title.str()]->Fill( pullGPX_tr_mono );
00496         title.str("");
00497         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
00498         hPullGP_Y_tr_mono[title.str()]->Fill( pullGPY_tr_mono );
00499         title.str("");
00500         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
00501         hPullGP_Z_tr_mono[title.str()]->Fill( pullGPZ_tr_mono );
00502 
00503         double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x())/sqrt(monoTsosGPEr.cxx());
00504         double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y())/sqrt(monoTsosGPEr.cyy());
00505         double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z())/sqrt(monoTsosGPEr.czz());
00506         //double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x());
00507         //double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y());
00508         //double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z());
00509 
00510         title.str("");
00511         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
00512         hPullGP_X_ts_mono[title.str()]->Fill( pullGPX_ts_mono );
00513         title.str("");
00514         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
00515         hPullGP_Y_ts_mono[title.str()]->Fill( pullGPY_ts_mono );
00516         title.str("");
00517         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
00518         hPullGP_Z_ts_mono[title.str()]->Fill( pullGPZ_ts_mono );
00519 
00520         double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x())/sqrt(monoTsosGMEr.cxx());
00521         double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y())/sqrt(monoTsosGMEr.cyy());
00522         double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z())/sqrt(monoTsosGMEr.czz());
00523         //double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x());
00524         //double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y());
00525         //double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z());
00526 
00527         title.str("");
00528         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
00529         hPullGM_X_ts_mono[title.str()]->Fill( pullGMX_ts_mono );
00530         title.str("");
00531         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
00532         hPullGM_Y_ts_mono[title.str()]->Fill( pullGMY_ts_mono );
00533         title.str("");
00534         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
00535         hPullGM_Z_ts_mono[title.str()]->Fill( pullGMZ_ts_mono );
00536 
00537         //stereo
00538         LogTrace("TestSmoothHits") << "STEREO HIT" << std::endl;
00539         auto s = dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->stereoHit();
00540         CTTRHp tStereoHit = 
00541           theBuilder->build(&s);
00542         if (tStereoHit==0) continue;
00543         vector<PSimHit> assStereoSimHits = hitAssociator->associateHit(*tStereoHit->hit());
00544         if (assStereoSimHits.size()==0) continue;
00545         const PSimHit sStereoHit = *(assSimHits.begin());
00546         const Surface * stereoSurf = &( tStereoHit->det()->surface() );
00547         if (stereoSurf==0) continue;
00548         TSOS stereoState = thePropagator->propagate(lastState,*stereoSurf);
00549         if (stereoState.isValid()==0) continue;
00550 
00551         LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
00552         GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
00553         LocalPoint stereoShitLPos = sStereoHit.localPosition();
00554         GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
00555 
00556         GlobalVector stereoTsosGMom = stereoState.globalMomentum();
00557         GlobalError  stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00558         GlobalPoint  stereoTsosGPos = stereoState.globalPosition();
00559         GlobalError  stereoTsosGPEr = stereoState.cartesianError().position();
00560 
00561         GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
00562         GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
00563 
00564         double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/sqrt(stereoRhitGPEr.cxx());
00565         double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/sqrt(stereoRhitGPEr.cyy());
00566         double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/sqrt(stereoRhitGPEr.czz());
00567         //double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x());
00568         //double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y());
00569         //double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z());
00570 
00571         title.str("");
00572         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
00573         hPullGP_X_rs_stereo[title.str()]->Fill( pullGPX_rs_stereo );
00574         title.str("");
00575         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
00576         hPullGP_Y_rs_stereo[title.str()]->Fill( pullGPY_rs_stereo );
00577         title.str("");
00578         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
00579         hPullGP_Z_rs_stereo[title.str()]->Fill( pullGPZ_rs_stereo );
00580 
00581         double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x())/sqrt(stereoTsosGPEr.cxx()+stereoRhitGPEr.cxx());
00582         double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y())/sqrt(stereoTsosGPEr.cyy()+stereoRhitGPEr.cyy());
00583         double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z())/sqrt(stereoTsosGPEr.czz()+stereoRhitGPEr.czz());
00584         //double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x());
00585         //double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y());
00586         //double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z());
00587 
00588         title.str("");
00589         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
00590         hPullGP_X_tr_stereo[title.str()]->Fill( pullGPX_tr_stereo );
00591         title.str("");
00592         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
00593         hPullGP_Y_tr_stereo[title.str()]->Fill( pullGPY_tr_stereo );
00594         title.str("");
00595         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
00596         hPullGP_Z_tr_stereo[title.str()]->Fill( pullGPZ_tr_stereo );
00597 
00598         double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x())/sqrt(stereoTsosGPEr.cxx());
00599         double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y())/sqrt(stereoTsosGPEr.cyy());
00600         double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z())/sqrt(stereoTsosGPEr.czz());
00601         //double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x());
00602         //double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y());
00603         //double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z());
00604 
00605         title.str("");
00606         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
00607         hPullGP_X_ts_stereo[title.str()]->Fill( pullGPX_ts_stereo );
00608         title.str("");
00609         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00610         hPullGP_Y_ts_stereo[title.str()]->Fill( pullGPY_ts_stereo );
00611         title.str("");
00612         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00613         hPullGP_Z_ts_stereo[title.str()]->Fill( pullGPZ_ts_stereo );
00614 
00615         double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x())/sqrt(stereoTsosGMEr.cxx());
00616         double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y())/sqrt(stereoTsosGMEr.cyy());
00617         double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z())/sqrt(stereoTsosGMEr.czz());
00618         //double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x());
00619         //double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y());
00620         //double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z());
00621 
00622         title.str("");
00623         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
00624         hPullGM_X_ts_stereo[title.str()]->Fill( pullGMX_ts_stereo );
00625         title.str("");
00626         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00627         hPullGM_Y_ts_stereo[title.str()]->Fill( pullGMY_ts_stereo );
00628         title.str("");
00629         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00630         hPullGM_Z_ts_stereo[title.str()]->Fill( pullGMZ_ts_stereo );
00631       }    
00632       lastState = updatedState;
00633       //#endif
00634     }
00635   }
00636   delete hitAssociator;
00637   LogTrace("TestSmoothHits") << "end of event" << std::endl;
00638 }
00639 
00640 void TestSmoothHits::endJob() {
00641   //file->Write();
00642   TDirectory * chi2i = file->mkdir("Chi2_Increment");
00643 
00644   TDirectory * gp_ts = file->mkdir("GP_TSOS-SimHit");
00645   TDirectory * gm_ts = file->mkdir("GM_TSOS-SimHit");
00646   TDirectory * gp_tr = file->mkdir("GP_TSOS-RecHit");
00647   TDirectory * gp_rs = file->mkdir("GP_RecHit-SimHit");
00648 
00649   TDirectory * gp_tsx = gp_ts->mkdir("X");
00650   TDirectory * gp_tsy = gp_ts->mkdir("Y");
00651   TDirectory * gp_tsz = gp_ts->mkdir("Z");
00652   TDirectory * gm_tsx = gm_ts->mkdir("X");
00653   TDirectory * gm_tsy = gm_ts->mkdir("Y");
00654   TDirectory * gm_tsz = gm_ts->mkdir("Z");
00655   TDirectory * gp_trx = gp_tr->mkdir("X");
00656   TDirectory * gp_try = gp_tr->mkdir("Y");
00657   TDirectory * gp_trz = gp_tr->mkdir("Z");
00658   TDirectory * gp_rsx = gp_rs->mkdir("X");
00659   TDirectory * gp_rsy = gp_rs->mkdir("Y");
00660   TDirectory * gp_rsz = gp_rs->mkdir("Z");
00661 
00662   TDirectory * gp_tsx_mono = gp_ts->mkdir("MONOX");
00663   TDirectory * gp_tsy_mono = gp_ts->mkdir("MONOY");
00664   TDirectory * gp_tsz_mono = gp_ts->mkdir("MONOZ");
00665   TDirectory * gm_tsx_mono = gm_ts->mkdir("MONOX");
00666   TDirectory * gm_tsy_mono = gm_ts->mkdir("MONOY");
00667   TDirectory * gm_tsz_mono = gm_ts->mkdir("MONOZ");
00668   TDirectory * gp_trx_mono = gp_tr->mkdir("MONOX");
00669   TDirectory * gp_try_mono = gp_tr->mkdir("MONOY");
00670   TDirectory * gp_trz_mono = gp_tr->mkdir("MONOZ");
00671   TDirectory * gp_rsx_mono = gp_rs->mkdir("MONOX");
00672   TDirectory * gp_rsy_mono = gp_rs->mkdir("MONOY");
00673   TDirectory * gp_rsz_mono = gp_rs->mkdir("MONOZ");
00674 
00675   TDirectory * gp_tsx_stereo = gp_ts->mkdir("STEREOX");
00676   TDirectory * gp_tsy_stereo = gp_ts->mkdir("STEREOY");
00677   TDirectory * gp_tsz_stereo = gp_ts->mkdir("STEREOZ");
00678   TDirectory * gm_tsx_stereo = gm_ts->mkdir("STEREOX");
00679   TDirectory * gm_tsy_stereo = gm_ts->mkdir("STEREOY");
00680   TDirectory * gm_tsz_stereo = gm_ts->mkdir("STEREOZ");
00681   TDirectory * gp_trx_stereo = gp_tr->mkdir("STEREOX");
00682   TDirectory * gp_try_stereo = gp_tr->mkdir("STEREOY");
00683   TDirectory * gp_trz_stereo = gp_tr->mkdir("STEREOZ");
00684   TDirectory * gp_rsx_stereo = gp_rs->mkdir("STEREOX");
00685   TDirectory * gp_rsy_stereo = gp_rs->mkdir("STEREOY");
00686   TDirectory * gp_rsz_stereo = gp_rs->mkdir("STEREOZ");
00687 
00688   chi2i->cd();
00689   hTotChi2Increment->Write();
00690   hChi2_vs_Process->Write();
00691   hChi2_vs_clsize->Write();
00692   for (int i=0; i!=6; i++)
00693     for (int j=0; j!=9; j++){
00694       if (i==0 && j>2) break;
00695       if (i==1 && j>1) break;
00696       if (i==2 && j>3) break;
00697       if (i==3 && j>2) break;
00698       if (i==4 && j>5) break;
00699       if (i==5 && j>8) break;
00700       chi2i->cd();
00701       title.str("");
00702       title << "Chi2Increment_" << i+1 << "-" << j+1;
00703       hChi2Increment[title.str()]->Write();
00704 
00705       gp_ts->cd();
00706       gp_tsx->cd();
00707       title.str("");
00708       title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
00709       hPullGP_X_ts[title.str()]->Write();
00710       gp_tsy->cd();
00711       title.str("");
00712       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
00713       hPullGP_Y_ts[title.str()]->Write();
00714       gp_tsz->cd();
00715       title.str("");
00716       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
00717       hPullGP_Z_ts[title.str()]->Write();
00718 
00719       gm_ts->cd();
00720       gm_tsx->cd();
00721       title.str("");
00722       title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
00723       hPullGM_X_ts[title.str()]->Write();
00724       gm_tsy->cd();
00725       title.str("");
00726       title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
00727       hPullGM_Y_ts[title.str()]->Write();
00728       gm_tsz->cd();
00729       title.str("");
00730       title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
00731       hPullGM_Z_ts[title.str()]->Write();
00732 
00733       gp_tr->cd();
00734       gp_trx->cd();
00735       title.str("");
00736       title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
00737       hPullGP_X_tr[title.str()]->Write();
00738       gp_try->cd();
00739       title.str("");
00740       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
00741       hPullGP_Y_tr[title.str()]->Write();
00742       gp_trz->cd();
00743       title.str("");
00744       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
00745       hPullGP_Z_tr[title.str()]->Write();
00746 
00747       gp_rs->cd();
00748       gp_rsx->cd();
00749       title.str("");
00750       title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
00751       hPullGP_X_rs[title.str()]->Write();
00752       gp_rsy->cd();
00753       title.str("");
00754       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
00755       hPullGP_Y_rs[title.str()]->Write();
00756       gp_rsz->cd();
00757       title.str("");
00758       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
00759       hPullGP_Z_rs[title.str()]->Write();
00760 
00761       if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
00762         //mono
00763         gp_ts->cd();
00764         gp_tsx_mono->cd();
00765         title.str("");
00766         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
00767         hPullGP_X_ts_mono[title.str()]->Write();
00768         gp_tsy_mono->cd();
00769         title.str("");
00770         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00771         hPullGP_Y_ts_mono[title.str()]->Write();
00772         gp_tsz_mono->cd();
00773         title.str("");
00774         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00775         hPullGP_Z_ts_mono[title.str()]->Write();
00776 
00777         gm_ts->cd();
00778         gm_tsx_mono->cd();
00779         title.str("");
00780         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
00781         hPullGM_X_ts_mono[title.str()]->Write();
00782         gm_tsy_mono->cd();
00783         title.str("");
00784         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00785         hPullGM_Y_ts_mono[title.str()]->Write();
00786         gm_tsz_mono->cd();
00787         title.str("");
00788         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00789         hPullGM_Z_ts_mono[title.str()]->Write();
00790 
00791         gp_tr->cd();
00792         gp_trx_mono->cd();
00793         title.str("");
00794         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
00795         hPullGP_X_tr_mono[title.str()]->Write();
00796         gp_try_mono->cd();
00797         title.str("");
00798         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
00799         hPullGP_Y_tr_mono[title.str()]->Write();
00800         gp_trz_mono->cd();
00801         title.str("");
00802         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
00803         hPullGP_Z_tr_mono[title.str()]->Write();
00804 
00805         gp_rs->cd();
00806         gp_rsx_mono->cd();
00807         title.str("");
00808         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
00809         hPullGP_X_rs_mono[title.str()]->Write();
00810         gp_rsy_mono->cd();
00811         title.str("");
00812         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
00813         hPullGP_Y_rs_mono[title.str()]->Write();
00814         gp_rsz_mono->cd();
00815         title.str("");
00816         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
00817         hPullGP_Z_rs_mono[title.str()]->Write();
00818 
00819         //stereo
00820         gp_ts->cd();
00821         gp_tsx_stereo->cd();
00822         title.str("");
00823         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00824         hPullGP_X_ts_stereo[title.str()]->Write();
00825         gp_tsy_stereo->cd();
00826         title.str("");
00827         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00828         hPullGP_Y_ts_stereo[title.str()]->Write();
00829         gp_tsz_stereo->cd();
00830         title.str("");
00831         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00832         hPullGP_Z_ts_stereo[title.str()]->Write();
00833 
00834         gm_ts->cd();
00835         gm_tsx_stereo->cd();
00836         title.str("");
00837         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00838         hPullGM_X_ts_stereo[title.str()]->Write();
00839         gm_tsy_stereo->cd();
00840         title.str("");
00841         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00842         hPullGM_Y_ts_stereo[title.str()]->Write();
00843         gm_tsz_stereo->cd();
00844         title.str("");
00845         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00846         hPullGM_Z_ts_stereo[title.str()]->Write();
00847 
00848         gp_tr->cd();
00849         gp_trx_stereo->cd();
00850         title.str("");
00851         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
00852         hPullGP_X_tr_stereo[title.str()]->Write();
00853         gp_try_stereo->cd();
00854         title.str("");
00855         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
00856         hPullGP_Y_tr_stereo[title.str()]->Write();
00857         gp_trz_stereo->cd();
00858         title.str("");
00859         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
00860         hPullGP_Z_tr_stereo[title.str()]->Write();
00861 
00862         gp_rs->cd();
00863         gp_rsx_stereo->cd();
00864         title.str("");
00865         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
00866         hPullGP_X_rs_stereo[title.str()]->Write();
00867         gp_rsy_stereo->cd();
00868         title.str("");
00869         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
00870         hPullGP_Y_rs_stereo[title.str()]->Write();
00871         gp_rsz_stereo->cd();
00872         title.str("");
00873         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
00874         hPullGP_Z_rs_stereo[title.str()]->Write();
00875       }
00876     }
00877 
00878   file->Close();
00879 }
00880 
00881 //needed by to do the residual for matched hits
00882 //taken from SiStripTrackingRecHitsValid.cc
00883 std::pair<LocalPoint,LocalVector> 
00884 TestSmoothHits::projectHit( const PSimHit& hit, const StripGeomDetUnit* stripDet, const BoundPlane& plane) 
00885 { 
00886    const StripTopology& topol = stripDet->specificTopology();
00887    GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
00888    LocalPoint localHit = plane.toLocal(globalpos);
00889    //track direction
00890    LocalVector locdir=hit.localDirection();
00891    //rotate track in new frame
00892    
00893    GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
00894    LocalVector dir=plane.toLocal(globaldir);
00895    float scale = -localHit.z() / dir.z();
00896    
00897    LocalPoint projectedPos = localHit + scale*dir;
00898    
00899    float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
00900    
00901    LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame
00902    
00903    LocalVector localStripDir( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
00904    
00905    return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
00906 }
00907 
00908 #include "FWCore/Framework/interface/ModuleFactory.h"
00909 #include "FWCore/Framework/interface/MakerMacros.h"
00910 
00911 DEFINE_FWK_MODULE(TestSmoothHits);
/data/refman/pasoursint/CMSSW_5_3_1/src/RecoTracker/DebugTools/plugins/TestSmoothHits.cc
