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     TrajectoryStateTransform transformer;
00212     
00213     DetId  detId(state.detId());
00214     TrajectoryStateOnSurface theTSOS=
00215       transformer.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         CTTRHp tMonoHit = 
00445           theBuilder->build(dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->monoHit());
00446         if (tMonoHit==0) continue;
00447         vector<PSimHit> assMonoSimHits = hitAssociator->associateHit(*tMonoHit->hit());
00448         if (assMonoSimHits.size()==0) continue;
00449         const PSimHit sMonoHit = *(assSimHits.begin());
00450         const Surface * monoSurf = &( tMonoHit->det()->surface() );
00451         if (monoSurf==0) continue;
00452         TSOS monoState = thePropagator->propagate(lastState,*monoSurf);
00453         if (monoState.isValid()==0) continue;
00454 
00455         LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
00456         GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
00457         LocalPoint monoShitLPos = sMonoHit.localPosition();
00458         GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
00459 
00460         GlobalVector monoTsosGMom = monoState.globalMomentum();
00461         GlobalError  monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00462         GlobalPoint  monoTsosGPos = monoState.globalPosition();
00463         GlobalError  monoTsosGPEr = monoState.cartesianError().position();
00464 
00465         GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
00466         GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
00467 
00468         double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x())/sqrt(monoRhitGPEr.cxx());
00469         double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y())/sqrt(monoRhitGPEr.cyy());
00470         double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z())/sqrt(monoRhitGPEr.czz());
00471         //double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x());
00472         //double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y());
00473         //double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z());
00474 
00475         title.str("");
00476         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
00477         hPullGP_X_rs_mono[title.str()]->Fill( pullGPX_rs_mono );
00478         title.str("");
00479         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
00480         hPullGP_Y_rs_mono[title.str()]->Fill( pullGPY_rs_mono );
00481         title.str("");
00482         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
00483         hPullGP_Z_rs_mono[title.str()]->Fill( pullGPZ_rs_mono );
00484 
00485         double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x())/sqrt(monoTsosGPEr.cxx()+monoRhitGPEr.cxx());
00486         double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y())/sqrt(monoTsosGPEr.cyy()+monoRhitGPEr.cyy());
00487         double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z())/sqrt(monoTsosGPEr.czz()+monoRhitGPEr.czz());
00488         //double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x());
00489         //double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y());
00490         //double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z());
00491 
00492         title.str("");
00493         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
00494         hPullGP_X_tr_mono[title.str()]->Fill( pullGPX_tr_mono );
00495         title.str("");
00496         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
00497         hPullGP_Y_tr_mono[title.str()]->Fill( pullGPY_tr_mono );
00498         title.str("");
00499         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
00500         hPullGP_Z_tr_mono[title.str()]->Fill( pullGPZ_tr_mono );
00501 
00502         double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x())/sqrt(monoTsosGPEr.cxx());
00503         double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y())/sqrt(monoTsosGPEr.cyy());
00504         double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z())/sqrt(monoTsosGPEr.czz());
00505         //double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x());
00506         //double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y());
00507         //double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z());
00508 
00509         title.str("");
00510         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
00511         hPullGP_X_ts_mono[title.str()]->Fill( pullGPX_ts_mono );
00512         title.str("");
00513         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
00514         hPullGP_Y_ts_mono[title.str()]->Fill( pullGPY_ts_mono );
00515         title.str("");
00516         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
00517         hPullGP_Z_ts_mono[title.str()]->Fill( pullGPZ_ts_mono );
00518 
00519         double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x())/sqrt(monoTsosGMEr.cxx());
00520         double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y())/sqrt(monoTsosGMEr.cyy());
00521         double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z())/sqrt(monoTsosGMEr.czz());
00522         //double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x());
00523         //double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y());
00524         //double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z());
00525 
00526         title.str("");
00527         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
00528         hPullGM_X_ts_mono[title.str()]->Fill( pullGMX_ts_mono );
00529         title.str("");
00530         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
00531         hPullGM_Y_ts_mono[title.str()]->Fill( pullGMY_ts_mono );
00532         title.str("");
00533         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
00534         hPullGM_Z_ts_mono[title.str()]->Fill( pullGMZ_ts_mono );
00535 
00536         //stereo
00537         LogTrace("TestSmoothHits") << "STEREO HIT" << std::endl;
00538         CTTRHp tStereoHit = 
00539           theBuilder->build(dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->stereoHit());
00540         if (tStereoHit==0) continue;
00541         vector<PSimHit> assStereoSimHits = hitAssociator->associateHit(*tStereoHit->hit());
00542         if (assStereoSimHits.size()==0) continue;
00543         const PSimHit sStereoHit = *(assSimHits.begin());
00544         const Surface * stereoSurf = &( tStereoHit->det()->surface() );
00545         if (stereoSurf==0) continue;
00546         TSOS stereoState = thePropagator->propagate(lastState,*stereoSurf);
00547         if (stereoState.isValid()==0) continue;
00548 
00549         LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
00550         GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
00551         LocalPoint stereoShitLPos = sStereoHit.localPosition();
00552         GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
00553 
00554         GlobalVector stereoTsosGMom = stereoState.globalMomentum();
00555         GlobalError  stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00556         GlobalPoint  stereoTsosGPos = stereoState.globalPosition();
00557         GlobalError  stereoTsosGPEr = stereoState.cartesianError().position();
00558 
00559         GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
00560         GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
00561 
00562         double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/sqrt(stereoRhitGPEr.cxx());
00563         double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/sqrt(stereoRhitGPEr.cyy());
00564         double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/sqrt(stereoRhitGPEr.czz());
00565         //double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x());
00566         //double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y());
00567         //double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z());
00568 
00569         title.str("");
00570         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
00571         hPullGP_X_rs_stereo[title.str()]->Fill( pullGPX_rs_stereo );
00572         title.str("");
00573         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
00574         hPullGP_Y_rs_stereo[title.str()]->Fill( pullGPY_rs_stereo );
00575         title.str("");
00576         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
00577         hPullGP_Z_rs_stereo[title.str()]->Fill( pullGPZ_rs_stereo );
00578 
00579         double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x())/sqrt(stereoTsosGPEr.cxx()+stereoRhitGPEr.cxx());
00580         double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y())/sqrt(stereoTsosGPEr.cyy()+stereoRhitGPEr.cyy());
00581         double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z())/sqrt(stereoTsosGPEr.czz()+stereoRhitGPEr.czz());
00582         //double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x());
00583         //double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y());
00584         //double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z());
00585 
00586         title.str("");
00587         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
00588         hPullGP_X_tr_stereo[title.str()]->Fill( pullGPX_tr_stereo );
00589         title.str("");
00590         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
00591         hPullGP_Y_tr_stereo[title.str()]->Fill( pullGPY_tr_stereo );
00592         title.str("");
00593         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
00594         hPullGP_Z_tr_stereo[title.str()]->Fill( pullGPZ_tr_stereo );
00595 
00596         double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x())/sqrt(stereoTsosGPEr.cxx());
00597         double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y())/sqrt(stereoTsosGPEr.cyy());
00598         double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z())/sqrt(stereoTsosGPEr.czz());
00599         //double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x());
00600         //double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y());
00601         //double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z());
00602 
00603         title.str("");
00604         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
00605         hPullGP_X_ts_stereo[title.str()]->Fill( pullGPX_ts_stereo );
00606         title.str("");
00607         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00608         hPullGP_Y_ts_stereo[title.str()]->Fill( pullGPY_ts_stereo );
00609         title.str("");
00610         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00611         hPullGP_Z_ts_stereo[title.str()]->Fill( pullGPZ_ts_stereo );
00612 
00613         double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x())/sqrt(stereoTsosGMEr.cxx());
00614         double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y())/sqrt(stereoTsosGMEr.cyy());
00615         double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z())/sqrt(stereoTsosGMEr.czz());
00616         //double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x());
00617         //double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y());
00618         //double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z());
00619 
00620         title.str("");
00621         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
00622         hPullGM_X_ts_stereo[title.str()]->Fill( pullGMX_ts_stereo );
00623         title.str("");
00624         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00625         hPullGM_Y_ts_stereo[title.str()]->Fill( pullGMY_ts_stereo );
00626         title.str("");
00627         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00628         hPullGM_Z_ts_stereo[title.str()]->Fill( pullGMZ_ts_stereo );
00629       }    
00630       lastState = updatedState;
00631       //#endif
00632     }
00633   }
00634   delete hitAssociator;
00635   LogTrace("TestSmoothHits") << "end of event" << std::endl;
00636 }
00637 
00638 void TestSmoothHits::endJob() {
00639   //file->Write();
00640   TDirectory * chi2i = file->mkdir("Chi2_Increment");
00641 
00642   TDirectory * gp_ts = file->mkdir("GP_TSOS-SimHit");
00643   TDirectory * gm_ts = file->mkdir("GM_TSOS-SimHit");
00644   TDirectory * gp_tr = file->mkdir("GP_TSOS-RecHit");
00645   TDirectory * gp_rs = file->mkdir("GP_RecHit-SimHit");
00646 
00647   TDirectory * gp_tsx = gp_ts->mkdir("X");
00648   TDirectory * gp_tsy = gp_ts->mkdir("Y");
00649   TDirectory * gp_tsz = gp_ts->mkdir("Z");
00650   TDirectory * gm_tsx = gm_ts->mkdir("X");
00651   TDirectory * gm_tsy = gm_ts->mkdir("Y");
00652   TDirectory * gm_tsz = gm_ts->mkdir("Z");
00653   TDirectory * gp_trx = gp_tr->mkdir("X");
00654   TDirectory * gp_try = gp_tr->mkdir("Y");
00655   TDirectory * gp_trz = gp_tr->mkdir("Z");
00656   TDirectory * gp_rsx = gp_rs->mkdir("X");
00657   TDirectory * gp_rsy = gp_rs->mkdir("Y");
00658   TDirectory * gp_rsz = gp_rs->mkdir("Z");
00659 
00660   TDirectory * gp_tsx_mono = gp_ts->mkdir("MONOX");
00661   TDirectory * gp_tsy_mono = gp_ts->mkdir("MONOY");
00662   TDirectory * gp_tsz_mono = gp_ts->mkdir("MONOZ");
00663   TDirectory * gm_tsx_mono = gm_ts->mkdir("MONOX");
00664   TDirectory * gm_tsy_mono = gm_ts->mkdir("MONOY");
00665   TDirectory * gm_tsz_mono = gm_ts->mkdir("MONOZ");
00666   TDirectory * gp_trx_mono = gp_tr->mkdir("MONOX");
00667   TDirectory * gp_try_mono = gp_tr->mkdir("MONOY");
00668   TDirectory * gp_trz_mono = gp_tr->mkdir("MONOZ");
00669   TDirectory * gp_rsx_mono = gp_rs->mkdir("MONOX");
00670   TDirectory * gp_rsy_mono = gp_rs->mkdir("MONOY");
00671   TDirectory * gp_rsz_mono = gp_rs->mkdir("MONOZ");
00672 
00673   TDirectory * gp_tsx_stereo = gp_ts->mkdir("STEREOX");
00674   TDirectory * gp_tsy_stereo = gp_ts->mkdir("STEREOY");
00675   TDirectory * gp_tsz_stereo = gp_ts->mkdir("STEREOZ");
00676   TDirectory * gm_tsx_stereo = gm_ts->mkdir("STEREOX");
00677   TDirectory * gm_tsy_stereo = gm_ts->mkdir("STEREOY");
00678   TDirectory * gm_tsz_stereo = gm_ts->mkdir("STEREOZ");
00679   TDirectory * gp_trx_stereo = gp_tr->mkdir("STEREOX");
00680   TDirectory * gp_try_stereo = gp_tr->mkdir("STEREOY");
00681   TDirectory * gp_trz_stereo = gp_tr->mkdir("STEREOZ");
00682   TDirectory * gp_rsx_stereo = gp_rs->mkdir("STEREOX");
00683   TDirectory * gp_rsy_stereo = gp_rs->mkdir("STEREOY");
00684   TDirectory * gp_rsz_stereo = gp_rs->mkdir("STEREOZ");
00685 
00686   chi2i->cd();
00687   hTotChi2Increment->Write();
00688   hChi2_vs_Process->Write();
00689   hChi2_vs_clsize->Write();
00690   for (int i=0; i!=6; i++)
00691     for (int j=0; j!=9; j++){
00692       if (i==0 && j>2) break;
00693       if (i==1 && j>1) break;
00694       if (i==2 && j>3) break;
00695       if (i==3 && j>2) break;
00696       if (i==4 && j>5) break;
00697       if (i==5 && j>8) break;
00698       chi2i->cd();
00699       title.str("");
00700       title << "Chi2Increment_" << i+1 << "-" << j+1;
00701       hChi2Increment[title.str()]->Write();
00702 
00703       gp_ts->cd();
00704       gp_tsx->cd();
00705       title.str("");
00706       title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
00707       hPullGP_X_ts[title.str()]->Write();
00708       gp_tsy->cd();
00709       title.str("");
00710       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
00711       hPullGP_Y_ts[title.str()]->Write();
00712       gp_tsz->cd();
00713       title.str("");
00714       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
00715       hPullGP_Z_ts[title.str()]->Write();
00716 
00717       gm_ts->cd();
00718       gm_tsx->cd();
00719       title.str("");
00720       title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
00721       hPullGM_X_ts[title.str()]->Write();
00722       gm_tsy->cd();
00723       title.str("");
00724       title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
00725       hPullGM_Y_ts[title.str()]->Write();
00726       gm_tsz->cd();
00727       title.str("");
00728       title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
00729       hPullGM_Z_ts[title.str()]->Write();
00730 
00731       gp_tr->cd();
00732       gp_trx->cd();
00733       title.str("");
00734       title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
00735       hPullGP_X_tr[title.str()]->Write();
00736       gp_try->cd();
00737       title.str("");
00738       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
00739       hPullGP_Y_tr[title.str()]->Write();
00740       gp_trz->cd();
00741       title.str("");
00742       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
00743       hPullGP_Z_tr[title.str()]->Write();
00744 
00745       gp_rs->cd();
00746       gp_rsx->cd();
00747       title.str("");
00748       title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
00749       hPullGP_X_rs[title.str()]->Write();
00750       gp_rsy->cd();
00751       title.str("");
00752       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
00753       hPullGP_Y_rs[title.str()]->Write();
00754       gp_rsz->cd();
00755       title.str("");
00756       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
00757       hPullGP_Z_rs[title.str()]->Write();
00758 
00759       if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
00760         //mono
00761         gp_ts->cd();
00762         gp_tsx_mono->cd();
00763         title.str("");
00764         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
00765         hPullGP_X_ts_mono[title.str()]->Write();
00766         gp_tsy_mono->cd();
00767         title.str("");
00768         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00769         hPullGP_Y_ts_mono[title.str()]->Write();
00770         gp_tsz_mono->cd();
00771         title.str("");
00772         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00773         hPullGP_Z_ts_mono[title.str()]->Write();
00774 
00775         gm_ts->cd();
00776         gm_tsx_mono->cd();
00777         title.str("");
00778         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
00779         hPullGM_X_ts_mono[title.str()]->Write();
00780         gm_tsy_mono->cd();
00781         title.str("");
00782         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00783         hPullGM_Y_ts_mono[title.str()]->Write();
00784         gm_tsz_mono->cd();
00785         title.str("");
00786         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00787         hPullGM_Z_ts_mono[title.str()]->Write();
00788 
00789         gp_tr->cd();
00790         gp_trx_mono->cd();
00791         title.str("");
00792         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
00793         hPullGP_X_tr_mono[title.str()]->Write();
00794         gp_try_mono->cd();
00795         title.str("");
00796         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
00797         hPullGP_Y_tr_mono[title.str()]->Write();
00798         gp_trz_mono->cd();
00799         title.str("");
00800         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
00801         hPullGP_Z_tr_mono[title.str()]->Write();
00802 
00803         gp_rs->cd();
00804         gp_rsx_mono->cd();
00805         title.str("");
00806         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
00807         hPullGP_X_rs_mono[title.str()]->Write();
00808         gp_rsy_mono->cd();
00809         title.str("");
00810         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
00811         hPullGP_Y_rs_mono[title.str()]->Write();
00812         gp_rsz_mono->cd();
00813         title.str("");
00814         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
00815         hPullGP_Z_rs_mono[title.str()]->Write();
00816 
00817         //stereo
00818         gp_ts->cd();
00819         gp_tsx_stereo->cd();
00820         title.str("");
00821         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00822         hPullGP_X_ts_stereo[title.str()]->Write();
00823         gp_tsy_stereo->cd();
00824         title.str("");
00825         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00826         hPullGP_Y_ts_stereo[title.str()]->Write();
00827         gp_tsz_stereo->cd();
00828         title.str("");
00829         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00830         hPullGP_Z_ts_stereo[title.str()]->Write();
00831 
00832         gm_ts->cd();
00833         gm_tsx_stereo->cd();
00834         title.str("");
00835         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00836         hPullGM_X_ts_stereo[title.str()]->Write();
00837         gm_tsy_stereo->cd();
00838         title.str("");
00839         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00840         hPullGM_Y_ts_stereo[title.str()]->Write();
00841         gm_tsz_stereo->cd();
00842         title.str("");
00843         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00844         hPullGM_Z_ts_stereo[title.str()]->Write();
00845 
00846         gp_tr->cd();
00847         gp_trx_stereo->cd();
00848         title.str("");
00849         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
00850         hPullGP_X_tr_stereo[title.str()]->Write();
00851         gp_try_stereo->cd();
00852         title.str("");
00853         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
00854         hPullGP_Y_tr_stereo[title.str()]->Write();
00855         gp_trz_stereo->cd();
00856         title.str("");
00857         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
00858         hPullGP_Z_tr_stereo[title.str()]->Write();
00859 
00860         gp_rs->cd();
00861         gp_rsx_stereo->cd();
00862         title.str("");
00863         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
00864         hPullGP_X_rs_stereo[title.str()]->Write();
00865         gp_rsy_stereo->cd();
00866         title.str("");
00867         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
00868         hPullGP_Y_rs_stereo[title.str()]->Write();
00869         gp_rsz_stereo->cd();
00870         title.str("");
00871         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
00872         hPullGP_Z_rs_stereo[title.str()]->Write();
00873       }
00874     }
00875 
00876   file->Close();
00877 }
00878 
00879 //needed by to do the residual for matched hits
00880 //taken from SiStripTrackingRecHitsValid.cc
00881 std::pair<LocalPoint,LocalVector> 
00882 TestSmoothHits::projectHit( const PSimHit& hit, const StripGeomDetUnit* stripDet, const BoundPlane& plane) 
00883 { 
00884    const StripTopology& topol = stripDet->specificTopology();
00885    GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
00886    LocalPoint localHit = plane.toLocal(globalpos);
00887    //track direction
00888    LocalVector locdir=hit.localDirection();
00889    //rotate track in new frame
00890    
00891    GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
00892    LocalVector dir=plane.toLocal(globaldir);
00893    float scale = -localHit.z() / dir.z();
00894    
00895    LocalPoint projectedPos = localHit + scale*dir;
00896    
00897    float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
00898    
00899    LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame
00900    
00901    LocalVector localStripDir( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
00902    
00903    return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
00904 }
00905 
00906 #include "FWCore/Framework/interface/ModuleFactory.h"
00907 #include "FWCore/Framework/interface/MakerMacros.h"
00908 
00909 DEFINE_FWK_MODULE(TestSmoothHits);
/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoTracker/DebugTools/plugins/TestSmoothHits.cc
