CMS 3D CMS Logo
RecoTracker » DebugTools » plugins

TestSmoothHits.cc

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