CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/RecoTracker/DebugTools/plugins/TestHits.cc

Go to the documentation of this file.
00001 #include "RecoTracker/DebugTools/interface/TestHits.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 #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
00014 
00015 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
00016 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00017 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00018 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00019 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00020 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00021 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00022 
00023 typedef TrajectoryStateOnSurface TSOS;
00024 typedef TransientTrackingRecHit::ConstRecHitPointer CTTRHp;
00025 using namespace std;
00026 using namespace edm;
00027 
00028 TestHits::TestHits(const edm::ParameterSet& iConfig):
00029   conf_(iConfig){
00030   LogTrace("TestHits") << conf_<< std::endl;
00031   propagatorName = conf_.getParameter<std::string>("Propagator");   
00032   builderName = conf_.getParameter<std::string>("TTRHBuilder");   
00033   srcName = conf_.getParameter<std::string>("src");   
00034   fname = conf_.getParameter<std::string>("Fitter");
00035   mineta = conf_.getParameter<double>("mineta");
00036   maxeta = conf_.getParameter<double>("maxeta");
00037 }
00038 
00039 TestHits::~TestHits(){}
00040 
00041 void TestHits::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  
00050   file = new TFile("testhits.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   hProcess_vs_Chi2  = new TH2F("Process_vs_Chi2","Process_vs_Chi2",1000,0,100,17,-0.5,16.5);  
00188   hClsize_vs_Chi2  = new TH2F("Clsize_vs_Chi2","Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00189 }
00190 
00191 
00192 void TestHits::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00193 {
00194   LogTrace("TestHits") << "\nnew event";
00195 
00196   iEvent.getByLabel(srcName,theTCCollection ); 
00197   hitAssociator = new TrackerHitAssociator(iEvent);
00198 
00199   TrajectoryStateCombiner combiner;
00200 
00201   for (TrackCandidateCollection::const_iterator i=theTCCollection->begin(); i!=theTCCollection->end();i++){
00202 
00203     LogTrace("TestHits") << "\n*****************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     
00211     
00212     DetId  detId(state.detId());
00213     TrajectoryStateOnSurface theTSOS=
00214       trajectoryStateTransform::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     std::vector<Trajectory> result = fit->fit(theTC->seed(), hits, theTSOS);
00228     if (result.size()==0) continue;
00229     std::vector<TrajectoryMeasurement> vtm = result[0].measurements();
00230     double tchi2 = 0;
00231 
00232     TSOS lastState = theTSOS;
00233     for (std::vector<TrajectoryMeasurement>::iterator tm=vtm.begin(); tm!=vtm.end();tm++){
00234 
00235       TransientTrackingRecHit::ConstRecHitPointer rhit = tm->recHit();
00236       if ((rhit)->isValid()==0&&rhit->det()!=0) continue;
00237       LogTrace("TestHits") << "*****************new hit*****************" ;
00238 
00239       int subdetId = rhit->det()->geographicalId().subdetId();
00240       int layerId  = 0;
00241       DetId id = rhit->det()->geographicalId();
00242       if (id.subdetId()==3) layerId = ((TIBDetId)(id)).layer();
00243       if (id.subdetId()==5) layerId = ((TOBDetId)(id)).layer();
00244       if (id.subdetId()==1) layerId = ((PXBDetId)(id)).layer();
00245       if (id.subdetId()==4) layerId = ((TIDDetId)(id)).wheel();
00246       if (id.subdetId()==6) layerId = ((TECDetId)(id)).wheel();
00247       if (id.subdetId()==2) layerId = ((PXFDetId)(id)).disk();
00248       LogTrace("TestHits") << "subdetId=" << subdetId << " layerId=" << layerId ;
00249 
00250       double delta = 99999;
00251       LocalPoint rhitLPv = rhit->localPosition();
00252 
00253       std::vector<PSimHit> assSimHits = hitAssociator->associateHit(*(rhit->hit()));
00254       if (assSimHits.size()==0) continue;
00255       PSimHit shit;
00256       for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
00257         if ((m->localPosition()-rhitLPv).mag()<delta) {
00258           shit=*m;
00259           delta = (m->localPosition()-rhitLPv).mag();
00260         }
00261       }
00262 
00263       TSOS currentState = tm->forwardPredictedState();
00264       if (tm->backwardPredictedState().isValid()) 
00265         currentState = combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
00266       TSOS updatedState = tm->updatedState();
00267       tchi2+=tm->estimate();
00268 
00269       //plot chi2 increment
00270       double chi2increment = tm->estimate();
00271       LogTrace("TestHits") << "tm->estimate()=" << tm->estimate();
00272       title.str("");
00273       title << "Chi2Increment_" << subdetId << "-" << layerId;
00274       hChi2Increment[title.str()]->Fill( chi2increment );
00275       hTotChi2Increment->Fill( chi2increment );
00276       hProcess_vs_Chi2->Fill( chi2increment, shit.processType() );
00277       if (dynamic_cast<const SiPixelRecHit*>(rhit->hit()))      
00278         hClsize_vs_Chi2->Fill( chi2increment, ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size() );
00279       if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit()))    
00280         hClsize_vs_Chi2->Fill( chi2increment, ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() );
00281      
00282       //test hits
00283       const Surface * surf = &( (rhit)->det()->surface() );
00284       LocalVector shitLMom;
00285       LocalPoint shitLPos;
00286       if (dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())) {
00287         double rechitmatchedx = rhit->localPosition().x();
00288         double rechitmatchedy = rhit->localPosition().y();
00289         double mindist = 999999;
00290         double distx, disty;
00291         std::pair<LocalPoint,LocalVector> closestPair;
00292         const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) ((const GluedGeomDet *)(rhit)->det())->stereoDet();
00293         const BoundPlane& plane = (rhit)->det()->surface();
00294         for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++) {
00295           //project simhit;
00296           std::pair<LocalPoint,LocalVector> hitPair = projectHit((*m),stripDet,plane);
00297           distx = fabs(rechitmatchedx - hitPair.first.x());
00298           disty = fabs(rechitmatchedy - hitPair.first.y());
00299           double dist = distx*distx+disty*disty;
00300           if(sqrt(dist)<mindist){
00301             mindist = dist;
00302             closestPair = hitPair;
00303           }
00304         }
00305         shitLPos = closestPair.first;
00306         shitLMom = closestPair.second;
00307       } else {
00308         shitLPos = shit.localPosition();        
00309         shitLMom = shit.momentumAtEntry();
00310       }
00311       GlobalVector shitGMom = surf->toGlobal(shitLMom);
00312       GlobalPoint shitGPos = surf->toGlobal(shitLPos);
00313 
00314       GlobalVector tsosGMom = currentState.globalMomentum();
00315       GlobalError  tsosGMEr(currentState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00316       GlobalPoint  tsosGPos = currentState.globalPosition();
00317       GlobalError  tsosGPEr = currentState.cartesianError().position();
00318 
00319       GlobalPoint rhitGPos = (rhit)->globalPosition();
00320       GlobalError rhitGPEr = (rhit)->globalPositionError();
00321 
00322       double pullGPX_rs = (rhitGPos.x()-shitGPos.x())/sqrt(rhitGPEr.cxx());
00323       double pullGPY_rs = (rhitGPos.y()-shitGPos.y())/sqrt(rhitGPEr.cyy());
00324       double pullGPZ_rs = (rhitGPos.z()-shitGPos.z())/sqrt(rhitGPEr.czz());
00325       //double pullGPX_rs = (rhitGPos.x()-shitGPos.x());
00326       //double pullGPY_rs = (rhitGPos.y()-shitGPos.y());
00327       //double pullGPZ_rs = (rhitGPos.z()-shitGPos.z());
00328       
00329       LogTrace("TestHits") << "rs" << std::endl;
00330       LogVerbatim("TestHits") << "assSimHits.size()=" << assSimHits.size() ;
00331       LogVerbatim("TestHits") << "tsos globalPos   =" << tsosGPos ;
00332       LogVerbatim("TestHits") << "sim hit globalPos=" << shitGPos ;
00333       LogVerbatim("TestHits") << "rec hit globalPos=" << rhitGPos ;
00334       LogVerbatim("TestHits") << "geographicalId   =" << rhit->det()->geographicalId().rawId() ;
00335       LogVerbatim("TestHits") << "surface position =" << surf->position() ;
00336 
00337       title.str("");
00338       title << "PullGP_X_" << subdetId << "-" << layerId << "_rs";
00339       hPullGP_X_rs[title.str()]->Fill( pullGPX_rs );
00340       title.str("");
00341       title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs";
00342       hPullGP_Y_rs[title.str()]->Fill( pullGPY_rs );
00343       title.str("");
00344       title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs";
00345       hPullGP_Z_rs[title.str()]->Fill( pullGPZ_rs );
00346 
00347       double pullGPX_tr = (tsosGPos.x()-rhitGPos.x())/sqrt(tsosGPEr.cxx()+rhitGPEr.cxx());
00348       double pullGPY_tr = (tsosGPos.y()-rhitGPos.y())/sqrt(tsosGPEr.cyy()+rhitGPEr.cyy());
00349       double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z())/sqrt(tsosGPEr.czz()+rhitGPEr.czz());
00350       //double pullGPX_tr = (tsosGPos.x()-rhitGPos.x());
00351       //double pullGPY_tr = (tsosGPos.y()-rhitGPos.y());
00352       //double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z());
00353 
00354       LogTrace("TestHits") << "tr" << std::endl;
00355 
00356       title.str("");
00357       title << "PullGP_X_" << subdetId << "-" << layerId << "_tr";
00358       hPullGP_X_tr[title.str()]->Fill( pullGPX_tr );
00359       title.str("");
00360       title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr";
00361       hPullGP_Y_tr[title.str()]->Fill( pullGPY_tr );
00362       title.str("");
00363       title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr";
00364       hPullGP_Z_tr[title.str()]->Fill( pullGPZ_tr );
00365 
00366       double pullGPX_ts = (tsosGPos.x()-shitGPos.x())/sqrt(tsosGPEr.cxx());
00367       double pullGPY_ts = (tsosGPos.y()-shitGPos.y())/sqrt(tsosGPEr.cyy());
00368       double pullGPZ_ts = (tsosGPos.z()-shitGPos.z())/sqrt(tsosGPEr.czz());
00369       //double pullGPX_ts = (tsosGPos.x()-shitGPos.x());
00370       //double pullGPY_ts = (tsosGPos.y()-shitGPos.y());
00371       //double pullGPZ_ts = (tsosGPos.z()-shitGPos.z());
00372 
00373       LogTrace("TestHits") << "ts1" << std::endl;
00374 
00375       title.str("");
00376       title << "PullGP_X_" << subdetId << "-" << layerId << "_ts";
00377       hPullGP_X_ts[title.str()]->Fill( pullGPX_ts );
00378       title.str("");
00379       title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts";
00380       hPullGP_Y_ts[title.str()]->Fill( pullGPY_ts );
00381       title.str("");
00382       title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts";
00383       hPullGP_Z_ts[title.str()]->Fill( pullGPZ_ts );
00384 
00385       double pullGMX_ts = (tsosGMom.x()-shitGMom.x())/sqrt(tsosGMEr.cxx());
00386       double pullGMY_ts = (tsosGMom.y()-shitGMom.y())/sqrt(tsosGMEr.cyy());
00387       double pullGMZ_ts = (tsosGMom.z()-shitGMom.z())/sqrt(tsosGMEr.czz());
00388       //double pullGMX_ts = (tsosGMom.x()-shitGMom.x());
00389       //double pullGMY_ts = (tsosGMom.y()-shitGMom.y());
00390       //double pullGMZ_ts = (tsosGMom.z()-shitGMom.z());
00391 
00392       LogTrace("TestHits") << "ts2" << std::endl;
00393 
00394       title.str("");
00395       title << "PullGM_X_" << subdetId << "-" << layerId << "_ts";
00396       hPullGM_X_ts[title.str()]->Fill( pullGMX_ts );
00397       title.str("");
00398       title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts";
00399       hPullGM_Y_ts[title.str()]->Fill( pullGMY_ts );
00400       title.str("");
00401       title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts";
00402       hPullGM_Z_ts[title.str()]->Fill( pullGMZ_ts );
00403 
00404       if (dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())) {
00405         //mono
00406         LogTrace("TestHits") << "MONO HIT" << std::endl;
00407         auto m = dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->monoHit();
00408         CTTRHp tMonoHit = theBuilder->build(&m);
00409         if (tMonoHit==0) continue;
00410         vector<PSimHit> assMonoSimHits = hitAssociator->associateHit(*tMonoHit->hit());
00411         if (assMonoSimHits.size()==0) continue;
00412         const PSimHit sMonoHit = *(assSimHits.begin());
00413         const Surface * monoSurf = &( tMonoHit->det()->surface() );
00414         if (monoSurf==0) continue;
00415         TSOS monoState = thePropagator->propagate(lastState,*monoSurf);
00416         if (monoState.isValid()==0) continue;
00417 
00418         LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
00419         GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
00420         LocalPoint monoShitLPos = sMonoHit.localPosition();
00421         GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
00422 
00423         GlobalVector monoTsosGMom = monoState.globalMomentum();
00424         GlobalError  monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00425         GlobalPoint  monoTsosGPos = monoState.globalPosition();
00426         GlobalError  monoTsosGPEr = monoState.cartesianError().position();
00427 
00428         GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
00429         GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
00430 
00431         double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x())/sqrt(monoRhitGPEr.cxx());
00432         double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y())/sqrt(monoRhitGPEr.cyy());
00433         double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z())/sqrt(monoRhitGPEr.czz());
00434         //double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x());
00435         //double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y());
00436         //double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z());
00437 
00438         title.str("");
00439         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
00440         hPullGP_X_rs_mono[title.str()]->Fill( pullGPX_rs_mono );
00441         title.str("");
00442         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
00443         hPullGP_Y_rs_mono[title.str()]->Fill( pullGPY_rs_mono );
00444         title.str("");
00445         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
00446         hPullGP_Z_rs_mono[title.str()]->Fill( pullGPZ_rs_mono );
00447 
00448         double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x())/sqrt(monoTsosGPEr.cxx()+monoRhitGPEr.cxx());
00449         double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y())/sqrt(monoTsosGPEr.cyy()+monoRhitGPEr.cyy());
00450         double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z())/sqrt(monoTsosGPEr.czz()+monoRhitGPEr.czz());
00451         //double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x());
00452         //double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y());
00453         //double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z());
00454 
00455         title.str("");
00456         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
00457         hPullGP_X_tr_mono[title.str()]->Fill( pullGPX_tr_mono );
00458         title.str("");
00459         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
00460         hPullGP_Y_tr_mono[title.str()]->Fill( pullGPY_tr_mono );
00461         title.str("");
00462         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
00463         hPullGP_Z_tr_mono[title.str()]->Fill( pullGPZ_tr_mono );
00464 
00465         double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x())/sqrt(monoTsosGPEr.cxx());
00466         double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y())/sqrt(monoTsosGPEr.cyy());
00467         double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z())/sqrt(monoTsosGPEr.czz());
00468         //double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x());
00469         //double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y());
00470         //double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z());
00471 
00472         title.str("");
00473         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
00474         hPullGP_X_ts_mono[title.str()]->Fill( pullGPX_ts_mono );
00475         title.str("");
00476         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
00477         hPullGP_Y_ts_mono[title.str()]->Fill( pullGPY_ts_mono );
00478         title.str("");
00479         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
00480         hPullGP_Z_ts_mono[title.str()]->Fill( pullGPZ_ts_mono );
00481 
00482         double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x())/sqrt(monoTsosGMEr.cxx());
00483         double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y())/sqrt(monoTsosGMEr.cyy());
00484         double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z())/sqrt(monoTsosGMEr.czz());
00485         //double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x());
00486         //double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y());
00487         //double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z());
00488 
00489         title.str("");
00490         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
00491         hPullGM_X_ts_mono[title.str()]->Fill( pullGMX_ts_mono );
00492         title.str("");
00493         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
00494         hPullGM_Y_ts_mono[title.str()]->Fill( pullGMY_ts_mono );
00495         title.str("");
00496         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
00497         hPullGM_Z_ts_mono[title.str()]->Fill( pullGMZ_ts_mono );
00498 
00499         //stereo
00500         LogTrace("TestHits") << "STEREO HIT" << std::endl;
00501         auto s = dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->stereoHit();
00502         CTTRHp tStereoHit = 
00503           theBuilder->build(&s);
00504         if (tStereoHit==0) continue;
00505         vector<PSimHit> assStereoSimHits = hitAssociator->associateHit(*tStereoHit->hit());
00506         if (assStereoSimHits.size()==0) continue;
00507         const PSimHit sStereoHit = *(assSimHits.begin());
00508         const Surface * stereoSurf = &( tStereoHit->det()->surface() );
00509         if (stereoSurf==0) continue;
00510         TSOS stereoState = thePropagator->propagate(lastState,*stereoSurf);
00511         if (stereoState.isValid()==0) continue;
00512 
00513         LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
00514         GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
00515         LocalPoint stereoShitLPos = sStereoHit.localPosition();
00516         GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
00517 
00518         GlobalVector stereoTsosGMom = stereoState.globalMomentum();
00519         GlobalError  stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00520         GlobalPoint  stereoTsosGPos = stereoState.globalPosition();
00521         GlobalError  stereoTsosGPEr = stereoState.cartesianError().position();
00522 
00523         GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
00524         GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
00525 
00526         double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/sqrt(stereoRhitGPEr.cxx());
00527         double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/sqrt(stereoRhitGPEr.cyy());
00528         double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/sqrt(stereoRhitGPEr.czz());
00529         //double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x());
00530         //double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y());
00531         //double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z());
00532 
00533         title.str("");
00534         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
00535         hPullGP_X_rs_stereo[title.str()]->Fill( pullGPX_rs_stereo );
00536         title.str("");
00537         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
00538         hPullGP_Y_rs_stereo[title.str()]->Fill( pullGPY_rs_stereo );
00539         title.str("");
00540         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
00541         hPullGP_Z_rs_stereo[title.str()]->Fill( pullGPZ_rs_stereo );
00542 
00543         double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x())/sqrt(stereoTsosGPEr.cxx()+stereoRhitGPEr.cxx());
00544         double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y())/sqrt(stereoTsosGPEr.cyy()+stereoRhitGPEr.cyy());
00545         double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z())/sqrt(stereoTsosGPEr.czz()+stereoRhitGPEr.czz());
00546         //double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x());
00547         //double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y());
00548         //double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z());
00549 
00550         title.str("");
00551         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
00552         hPullGP_X_tr_stereo[title.str()]->Fill( pullGPX_tr_stereo );
00553         title.str("");
00554         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
00555         hPullGP_Y_tr_stereo[title.str()]->Fill( pullGPY_tr_stereo );
00556         title.str("");
00557         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
00558         hPullGP_Z_tr_stereo[title.str()]->Fill( pullGPZ_tr_stereo );
00559 
00560         double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x())/sqrt(stereoTsosGPEr.cxx());
00561         double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y())/sqrt(stereoTsosGPEr.cyy());
00562         double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z())/sqrt(stereoTsosGPEr.czz());
00563         //double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x());
00564         //double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y());
00565         //double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z());
00566 
00567         title.str("");
00568         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
00569         hPullGP_X_ts_stereo[title.str()]->Fill( pullGPX_ts_stereo );
00570         title.str("");
00571         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00572         hPullGP_Y_ts_stereo[title.str()]->Fill( pullGPY_ts_stereo );
00573         title.str("");
00574         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00575         hPullGP_Z_ts_stereo[title.str()]->Fill( pullGPZ_ts_stereo );
00576 
00577         double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x())/sqrt(stereoTsosGMEr.cxx());
00578         double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y())/sqrt(stereoTsosGMEr.cyy());
00579         double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z())/sqrt(stereoTsosGMEr.czz());
00580         //double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x());
00581         //double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y());
00582         //double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z());
00583 
00584         title.str("");
00585         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
00586         hPullGM_X_ts_stereo[title.str()]->Fill( pullGMX_ts_stereo );
00587         title.str("");
00588         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00589         hPullGM_Y_ts_stereo[title.str()]->Fill( pullGMY_ts_stereo );
00590         title.str("");
00591         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00592         hPullGM_Z_ts_stereo[title.str()]->Fill( pullGMZ_ts_stereo );
00593       }    
00594       lastState = updatedState;
00595     }
00596     LogTrace("TestHits") << "traj chi2="  << tchi2 ;
00597     LogTrace("TestHits") << "track chi2=" << result[0].chiSquared() ;
00598   }
00599   delete hitAssociator;
00600   LogTrace("TestHits") << "end of event" << std::endl;
00601 }
00602 //     TSOS lastState = theTSOS;
00603 //     for (std::vector<TrajectoryMeasurement>::iterator tm=vtm.begin(); tm!=vtm.end();tm++){
00604 
00605 //       TransientTrackingRecHit::ConstRecHitPointer rhit = tm->recHit();
00606 //       if ((rhit)->isValid()==0&&rhit->det()!=0) continue;
00607     
00608 // //     //test hits
00609 // //     TSOS lastState, currentState, updatedState;
00610 // //     updatedState=theTSOS;
00611 // //     for (TransientTrackingRecHit::RecHitContainer::iterator rhit=hits.begin()+1;rhit!=hits.end();++rhit){
00612 
00613 //       lastState=updatedState;
00614 
00615 //       LogTrace("TestHits") << "new hit" << std::endl;
00616 
00617 //       if ((*rhit)->isValid()==0) continue;
00618 
00619 //       int subdetId = (*rhit)->det()->geographicalId().subdetId();
00620 //       int layerId  = 0;
00621 //       DetId id = (*rhit)->det()->geographicalId();
00622 //       if (id.subdetId()==3) layerId = ((TIBDetId)(id)).layer();
00623 //       if (id.subdetId()==5) layerId = ((TOBDetId)(id)).layer();
00624 //       if (id.subdetId()==1) layerId = ((PXBDetId)(id)).layer();
00625 //       if (id.subdetId()==4) layerId = ((TIDDetId)(id)).wheel();
00626 //       if (id.subdetId()==6) layerId = ((TECDetId)(id)).wheel();
00627 //       if (id.subdetId()==2) layerId = ((PXFDetId)(id)).disk();
00628 //       const Surface * surf = &( (*rhit)->det()->surface() );
00629 //       currentState=thePropagator->propagate(lastState,*surf);        
00630 //       if (currentState.isValid()==0) continue;
00631 //       updatedState=theUpdator->update(currentState,**rhit);
00632 
00633 //       double delta = 99999;
00634 //       LocalPoint rhitLP = rhit->localPosition();
00635 
00636 //       std::vector<PSimHit> assSimHits = hitAssociator->associateHit(*(*rhit)->hit());
00637 //       if (assSimHits.size()==0) continue;
00638 //       PSimHit shit;
00639 //       for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
00640 //      if ((*m-rhitLP).mag()<delta) {
00641 //        shit=*m;
00642 //        delta = (*m-rhitLP).mag();
00643 //      }
00644 //       }
00645 
00646 //       LocalVector shitLMom = shit.momentumAtEntry();
00647 //       GlobalVector shitGMom = surf->toGlobal(shitLMom);
00648 //       LocalPoint shitLPos = shit.localPosition();
00649 //       GlobalPoint shitGPos = surf->toGlobal(shitLPos);
00650 
00651 //       GlobalVector tsosGMom = currentState.globalMomentum();
00652 //       GlobalError  tsosGMEr(currentState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00653 //       GlobalPoint  tsosGPos = currentState.globalPosition();
00654 //       GlobalError  tsosGPEr = currentState.cartesianError().position();
00655 
00656 //       GlobalPoint rhitGPos = (*rhit)->globalPosition();
00657 //       GlobalError rhitGPEr = (*rhit)->globalPositionError();
00658 
00659 //       double pullGPX_rs = (rhitGPos.x()-shitGPos.x())/sqrt(rhitGPEr.cxx());
00660 //       double pullGPY_rs = (rhitGPos.y()-shitGPos.y())/sqrt(rhitGPEr.cyy());
00661 //       double pullGPZ_rs = (rhitGPos.z()-shitGPos.z())/sqrt(rhitGPEr.czz());
00662 // //       double pullGPX_rs = (rhitGPos.x()-shitGPos.x())/*/sqrt(rhitGPEr.cxx())*/;
00663 // //       double pullGPY_rs = (rhitGPos.y()-shitGPos.y())/*/sqrt(rhitGPEr.cyy())*/;
00664 // //       double pullGPZ_rs = (rhitGPos.z()-shitGPos.z())/*/sqrt(rhitGPEr.czz())*/;
00665 
00666 //       //plot chi2 increment
00667 //       MeasurementExtractor me(currentState);
00668 //       double chi2increment = computeChi2Increment(me,*rhit);
00669 //       LogTrace("TestHits") << "chi2increment=" << chi2increment << std::endl;
00670 //       title.str("");
00671 //       title << "Chi2Increment_" << subdetId << "-" << layerId;
00672 //       hChi2Increment[title.str()]->Fill( chi2increment );
00673 //       hTotChi2Increment->Fill( chi2increment );
00674 //       hChi2_vs_Process->Fill( chi2increment, shit.processType() );
00675 //       if (dynamic_cast<const SiPixelRecHit*>((*rhit)->hit()))        
00676 //      hChi2_vs_clsize->Fill( chi2increment, ((const SiPixelRecHit*)(*rhit)->hit())->cluster()->size() );
00677 //       if (dynamic_cast<const SiStripRecHit2D*>((*rhit)->hit()))      
00678 //      hChi2_vs_clsize->Fill( chi2increment, ((const SiStripRecHit2D*)(*rhit)->hit())->cluster()->amplitudes().size() );
00679       
00680 //       LogTrace("TestHits") << "rs" << std::endl;
00681 
00682 //       title.str("");
00683 //       title << "PullGP_X_" << subdetId << "-" << layerId << "_rs";
00684 //       hPullGP_X_rs[title.str()]->Fill( pullGPX_rs );
00685 //       title.str("");
00686 //       title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs";
00687 //       hPullGP_Y_rs[title.str()]->Fill( pullGPY_rs );
00688 //       title.str("");
00689 //       title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs";
00690 //       hPullGP_Z_rs[title.str()]->Fill( pullGPZ_rs );
00691 
00692 //       double pullGPX_tr = (tsosGPos.x()-rhitGPos.x())/sqrt(tsosGPEr.cxx()+rhitGPEr.cxx());
00693 //       double pullGPY_tr = (tsosGPos.y()-rhitGPos.y())/sqrt(tsosGPEr.cyy()+rhitGPEr.cyy());
00694 //       double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z())/sqrt(tsosGPEr.czz()+rhitGPEr.czz());
00695 // //       double pullGPX_tr = (tsosGPos.x()-rhitGPos.x())/*/sqrt(tsosGPEr.cxx()+rhitGPEr.cxx())*/;
00696 // //       double pullGPY_tr = (tsosGPos.y()-rhitGPos.y())/*/sqrt(tsosGPEr.cyy()+rhitGPEr.cyy())*/;
00697 // //       double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z())/*/sqrt(tsosGPEr.czz()+rhitGPEr.czz())*/;
00698 
00699 //       LogTrace("TestHits") << "tr" << std::endl;
00700 
00701 //       title.str("");
00702 //       title << "PullGP_X_" << subdetId << "-" << layerId << "_tr";
00703 //       hPullGP_X_tr[title.str()]->Fill( pullGPX_tr );
00704 //       title.str("");
00705 //       title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr";
00706 //       hPullGP_Y_tr[title.str()]->Fill( pullGPY_tr );
00707 //       title.str("");
00708 //       title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr";
00709 //       hPullGP_Z_tr[title.str()]->Fill( pullGPZ_tr );
00710 
00711 //       double pullGPX_ts = (tsosGPos.x()-shitGPos.x())/sqrt(tsosGPEr.cxx());
00712 //       double pullGPY_ts = (tsosGPos.y()-shitGPos.y())/sqrt(tsosGPEr.cyy());
00713 //       double pullGPZ_ts = (tsosGPos.z()-shitGPos.z())/sqrt(tsosGPEr.czz());
00714 // //       double pullGPX_ts = (tsosGPos.x()-shitGPos.x())/*/sqrt(tsosGPEr.cxx())*/;
00715 // //       double pullGPY_ts = (tsosGPos.y()-shitGPos.y())/*/sqrt(tsosGPEr.cyy())*/;
00716 // //       double pullGPZ_ts = (tsosGPos.z()-shitGPos.z())/*/sqrt(tsosGPEr.czz())*/;
00717 
00718 //       LogTrace("TestHits") << "ts1" << std::endl;
00719 
00720 //       title.str("");
00721 //       title << "PullGP_X_" << subdetId << "-" << layerId << "_ts";
00722 //       hPullGP_X_ts[title.str()]->Fill( pullGPX_ts );
00723 //       title.str("");
00724 //       title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts";
00725 //       hPullGP_Y_ts[title.str()]->Fill( pullGPY_ts );
00726 //       title.str("");
00727 //       title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts";
00728 //       hPullGP_Z_ts[title.str()]->Fill( pullGPZ_ts );
00729 
00730 //       double pullGMX_ts = (tsosGMom.x()-shitGMom.x())/sqrt(tsosGMEr.cxx());
00731 //       double pullGMY_ts = (tsosGMom.y()-shitGMom.y())/sqrt(tsosGMEr.cyy());
00732 //       double pullGMZ_ts = (tsosGMom.z()-shitGMom.z())/sqrt(tsosGMEr.czz());
00733 // //       double pullGMX_ts = (tsosGMom.x()-shitGMom.x())/*/sqrt(tsosGMEr.cxx())*/;
00734 // //       double pullGMY_ts = (tsosGMom.y()-shitGMom.y())/*/sqrt(tsosGMEr.cyy())*/;
00735 // //       double pullGMZ_ts = (tsosGMom.z()-shitGMom.z())/*/sqrt(tsosGMEr.czz())*/;
00736 
00737 //       LogTrace("TestHits") << "ts2" << std::endl;
00738 
00739 //       title.str("");
00740 //       title << "PullGM_X_" << subdetId << "-" << layerId << "_ts";
00741 //       hPullGM_X_ts[title.str()]->Fill( pullGMX_ts );
00742 //       title.str("");
00743 //       title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts";
00744 //       hPullGM_Y_ts[title.str()]->Fill( pullGMY_ts );
00745 //       title.str("");
00746 //       title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts";
00747 //       hPullGM_Z_ts[title.str()]->Fill( pullGMZ_ts );
00748 
00749 //       if (dynamic_cast<const SiStripMatchedRecHit2D*>((*rhit)->hit())) {
00750 //      //mono
00751 //      LogTrace("TestHits") << "MONO HIT" << std::endl;
00752 //      CTTRHp tMonoHit = 
00753 //        theBuilder->build(dynamic_cast<const SiStripMatchedRecHit2D*>((*rhit)->hit())->monoHit());
00754 //      if (tMonoHit==0) continue;
00755 //      vector<PSimHit> assMonoSimHits = hitAssociator->associateHit(*tMonoHit->hit());
00756 //      if (assMonoSimHits.size()==0) continue;
00757 //      const PSimHit sMonoHit = *(assSimHits.begin());
00758 //      const Surface * monoSurf = &( tMonoHit->det()->surface() );
00759 //      if (monoSurf==0) continue;
00760 //      TSOS monoState = thePropagator->propagate(lastState,*monoSurf);
00761 //      if (monoState.isValid()==0) continue;
00762 
00763 //      LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
00764 //      GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
00765 //      LocalPoint monoShitLPos = sMonoHit.localPosition();
00766 //      GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
00767 
00768 //      GlobalVector monoTsosGMom = monoState.globalMomentum();
00769 //      GlobalError  monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00770 //      GlobalPoint  monoTsosGPos = monoState.globalPosition();
00771 //      GlobalError  monoTsosGPEr = monoState.cartesianError().position();
00772 
00773 //      GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
00774 //      GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
00775 
00776 //      double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x())/sqrt(monoRhitGPEr.cxx());
00777 //      double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y())/sqrt(monoRhitGPEr.cyy());
00778 //      double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z())/sqrt(monoRhitGPEr.czz());
00779 // //   double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x())/*/sqrt(monoRhitGPEr.cxx())*/;
00780 // //   double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y())/*/sqrt(monoRhitGPEr.cyy())*/;
00781 // //   double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z())/*/sqrt(monoRhitGPEr.czz())*/;
00782 
00783 //      title.str("");
00784 //      title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
00785 //      hPullGP_X_rs_mono[title.str()]->Fill( pullGPX_rs_mono );
00786 //      title.str("");
00787 //      title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
00788 //      hPullGP_Y_rs_mono[title.str()]->Fill( pullGPY_rs_mono );
00789 //      title.str("");
00790 //      title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
00791 //      hPullGP_Z_rs_mono[title.str()]->Fill( pullGPZ_rs_mono );
00792 
00793 //      double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x())/sqrt(monoTsosGPEr.cxx()+monoRhitGPEr.cxx());
00794 //      double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y())/sqrt(monoTsosGPEr.cyy()+monoRhitGPEr.cyy());
00795 //      double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z())/sqrt(monoTsosGPEr.czz()+monoRhitGPEr.czz());
00796 // //   double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x())/*/sqrt(monoTsosGPEr.cxx()+monoRhitGPEr.cxx())*/;
00797 // //   double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y())/*/sqrt(monoTsosGPEr.cyy()+monoRhitGPEr.cyy())*/;
00798 // //   double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z())/*/sqrt(monoTsosGPEr.czz()+monoRhitGPEr.czz())*/;
00799 
00800 //      title.str("");
00801 //      title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
00802 //      hPullGP_X_tr_mono[title.str()]->Fill( pullGPX_tr_mono );
00803 //      title.str("");
00804 //      title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
00805 //      hPullGP_Y_tr_mono[title.str()]->Fill( pullGPY_tr_mono );
00806 //      title.str("");
00807 //      title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
00808 //      hPullGP_Z_tr_mono[title.str()]->Fill( pullGPZ_tr_mono );
00809 
00810 //      double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x())/sqrt(monoTsosGPEr.cxx());
00811 //      double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y())/sqrt(monoTsosGPEr.cyy());
00812 //      double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z())/sqrt(monoTsosGPEr.czz());
00813 // //   double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x())/*/sqrt(monoTsosGPEr.cxx())*/;
00814 // //   double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y())/*/sqrt(monoTsosGPEr.cyy())*/;
00815 // //   double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z())/*/sqrt(monoTsosGPEr.czz())*/;
00816 
00817 //      title.str("");
00818 //      title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
00819 //      hPullGP_X_ts_mono[title.str()]->Fill( pullGPX_ts_mono );
00820 //      title.str("");
00821 //      title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
00822 //      hPullGP_Y_ts_mono[title.str()]->Fill( pullGPY_ts_mono );
00823 //      title.str("");
00824 //      title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
00825 //      hPullGP_Z_ts_mono[title.str()]->Fill( pullGPZ_ts_mono );
00826 
00827 //      double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x())/sqrt(monoTsosGMEr.cxx());
00828 //      double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y())/sqrt(monoTsosGMEr.cyy());
00829 //      double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z())/sqrt(monoTsosGMEr.czz());
00830 // //   double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x())/*/sqrt(monoTsosGMEr.cxx())*/;
00831 // //   double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y())/*/sqrt(monoTsosGMEr.cyy())*/;
00832 // //   double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z())/*/sqrt(monoTsosGMEr.czz())*/;
00833 
00834 //      title.str("");
00835 //      title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
00836 //      hPullGM_X_ts_mono[title.str()]->Fill( pullGMX_ts_mono );
00837 //      title.str("");
00838 //      title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
00839 //      hPullGM_Y_ts_mono[title.str()]->Fill( pullGMY_ts_mono );
00840 //      title.str("");
00841 //      title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
00842 //      hPullGM_Z_ts_mono[title.str()]->Fill( pullGMZ_ts_mono );
00843 
00844 //      //stereo
00845 //      LogTrace("TestHits") << "STEREO HIT" << std::endl;
00846 //      CTTRHp tStereoHit = 
00847 //        theBuilder->build(dynamic_cast<const SiStripMatchedRecHit2D*>((*rhit)->hit())->stereoHit());
00848 //      if (tStereoHit==0) continue;
00849 //      vector<PSimHit> assStereoSimHits = hitAssociator->associateHit(*tStereoHit->hit());
00850 //      if (assStereoSimHits.size()==0) continue;
00851 //      const PSimHit sStereoHit = *(assSimHits.begin());
00852 //      const Surface * stereoSurf = &( tStereoHit->det()->surface() );
00853 //      if (stereoSurf==0) continue;
00854 //      TSOS stereoState = thePropagator->propagate(lastState,*stereoSurf);
00855 //      if (stereoState.isValid()==0) continue;
00856 
00857 //      LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
00858 //      GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
00859 //      LocalPoint stereoShitLPos = sStereoHit.localPosition();
00860 //      GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
00861 
00862 //      GlobalVector stereoTsosGMom = stereoState.globalMomentum();
00863 //      GlobalError  stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00864 //      GlobalPoint  stereoTsosGPos = stereoState.globalPosition();
00865 //      GlobalError  stereoTsosGPEr = stereoState.cartesianError().position();
00866 
00867 //      GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
00868 //      GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
00869 
00870 //      double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/sqrt(stereoRhitGPEr.cxx());
00871 //      double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/sqrt(stereoRhitGPEr.cyy());
00872 //      double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/sqrt(stereoRhitGPEr.czz());
00873 // //   double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/*/sqrt(stereoRhitGPEr.cxx())*/;
00874 // //   double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/*/sqrt(stereoRhitGPEr.cyy())*/;
00875 // //   double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/*/sqrt(stereoRhitGPEr.czz())*/;
00876 
00877 //      title.str("");
00878 //      title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
00879 //      hPullGP_X_rs_stereo[title.str()]->Fill( pullGPX_rs_stereo );
00880 //      title.str("");
00881 //      title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
00882 //      hPullGP_Y_rs_stereo[title.str()]->Fill( pullGPY_rs_stereo );
00883 //      title.str("");
00884 //      title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
00885 //      hPullGP_Z_rs_stereo[title.str()]->Fill( pullGPZ_rs_stereo );
00886 
00887 //      double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x())/sqrt(stereoTsosGPEr.cxx()+stereoRhitGPEr.cxx());
00888 //      double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y())/sqrt(stereoTsosGPEr.cyy()+stereoRhitGPEr.cyy());
00889 //      double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z())/sqrt(stereoTsosGPEr.czz()+stereoRhitGPEr.czz());
00890 // //   double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x())/*/sqrt(stereoTsosGPEr.cxx()+stereoRhitGPEr.cxx())*/;
00891 // //   double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y())/*/sqrt(stereoTsosGPEr.cyy()+stereoRhitGPEr.cyy())*/;
00892 // //   double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z())/*/sqrt(stereoTsosGPEr.czz()+stereoRhitGPEr.czz())*/;
00893 
00894 //      title.str("");
00895 //      title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
00896 //      hPullGP_X_tr_stereo[title.str()]->Fill( pullGPX_tr_stereo );
00897 //      title.str("");
00898 //      title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
00899 //      hPullGP_Y_tr_stereo[title.str()]->Fill( pullGPY_tr_stereo );
00900 //      title.str("");
00901 //      title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
00902 //      hPullGP_Z_tr_stereo[title.str()]->Fill( pullGPZ_tr_stereo );
00903 
00904 //      double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x())/sqrt(stereoTsosGPEr.cxx());
00905 //      double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y())/sqrt(stereoTsosGPEr.cyy());
00906 //      double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z())/sqrt(stereoTsosGPEr.czz());
00907 // //   double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x())/*/sqrt(stereoTsosGPEr.cxx())*/;
00908 // //   double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y())/*/sqrt(stereoTsosGPEr.cyy())*/;
00909 // //   double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z())/*/sqrt(stereoTsosGPEr.czz())*/;
00910 
00911 //      title.str("");
00912 //      title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
00913 //      hPullGP_X_ts_stereo[title.str()]->Fill( pullGPX_ts_stereo );
00914 //      title.str("");
00915 //      title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00916 //      hPullGP_Y_ts_stereo[title.str()]->Fill( pullGPY_ts_stereo );
00917 //      title.str("");
00918 //      title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00919 //      hPullGP_Z_ts_stereo[title.str()]->Fill( pullGPZ_ts_stereo );
00920 
00921 //      double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x())/sqrt(stereoTsosGMEr.cxx());
00922 //      double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y())/sqrt(stereoTsosGMEr.cyy());
00923 //      double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z())/sqrt(stereoTsosGMEr.czz());
00924 // //   double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x())/*/sqrt(stereoTsosGMEr.cxx())*/;
00925 // //   double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y())/*/sqrt(stereoTsosGMEr.cyy())*/;
00926 // //   double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z())/*/sqrt(stereoTsosGMEr.czz())*/;
00927 
00928 //      title.str("");
00929 //      title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
00930 //      hPullGM_X_ts_stereo[title.str()]->Fill( pullGMX_ts_stereo );
00931 //      title.str("");
00932 //      title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00933 //      hPullGM_Y_ts_stereo[title.str()]->Fill( pullGMY_ts_stereo );
00934 //      title.str("");
00935 //      title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00936 //      hPullGM_Z_ts_stereo[title.str()]->Fill( pullGMZ_ts_stereo );
00937 //       }
00938 //     }
00939 //   }
00940 //   delete hitAssociator;
00941 //   LogTrace("TestHits") << "end of event" << std::endl;
00942 // }
00943 
00944 void TestHits::endJob() {
00945   //file->Write();
00946   TDirectory * chi2i = file->mkdir("Chi2_Increment");
00947 
00948   TDirectory * gp_ts = file->mkdir("GP_TSOS-SimHit");
00949   TDirectory * gm_ts = file->mkdir("GM_TSOS-SimHit");
00950   TDirectory * gp_tr = file->mkdir("GP_TSOS-RecHit");
00951   TDirectory * gp_rs = file->mkdir("GP_RecHit-SimHit");
00952 
00953   TDirectory * gp_tsx = gp_ts->mkdir("X");
00954   TDirectory * gp_tsy = gp_ts->mkdir("Y");
00955   TDirectory * gp_tsz = gp_ts->mkdir("Z");
00956   TDirectory * gm_tsx = gm_ts->mkdir("X");
00957   TDirectory * gm_tsy = gm_ts->mkdir("Y");
00958   TDirectory * gm_tsz = gm_ts->mkdir("Z");
00959   TDirectory * gp_trx = gp_tr->mkdir("X");
00960   TDirectory * gp_try = gp_tr->mkdir("Y");
00961   TDirectory * gp_trz = gp_tr->mkdir("Z");
00962   TDirectory * gp_rsx = gp_rs->mkdir("X");
00963   TDirectory * gp_rsy = gp_rs->mkdir("Y");
00964   TDirectory * gp_rsz = gp_rs->mkdir("Z");
00965 
00966   TDirectory * gp_tsx_mono = gp_ts->mkdir("MONOX");
00967   TDirectory * gp_tsy_mono = gp_ts->mkdir("MONOY");
00968   TDirectory * gp_tsz_mono = gp_ts->mkdir("MONOZ");
00969   TDirectory * gm_tsx_mono = gm_ts->mkdir("MONOX");
00970   TDirectory * gm_tsy_mono = gm_ts->mkdir("MONOY");
00971   TDirectory * gm_tsz_mono = gm_ts->mkdir("MONOZ");
00972   TDirectory * gp_trx_mono = gp_tr->mkdir("MONOX");
00973   TDirectory * gp_try_mono = gp_tr->mkdir("MONOY");
00974   TDirectory * gp_trz_mono = gp_tr->mkdir("MONOZ");
00975   TDirectory * gp_rsx_mono = gp_rs->mkdir("MONOX");
00976   TDirectory * gp_rsy_mono = gp_rs->mkdir("MONOY");
00977   TDirectory * gp_rsz_mono = gp_rs->mkdir("MONOZ");
00978 
00979   TDirectory * gp_tsx_stereo = gp_ts->mkdir("STEREOX");
00980   TDirectory * gp_tsy_stereo = gp_ts->mkdir("STEREOY");
00981   TDirectory * gp_tsz_stereo = gp_ts->mkdir("STEREOZ");
00982   TDirectory * gm_tsx_stereo = gm_ts->mkdir("STEREOX");
00983   TDirectory * gm_tsy_stereo = gm_ts->mkdir("STEREOY");
00984   TDirectory * gm_tsz_stereo = gm_ts->mkdir("STEREOZ");
00985   TDirectory * gp_trx_stereo = gp_tr->mkdir("STEREOX");
00986   TDirectory * gp_try_stereo = gp_tr->mkdir("STEREOY");
00987   TDirectory * gp_trz_stereo = gp_tr->mkdir("STEREOZ");
00988   TDirectory * gp_rsx_stereo = gp_rs->mkdir("STEREOX");
00989   TDirectory * gp_rsy_stereo = gp_rs->mkdir("STEREOY");
00990   TDirectory * gp_rsz_stereo = gp_rs->mkdir("STEREOZ");
00991 
00992   chi2i->cd();
00993   hTotChi2Increment->Write();
00994   hProcess_vs_Chi2->Write();
00995   hClsize_vs_Chi2->Write();
00996   for (int i=0; i!=6; i++)
00997     for (int j=0; j!=9; j++){
00998       if (i==0 && j>2) break;
00999       if (i==1 && j>1) break;
01000       if (i==2 && j>3) break;
01001       if (i==3 && j>2) break;
01002       if (i==4 && j>5) break;
01003       if (i==5 && j>8) break;
01004       chi2i->cd();
01005       title.str("");
01006       title << "Chi2Increment_" << i+1 << "-" << j+1;
01007       hChi2Increment[title.str()]->Write();
01008 
01009       gp_ts->cd();
01010       gp_tsx->cd();
01011       title.str("");
01012       title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
01013       hPullGP_X_ts[title.str()]->Write();
01014       gp_tsy->cd();
01015       title.str("");
01016       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
01017       hPullGP_Y_ts[title.str()]->Write();
01018       gp_tsz->cd();
01019       title.str("");
01020       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
01021       hPullGP_Z_ts[title.str()]->Write();
01022 
01023       gm_ts->cd();
01024       gm_tsx->cd();
01025       title.str("");
01026       title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
01027       hPullGM_X_ts[title.str()]->Write();
01028       gm_tsy->cd();
01029       title.str("");
01030       title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
01031       hPullGM_Y_ts[title.str()]->Write();
01032       gm_tsz->cd();
01033       title.str("");
01034       title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
01035       hPullGM_Z_ts[title.str()]->Write();
01036 
01037       gp_tr->cd();
01038       gp_trx->cd();
01039       title.str("");
01040       title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
01041       hPullGP_X_tr[title.str()]->Write();
01042       gp_try->cd();
01043       title.str("");
01044       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
01045       hPullGP_Y_tr[title.str()]->Write();
01046       gp_trz->cd();
01047       title.str("");
01048       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
01049       hPullGP_Z_tr[title.str()]->Write();
01050 
01051       gp_rs->cd();
01052       gp_rsx->cd();
01053       title.str("");
01054       title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
01055       hPullGP_X_rs[title.str()]->Write();
01056       gp_rsy->cd();
01057       title.str("");
01058       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
01059       hPullGP_Y_rs[title.str()]->Write();
01060       gp_rsz->cd();
01061       title.str("");
01062       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
01063       hPullGP_Z_rs[title.str()]->Write();
01064 
01065       if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
01066         //mono
01067         gp_ts->cd();
01068         gp_tsx_mono->cd();
01069         title.str("");
01070         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
01071         hPullGP_X_ts_mono[title.str()]->Write();
01072         gp_tsy_mono->cd();
01073         title.str("");
01074         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
01075         hPullGP_Y_ts_mono[title.str()]->Write();
01076         gp_tsz_mono->cd();
01077         title.str("");
01078         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
01079         hPullGP_Z_ts_mono[title.str()]->Write();
01080 
01081         gm_ts->cd();
01082         gm_tsx_mono->cd();
01083         title.str("");
01084         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
01085         hPullGM_X_ts_mono[title.str()]->Write();
01086         gm_tsy_mono->cd();
01087         title.str("");
01088         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
01089         hPullGM_Y_ts_mono[title.str()]->Write();
01090         gm_tsz_mono->cd();
01091         title.str("");
01092         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
01093         hPullGM_Z_ts_mono[title.str()]->Write();
01094 
01095         gp_tr->cd();
01096         gp_trx_mono->cd();
01097         title.str("");
01098         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
01099         hPullGP_X_tr_mono[title.str()]->Write();
01100         gp_try_mono->cd();
01101         title.str("");
01102         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
01103         hPullGP_Y_tr_mono[title.str()]->Write();
01104         gp_trz_mono->cd();
01105         title.str("");
01106         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
01107         hPullGP_Z_tr_mono[title.str()]->Write();
01108 
01109         gp_rs->cd();
01110         gp_rsx_mono->cd();
01111         title.str("");
01112         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
01113         hPullGP_X_rs_mono[title.str()]->Write();
01114         gp_rsy_mono->cd();
01115         title.str("");
01116         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
01117         hPullGP_Y_rs_mono[title.str()]->Write();
01118         gp_rsz_mono->cd();
01119         title.str("");
01120         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
01121         hPullGP_Z_rs_mono[title.str()]->Write();
01122 
01123         //stereo
01124         gp_ts->cd();
01125         gp_tsx_stereo->cd();
01126         title.str("");
01127         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
01128         hPullGP_X_ts_stereo[title.str()]->Write();
01129         gp_tsy_stereo->cd();
01130         title.str("");
01131         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
01132         hPullGP_Y_ts_stereo[title.str()]->Write();
01133         gp_tsz_stereo->cd();
01134         title.str("");
01135         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
01136         hPullGP_Z_ts_stereo[title.str()]->Write();
01137 
01138         gm_ts->cd();
01139         gm_tsx_stereo->cd();
01140         title.str("");
01141         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
01142         hPullGM_X_ts_stereo[title.str()]->Write();
01143         gm_tsy_stereo->cd();
01144         title.str("");
01145         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
01146         hPullGM_Y_ts_stereo[title.str()]->Write();
01147         gm_tsz_stereo->cd();
01148         title.str("");
01149         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
01150         hPullGM_Z_ts_stereo[title.str()]->Write();
01151 
01152         gp_tr->cd();
01153         gp_trx_stereo->cd();
01154         title.str("");
01155         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
01156         hPullGP_X_tr_stereo[title.str()]->Write();
01157         gp_try_stereo->cd();
01158         title.str("");
01159         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
01160         hPullGP_Y_tr_stereo[title.str()]->Write();
01161         gp_trz_stereo->cd();
01162         title.str("");
01163         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
01164         hPullGP_Z_tr_stereo[title.str()]->Write();
01165 
01166         gp_rs->cd();
01167         gp_rsx_stereo->cd();
01168         title.str("");
01169         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
01170         hPullGP_X_rs_stereo[title.str()]->Write();
01171         gp_rsy_stereo->cd();
01172         title.str("");
01173         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
01174         hPullGP_Y_rs_stereo[title.str()]->Write();
01175         gp_rsz_stereo->cd();
01176         title.str("");
01177         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
01178         hPullGP_Z_rs_stereo[title.str()]->Write();
01179       }
01180     }
01181 
01182   file->Close();
01183 }
01184 
01185 //needed by to do the residual for matched hits
01186 //taken from SiStripTrackingRecHitsValid.cc
01187 std::pair<LocalPoint,LocalVector> 
01188 TestHits::projectHit( const PSimHit& hit, const StripGeomDetUnit* stripDet, const BoundPlane& plane) 
01189 { 
01190    const StripTopology& topol = stripDet->specificTopology();
01191    GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
01192    LocalPoint localHit = plane.toLocal(globalpos);
01193    //track direction
01194    LocalVector locdir=hit.localDirection();
01195    //rotate track in new frame
01196    
01197    GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
01198    LocalVector dir=plane.toLocal(globaldir);
01199    float scale = -localHit.z() / dir.z();
01200    
01201    LocalPoint projectedPos = localHit + scale*dir;
01202    
01203    float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
01204    
01205    LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame
01206    
01207    LocalVector localStripDir( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
01208    
01209    return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
01210 }
01211 
01212 #include "FWCore/Framework/interface/ModuleFactory.h"
01213 #include "FWCore/Framework/interface/MakerMacros.h"
01214 
01215 DEFINE_FWK_MODULE(TestHits);