CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoTracker/DebugTools/plugins/TestTrackHits.cc

Go to the documentation of this file.
00001 #include "RecoTracker/DebugTools/interface/TestTrackHits.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00004 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00005 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00006 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
00007 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00008 #include <TDirectory.h>
00009 
00010 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
00011 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00012 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00013 
00014 typedef TrajectoryStateOnSurface TSOS;
00015 typedef TransientTrackingRecHit::ConstRecHitPointer CTTRHp;
00016 using namespace std;
00017 using namespace edm;
00018 
00019 TestTrackHits::TestTrackHits(const edm::ParameterSet& iConfig):
00020   conf_(iConfig) {
00021   LogTrace("TestTrackHits") << conf_;
00022   propagatorName = conf_.getParameter<std::string>("Propagator");   
00023   builderName = conf_.getParameter<std::string>("TTRHBuilder");   
00024   srcName = conf_.getParameter<std::string>("src");   
00025   tpName = conf_.getParameter<std::string>("tpname");   
00026   updatorName = conf_.getParameter<std::string>("updator");
00027   out = conf_.getParameter<std::string>("out");
00028 //   ParameterSet cuts = conf_.getParameter<ParameterSet>("RecoTracksCuts"); 
00029 //   selectRecoTracks = RecoTrackSelector(cuts.getParameter<double>("ptMin"),
00030 //                                     cuts.getParameter<double>("minRapidity"),
00031 //                                     cuts.getParameter<double>("maxRapidity"),
00032 //                                     cuts.getParameter<double>("tip"),
00033 //                                     cuts.getParameter<double>("lip"),
00034 //                                     cuts.getParameter<int>("minHit"),
00035 //                                     cuts.getParameter<double>("maxChi2"));
00036 }
00037 
00038 TestTrackHits::~TestTrackHits(){}
00039 
00040 void TestTrackHits::beginRun(edm::Run & run, const edm::EventSetup& iSetup)
00041 {
00042   iSetup.get<TrackerDigiGeometryRecord>().get(theG);
00043   iSetup.get<IdealMagneticFieldRecord>().get(theMF);  
00044   iSetup.get<TrackingComponentsRecord>().get(propagatorName,thePropagator);
00045   iSetup.get<TransientRecHitRecord>().get(builderName,theBuilder);
00046   iSetup.get<TrackingComponentsRecord>().get(updatorName,theUpdator);
00047   iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits",trackAssociator);
00048 
00049   file = new TFile(out.c_str(),"recreate");
00050   for (int i=0; i!=6; i++)
00051     for (int j=0; j!=9; j++){
00052       if (i==0 && j>2) break;
00053       if (i==1 && j>1) break;
00054       if (i==2 && j>3) break;
00055       if (i==3 && j>2) break;
00056       if (i==4 && j>5) break;
00057       if (i==5 && j>8) break;
00058 
00059       title.str("");
00060       title << "Chi2Increment_" << i+1 << "-" << j+1 ;
00061       hChi2Increment[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00062       title.str("");
00063       title << "Chi2IncrementVsEta_" << i+1 << "-" << j+1 ;
00064       hChi2IncrementVsEta[title.str()] = new TH2F(title.str().c_str(),title.str().c_str(),50,-2.5,2.5,1000,0,100);
00065       title.str("");
00066       title << "Chi2GoodHit_" << i+1 << "-" << j+1 ;
00067       hChi2GoodHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00068       title.str("");
00069       title << "Chi2BadHit_" << i+1 << "-" << j+1 ;
00070       hChi2BadHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00071       title.str("");
00072       title << "Chi2DeltaHit_" << i+1 << "-" << j+1 ;
00073       hChi2DeltaHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00074       title.str("");
00075       title << "Chi2NSharedHit_" << i+1 << "-" << j+1 ;
00076       hChi2NSharedHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00077       title.str("");
00078       title << "Chi2SharedHit_" << i+1 << "-" << j+1 ;
00079       hChi2SharedHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00080 
00081       title.str("");
00082       title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
00083       hPullGP_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00084       title.str("");
00085       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
00086       hPullGP_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00087       title.str("");
00088       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
00089       hPullGP_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00090 
00091       title.str("");
00092       title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
00093       hPullGM_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00094       title.str("");
00095       title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
00096       hPullGM_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00097       title.str("");
00098       title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
00099       hPullGM_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00100 
00101       title.str("");
00102       title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
00103       hPullGP_X_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00104       title.str("");
00105       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
00106       hPullGP_Y_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00107       title.str("");
00108       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
00109       hPullGP_Z_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00110 
00111       title.str("");
00112       title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
00113       hPullGP_X_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00114       title.str("");
00115       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
00116       hPullGP_Y_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00117       title.str("");
00118       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
00119       hPullGP_Z_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00120 
00121       if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
00122         //mono
00123         title.str("");
00124         title << "Chi2Increment_mono_" << i+1 << "-" << j+1 ;
00125         hChi2Increment_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00126 
00127         title.str("");
00128         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
00129         hPullGP_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00130         title.str("");
00131         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00132         hPullGP_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00133         title.str("");
00134         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00135         hPullGP_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00136 
00137         title.str("");
00138         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
00139         hPullGM_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00140         title.str("");
00141         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00142         hPullGM_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00143         title.str("");
00144         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00145         hPullGM_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00146 
00147         title.str("");
00148         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
00149         hPullGP_X_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00150         title.str("");
00151         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
00152         hPullGP_Y_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00153         title.str("");
00154         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
00155         hPullGP_Z_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00156 
00157         title.str("");
00158         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
00159         hPullGP_X_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00160         title.str("");
00161         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
00162         hPullGP_Y_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00163         title.str("");
00164         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
00165         hPullGP_Z_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00166 
00167         //stereo
00168         title.str("");
00169         title << "Chi2Increment_stereo_" << i+1 << "-" << j+1 ;
00170         hChi2Increment_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00171 
00172         title.str("");
00173         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00174         hPullGP_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00175         title.str("");
00176         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00177         hPullGP_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00178         title.str("");
00179         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00180         hPullGP_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00181 
00182         title.str("");
00183         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00184         hPullGM_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00185         title.str("");
00186         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00187         hPullGM_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00188         title.str("");
00189         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00190         hPullGM_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00191 
00192         title.str("");
00193         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
00194         hPullGP_X_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00195         title.str("");
00196         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
00197         hPullGP_Y_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00198         title.str("");
00199         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
00200         hPullGP_Z_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00201 
00202         title.str("");
00203         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
00204         hPullGP_X_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00205         title.str("");
00206         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
00207         hPullGP_Y_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00208         title.str("");
00209         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
00210         hPullGP_Z_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00211       }
00212     }
00213   hTotChi2Increment = new TH1F("TotChi2Increment","TotChi2Increment",1000,0,100);
00214   hTotChi2GoodHit = new TH1F("TotChi2GoodHit","TotChi2GoodHit",1000,0,100);
00215   hTotChi2BadHit = new TH1F("TotChi2BadHit","TotChi2BadHit",1000,0,100);
00216   hTotChi2DeltaHit = new TH1F("TotChi2DeltaHit","TotChi2DeltaHit",1000,0,100);
00217   hTotChi2NSharedHit = new TH1F("TotChi2NSharedHit","TotChi2NSharedHit",1000,0,100);
00218   hTotChi2SharedHit = new TH1F("TotChi2SharedHit","TotChi2SharedHit",1000,0,100);
00219   hProcess_vs_Chi2  = new TH2F("Process_vs_Chi2","Process_vs_Chi2",1000,0,100,17,-0.5,16.5);  
00220   hClsize_vs_Chi2   = new TH2F("Clsize_vs_Chi2","Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00221   hPixClsize_vs_Chi2= new TH2F("PixClsize_vs_Chi2","PixClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00222   hPrjClsize_vs_Chi2= new TH2F("PrjClsize_vs_Chi2","PrjClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00223   hSt1Clsize_vs_Chi2= new TH2F("St1Clsize_vs_Chi2","St1Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00224   hSt2Clsize_vs_Chi2= new TH2F("St2Clsize_vs_Chi2","St2Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00225   hGoodHit_vs_Chi2  = new TH2F("GoodHit_vs_Chi2","GoodHit_vs_Chi2",10000,0,1000,2,-0.5,1.5);
00226   hClusterSize = new TH1F("ClusterSize","ClusterSize",40,-0.5,39.5);
00227   hPixClusterSize = new TH1F("PixClusterSize","PixClusterSize",40,-0.5,39.5);
00228   hPrjClusterSize = new TH1F("PrjClusterSize","PrjClusterSize",40,-0.5,39.5);
00229   hSt1ClusterSize = new TH1F("St1ClusterSize","St1ClusterSize",40,-0.5,39.5);
00230   hSt2ClusterSize = new TH1F("St2ClusterSize","St2ClusterSize",40,-0.5,39.5);
00231   hSimHitVecSize = new TH1F("hSimHitVecSize","hSimHitVecSize",40,-0.5,39.5);
00232   hPixSimHitVecSize = new TH1F("PixSimHitVecSize","PixSimHitVecSize",40,-0.5,39.5);
00233   hPrjSimHitVecSize = new TH1F("PrjSimHitVecSize","PrjSimHitVecSize",40,-0.5,39.5);
00234   hSt1SimHitVecSize = new TH1F("St1SimHitVecSize","St1SimHitVecSize",40,-0.5,39.5);
00235   hSt2SimHitVecSize = new TH1F("St2SimHitVecSize","St2SimHitVecSize",40,-0.5,39.5);
00236   goodbadmerged = new TH1F("goodbadmerged","goodbadmerged",5,0.5,5.5);
00237   energyLossRatio = new TH1F("energyLossRatio","energyLossRatio",100,0,1);
00238   mergedPull = new TH1F("mergedPull","mergedPull",200,0,20);
00239   probXgood = new TH1F("probXgood","probXgood",110,0,1.1);
00240   probXbad = new TH1F("probXbad","probXbad",110,0,1.1);
00241   probXdelta = new TH1F("probXdelta","probXdelta",110,0,1.1);
00242   probXshared = new TH1F("probXshared","probXshared",110,0,1.1);
00243   probXnoshare = new TH1F("probXnoshare","probXnoshare",110,0,1.1);
00244   probYgood = new TH1F("probYgood","probYgood",110,0,1.1);
00245   probYbad = new TH1F("probYbad","probYbad",110,0,1.1);
00246   probYdelta = new TH1F("probYdelta","probYdelta",110,0,1.1);
00247   probYshared = new TH1F("probYshared","probYshared",110,0,1.1);
00248   probYnoshare = new TH1F("probYnoshare","probYnoshare",110,0,1.1);
00249 }
00250 
00251 void TestTrackHits::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00252 {
00253   //Retrieve tracker topology from geometry
00254   edm::ESHandle<TrackerTopology> tTopo;
00255   iSetup.get<IdealGeometryRecord>().get(tTopo);
00256 
00257 
00258   LogDebug("TestTrackHits") << "new event" ;
00259   iEvent.getByLabel(srcName,trajCollectionHandle);
00260   iEvent.getByLabel(srcName,trackCollectionHandle);
00261   iEvent.getByLabel(srcName,trajTrackAssociationCollectionHandle);
00262   iEvent.getByLabel(tpName,trackingParticleCollectionHandle);
00263   edm::Handle<reco::BeamSpot> beamSpot;
00264   iEvent.getByLabel("offlineBeamSpot",beamSpot); 
00265 
00266   LogTrace("TestTrackHits") << "Tr collection size=" << trackCollectionHandle->size();
00267   LogTrace("TestTrackHits") << "TP collection size=" << trackingParticleCollectionHandle->size();
00268 
00269   hitAssociator = new TrackerHitAssociator(iEvent);
00270   
00271   reco::RecoToSimCollection recSimColl=trackAssociator->associateRecoToSim(trackCollectionHandle,
00272                                                                            trackingParticleCollectionHandle,
00273                                                                            &iEvent,&iSetup);
00274   
00275   TrajectoryStateCombiner combiner;
00276 
00277   int evtHits = 0;
00278   int i=0;
00279   int yy=0;
00280   int yyy=0;
00281   for(std::vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it!=trajCollectionHandle->end();it++){
00282     
00283     LogTrace("TestTrackHits") << "\n*****************new trajectory********************" ;
00284     double tchi2 = 0;
00285 
00286     std::vector<TrajectoryMeasurement> tmColl = it->measurements();
00287 
00288     edm::Ref<std::vector<Trajectory> >  traj(trajCollectionHandle, i);
00289     reco::TrackRef tmptrack = (*trajTrackAssociationCollectionHandle.product())[traj];
00290     edm::RefToBase<reco::Track> track(tmptrack);
00291     //     if ( !selectRecoTracks( *track,beamSpot.product() ) ) {
00292     //       LogTrace("TestTrackHits") << "track does not pass quality cuts: skippingtrack #" << ++yyy;
00293     //       i++;
00294     //       continue;
00295     //     }
00296 
00297     std::vector<std::pair<TrackingParticleRef, double> > tP;
00298     if(recSimColl.find(track) != recSimColl.end()){
00299       tP = recSimColl[track];
00300       if (tP.size()!=0) {
00301         edm::LogVerbatim("TestTrackHits") << "reco::Track #" << ++yyy << " with pt=" << track->pt() 
00302                                           << " associated with quality:" << tP.begin()->second <<" good track #" << ++yy << " has hits:" << track->numberOfValidHits() << "\n";
00303       }
00304     }else{
00305       edm::LogVerbatim("TestTrackHits") << "reco::Track #" << ++yyy << " with pt=" << track->pt()
00306                                          << " NOT associated to any TrackingParticle" << "\n";
00307       i++;
00308       continue;
00309     }
00310     //     if(recSimColl.find(track) != recSimColl.end()) {
00311     //       tP = recSimColl[track];
00312     //     } else {
00313     //       LogTrace("TestTrackHits") << "fake track: skipping track " << ++yyy;
00314     //       continue;//skip fake tracks
00315     //     }
00316     //     if (tP.size()==0) {
00317     //       LogTrace("TestTrackHits") << "fake track: skipping track " << ++yyy;
00318     //       continue;
00319     //     }
00320     TrackingParticleRef tp = tP.begin()->first;
00321     LogTrace("TestTrackHits") << "a tp is associated with fraction=" << tP.begin()->second;
00322     //LogTrace("TestTrackHits") << "last tp is associated with fraction=" << (tP.end()-1)->second;
00323     std::vector<unsigned int> tpids;
00324     for (TrackingParticle::g4t_iterator g4T=tp->g4Track_begin(); g4T!=tp->g4Track_end(); ++g4T) {
00325       LogTrace("TestTrackHits") << "tp id=" << g4T->trackId();
00326       tpids.push_back(g4T->trackId());
00327     }
00328 
00329     //LogTrace("TestTrackHits") << "Analyzing hits of track number " << ++yyy << " good track number " << ++yy;
00330     int pp = 0;
00331     for (std::vector<TrajectoryMeasurement>::iterator tm=tmColl.begin();tm!=tmColl.end();++tm){
00332       
00333       tchi2+=tm->estimate();
00334 
00335       LogTrace("TestTrackHits") << "+++++++++++++++++new hit+++++++++++++++++" ;
00336       CTTRHp rhit = tm->recHit();
00337       //TSOS state = tm->backwardPredictedState();
00338       //TSOS state = tm->forwardPredictedState();
00339       TSOS state = combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
00340 
00341       if (rhit->isValid()==0 && rhit->det()!=0) continue;
00342       evtHits++;
00343       LogTrace("TestTrackHits") << "valid hit #" << ++pp << "of hits=" << track->numberOfValidHits();
00344 
00345       int subdetId = rhit->det()->geographicalId().subdetId();
00346       DetId id = rhit->det()->geographicalId();
00347       int layerId  = tTopo->layer(id);
00348       LogTrace("TestTrackHits") << "subdetId=" << subdetId << " layerId=" << layerId ;
00349 
00350       const Surface * surf = rhit->surface();
00351       if (surf==0) continue;
00352 
00353       double energyLoss_ = 0.;
00354       unsigned int monoId = 0;
00355       std::vector<double> energyLossM;
00356       std::vector<double> energyLossS;
00357       std::vector<PSimHit> assSimHits = hitAssociator->associateHit(*(rhit)->hit());
00358       unsigned int  simhitvecsize = assSimHits.size();
00359       if (simhitvecsize==0) continue;
00360       PSimHit shit;
00361       std::vector<unsigned int> trackIds;
00362       energyLossS.clear();
00363       energyLossM.clear();
00364       for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
00365         unsigned int tId = m->trackId();
00366         if (find(trackIds.begin(),trackIds.end(),tId)==trackIds.end()) trackIds.push_back(tId);
00367         if (m->energyLoss()>energyLoss_) {
00368           shit=*m;
00369           energyLoss_ = m->energyLoss();
00370         }
00371         if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
00372           if (monoId==0) monoId = m->detUnitId();
00373           if (monoId == m->detUnitId()){
00374             energyLossM.push_back(m->energyLoss());
00375           }
00376           else {
00377             energyLossS.push_back(m->energyLoss());
00378           }
00379           //std::cout << "detUnitId="  << m->detUnitId() << " trackId=" << m->trackId() << " energyLoss=" << m->energyLoss() << std::endl;
00380         } else {
00381           energyLossM.push_back(m->energyLoss());
00382         }
00383       }
00384       //double delta = 99999;
00385       //LocalPoint rhitLPv = rhit->localPosition();
00386       //vector<PSimHit> assSimHits = hitAssociator->associateHit(*(rhit)->hit());
00387       //unsigned int  simhitvecsize = assSimHits.size();
00388       //if (simhitvecsize==0) continue;
00389       //PSimHit shit;
00390       //for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
00391       //  if ((m->localPosition()-rhitLPv).mag()<delta) {
00392       //    shit=*m;
00393       //    delta = (m->localPosition()-rhitLPv).mag();
00394       //  }
00395       //}
00396 
00397       //plot chi2 increment
00398       double chi2increment = tm->estimate();
00399       LogTrace("TestTrackHits") << "tm->estimate()=" << tm->estimate();
00400       title.str("");
00401       title << "Chi2Increment_" << subdetId << "-" << layerId;
00402       hChi2Increment[title.str()]->Fill( chi2increment );
00403       title.str("");
00404       title << "Chi2IncrementVsEta_" << subdetId << "-" << layerId;
00405       hChi2IncrementVsEta[title.str()]->Fill( track->eta(), chi2increment );
00406       hTotChi2Increment->Fill( chi2increment );
00407       hProcess_vs_Chi2->Fill( chi2increment, shit.processType() );
00408 
00409       int clustersize = 0;
00410       bool mergedhit = false;
00411       if (dynamic_cast<const SiPixelRecHit*>(rhit->hit())){     
00412         clustersize =  ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size() ;
00413         hPixClsize_vs_Chi2->Fill(chi2increment, clustersize);
00414         hPixClusterSize->Fill(clustersize);
00415         hPixSimHitVecSize->Fill(simhitvecsize);
00416         if (simhitvecsize>1) mergedhit = true;
00417       }
00418       else if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit())){
00419         clustersize =  ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() ;
00420         hSt1Clsize_vs_Chi2->Fill(chi2increment, clustersize);
00421         hSt1ClusterSize->Fill(clustersize);
00422         hSt1SimHitVecSize->Fill(simhitvecsize);
00423         if (simhitvecsize>1) mergedhit = true;
00424       }
00425       else if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())){
00426         int clsize1 = ((const SiStripMatchedRecHit2D*)(rhit->hit()))->monoCluster().amplitudes().size();
00427         int clsize2 =  ((const SiStripMatchedRecHit2D*)(rhit->hit()))->stereoCluster().amplitudes().size();
00428         if (clsize1>clsize2) clustersize = clsize1;
00429         else clustersize = clsize2;
00430         hSt2Clsize_vs_Chi2->Fill(chi2increment, clustersize);
00431         hSt2ClusterSize->Fill(clustersize);
00432         hSt2SimHitVecSize->Fill(simhitvecsize);
00433         if (simhitvecsize>2) mergedhit = true;
00434       }
00435       else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(rhit->hit())){
00436         clustersize =  ((const ProjectedSiStripRecHit2D*)(rhit->hit()))->originalHit().cluster()->amplitudes().size();
00437         hPrjClsize_vs_Chi2->Fill(chi2increment, clustersize);
00438         hPrjClusterSize->Fill(clustersize);
00439         hPrjSimHitVecSize->Fill(simhitvecsize);
00440         if (simhitvecsize>1) mergedhit = true;
00441       }
00442       hClsize_vs_Chi2->Fill( chi2increment, clustersize);
00443       hClusterSize->Fill(clustersize);
00444       hSimHitVecSize->Fill(simhitvecsize);
00445 
00446       //       if (dynamic_cast<const SiPixelRecHit*>(rhit->hit()))     
00447       //        hClsize_vs_Chi2->Fill( chi2increment, ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size() );
00448       //       if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit()))   
00449       //        hClsize_vs_Chi2->Fill( chi2increment, ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() );
00450 
00451       std::vector<SimHitIdpr> simTrackIds = hitAssociator->associateHitId(*(rhit)->hit());
00452       bool goodhit = false;
00453       for(size_t j=0; j<simTrackIds.size(); j++){
00454         LogTrace("TestTrackHits") << "hit id=" << simTrackIds[j].first;
00455         for (size_t jj=0; jj<tpids.size(); jj++){
00456           if (simTrackIds[j].first == tpids[jj]) goodhit = true;
00457           break;
00458         }
00459       }
00460       bool shared = true;
00461       bool ioniOnly = true;
00462       const SiPixelRecHit* pix = dynamic_cast<const SiPixelRecHit*>(rhit->hit());
00463       if (goodhit) {
00464         if (energyLossM.size()>1&&energyLossS.size()<=1) {
00465           sort(energyLossM.begin(),energyLossM.end(),greater<double>());
00466           energyLossRatio->Fill(energyLossM[1]/energyLossM[0]);
00467         }
00468         else if (energyLossS.size()>1&&energyLossM.size()<=1) {
00469           sort(energyLossS.begin(),energyLossS.end(),greater<double>());
00470           energyLossRatio->Fill(energyLossS[1]/energyLossS[0]);
00471         }
00472         else if (energyLossS.size()>1&&energyLossM.size()>1) {
00473           sort(energyLossM.begin(),energyLossM.end(),greater<double>());
00474           sort(energyLossS.begin(),energyLossS.end(),greater<double>());
00475           energyLossM[1]/energyLossM[0] > energyLossS[1]/energyLossS[0] 
00476             ? energyLossRatio->Fill(energyLossM[1]/energyLossM[0]) 
00477             : energyLossRatio->Fill(energyLossS[1]/energyLossS[0]);
00478         }
00479         if ( mergedhit ) {
00480           //not optimized for matched hits
00481           LogVerbatim("TestTrackHits") << "MERGED HIT" << std::endl;
00482           unsigned int idc = 0;
00483           for (size_t jj=0; jj<tpids.size(); jj++){
00484             idc += std::count(trackIds.begin(),trackIds.end(),tpids[jj]);
00485           }
00486           if (idc==trackIds.size()) {
00487             shared = false;
00488           }
00489           for(std::vector<PSimHit>::const_iterator m=assSimHits.begin()+1; m<assSimHits.end(); m++){
00490             if ((m->processType()!=7&&m->processType()!=8&&m->processType()!=9)&&abs(m->particleType())!=11){
00491               ioniOnly = false;
00492               break;
00493             }
00494           }
00495           if (ioniOnly&&!shared) {
00496             title.str("");
00497             title << "Chi2DeltaHit_" << subdetId << "-" << layerId;
00498             hChi2DeltaHit[title.str()]->Fill( chi2increment );
00499             hTotChi2DeltaHit->Fill( chi2increment );      
00500             if (pix) {
00501               probXdelta->Fill(pix->probabilityX());
00502               probYdelta->Fill(pix->probabilityY());
00503             }       
00504           } else if(!ioniOnly&&!shared) {
00505             title.str("");
00506             title << "Chi2NSharedHit_" << subdetId << "-" << layerId;
00507             hChi2NSharedHit[title.str()]->Fill( chi2increment );
00508             hTotChi2NSharedHit->Fill( chi2increment );    
00509             if (pix) {
00510               probXnoshare->Fill(pix->probabilityX());
00511               probYnoshare->Fill(pix->probabilityY());
00512             }       
00513           } else {
00514             title.str("");
00515             title << "Chi2SharedHit_" << subdetId << "-" << layerId;
00516             hChi2SharedHit[title.str()]->Fill( chi2increment );
00517             hTotChi2SharedHit->Fill( chi2increment );     
00518             if (pix) {
00519               probXshared->Fill(pix->probabilityX());
00520               probYshared->Fill(pix->probabilityY());
00521             }
00522           }
00523           
00524           for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++) {
00525             unsigned int tId = m->trackId();
00526             LogVerbatim("TestTrackHits") << "component with id=" << tId <<  " eLoss=" << m->energyLoss() << " pType=" <<  m->processType();
00527             if (find(tpids.begin(),tpids.end(),tId)==tpids.end()) continue;
00528             if (m->processType()==2) {
00529               GlobalPoint gpr = rhit->globalPosition();
00530               AlgebraicSymMatrix33 ger = rhit->globalPositionError().matrix();
00531               GlobalPoint gps = surf->toGlobal(m->localPosition());
00532               LogVerbatim("TestTrackHits") << gpr << " " << gps << " " << ger;
00533               ROOT::Math::SVector<double,3> delta;
00534               delta[0]=gpr.x()-gps.x();
00535               delta[1]=gpr.y()-gps.y();
00536               delta[2]=gpr.z()-gps.z();
00537               LogVerbatim("TestTrackHits") << delta << " " << ger ;
00538               double mpull = sqrt(delta[0]*delta[0]/ger[0][0]+delta[1]*delta[1]/ger[1][1]+delta[2]*delta[2]/ger[2][2]);
00539               LogVerbatim("TestTrackHits") << "hit pull=" << mpull;//ger.similarity(delta);
00540               mergedPull->Fill(mpull);
00541               break;
00542             }
00543           }
00544         } else {
00545           LogVerbatim("TestTrackHits") << "good hit" ;
00546           title.str("");
00547           title << "Chi2GoodHit_" << subdetId << "-" << layerId;
00548           hChi2GoodHit[title.str()]->Fill( chi2increment );
00549           hTotChi2GoodHit->Fill( chi2increment );         
00550           if (pix) {
00551             probXgood->Fill(pix->probabilityX());
00552             probYgood->Fill(pix->probabilityY());
00553           }
00554         }
00555       } else {
00556         LogVerbatim("TestTrackHits") << "bad hit" ;
00557         title.str("");
00558         title << "Chi2BadHit_" << subdetId << "-" << layerId;
00559         hChi2BadHit[title.str()]->Fill( chi2increment );
00560         hTotChi2BadHit->Fill( chi2increment );
00561         goodbadmerged->Fill(2);
00562         if (pix) {
00563           probXbad->Fill(pix->probabilityX());
00564           probYbad->Fill(pix->probabilityY());
00565         }
00566       }
00567       hGoodHit_vs_Chi2->Fill(chi2increment,goodhit);
00568 
00569       LocalVector shitLMom;
00570       LocalPoint shitLPos;
00571       if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
00572         if (simhitvecsize>2 && goodhit) { 
00573           if (ioniOnly&&!shared) goodbadmerged->Fill(3);
00574           else if(!ioniOnly&&!shared) goodbadmerged->Fill(4);
00575           else goodbadmerged->Fill(5);
00576         }
00577         else if (goodhit) goodbadmerged->Fill(1);
00578         double rechitmatchedx = rhit->localPosition().x();
00579         double rechitmatchedy = rhit->localPosition().y();
00580         double mindist = 999999;
00581         float distx, disty;
00582         std::pair<LocalPoint,LocalVector> closestPair;
00583         const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) ((const GluedGeomDet *)(rhit)->det())->stereoDet();
00584         const BoundPlane& plane = (rhit)->det()->surface();
00585         for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
00586           //project simhit;
00587           std::pair<LocalPoint,LocalVector> hitPair = projectHit((*m),stripDet,plane);
00588           distx = fabs(rechitmatchedx - hitPair.first.x());
00589           disty = fabs(rechitmatchedy - hitPair.first.y());
00590           double dist = distx*distx+disty*disty;
00591           if(sqrt(dist)<mindist){
00592             mindist = dist;
00593             closestPair = hitPair;
00594           }
00595         }
00596         shitLPos = closestPair.first;
00597         shitLMom = closestPair.second;
00598       } else {
00599         if (simhitvecsize>1 && goodhit) { 
00600           if (ioniOnly&&!shared) goodbadmerged->Fill(3);
00601           else if(!ioniOnly&&!shared) goodbadmerged->Fill(4);
00602           else goodbadmerged->Fill(5);
00603         }
00604         else if (goodhit) goodbadmerged->Fill(1);
00605         shitLPos = shit.localPosition();        
00606         shitLMom = shit.momentumAtEntry();
00607       }
00608       GlobalVector shitGMom = surf->toGlobal(shitLMom);
00609       GlobalPoint shitGPos = surf->toGlobal(shitLPos);
00610 
00611       GlobalVector tsosGMom = state.globalMomentum();
00612       GlobalError  tsosGMEr(state.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00613       GlobalPoint  tsosGPos = state.globalPosition();
00614       GlobalError  tsosGPEr = state.cartesianError().position();
00615 
00616       GlobalPoint rhitGPos = (rhit)->globalPosition();
00617       GlobalError rhitGPEr = (rhit)->globalPositionError();
00618 
00619       LogVerbatim("TestTrackHits") << "assSimHits.size()=" << assSimHits.size() ;
00620       LogVerbatim("TestTrackHits") << "tsos globalPos   =" << tsosGPos ;
00621       LogVerbatim("TestTrackHits") << "sim hit globalPos=" << shitGPos ;
00622       LogVerbatim("TestTrackHits") << "rec hit globalPos=" << rhitGPos ;
00623       LogVerbatim("TestTrackHits") << "geographicalId   =" << rhit->det()->geographicalId().rawId() ;
00624       LogVerbatim("TestTrackHits") << "surface position =" << surf->position() ;
00625 
00626 # if 0
00627       if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->geographicalId()=" 
00628                                                      << rhit->detUnit()->geographicalId().rawId() ;
00629       LogTrace("TestTrackHits") << "rhit->det()->surface().position()=" 
00630                                 << rhit->det()->surface().position() ;
00631       if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().position()="  
00632                                                      << rhit->detUnit()->surface().position() ;
00633       LogTrace("TestTrackHits") << "rhit->det()->position()=" << rhit->det()->position() ;
00634       if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->position()="  
00635                                                      << rhit->detUnit()->position() ;
00636       LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().length()=" 
00637                                 << rhit->det()->surface().bounds().length() ;
00638       if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().length()="  
00639                                                      << rhit->detUnit()->surface().bounds().length() ;
00640       LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().width()=" 
00641                                 << rhit->det()->surface().bounds().width() ;
00642       if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().width()="  
00643                                                      << rhit->detUnit()->surface().bounds().width() ;
00644       LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().thickness()=" 
00645                                 << rhit->det()->surface().bounds().thickness() ;
00646       if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().thickness()="  
00647                                                      << rhit->detUnit()->surface().bounds().thickness() ;
00648 #endif      
00649 
00650       double pullGPX_rs = (rhitGPos.x()-shitGPos.x())/sqrt(rhitGPEr.cxx());
00651       double pullGPY_rs = (rhitGPos.y()-shitGPos.y())/sqrt(rhitGPEr.cyy());
00652       double pullGPZ_rs = (rhitGPos.z()-shitGPos.z())/sqrt(rhitGPEr.czz());
00653       //double pullGPX_rs = (rhitGPos.x()-shitGPos.x());
00654       //double pullGPY_rs = (rhitGPos.y()-shitGPos.y());
00655       //double pullGPZ_rs = (rhitGPos.z()-shitGPos.z());
00656 
00657       LogTrace("TestTrackHits") << "rs" ;
00658 
00659       title.str("");
00660       title << "PullGP_X_" << subdetId << "-" << layerId << "_rs";
00661       hPullGP_X_rs[title.str()]->Fill( pullGPX_rs );
00662       title.str("");
00663       title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs";
00664       hPullGP_Y_rs[title.str()]->Fill( pullGPY_rs );
00665       title.str("");
00666       title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs";
00667       hPullGP_Z_rs[title.str()]->Fill( pullGPZ_rs );
00668 
00669       double pullGPX_tr = (tsosGPos.x()-rhitGPos.x())/sqrt(tsosGPEr.cxx()+rhitGPEr.cxx());
00670       double pullGPY_tr = (tsosGPos.y()-rhitGPos.y())/sqrt(tsosGPEr.cyy()+rhitGPEr.cyy());
00671       double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z())/sqrt(tsosGPEr.czz()+rhitGPEr.czz());
00672       //double pullGPX_tr = (tsosGPos.x()-rhitGPos.x());
00673       //double pullGPY_tr = (tsosGPos.y()-rhitGPos.y());
00674       //double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z());
00675 
00676       LogTrace("TestTrackHits") << "tr" ;
00677 
00678       title.str("");
00679       title << "PullGP_X_" << subdetId << "-" << layerId << "_tr";
00680       hPullGP_X_tr[title.str()]->Fill( pullGPX_tr );
00681       title.str("");
00682       title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr";
00683       hPullGP_Y_tr[title.str()]->Fill( pullGPY_tr );
00684       title.str("");
00685       title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr";
00686       hPullGP_Z_tr[title.str()]->Fill( pullGPZ_tr );
00687 
00688       double pullGPX_ts = (tsosGPos.x()-shitGPos.x())/sqrt(tsosGPEr.cxx());
00689       double pullGPY_ts = (tsosGPos.y()-shitGPos.y())/sqrt(tsosGPEr.cyy());
00690       double pullGPZ_ts = (tsosGPos.z()-shitGPos.z())/sqrt(tsosGPEr.czz());
00691       //double pullGPX_ts = (tsosGPos.x()-shitGPos.x());
00692       //double pullGPY_ts = (tsosGPos.y()-shitGPos.y());
00693       //double pullGPZ_ts = (tsosGPos.z()-shitGPos.z());
00694 
00695       LogTrace("TestTrackHits") << "ts1" ;
00696 
00697       title.str("");
00698       title << "PullGP_X_" << subdetId << "-" << layerId << "_ts";
00699       hPullGP_X_ts[title.str()]->Fill( pullGPX_ts );
00700       title.str("");
00701       title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts";
00702       hPullGP_Y_ts[title.str()]->Fill( pullGPY_ts );
00703       title.str("");
00704       title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts";
00705       hPullGP_Z_ts[title.str()]->Fill( pullGPZ_ts );
00706 
00707       double pullGMX_ts = (tsosGMom.x()-shitGMom.x())/sqrt(tsosGMEr.cxx());
00708       double pullGMY_ts = (tsosGMom.y()-shitGMom.y())/sqrt(tsosGMEr.cyy());
00709       double pullGMZ_ts = (tsosGMom.z()-shitGMom.z())/sqrt(tsosGMEr.czz());
00710       //double pullGMX_ts = (tsosGMom.x()-shitGMom.x());
00711       //double pullGMY_ts = (tsosGMom.y()-shitGMom.y());
00712       //double pullGMZ_ts = (tsosGMom.z()-shitGMom.z());
00713 
00714       LogTrace("TestTrackHits") << "ts2" ;
00715 
00716       title.str("");
00717       title << "PullGM_X_" << subdetId << "-" << layerId << "_ts";
00718       hPullGM_X_ts[title.str()]->Fill( pullGMX_ts );
00719       title.str("");
00720       title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts";
00721       hPullGM_Y_ts[title.str()]->Fill( pullGMY_ts );
00722       title.str("");
00723       title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts";
00724       hPullGM_Z_ts[title.str()]->Fill( pullGMZ_ts );
00725       if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
00726         Propagator* thePropagatorAnyDir = new PropagatorWithMaterial(anyDirection,0.105,theMF.product(),1.6);
00727         //mono
00728         LogTrace("TestTrackHits") << "MONO HIT" ;
00729         auto m = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->monoHit();
00730         CTTRHp tMonoHit = theBuilder->build(&m);
00731         if (tMonoHit==0) continue;
00732         vector<PSimHit> assMonoSimHits = hitAssociator->associateHit(*tMonoHit->hit());
00733         if (assMonoSimHits.size()==0) continue;
00734         const PSimHit sMonoHit = *(assMonoSimHits.begin());
00735         const Surface * monoSurf = &( tMonoHit->det()->surface() );
00736         if (monoSurf==0) continue;
00737         TSOS monoState = thePropagatorAnyDir->propagate(state,*monoSurf);
00738         if (monoState.isValid()==0) continue;
00739 
00740         LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
00741         GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
00742         LocalPoint monoShitLPos = sMonoHit.localPosition();
00743         GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
00744 
00745         //LogTrace("TestTrackHits") << "assMonoSimHits.size()=" << assMonoSimHits.size() ;
00746         //LogTrace("TestTrackHits") << "mono shit=" << monoShitGPos ;
00747 
00748         GlobalVector monoTsosGMom = monoState.globalMomentum();
00749         GlobalError  monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00750         GlobalPoint  monoTsosGPos = monoState.globalPosition();
00751         GlobalError  monoTsosGPEr = monoState.cartesianError().position();
00752 
00753         GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
00754         GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
00755 
00756         double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x())/sqrt(monoRhitGPEr.cxx());
00757         double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y())/sqrt(monoRhitGPEr.cyy());
00758         double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z())/sqrt(monoRhitGPEr.czz());
00759         //double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x());
00760         //double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y());
00761         //double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z());
00762 
00763         MeasurementExtractor meMo(monoState);
00764         double chi2mono = computeChi2Increment(meMo,tMonoHit);
00765 
00766         title.str("");
00767         title << "Chi2Increment_mono_" << subdetId << "-" << layerId ;
00768         hChi2Increment_mono[title.str()]->Fill(chi2mono);
00769 
00770         title.str("");
00771         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
00772         hPullGP_X_rs_mono[title.str()]->Fill( pullGPX_rs_mono );
00773         title.str("");
00774         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
00775         hPullGP_Y_rs_mono[title.str()]->Fill( pullGPY_rs_mono );
00776         title.str("");
00777         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
00778         hPullGP_Z_rs_mono[title.str()]->Fill( pullGPZ_rs_mono );
00779 
00780         double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x())/sqrt(monoTsosGPEr.cxx()+monoRhitGPEr.cxx());
00781         double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y())/sqrt(monoTsosGPEr.cyy()+monoRhitGPEr.cyy());
00782         double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z())/sqrt(monoTsosGPEr.czz()+monoRhitGPEr.czz());
00783         //double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x());
00784         //double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y());
00785         //double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z());
00786 
00787         title.str("");
00788         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
00789         hPullGP_X_tr_mono[title.str()]->Fill( pullGPX_tr_mono );
00790         title.str("");
00791         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
00792         hPullGP_Y_tr_mono[title.str()]->Fill( pullGPY_tr_mono );
00793         title.str("");
00794         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
00795         hPullGP_Z_tr_mono[title.str()]->Fill( pullGPZ_tr_mono );
00796 
00797         double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x())/sqrt(monoTsosGPEr.cxx());
00798         double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y())/sqrt(monoTsosGPEr.cyy());
00799         double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z())/sqrt(monoTsosGPEr.czz());
00800         //double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x());
00801         //double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y());
00802         //double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z());
00803 
00804         title.str("");
00805         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
00806         hPullGP_X_ts_mono[title.str()]->Fill( pullGPX_ts_mono );
00807         title.str("");
00808         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
00809         hPullGP_Y_ts_mono[title.str()]->Fill( pullGPY_ts_mono );
00810         title.str("");
00811         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
00812         hPullGP_Z_ts_mono[title.str()]->Fill( pullGPZ_ts_mono );
00813 
00814         double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x())/sqrt(monoTsosGMEr.cxx());
00815         double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y())/sqrt(monoTsosGMEr.cyy());
00816         double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z())/sqrt(monoTsosGMEr.czz());
00817         //double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x());
00818         //double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y());
00819         //double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z());
00820 
00821         title.str("");
00822         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
00823         hPullGM_X_ts_mono[title.str()]->Fill( pullGMX_ts_mono );
00824         title.str("");
00825         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
00826         hPullGM_Y_ts_mono[title.str()]->Fill( pullGMY_ts_mono );
00827         title.str("");
00828         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
00829         hPullGM_Z_ts_mono[title.str()]->Fill( pullGMZ_ts_mono );
00830 
00831         //stereo
00832         LogTrace("TestTrackHits") << "STEREO HIT" ;
00833         auto s = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->stereoHit();
00834         CTTRHp tStereoHit = theBuilder->build(&s);
00835         if (tStereoHit==0) continue;
00836         vector<PSimHit> assStereoSimHits = hitAssociator->associateHit(*tStereoHit->hit());
00837         if (assStereoSimHits.size()==0) continue;
00838         const PSimHit sStereoHit = *(assStereoSimHits.begin());
00839         const Surface * stereoSurf = &( tStereoHit->det()->surface() );
00840         if (stereoSurf==0) continue;
00841         TSOS stereoState = thePropagatorAnyDir->propagate(state,*stereoSurf);
00842         if (stereoState.isValid()==0) continue;
00843 
00844         LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
00845         GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
00846         LocalPoint stereoShitLPos = sStereoHit.localPosition();
00847         GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
00848 
00849         //LogTrace("TestTrackHits") << "assStereoSimHits.size()=" << assStereoSimHits.size() ;
00850         //LogTrace("TestTrackHits") << "stereo shit=" << stereoShitGPos ;
00851 
00852         GlobalVector stereoTsosGMom = stereoState.globalMomentum();
00853         GlobalError  stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00854         GlobalPoint  stereoTsosGPos = stereoState.globalPosition();
00855         GlobalError  stereoTsosGPEr = stereoState.cartesianError().position();
00856 
00857         GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
00858         GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
00859 
00860         MeasurementExtractor meSt(stereoState);
00861         double chi2stereo = computeChi2Increment(meSt,tStereoHit);
00862 
00863         title.str("");
00864         title << "Chi2Increment_stereo_" << subdetId << "-" << layerId ;
00865         hChi2Increment_stereo[title.str()]->Fill(chi2stereo);
00866 
00867         double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/sqrt(stereoRhitGPEr.cxx());
00868         double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/sqrt(stereoRhitGPEr.cyy());
00869         double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/sqrt(stereoRhitGPEr.czz());
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 
00874         title.str("");
00875         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
00876         hPullGP_X_rs_stereo[title.str()]->Fill( pullGPX_rs_stereo );
00877         title.str("");
00878         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
00879         hPullGP_Y_rs_stereo[title.str()]->Fill( pullGPY_rs_stereo );
00880         title.str("");
00881         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
00882         hPullGP_Z_rs_stereo[title.str()]->Fill( pullGPZ_rs_stereo );
00883 
00884         double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x())/sqrt(stereoTsosGPEr.cxx()+stereoRhitGPEr.cxx());
00885         double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y())/sqrt(stereoTsosGPEr.cyy()+stereoRhitGPEr.cyy());
00886         double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z())/sqrt(stereoTsosGPEr.czz()+stereoRhitGPEr.czz());
00887         //double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x());
00888         //double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y());
00889         //double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z());
00890 
00891         title.str("");
00892         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
00893         hPullGP_X_tr_stereo[title.str()]->Fill( pullGPX_tr_stereo );
00894         title.str("");
00895         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
00896         hPullGP_Y_tr_stereo[title.str()]->Fill( pullGPY_tr_stereo );
00897         title.str("");
00898         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
00899         hPullGP_Z_tr_stereo[title.str()]->Fill( pullGPZ_tr_stereo );
00900 
00901         double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x())/sqrt(stereoTsosGPEr.cxx());
00902         double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y())/sqrt(stereoTsosGPEr.cyy());
00903         double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z())/sqrt(stereoTsosGPEr.czz());
00904         //double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x());
00905         //double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y());
00906         //double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z());
00907 
00908         title.str("");
00909         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
00910         hPullGP_X_ts_stereo[title.str()]->Fill( pullGPX_ts_stereo );
00911         title.str("");
00912         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00913         hPullGP_Y_ts_stereo[title.str()]->Fill( pullGPY_ts_stereo );
00914         title.str("");
00915         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00916         hPullGP_Z_ts_stereo[title.str()]->Fill( pullGPZ_ts_stereo );
00917 
00918         double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x())/sqrt(stereoTsosGMEr.cxx());
00919         double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y())/sqrt(stereoTsosGMEr.cyy());
00920         double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z())/sqrt(stereoTsosGMEr.czz());
00921         //double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x());
00922         //double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y());
00923         //double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z());
00924 
00925         title.str("");
00926         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
00927         hPullGM_X_ts_stereo[title.str()]->Fill( pullGMX_ts_stereo );
00928         title.str("");
00929         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00930         hPullGM_Y_ts_stereo[title.str()]->Fill( pullGMY_ts_stereo );
00931         title.str("");
00932         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00933         hPullGM_Z_ts_stereo[title.str()]->Fill( pullGMZ_ts_stereo );
00934       }
00935     }
00936     LogTrace("TestTrackHits") << "traj chi2="  << tchi2 ;
00937     LogTrace("TestTrackHits") << "track chi2=" << track->chi2() ;
00938     i++;
00939   }
00940   LogTrace("TestTrackHits") << "end of event: processd hits=" << evtHits ;
00941   delete hitAssociator;
00942 }
00943 
00944 void TestTrackHits::endJob() {
00945   //file->Write();
00946   TDirectory * chi2d = file->mkdir("Chi2Increment");
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   chi2d->cd();
00993   hTotChi2Increment->Write();
00994   hTotChi2GoodHit->Write();
00995   hTotChi2BadHit->Write();
00996   hTotChi2DeltaHit->Write();
00997   hTotChi2NSharedHit->Write();
00998   hTotChi2SharedHit->Write();
00999   hProcess_vs_Chi2->Write();
01000   hClsize_vs_Chi2->Write();
01001   hPixClsize_vs_Chi2->Write();
01002   hPrjClsize_vs_Chi2->Write();
01003   hSt1Clsize_vs_Chi2->Write();
01004   hSt2Clsize_vs_Chi2->Write();
01005   hGoodHit_vs_Chi2->Write();
01006   hClusterSize->Write();
01007   hPixClusterSize->Write();
01008   hPrjClusterSize->Write();
01009   hSt1ClusterSize->Write();
01010   hSt2ClusterSize->Write();
01011   hSimHitVecSize->Write();
01012   hPixSimHitVecSize->Write();
01013   hPrjSimHitVecSize->Write();
01014   hSt1SimHitVecSize->Write();
01015   hSt2SimHitVecSize->Write();
01016   goodbadmerged->Write();
01017   energyLossRatio->Write();
01018   mergedPull->Write();
01019   probXgood->Write();
01020   probXbad->Write();
01021   probXdelta->Write();
01022   probXshared->Write();
01023   probXnoshare->Write();
01024   probYgood->Write();
01025   probYbad->Write();
01026   probYdelta->Write();
01027   probYshared->Write();
01028   probYnoshare->Write();
01029   for (int i=0; i!=6; i++)
01030     for (int j=0; j!=9; j++){
01031       if (i==0 && j>2) break;
01032       if (i==1 && j>1) break;
01033       if (i==2 && j>3) break;
01034       if (i==3 && j>2) break;
01035       if (i==4 && j>5) break;
01036       if (i==5 && j>8) break;
01037       chi2d->cd();
01038       title.str("");
01039       title << "Chi2Increment_" << i+1 << "-" << j+1 ;
01040       hChi2Increment[title.str()]->Write();
01041       title.str("");
01042       title << "Chi2IncrementVsEta_" << i+1 << "-" << j+1 ;
01043       hChi2IncrementVsEta[title.str()]->Write();
01044       title.str("");
01045       title << "Chi2GoodHit_" << i+1 << "-" << j+1 ;
01046       hChi2GoodHit[title.str()]->Write();
01047       title.str("");
01048       title << "Chi2BadHit_" << i+1 << "-" << j+1 ;
01049       hChi2BadHit[title.str()]->Write();
01050       title.str("");
01051       title << "Chi2DeltaHit_" << i+1 << "-" << j+1 ;
01052       hChi2DeltaHit[title.str()]->Write();
01053       title.str("");
01054       title << "Chi2NSharedHit_" << i+1 << "-" << j+1 ;
01055       hChi2NSharedHit[title.str()]->Write();
01056       title.str("");
01057       title << "Chi2SharedHit_" << i+1 << "-" << j+1 ;
01058       hChi2SharedHit[title.str()]->Write();
01059 
01060       gp_ts->cd();
01061       gp_tsx->cd();
01062       title.str("");
01063       title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
01064       hPullGP_X_ts[title.str()]->Write();
01065       gp_tsy->cd();
01066       title.str("");
01067       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
01068       hPullGP_Y_ts[title.str()]->Write();
01069       gp_tsz->cd();
01070       title.str("");
01071       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
01072       hPullGP_Z_ts[title.str()]->Write();
01073 
01074       gm_ts->cd();
01075       gm_tsx->cd();
01076       title.str("");
01077       title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
01078       hPullGM_X_ts[title.str()]->Write();
01079       gm_tsy->cd();
01080       title.str("");
01081       title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
01082       hPullGM_Y_ts[title.str()]->Write();
01083       gm_tsz->cd();
01084       title.str("");
01085       title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
01086       hPullGM_Z_ts[title.str()]->Write();
01087 
01088       gp_tr->cd();
01089       gp_trx->cd();
01090       title.str("");
01091       title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
01092       hPullGP_X_tr[title.str()]->Write();
01093       gp_try->cd();
01094       title.str("");
01095       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
01096       hPullGP_Y_tr[title.str()]->Write();
01097       gp_trz->cd();
01098       title.str("");
01099       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
01100       hPullGP_Z_tr[title.str()]->Write();
01101 
01102       gp_rs->cd();
01103       gp_rsx->cd();
01104       title.str("");
01105       title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
01106       hPullGP_X_rs[title.str()]->Write();
01107       gp_rsy->cd();
01108       title.str("");
01109       title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
01110       hPullGP_Y_rs[title.str()]->Write();
01111       gp_rsz->cd();
01112       title.str("");
01113       title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
01114       hPullGP_Z_rs[title.str()]->Write();
01115 
01116       if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
01117         chi2d->cd();
01118         title.str("");
01119         title << "Chi2Increment_mono_" << i+1 << "-" << j+1 ;
01120         hChi2Increment_mono[title.str()]->Write();
01121         title.str("");
01122         title << "Chi2Increment_stereo_" << i+1 << "-" << j+1 ;
01123         hChi2Increment_stereo[title.str()]->Write();
01124         //mono
01125         gp_ts->cd();
01126         gp_tsx_mono->cd();
01127         title.str("");
01128         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
01129         hPullGP_X_ts_mono[title.str()]->Write();
01130         gp_tsy_mono->cd();
01131         title.str("");
01132         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
01133         hPullGP_Y_ts_mono[title.str()]->Write();
01134         gp_tsz_mono->cd();
01135         title.str("");
01136         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
01137         hPullGP_Z_ts_mono[title.str()]->Write();
01138 
01139         gm_ts->cd();
01140         gm_tsx_mono->cd();
01141         title.str("");
01142         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
01143         hPullGM_X_ts_mono[title.str()]->Write();
01144         gm_tsy_mono->cd();
01145         title.str("");
01146         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
01147         hPullGM_Y_ts_mono[title.str()]->Write();
01148         gm_tsz_mono->cd();
01149         title.str("");
01150         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
01151         hPullGM_Z_ts_mono[title.str()]->Write();
01152 
01153         gp_tr->cd();
01154         gp_trx_mono->cd();
01155         title.str("");
01156         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
01157         hPullGP_X_tr_mono[title.str()]->Write();
01158         gp_try_mono->cd();
01159         title.str("");
01160         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
01161         hPullGP_Y_tr_mono[title.str()]->Write();
01162         gp_trz_mono->cd();
01163         title.str("");
01164         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
01165         hPullGP_Z_tr_mono[title.str()]->Write();
01166 
01167         gp_rs->cd();
01168         gp_rsx_mono->cd();
01169         title.str("");
01170         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
01171         hPullGP_X_rs_mono[title.str()]->Write();
01172         gp_rsy_mono->cd();
01173         title.str("");
01174         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
01175         hPullGP_Y_rs_mono[title.str()]->Write();
01176         gp_rsz_mono->cd();
01177         title.str("");
01178         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
01179         hPullGP_Z_rs_mono[title.str()]->Write();
01180 
01181         //stereo
01182         gp_ts->cd();
01183         gp_tsx_stereo->cd();
01184         title.str("");
01185         title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
01186         hPullGP_X_ts_stereo[title.str()]->Write();
01187         gp_tsy_stereo->cd();
01188         title.str("");
01189         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
01190         hPullGP_Y_ts_stereo[title.str()]->Write();
01191         gp_tsz_stereo->cd();
01192         title.str("");
01193         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
01194         hPullGP_Z_ts_stereo[title.str()]->Write();
01195 
01196         gm_ts->cd();
01197         gm_tsx_stereo->cd();
01198         title.str("");
01199         title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
01200         hPullGM_X_ts_stereo[title.str()]->Write();
01201         gm_tsy_stereo->cd();
01202         title.str("");
01203         title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
01204         hPullGM_Y_ts_stereo[title.str()]->Write();
01205         gm_tsz_stereo->cd();
01206         title.str("");
01207         title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
01208         hPullGM_Z_ts_stereo[title.str()]->Write();
01209 
01210         gp_tr->cd();
01211         gp_trx_stereo->cd();
01212         title.str("");
01213         title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
01214         hPullGP_X_tr_stereo[title.str()]->Write();
01215         gp_try_stereo->cd();
01216         title.str("");
01217         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
01218         hPullGP_Y_tr_stereo[title.str()]->Write();
01219         gp_trz_stereo->cd();
01220         title.str("");
01221         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
01222         hPullGP_Z_tr_stereo[title.str()]->Write();
01223 
01224         gp_rs->cd();
01225         gp_rsx_stereo->cd();
01226         title.str("");
01227         title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
01228         hPullGP_X_rs_stereo[title.str()]->Write();
01229         gp_rsy_stereo->cd();
01230         title.str("");
01231         title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
01232         hPullGP_Y_rs_stereo[title.str()]->Write();
01233         gp_rsz_stereo->cd();
01234         title.str("");
01235         title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
01236         hPullGP_Z_rs_stereo[title.str()]->Write();
01237       }
01238     }
01239 
01240   file->Close();
01241 }
01242 
01243 
01244 //needed by to do the residual for matched hits
01245 //taken from SiStripTrackingRecHitsValid.cc
01246 std::pair<LocalPoint,LocalVector> 
01247 TestTrackHits::projectHit( const PSimHit& hit, const StripGeomDetUnit* stripDet, const BoundPlane& plane) 
01248 { 
01249   const StripTopology& topol = stripDet->specificTopology();
01250   GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
01251   LocalPoint localHit = plane.toLocal(globalpos);
01252   //track direction
01253   LocalVector locdir=hit.localDirection();
01254   //rotate track in new frame
01255    
01256   GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
01257   LocalVector dir=plane.toLocal(globaldir);
01258   float scale = -localHit.z() / dir.z();
01259    
01260   LocalPoint projectedPos = localHit + scale*dir;
01261    
01262   float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
01263    
01264   LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame
01265    
01266   LocalVector localStripDir( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
01267    
01268   return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
01269 }
01270 
01271 template<unsigned int D> 
01272 double TestTrackHits::computeChi2Increment(MeasurementExtractor me, 
01273                                            TransientTrackingRecHit::ConstRecHitPointer rhit) {
01274   typedef typename AlgebraicROOTObject<D>::Vector VecD;
01275   typedef typename AlgebraicROOTObject<D,D>::SymMatrix SMatDD;
01276   VecD r = asSVector<D>(rhit->parameters()) - me.measuredParameters<D>(*rhit);
01277   
01278   SMatDD R = asSMatrix<D>(rhit->parametersError()) + me.measuredError<D>(*rhit);
01279   R.Invert();
01280   return ROOT::Math::Similarity(r,R) ;
01281 }
01282 
01283 #include "FWCore/Framework/interface/ModuleFactory.h"
01284 #include "FWCore/Framework/interface/MakerMacros.h"
01285 DEFINE_FWK_MODULE(TestTrackHits);