CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TestTrackHits.cc
Go to the documentation of this file.
7 #include <TDirectory.h>
8 
12 
15 using namespace std;
16 using namespace edm;
17 
19  trackerHitAssociatorConfig_(consumesCollector()) {
20  LogTrace("TestTrackHits") << iConfig;
21  propagatorName = iConfig.getParameter<std::string>("Propagator");
22  builderName = iConfig.getParameter<std::string>("TTRHBuilder");
23  srcName = iConfig.getParameter<std::string>("src");
24  tpName = iConfig.getParameter<std::string>("tpname");
25  updatorName = iConfig.getParameter<std::string>("updator");
26  out = iConfig.getParameter<std::string>("out");
27 // ParameterSet cuts = iConfig.getParameter<ParameterSet>("RecoTracksCuts");
28 // selectRecoTracks = RecoTrackSelector(cuts.getParameter<double>("ptMin"),
29 // cuts.getParameter<double>("minRapidity"),
30 // cuts.getParameter<double>("maxRapidity"),
31 // cuts.getParameter<double>("tip"),
32 // cuts.getParameter<double>("lip"),
33 // cuts.getParameter<int>("minHit"),
34 // cuts.getParameter<double>("maxChi2"));
35 }
36 
38 
40 {
41  iSetup.get<TrackerDigiGeometryRecord>().get(theG);
42  iSetup.get<IdealMagneticFieldRecord>().get(theMF);
46 
47  file = new TFile(out.c_str(),"recreate");
48  for (int i=0; i!=6; i++)
49  for (int j=0; j!=9; j++){
50  if (i==0 && j>2) break;
51  if (i==1 && j>1) break;
52  if (i==2 && j>3) break;
53  if (i==3 && j>2) break;
54  if (i==4 && j>5) break;
55  if (i==5 && j>8) break;
56 
57  title.str("");
58  title << "Chi2Increment_" << i+1 << "-" << j+1 ;
59  hChi2Increment[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
60  title.str("");
61  title << "Chi2IncrementVsEta_" << i+1 << "-" << j+1 ;
62  hChi2IncrementVsEta[title.str()] = new TH2F(title.str().c_str(),title.str().c_str(),50,-2.5,2.5,1000,0,100);
63  title.str("");
64  title << "Chi2GoodHit_" << i+1 << "-" << j+1 ;
65  hChi2GoodHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
66  title.str("");
67  title << "Chi2BadHit_" << i+1 << "-" << j+1 ;
68  hChi2BadHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
69  title.str("");
70  title << "Chi2DeltaHit_" << i+1 << "-" << j+1 ;
71  hChi2DeltaHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
72  title.str("");
73  title << "Chi2NSharedHit_" << i+1 << "-" << j+1 ;
74  hChi2NSharedHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
75  title.str("");
76  title << "Chi2SharedHit_" << i+1 << "-" << j+1 ;
77  hChi2SharedHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
78 
79  title.str("");
80  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
81  hPullGP_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
82  title.str("");
83  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
84  hPullGP_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
85  title.str("");
86  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
87  hPullGP_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
88 
89  title.str("");
90  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
91  hPullGM_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
92  title.str("");
93  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
94  hPullGM_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
95  title.str("");
96  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
97  hPullGM_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
98 
99  title.str("");
100  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
101  hPullGP_X_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
102  title.str("");
103  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
104  hPullGP_Y_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
105  title.str("");
106  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
107  hPullGP_Z_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
108 
109  title.str("");
110  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
111  hPullGP_X_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
112  title.str("");
113  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
114  hPullGP_Y_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
115  title.str("");
116  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
117  hPullGP_Z_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
118 
119  if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
120  //mono
121  title.str("");
122  title << "Chi2Increment_mono_" << i+1 << "-" << j+1 ;
123  hChi2Increment_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
124 
125  title.str("");
126  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
127  hPullGP_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
128  title.str("");
129  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
130  hPullGP_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
131  title.str("");
132  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
133  hPullGP_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
134 
135  title.str("");
136  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
137  hPullGM_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
138  title.str("");
139  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
140  hPullGM_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
141  title.str("");
142  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
143  hPullGM_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
144 
145  title.str("");
146  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
147  hPullGP_X_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
148  title.str("");
149  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
150  hPullGP_Y_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
151  title.str("");
152  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
153  hPullGP_Z_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
154 
155  title.str("");
156  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
157  hPullGP_X_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
158  title.str("");
159  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
160  hPullGP_Y_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
161  title.str("");
162  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
163  hPullGP_Z_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
164 
165  //stereo
166  title.str("");
167  title << "Chi2Increment_stereo_" << i+1 << "-" << j+1 ;
168  hChi2Increment_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
169 
170  title.str("");
171  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
172  hPullGP_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
173  title.str("");
174  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
175  hPullGP_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
176  title.str("");
177  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
178  hPullGP_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
179 
180  title.str("");
181  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
182  hPullGM_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
183  title.str("");
184  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
185  hPullGM_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
186  title.str("");
187  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
188  hPullGM_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
189 
190  title.str("");
191  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
192  hPullGP_X_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
193  title.str("");
194  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
195  hPullGP_Y_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
196  title.str("");
197  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
198  hPullGP_Z_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
199 
200  title.str("");
201  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
202  hPullGP_X_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
203  title.str("");
204  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
205  hPullGP_Y_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
206  title.str("");
207  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
208  hPullGP_Z_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
209  }
210  }
211  hTotChi2Increment = new TH1F("TotChi2Increment","TotChi2Increment",1000,0,100);
212  hTotChi2GoodHit = new TH1F("TotChi2GoodHit","TotChi2GoodHit",1000,0,100);
213  hTotChi2BadHit = new TH1F("TotChi2BadHit","TotChi2BadHit",1000,0,100);
214  hTotChi2DeltaHit = new TH1F("TotChi2DeltaHit","TotChi2DeltaHit",1000,0,100);
215  hTotChi2NSharedHit = new TH1F("TotChi2NSharedHit","TotChi2NSharedHit",1000,0,100);
216  hTotChi2SharedHit = new TH1F("TotChi2SharedHit","TotChi2SharedHit",1000,0,100);
217  hProcess_vs_Chi2 = new TH2F("Process_vs_Chi2","Process_vs_Chi2",1000,0,100,17,-0.5,16.5);
218  hClsize_vs_Chi2 = new TH2F("Clsize_vs_Chi2","Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
219  hPixClsize_vs_Chi2= new TH2F("PixClsize_vs_Chi2","PixClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
220  hPrjClsize_vs_Chi2= new TH2F("PrjClsize_vs_Chi2","PrjClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
221  hSt1Clsize_vs_Chi2= new TH2F("St1Clsize_vs_Chi2","St1Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
222  hSt2Clsize_vs_Chi2= new TH2F("St2Clsize_vs_Chi2","St2Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
223  hGoodHit_vs_Chi2 = new TH2F("GoodHit_vs_Chi2","GoodHit_vs_Chi2",10000,0,1000,2,-0.5,1.5);
224  hClusterSize = new TH1F("ClusterSize","ClusterSize",40,-0.5,39.5);
225  hPixClusterSize = new TH1F("PixClusterSize","PixClusterSize",40,-0.5,39.5);
226  hPrjClusterSize = new TH1F("PrjClusterSize","PrjClusterSize",40,-0.5,39.5);
227  hSt1ClusterSize = new TH1F("St1ClusterSize","St1ClusterSize",40,-0.5,39.5);
228  hSt2ClusterSize = new TH1F("St2ClusterSize","St2ClusterSize",40,-0.5,39.5);
229  hSimHitVecSize = new TH1F("hSimHitVecSize","hSimHitVecSize",40,-0.5,39.5);
230  hPixSimHitVecSize = new TH1F("PixSimHitVecSize","PixSimHitVecSize",40,-0.5,39.5);
231  hPrjSimHitVecSize = new TH1F("PrjSimHitVecSize","PrjSimHitVecSize",40,-0.5,39.5);
232  hSt1SimHitVecSize = new TH1F("St1SimHitVecSize","St1SimHitVecSize",40,-0.5,39.5);
233  hSt2SimHitVecSize = new TH1F("St2SimHitVecSize","St2SimHitVecSize",40,-0.5,39.5);
234  goodbadmerged = new TH1F("goodbadmerged","goodbadmerged",5,0.5,5.5);
235  energyLossRatio = new TH1F("energyLossRatio","energyLossRatio",100,0,1);
236  mergedPull = new TH1F("mergedPull","mergedPull",200,0,20);
237  probXgood = new TH1F("probXgood","probXgood",110,0,1.1);
238  probXbad = new TH1F("probXbad","probXbad",110,0,1.1);
239  probXdelta = new TH1F("probXdelta","probXdelta",110,0,1.1);
240  probXshared = new TH1F("probXshared","probXshared",110,0,1.1);
241  probXnoshare = new TH1F("probXnoshare","probXnoshare",110,0,1.1);
242  probYgood = new TH1F("probYgood","probYgood",110,0,1.1);
243  probYbad = new TH1F("probYbad","probYbad",110,0,1.1);
244  probYdelta = new TH1F("probYdelta","probYdelta",110,0,1.1);
245  probYshared = new TH1F("probYshared","probYshared",110,0,1.1);
246  probYnoshare = new TH1F("probYnoshare","probYnoshare",110,0,1.1);
247 }
248 
250 {
251  //Retrieve tracker topology from geometry
253  iSetup.get<TrackerTopologyRcd>().get(tTopo);
254 
255 
256  LogDebug("TestTrackHits") << "new event" ;
262  iEvent.getByLabel("offlineBeamSpot",beamSpot);
263 
264  LogTrace("TestTrackHits") << "Tr collection size=" << trackCollectionHandle->size();
265  LogTrace("TestTrackHits") << "TP collection size=" << trackingParticleCollectionHandle->size();
266 
267  iEvent.getByLabel("trackAssociatorByHits",trackAssociator);
268 
269 
270  TrackerHitAssociator hitAssociator(iEvent, trackerHitAssociatorConfig_);
271 
272  reco::RecoToSimCollection recSimColl=trackAssociator->associateRecoToSim(trackCollectionHandle,
274 
276 
277  int evtHits = 0;
278  int i=0;
279  int yy=0;
280  int yyy=0;
281  for(std::vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it!=trajCollectionHandle->end();it++){
282 
283  LogTrace("TestTrackHits") << "\n*****************new trajectory********************" ;
284  double tchi2 = 0;
285 
286  std::vector<TrajectoryMeasurement> tmColl = it->measurements();
287 
289  reco::TrackRef tmptrack = (*trajTrackAssociationCollectionHandle.product())[traj];
290  edm::RefToBase<reco::Track> track(tmptrack);
291  // if ( !selectRecoTracks( *track,beamSpot.product() ) ) {
292  // LogTrace("TestTrackHits") << "track does not pass quality cuts: skippingtrack #" << ++yyy;
293  // i++;
294  // continue;
295  // }
296 
297  std::vector<std::pair<TrackingParticleRef, double> > tP;
298  if(recSimColl.find(track) != recSimColl.end()){
299  tP = recSimColl[track];
300  if (tP.size()!=0) {
301  edm::LogVerbatim("TestTrackHits") << "reco::Track #" << ++yyy << " with pt=" << track->pt()
302  << " associated with quality:" << tP.begin()->second <<" good track #" << ++yy << " has hits:" << track->numberOfValidHits() << "\n";
303  }
304  }else{
305  edm::LogVerbatim("TestTrackHits") << "reco::Track #" << ++yyy << " with pt=" << track->pt()
306  << " NOT associated to any TrackingParticle" << "\n";
307  i++;
308  continue;
309  }
310  // if(recSimColl.find(track) != recSimColl.end()) {
311  // tP = recSimColl[track];
312  // } else {
313  // LogTrace("TestTrackHits") << "fake track: skipping track " << ++yyy;
314  // continue;//skip fake tracks
315  // }
316  // if (tP.size()==0) {
317  // LogTrace("TestTrackHits") << "fake track: skipping track " << ++yyy;
318  // continue;
319  // }
320  TrackingParticleRef tp = tP.begin()->first;
321  LogTrace("TestTrackHits") << "a tp is associated with fraction=" << tP.begin()->second;
322  //LogTrace("TestTrackHits") << "last tp is associated with fraction=" << (tP.end()-1)->second;
323  std::vector<unsigned int> tpids;
324  for (TrackingParticle::g4t_iterator g4T=tp->g4Track_begin(); g4T!=tp->g4Track_end(); ++g4T) {
325  LogTrace("TestTrackHits") << "tp id=" << g4T->trackId();
326  tpids.push_back(g4T->trackId());
327  }
328 
329  //LogTrace("TestTrackHits") << "Analyzing hits of track number " << ++yyy << " good track number " << ++yy;
330  int pp = 0;
331  for (std::vector<TrajectoryMeasurement>::iterator tm=tmColl.begin();tm!=tmColl.end();++tm){
332 
333  tchi2+=tm->estimate();
334 
335  LogTrace("TestTrackHits") << "+++++++++++++++++new hit+++++++++++++++++" ;
336  CTTRHp rhit = tm->recHit();
337  //TSOS state = tm->backwardPredictedState();
338  //TSOS state = tm->forwardPredictedState();
339  TSOS state = combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
340 
341  if (rhit->isValid()==0 && rhit->det()!=0) continue;
342  evtHits++;
343  LogTrace("TestTrackHits") << "valid hit #" << ++pp << "of hits=" << track->numberOfValidHits();
344 
345  int subdetId = rhit->det()->geographicalId().subdetId();
346  DetId id = rhit->det()->geographicalId();
347  int layerId = tTopo->layer(id);
348  LogTrace("TestTrackHits") << "subdetId=" << subdetId << " layerId=" << layerId ;
349 
350  const Surface * surf = rhit->surface();
351  if (surf==0) continue;
352 
353  double energyLoss_ = 0.;
354  unsigned int monoId = 0;
355  std::vector<double> energyLossM;
356  std::vector<double> energyLossS;
357  std::vector<PSimHit> assSimHits = hitAssociator.associateHit(*(rhit)->hit());
358  unsigned int simhitvecsize = assSimHits.size();
359  if (simhitvecsize==0) continue;
360  PSimHit shit;
361  std::vector<unsigned int> trackIds;
362  energyLossS.clear();
363  energyLossM.clear();
364  for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
365  unsigned int tId = m->trackId();
366  if (find(trackIds.begin(),trackIds.end(),tId)==trackIds.end()) trackIds.push_back(tId);
367  if (m->energyLoss()>energyLoss_) {
368  shit=*m;
369  energyLoss_ = m->energyLoss();
370  }
371  if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
372  if (monoId==0) monoId = m->detUnitId();
373  if (monoId == m->detUnitId()){
374  energyLossM.push_back(m->energyLoss());
375  }
376  else {
377  energyLossS.push_back(m->energyLoss());
378  }
379  //std::cout << "detUnitId=" << m->detUnitId() << " trackId=" << m->trackId() << " energyLoss=" << m->energyLoss() << std::endl;
380  } else {
381  energyLossM.push_back(m->energyLoss());
382  }
383  }
384  //double delta = 99999;
385  //LocalPoint rhitLPv = rhit->localPosition();
386  //vector<PSimHit> assSimHits = hitAssociator.associateHit(*(rhit)->hit());
387  //unsigned int simhitvecsize = assSimHits.size();
388  //if (simhitvecsize==0) continue;
389  //PSimHit shit;
390  //for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
391  // if ((m->localPosition()-rhitLPv).mag()<delta) {
392  // shit=*m;
393  // delta = (m->localPosition()-rhitLPv).mag();
394  // }
395  //}
396 
397  //plot chi2 increment
398  double chi2increment = tm->estimate();
399  LogTrace("TestTrackHits") << "tm->estimate()=" << tm->estimate();
400  title.str("");
401  title << "Chi2Increment_" << subdetId << "-" << layerId;
402  hChi2Increment[title.str()]->Fill( chi2increment );
403  title.str("");
404  title << "Chi2IncrementVsEta_" << subdetId << "-" << layerId;
405  hChi2IncrementVsEta[title.str()]->Fill( track->eta(), chi2increment );
406  hTotChi2Increment->Fill( chi2increment );
407  hProcess_vs_Chi2->Fill( chi2increment, shit.processType() );
408 
409  int clustersize = 0;
410  bool mergedhit = false;
411  if (dynamic_cast<const SiPixelRecHit*>(rhit->hit())){
412  clustersize = ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size() ;
413  hPixClsize_vs_Chi2->Fill(chi2increment, clustersize);
414  hPixClusterSize->Fill(clustersize);
415  hPixSimHitVecSize->Fill(simhitvecsize);
416  if (simhitvecsize>1) mergedhit = true;
417  }
418  else if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit())){
419  clustersize = ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() ;
420  hSt1Clsize_vs_Chi2->Fill(chi2increment, clustersize);
421  hSt1ClusterSize->Fill(clustersize);
422  hSt1SimHitVecSize->Fill(simhitvecsize);
423  if (simhitvecsize>1) mergedhit = true;
424  }
425  else if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())){
426  int clsize1 = ((const SiStripMatchedRecHit2D*)(rhit->hit()))->monoCluster().amplitudes().size();
427  int clsize2 = ((const SiStripMatchedRecHit2D*)(rhit->hit()))->stereoCluster().amplitudes().size();
428  if (clsize1>clsize2) clustersize = clsize1;
429  else clustersize = clsize2;
430  hSt2Clsize_vs_Chi2->Fill(chi2increment, clustersize);
431  hSt2ClusterSize->Fill(clustersize);
432  hSt2SimHitVecSize->Fill(simhitvecsize);
433  if (simhitvecsize>2) mergedhit = true;
434  }
435  else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(rhit->hit())){
436  clustersize = ((const ProjectedSiStripRecHit2D*)(rhit->hit()))->originalHit().cluster()->amplitudes().size();
437  hPrjClsize_vs_Chi2->Fill(chi2increment, clustersize);
438  hPrjClusterSize->Fill(clustersize);
439  hPrjSimHitVecSize->Fill(simhitvecsize);
440  if (simhitvecsize>1) mergedhit = true;
441  }
442  hClsize_vs_Chi2->Fill( chi2increment, clustersize);
443  hClusterSize->Fill(clustersize);
444  hSimHitVecSize->Fill(simhitvecsize);
445 
446  // if (dynamic_cast<const SiPixelRecHit*>(rhit->hit()))
447  // hClsize_vs_Chi2->Fill( chi2increment, ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size() );
448  // if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit()))
449  // hClsize_vs_Chi2->Fill( chi2increment, ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() );
450 
451  std::vector<SimHitIdpr> simTrackIds = hitAssociator.associateHitId(*(rhit)->hit());
452  bool goodhit = false;
453  for(size_t j=0; j<simTrackIds.size(); j++){
454  LogTrace("TestTrackHits") << "hit id=" << simTrackIds[j].first;
455  for (size_t jj=0; jj<tpids.size(); jj++){
456  if (simTrackIds[j].first == tpids[jj]) goodhit = true;
457  break;
458  }
459  }
460  bool shared = true;
461  bool ioniOnly = true;
462  const SiPixelRecHit* pix = dynamic_cast<const SiPixelRecHit*>(rhit->hit());
463  if (goodhit) {
464  if (energyLossM.size()>1&&energyLossS.size()<=1) {
465  sort(energyLossM.begin(),energyLossM.end(),greater<double>());
466  energyLossRatio->Fill(energyLossM[1]/energyLossM[0]);
467  }
468  else if (energyLossS.size()>1&&energyLossM.size()<=1) {
469  sort(energyLossS.begin(),energyLossS.end(),greater<double>());
470  energyLossRatio->Fill(energyLossS[1]/energyLossS[0]);
471  }
472  else if (energyLossS.size()>1&&energyLossM.size()>1) {
473  sort(energyLossM.begin(),energyLossM.end(),greater<double>());
474  sort(energyLossS.begin(),energyLossS.end(),greater<double>());
475  energyLossM[1]/energyLossM[0] > energyLossS[1]/energyLossS[0]
476  ? energyLossRatio->Fill(energyLossM[1]/energyLossM[0])
477  : energyLossRatio->Fill(energyLossS[1]/energyLossS[0]);
478  }
479  if ( mergedhit ) {
480  //not optimized for matched hits
481  LogVerbatim("TestTrackHits") << "MERGED HIT" << std::endl;
482  unsigned int idc = 0;
483  for (size_t jj=0; jj<tpids.size(); jj++){
484  idc += std::count(trackIds.begin(),trackIds.end(),tpids[jj]);
485  }
486  if (idc==trackIds.size()) {
487  shared = false;
488  }
489  for(std::vector<PSimHit>::const_iterator m=assSimHits.begin()+1; m<assSimHits.end(); m++){
490  if ((m->processType()!=7&&m->processType()!=8&&m->processType()!=9)&&abs(m->particleType())!=11){
491  ioniOnly = false;
492  break;
493  }
494  }
495  if (ioniOnly&&!shared) {
496  title.str("");
497  title << "Chi2DeltaHit_" << subdetId << "-" << layerId;
498  hChi2DeltaHit[title.str()]->Fill( chi2increment );
499  hTotChi2DeltaHit->Fill( chi2increment );
500  if (pix) {
501  probXdelta->Fill(pix->probabilityX());
502  probYdelta->Fill(pix->probabilityY());
503  }
504  } else if(!ioniOnly&&!shared) {
505  title.str("");
506  title << "Chi2NSharedHit_" << subdetId << "-" << layerId;
507  hChi2NSharedHit[title.str()]->Fill( chi2increment );
508  hTotChi2NSharedHit->Fill( chi2increment );
509  if (pix) {
510  probXnoshare->Fill(pix->probabilityX());
511  probYnoshare->Fill(pix->probabilityY());
512  }
513  } else {
514  title.str("");
515  title << "Chi2SharedHit_" << subdetId << "-" << layerId;
516  hChi2SharedHit[title.str()]->Fill( chi2increment );
517  hTotChi2SharedHit->Fill( chi2increment );
518  if (pix) {
519  probXshared->Fill(pix->probabilityX());
520  probYshared->Fill(pix->probabilityY());
521  }
522  }
523 
524  for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++) {
525  unsigned int tId = m->trackId();
526  LogVerbatim("TestTrackHits") << "component with id=" << tId << " eLoss=" << m->energyLoss() << " pType=" << m->processType();
527  if (find(tpids.begin(),tpids.end(),tId)==tpids.end()) continue;
528  if (m->processType()==2) {
529  GlobalPoint gpr = rhit->globalPosition();
530  AlgebraicSymMatrix33 ger = rhit->globalPositionError().matrix();
531  GlobalPoint gps = surf->toGlobal(m->localPosition());
532  LogVerbatim("TestTrackHits") << gpr << " " << gps << " " << ger;
533  ROOT::Math::SVector<double,3> delta;
534  delta[0]=gpr.x()-gps.x();
535  delta[1]=gpr.y()-gps.y();
536  delta[2]=gpr.z()-gps.z();
537  LogVerbatim("TestTrackHits") << delta << " " << ger ;
538  double mpull = sqrt(delta[0]*delta[0]/ger[0][0]+delta[1]*delta[1]/ger[1][1]+delta[2]*delta[2]/ger[2][2]);
539  LogVerbatim("TestTrackHits") << "hit pull=" << mpull;//ger.similarity(delta);
540  mergedPull->Fill(mpull);
541  break;
542  }
543  }
544  } else {
545  LogVerbatim("TestTrackHits") << "good hit" ;
546  title.str("");
547  title << "Chi2GoodHit_" << subdetId << "-" << layerId;
548  hChi2GoodHit[title.str()]->Fill( chi2increment );
549  hTotChi2GoodHit->Fill( chi2increment );
550  if (pix) {
551  probXgood->Fill(pix->probabilityX());
552  probYgood->Fill(pix->probabilityY());
553  }
554  }
555  } else {
556  LogVerbatim("TestTrackHits") << "bad hit" ;
557  title.str("");
558  title << "Chi2BadHit_" << subdetId << "-" << layerId;
559  hChi2BadHit[title.str()]->Fill( chi2increment );
560  hTotChi2BadHit->Fill( chi2increment );
561  goodbadmerged->Fill(2);
562  if (pix) {
563  probXbad->Fill(pix->probabilityX());
564  probYbad->Fill(pix->probabilityY());
565  }
566  }
567  hGoodHit_vs_Chi2->Fill(chi2increment,goodhit);
568 
569  LocalVector shitLMom;
570  LocalPoint shitLPos;
571  if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
572  if (simhitvecsize>2 && goodhit) {
573  if (ioniOnly&&!shared) goodbadmerged->Fill(3);
574  else if(!ioniOnly&&!shared) goodbadmerged->Fill(4);
575  else goodbadmerged->Fill(5);
576  }
577  else if (goodhit) goodbadmerged->Fill(1);
578  double rechitmatchedx = rhit->localPosition().x();
579  double rechitmatchedy = rhit->localPosition().y();
580  double mindist = 999999;
581  float distx, disty;
582  std::pair<LocalPoint,LocalVector> closestPair;
583  const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) ((const GluedGeomDet *)(rhit)->det())->stereoDet();
584  const BoundPlane& plane = (rhit)->det()->surface();
585  for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
586  //project simhit;
587  std::pair<LocalPoint,LocalVector> hitPair = projectHit((*m),stripDet,plane);
588  distx = fabs(rechitmatchedx - hitPair.first.x());
589  disty = fabs(rechitmatchedy - hitPair.first.y());
590  double dist = distx*distx+disty*disty;
591  if(sqrt(dist)<mindist){
592  mindist = dist;
593  closestPair = hitPair;
594  }
595  }
596  shitLPos = closestPair.first;
597  shitLMom = closestPair.second;
598  } else {
599  if (simhitvecsize>1 && goodhit) {
600  if (ioniOnly&&!shared) goodbadmerged->Fill(3);
601  else if(!ioniOnly&&!shared) goodbadmerged->Fill(4);
602  else goodbadmerged->Fill(5);
603  }
604  else if (goodhit) goodbadmerged->Fill(1);
605  shitLPos = shit.localPosition();
606  shitLMom = shit.momentumAtEntry();
607  }
608  GlobalVector shitGMom = surf->toGlobal(shitLMom);
609  GlobalPoint shitGPos = surf->toGlobal(shitLPos);
610 
611  GlobalVector tsosGMom = state.globalMomentum();
612  GlobalError tsosGMEr(state.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
613  GlobalPoint tsosGPos = state.globalPosition();
614  GlobalError tsosGPEr = state.cartesianError().position();
615 
616  GlobalPoint rhitGPos = (rhit)->globalPosition();
617  GlobalError rhitGPEr = (rhit)->globalPositionError();
618 
619  LogVerbatim("TestTrackHits") << "assSimHits.size()=" << assSimHits.size() ;
620  LogVerbatim("TestTrackHits") << "tsos globalPos =" << tsosGPos ;
621  LogVerbatim("TestTrackHits") << "sim hit globalPos=" << shitGPos ;
622  LogVerbatim("TestTrackHits") << "rec hit globalPos=" << rhitGPos ;
623  LogVerbatim("TestTrackHits") << "geographicalId =" << rhit->det()->geographicalId().rawId() ;
624  LogVerbatim("TestTrackHits") << "surface position =" << surf->position() ;
625 
626 # if 0
627  if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->geographicalId()="
628  << rhit->detUnit()->geographicalId().rawId() ;
629  LogTrace("TestTrackHits") << "rhit->det()->surface().position()="
630  << rhit->det()->surface().position() ;
631  if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().position()="
632  << rhit->detUnit()->surface().position() ;
633  LogTrace("TestTrackHits") << "rhit->det()->position()=" << rhit->det()->position() ;
634  if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->position()="
635  << rhit->detUnit()->position() ;
636  LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().length()="
637  << rhit->det()->surface().bounds().length() ;
638  if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().length()="
639  << rhit->detUnit()->surface().bounds().length() ;
640  LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().width()="
641  << rhit->det()->surface().bounds().width() ;
642  if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().width()="
643  << rhit->detUnit()->surface().bounds().width() ;
644  LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().thickness()="
645  << rhit->det()->surface().bounds().thickness() ;
646  if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().thickness()="
647  << rhit->detUnit()->surface().bounds().thickness() ;
648 #endif
649 
650  double pullGPX_rs = (rhitGPos.x()-shitGPos.x())/sqrt(rhitGPEr.cxx());
651  double pullGPY_rs = (rhitGPos.y()-shitGPos.y())/sqrt(rhitGPEr.cyy());
652  double pullGPZ_rs = (rhitGPos.z()-shitGPos.z())/sqrt(rhitGPEr.czz());
653  //double pullGPX_rs = (rhitGPos.x()-shitGPos.x());
654  //double pullGPY_rs = (rhitGPos.y()-shitGPos.y());
655  //double pullGPZ_rs = (rhitGPos.z()-shitGPos.z());
656 
657  LogTrace("TestTrackHits") << "rs" ;
658 
659  title.str("");
660  title << "PullGP_X_" << subdetId << "-" << layerId << "_rs";
661  hPullGP_X_rs[title.str()]->Fill( pullGPX_rs );
662  title.str("");
663  title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs";
664  hPullGP_Y_rs[title.str()]->Fill( pullGPY_rs );
665  title.str("");
666  title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs";
667  hPullGP_Z_rs[title.str()]->Fill( pullGPZ_rs );
668 
669  double pullGPX_tr = (tsosGPos.x()-rhitGPos.x())/sqrt(tsosGPEr.cxx()+rhitGPEr.cxx());
670  double pullGPY_tr = (tsosGPos.y()-rhitGPos.y())/sqrt(tsosGPEr.cyy()+rhitGPEr.cyy());
671  double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z())/sqrt(tsosGPEr.czz()+rhitGPEr.czz());
672  //double pullGPX_tr = (tsosGPos.x()-rhitGPos.x());
673  //double pullGPY_tr = (tsosGPos.y()-rhitGPos.y());
674  //double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z());
675 
676  LogTrace("TestTrackHits") << "tr" ;
677 
678  title.str("");
679  title << "PullGP_X_" << subdetId << "-" << layerId << "_tr";
680  hPullGP_X_tr[title.str()]->Fill( pullGPX_tr );
681  title.str("");
682  title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr";
683  hPullGP_Y_tr[title.str()]->Fill( pullGPY_tr );
684  title.str("");
685  title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr";
686  hPullGP_Z_tr[title.str()]->Fill( pullGPZ_tr );
687 
688  double pullGPX_ts = (tsosGPos.x()-shitGPos.x())/sqrt(tsosGPEr.cxx());
689  double pullGPY_ts = (tsosGPos.y()-shitGPos.y())/sqrt(tsosGPEr.cyy());
690  double pullGPZ_ts = (tsosGPos.z()-shitGPos.z())/sqrt(tsosGPEr.czz());
691  //double pullGPX_ts = (tsosGPos.x()-shitGPos.x());
692  //double pullGPY_ts = (tsosGPos.y()-shitGPos.y());
693  //double pullGPZ_ts = (tsosGPos.z()-shitGPos.z());
694 
695  LogTrace("TestTrackHits") << "ts1" ;
696 
697  title.str("");
698  title << "PullGP_X_" << subdetId << "-" << layerId << "_ts";
699  hPullGP_X_ts[title.str()]->Fill( pullGPX_ts );
700  title.str("");
701  title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts";
702  hPullGP_Y_ts[title.str()]->Fill( pullGPY_ts );
703  title.str("");
704  title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts";
705  hPullGP_Z_ts[title.str()]->Fill( pullGPZ_ts );
706 
707  double pullGMX_ts = (tsosGMom.x()-shitGMom.x())/sqrt(tsosGMEr.cxx());
708  double pullGMY_ts = (tsosGMom.y()-shitGMom.y())/sqrt(tsosGMEr.cyy());
709  double pullGMZ_ts = (tsosGMom.z()-shitGMom.z())/sqrt(tsosGMEr.czz());
710  //double pullGMX_ts = (tsosGMom.x()-shitGMom.x());
711  //double pullGMY_ts = (tsosGMom.y()-shitGMom.y());
712  //double pullGMZ_ts = (tsosGMom.z()-shitGMom.z());
713 
714  LogTrace("TestTrackHits") << "ts2" ;
715 
716  title.str("");
717  title << "PullGM_X_" << subdetId << "-" << layerId << "_ts";
718  hPullGM_X_ts[title.str()]->Fill( pullGMX_ts );
719  title.str("");
720  title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts";
721  hPullGM_Y_ts[title.str()]->Fill( pullGMY_ts );
722  title.str("");
723  title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts";
724  hPullGM_Z_ts[title.str()]->Fill( pullGMZ_ts );
725  if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
726  Propagator* thePropagatorAnyDir = new PropagatorWithMaterial(anyDirection,0.105,theMF.product(),1.6);
727  //mono
728  LogTrace("TestTrackHits") << "MONO HIT" ;
729  auto m = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->monoHit();
730  CTTRHp tMonoHit = theBuilder->build(&m);
731  if (tMonoHit==0) continue;
732  vector<PSimHit> assMonoSimHits = hitAssociator.associateHit(*tMonoHit->hit());
733  if (assMonoSimHits.size()==0) continue;
734  const PSimHit sMonoHit = *(assMonoSimHits.begin());
735  const Surface * monoSurf = &( tMonoHit->det()->surface() );
736  if (monoSurf==0) continue;
737  TSOS monoState = thePropagatorAnyDir->propagate(state,*monoSurf);
738  if (monoState.isValid()==0) continue;
739 
740  LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
741  GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
742  LocalPoint monoShitLPos = sMonoHit.localPosition();
743  GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
744 
745  //LogTrace("TestTrackHits") << "assMonoSimHits.size()=" << assMonoSimHits.size() ;
746  //LogTrace("TestTrackHits") << "mono shit=" << monoShitGPos ;
747 
748  GlobalVector monoTsosGMom = monoState.globalMomentum();
749  GlobalError monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
750  GlobalPoint monoTsosGPos = monoState.globalPosition();
751  GlobalError monoTsosGPEr = monoState.cartesianError().position();
752 
753  GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
754  GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
755 
756  double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x())/sqrt(monoRhitGPEr.cxx());
757  double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y())/sqrt(monoRhitGPEr.cyy());
758  double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z())/sqrt(monoRhitGPEr.czz());
759  //double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x());
760  //double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y());
761  //double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z());
762 
763  MeasurementExtractor meMo(monoState);
764  double chi2mono = computeChi2Increment(meMo,tMonoHit);
765 
766  title.str("");
767  title << "Chi2Increment_mono_" << subdetId << "-" << layerId ;
768  hChi2Increment_mono[title.str()]->Fill(chi2mono);
769 
770  title.str("");
771  title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
772  hPullGP_X_rs_mono[title.str()]->Fill( pullGPX_rs_mono );
773  title.str("");
774  title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
775  hPullGP_Y_rs_mono[title.str()]->Fill( pullGPY_rs_mono );
776  title.str("");
777  title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
778  hPullGP_Z_rs_mono[title.str()]->Fill( pullGPZ_rs_mono );
779 
780  double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x())/sqrt(monoTsosGPEr.cxx()+monoRhitGPEr.cxx());
781  double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y())/sqrt(monoTsosGPEr.cyy()+monoRhitGPEr.cyy());
782  double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z())/sqrt(monoTsosGPEr.czz()+monoRhitGPEr.czz());
783  //double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x());
784  //double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y());
785  //double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z());
786 
787  title.str("");
788  title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
789  hPullGP_X_tr_mono[title.str()]->Fill( pullGPX_tr_mono );
790  title.str("");
791  title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
792  hPullGP_Y_tr_mono[title.str()]->Fill( pullGPY_tr_mono );
793  title.str("");
794  title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
795  hPullGP_Z_tr_mono[title.str()]->Fill( pullGPZ_tr_mono );
796 
797  double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x())/sqrt(monoTsosGPEr.cxx());
798  double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y())/sqrt(monoTsosGPEr.cyy());
799  double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z())/sqrt(monoTsosGPEr.czz());
800  //double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x());
801  //double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y());
802  //double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z());
803 
804  title.str("");
805  title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
806  hPullGP_X_ts_mono[title.str()]->Fill( pullGPX_ts_mono );
807  title.str("");
808  title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
809  hPullGP_Y_ts_mono[title.str()]->Fill( pullGPY_ts_mono );
810  title.str("");
811  title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
812  hPullGP_Z_ts_mono[title.str()]->Fill( pullGPZ_ts_mono );
813 
814  double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x())/sqrt(monoTsosGMEr.cxx());
815  double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y())/sqrt(monoTsosGMEr.cyy());
816  double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z())/sqrt(monoTsosGMEr.czz());
817  //double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x());
818  //double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y());
819  //double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z());
820 
821  title.str("");
822  title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
823  hPullGM_X_ts_mono[title.str()]->Fill( pullGMX_ts_mono );
824  title.str("");
825  title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
826  hPullGM_Y_ts_mono[title.str()]->Fill( pullGMY_ts_mono );
827  title.str("");
828  title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
829  hPullGM_Z_ts_mono[title.str()]->Fill( pullGMZ_ts_mono );
830 
831  //stereo
832  LogTrace("TestTrackHits") << "STEREO HIT" ;
833  auto s = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->stereoHit();
834  CTTRHp tStereoHit = theBuilder->build(&s);
835  if (tStereoHit==0) continue;
836  vector<PSimHit> assStereoSimHits = hitAssociator.associateHit(*tStereoHit->hit());
837  if (assStereoSimHits.size()==0) continue;
838  const PSimHit sStereoHit = *(assStereoSimHits.begin());
839  const Surface * stereoSurf = &( tStereoHit->det()->surface() );
840  if (stereoSurf==0) continue;
841  TSOS stereoState = thePropagatorAnyDir->propagate(state,*stereoSurf);
842  if (stereoState.isValid()==0) continue;
843 
844  LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
845  GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
846  LocalPoint stereoShitLPos = sStereoHit.localPosition();
847  GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
848 
849  //LogTrace("TestTrackHits") << "assStereoSimHits.size()=" << assStereoSimHits.size() ;
850  //LogTrace("TestTrackHits") << "stereo shit=" << stereoShitGPos ;
851 
852  GlobalVector stereoTsosGMom = stereoState.globalMomentum();
853  GlobalError stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
854  GlobalPoint stereoTsosGPos = stereoState.globalPosition();
855  GlobalError stereoTsosGPEr = stereoState.cartesianError().position();
856 
857  GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
858  GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
859 
860  MeasurementExtractor meSt(stereoState);
861  double chi2stereo = computeChi2Increment(meSt,tStereoHit);
862 
863  title.str("");
864  title << "Chi2Increment_stereo_" << subdetId << "-" << layerId ;
865  hChi2Increment_stereo[title.str()]->Fill(chi2stereo);
866 
867  double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/sqrt(stereoRhitGPEr.cxx());
868  double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/sqrt(stereoRhitGPEr.cyy());
869  double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/sqrt(stereoRhitGPEr.czz());
870  // double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/*/sqrt(stereoRhitGPEr.cxx())*/;
871  // double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/*/sqrt(stereoRhitGPEr.cyy())*/;
872  // double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/*/sqrt(stereoRhitGPEr.czz())*/;
873 
874  title.str("");
875  title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
876  hPullGP_X_rs_stereo[title.str()]->Fill( pullGPX_rs_stereo );
877  title.str("");
878  title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
879  hPullGP_Y_rs_stereo[title.str()]->Fill( pullGPY_rs_stereo );
880  title.str("");
881  title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
882  hPullGP_Z_rs_stereo[title.str()]->Fill( pullGPZ_rs_stereo );
883 
884  double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x())/sqrt(stereoTsosGPEr.cxx()+stereoRhitGPEr.cxx());
885  double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y())/sqrt(stereoTsosGPEr.cyy()+stereoRhitGPEr.cyy());
886  double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z())/sqrt(stereoTsosGPEr.czz()+stereoRhitGPEr.czz());
887  //double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x());
888  //double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y());
889  //double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z());
890 
891  title.str("");
892  title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
893  hPullGP_X_tr_stereo[title.str()]->Fill( pullGPX_tr_stereo );
894  title.str("");
895  title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
896  hPullGP_Y_tr_stereo[title.str()]->Fill( pullGPY_tr_stereo );
897  title.str("");
898  title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
899  hPullGP_Z_tr_stereo[title.str()]->Fill( pullGPZ_tr_stereo );
900 
901  double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x())/sqrt(stereoTsosGPEr.cxx());
902  double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y())/sqrt(stereoTsosGPEr.cyy());
903  double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z())/sqrt(stereoTsosGPEr.czz());
904  //double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x());
905  //double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y());
906  //double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z());
907 
908  title.str("");
909  title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
910  hPullGP_X_ts_stereo[title.str()]->Fill( pullGPX_ts_stereo );
911  title.str("");
912  title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
913  hPullGP_Y_ts_stereo[title.str()]->Fill( pullGPY_ts_stereo );
914  title.str("");
915  title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
916  hPullGP_Z_ts_stereo[title.str()]->Fill( pullGPZ_ts_stereo );
917 
918  double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x())/sqrt(stereoTsosGMEr.cxx());
919  double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y())/sqrt(stereoTsosGMEr.cyy());
920  double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z())/sqrt(stereoTsosGMEr.czz());
921  //double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x());
922  //double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y());
923  //double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z());
924 
925  title.str("");
926  title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
927  hPullGM_X_ts_stereo[title.str()]->Fill( pullGMX_ts_stereo );
928  title.str("");
929  title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
930  hPullGM_Y_ts_stereo[title.str()]->Fill( pullGMY_ts_stereo );
931  title.str("");
932  title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
933  hPullGM_Z_ts_stereo[title.str()]->Fill( pullGMZ_ts_stereo );
934  }
935  }
936  LogTrace("TestTrackHits") << "traj chi2=" << tchi2 ;
937  LogTrace("TestTrackHits") << "track chi2=" << track->chi2() ;
938  i++;
939  }
940  LogTrace("TestTrackHits") << "end of event: processd hits=" << evtHits ;
941 }
942 
944  //file->Write();
945  TDirectory * chi2d = file->mkdir("Chi2Increment");
946 
947  TDirectory * gp_ts = file->mkdir("GP_TSOS-SimHit");
948  TDirectory * gm_ts = file->mkdir("GM_TSOS-SimHit");
949  TDirectory * gp_tr = file->mkdir("GP_TSOS-RecHit");
950  TDirectory * gp_rs = file->mkdir("GP_RecHit-SimHit");
951 
952  TDirectory * gp_tsx = gp_ts->mkdir("X");
953  TDirectory * gp_tsy = gp_ts->mkdir("Y");
954  TDirectory * gp_tsz = gp_ts->mkdir("Z");
955  TDirectory * gm_tsx = gm_ts->mkdir("X");
956  TDirectory * gm_tsy = gm_ts->mkdir("Y");
957  TDirectory * gm_tsz = gm_ts->mkdir("Z");
958  TDirectory * gp_trx = gp_tr->mkdir("X");
959  TDirectory * gp_try = gp_tr->mkdir("Y");
960  TDirectory * gp_trz = gp_tr->mkdir("Z");
961  TDirectory * gp_rsx = gp_rs->mkdir("X");
962  TDirectory * gp_rsy = gp_rs->mkdir("Y");
963  TDirectory * gp_rsz = gp_rs->mkdir("Z");
964 
965  TDirectory * gp_tsx_mono = gp_ts->mkdir("MONOX");
966  TDirectory * gp_tsy_mono = gp_ts->mkdir("MONOY");
967  TDirectory * gp_tsz_mono = gp_ts->mkdir("MONOZ");
968  TDirectory * gm_tsx_mono = gm_ts->mkdir("MONOX");
969  TDirectory * gm_tsy_mono = gm_ts->mkdir("MONOY");
970  TDirectory * gm_tsz_mono = gm_ts->mkdir("MONOZ");
971  TDirectory * gp_trx_mono = gp_tr->mkdir("MONOX");
972  TDirectory * gp_try_mono = gp_tr->mkdir("MONOY");
973  TDirectory * gp_trz_mono = gp_tr->mkdir("MONOZ");
974  TDirectory * gp_rsx_mono = gp_rs->mkdir("MONOX");
975  TDirectory * gp_rsy_mono = gp_rs->mkdir("MONOY");
976  TDirectory * gp_rsz_mono = gp_rs->mkdir("MONOZ");
977 
978  TDirectory * gp_tsx_stereo = gp_ts->mkdir("STEREOX");
979  TDirectory * gp_tsy_stereo = gp_ts->mkdir("STEREOY");
980  TDirectory * gp_tsz_stereo = gp_ts->mkdir("STEREOZ");
981  TDirectory * gm_tsx_stereo = gm_ts->mkdir("STEREOX");
982  TDirectory * gm_tsy_stereo = gm_ts->mkdir("STEREOY");
983  TDirectory * gm_tsz_stereo = gm_ts->mkdir("STEREOZ");
984  TDirectory * gp_trx_stereo = gp_tr->mkdir("STEREOX");
985  TDirectory * gp_try_stereo = gp_tr->mkdir("STEREOY");
986  TDirectory * gp_trz_stereo = gp_tr->mkdir("STEREOZ");
987  TDirectory * gp_rsx_stereo = gp_rs->mkdir("STEREOX");
988  TDirectory * gp_rsy_stereo = gp_rs->mkdir("STEREOY");
989  TDirectory * gp_rsz_stereo = gp_rs->mkdir("STEREOZ");
990 
991  chi2d->cd();
992  hTotChi2Increment->Write();
993  hTotChi2GoodHit->Write();
994  hTotChi2BadHit->Write();
995  hTotChi2DeltaHit->Write();
996  hTotChi2NSharedHit->Write();
997  hTotChi2SharedHit->Write();
998  hProcess_vs_Chi2->Write();
999  hClsize_vs_Chi2->Write();
1000  hPixClsize_vs_Chi2->Write();
1001  hPrjClsize_vs_Chi2->Write();
1002  hSt1Clsize_vs_Chi2->Write();
1003  hSt2Clsize_vs_Chi2->Write();
1004  hGoodHit_vs_Chi2->Write();
1005  hClusterSize->Write();
1006  hPixClusterSize->Write();
1007  hPrjClusterSize->Write();
1008  hSt1ClusterSize->Write();
1009  hSt2ClusterSize->Write();
1010  hSimHitVecSize->Write();
1011  hPixSimHitVecSize->Write();
1012  hPrjSimHitVecSize->Write();
1013  hSt1SimHitVecSize->Write();
1014  hSt2SimHitVecSize->Write();
1015  goodbadmerged->Write();
1016  energyLossRatio->Write();
1017  mergedPull->Write();
1018  probXgood->Write();
1019  probXbad->Write();
1020  probXdelta->Write();
1021  probXshared->Write();
1022  probXnoshare->Write();
1023  probYgood->Write();
1024  probYbad->Write();
1025  probYdelta->Write();
1026  probYshared->Write();
1027  probYnoshare->Write();
1028  for (int i=0; i!=6; i++)
1029  for (int j=0; j!=9; j++){
1030  if (i==0 && j>2) break;
1031  if (i==1 && j>1) break;
1032  if (i==2 && j>3) break;
1033  if (i==3 && j>2) break;
1034  if (i==4 && j>5) break;
1035  if (i==5 && j>8) break;
1036  chi2d->cd();
1037  title.str("");
1038  title << "Chi2Increment_" << i+1 << "-" << j+1 ;
1039  hChi2Increment[title.str()]->Write();
1040  title.str("");
1041  title << "Chi2IncrementVsEta_" << i+1 << "-" << j+1 ;
1042  hChi2IncrementVsEta[title.str()]->Write();
1043  title.str("");
1044  title << "Chi2GoodHit_" << i+1 << "-" << j+1 ;
1045  hChi2GoodHit[title.str()]->Write();
1046  title.str("");
1047  title << "Chi2BadHit_" << i+1 << "-" << j+1 ;
1048  hChi2BadHit[title.str()]->Write();
1049  title.str("");
1050  title << "Chi2DeltaHit_" << i+1 << "-" << j+1 ;
1051  hChi2DeltaHit[title.str()]->Write();
1052  title.str("");
1053  title << "Chi2NSharedHit_" << i+1 << "-" << j+1 ;
1054  hChi2NSharedHit[title.str()]->Write();
1055  title.str("");
1056  title << "Chi2SharedHit_" << i+1 << "-" << j+1 ;
1057  hChi2SharedHit[title.str()]->Write();
1058 
1059  gp_ts->cd();
1060  gp_tsx->cd();
1061  title.str("");
1062  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
1063  hPullGP_X_ts[title.str()]->Write();
1064  gp_tsy->cd();
1065  title.str("");
1066  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
1067  hPullGP_Y_ts[title.str()]->Write();
1068  gp_tsz->cd();
1069  title.str("");
1070  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
1071  hPullGP_Z_ts[title.str()]->Write();
1072 
1073  gm_ts->cd();
1074  gm_tsx->cd();
1075  title.str("");
1076  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
1077  hPullGM_X_ts[title.str()]->Write();
1078  gm_tsy->cd();
1079  title.str("");
1080  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
1081  hPullGM_Y_ts[title.str()]->Write();
1082  gm_tsz->cd();
1083  title.str("");
1084  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
1085  hPullGM_Z_ts[title.str()]->Write();
1086 
1087  gp_tr->cd();
1088  gp_trx->cd();
1089  title.str("");
1090  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
1091  hPullGP_X_tr[title.str()]->Write();
1092  gp_try->cd();
1093  title.str("");
1094  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
1095  hPullGP_Y_tr[title.str()]->Write();
1096  gp_trz->cd();
1097  title.str("");
1098  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
1099  hPullGP_Z_tr[title.str()]->Write();
1100 
1101  gp_rs->cd();
1102  gp_rsx->cd();
1103  title.str("");
1104  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
1105  hPullGP_X_rs[title.str()]->Write();
1106  gp_rsy->cd();
1107  title.str("");
1108  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
1109  hPullGP_Y_rs[title.str()]->Write();
1110  gp_rsz->cd();
1111  title.str("");
1112  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
1113  hPullGP_Z_rs[title.str()]->Write();
1114 
1115  if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
1116  chi2d->cd();
1117  title.str("");
1118  title << "Chi2Increment_mono_" << i+1 << "-" << j+1 ;
1119  hChi2Increment_mono[title.str()]->Write();
1120  title.str("");
1121  title << "Chi2Increment_stereo_" << i+1 << "-" << j+1 ;
1122  hChi2Increment_stereo[title.str()]->Write();
1123  //mono
1124  gp_ts->cd();
1125  gp_tsx_mono->cd();
1126  title.str("");
1127  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
1128  hPullGP_X_ts_mono[title.str()]->Write();
1129  gp_tsy_mono->cd();
1130  title.str("");
1131  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
1132  hPullGP_Y_ts_mono[title.str()]->Write();
1133  gp_tsz_mono->cd();
1134  title.str("");
1135  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
1136  hPullGP_Z_ts_mono[title.str()]->Write();
1137 
1138  gm_ts->cd();
1139  gm_tsx_mono->cd();
1140  title.str("");
1141  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
1142  hPullGM_X_ts_mono[title.str()]->Write();
1143  gm_tsy_mono->cd();
1144  title.str("");
1145  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
1146  hPullGM_Y_ts_mono[title.str()]->Write();
1147  gm_tsz_mono->cd();
1148  title.str("");
1149  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
1150  hPullGM_Z_ts_mono[title.str()]->Write();
1151 
1152  gp_tr->cd();
1153  gp_trx_mono->cd();
1154  title.str("");
1155  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
1156  hPullGP_X_tr_mono[title.str()]->Write();
1157  gp_try_mono->cd();
1158  title.str("");
1159  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
1160  hPullGP_Y_tr_mono[title.str()]->Write();
1161  gp_trz_mono->cd();
1162  title.str("");
1163  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
1164  hPullGP_Z_tr_mono[title.str()]->Write();
1165 
1166  gp_rs->cd();
1167  gp_rsx_mono->cd();
1168  title.str("");
1169  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
1170  hPullGP_X_rs_mono[title.str()]->Write();
1171  gp_rsy_mono->cd();
1172  title.str("");
1173  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
1174  hPullGP_Y_rs_mono[title.str()]->Write();
1175  gp_rsz_mono->cd();
1176  title.str("");
1177  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
1178  hPullGP_Z_rs_mono[title.str()]->Write();
1179 
1180  //stereo
1181  gp_ts->cd();
1182  gp_tsx_stereo->cd();
1183  title.str("");
1184  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
1185  hPullGP_X_ts_stereo[title.str()]->Write();
1186  gp_tsy_stereo->cd();
1187  title.str("");
1188  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
1189  hPullGP_Y_ts_stereo[title.str()]->Write();
1190  gp_tsz_stereo->cd();
1191  title.str("");
1192  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
1193  hPullGP_Z_ts_stereo[title.str()]->Write();
1194 
1195  gm_ts->cd();
1196  gm_tsx_stereo->cd();
1197  title.str("");
1198  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
1199  hPullGM_X_ts_stereo[title.str()]->Write();
1200  gm_tsy_stereo->cd();
1201  title.str("");
1202  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
1203  hPullGM_Y_ts_stereo[title.str()]->Write();
1204  gm_tsz_stereo->cd();
1205  title.str("");
1206  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
1207  hPullGM_Z_ts_stereo[title.str()]->Write();
1208 
1209  gp_tr->cd();
1210  gp_trx_stereo->cd();
1211  title.str("");
1212  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
1213  hPullGP_X_tr_stereo[title.str()]->Write();
1214  gp_try_stereo->cd();
1215  title.str("");
1216  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
1217  hPullGP_Y_tr_stereo[title.str()]->Write();
1218  gp_trz_stereo->cd();
1219  title.str("");
1220  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
1221  hPullGP_Z_tr_stereo[title.str()]->Write();
1222 
1223  gp_rs->cd();
1224  gp_rsx_stereo->cd();
1225  title.str("");
1226  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
1227  hPullGP_X_rs_stereo[title.str()]->Write();
1228  gp_rsy_stereo->cd();
1229  title.str("");
1230  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
1231  hPullGP_Y_rs_stereo[title.str()]->Write();
1232  gp_rsz_stereo->cd();
1233  title.str("");
1234  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
1235  hPullGP_Z_rs_stereo[title.str()]->Write();
1236  }
1237  }
1238 
1239  file->Close();
1240 }
1241 
1242 
1243 //needed by to do the residual for matched hits
1244 //taken from SiStripTrackingRecHitsValid.cc
1245 std::pair<LocalPoint,LocalVector>
1246 TestTrackHits::projectHit( const PSimHit& hit, const StripGeomDetUnit* stripDet, const BoundPlane& plane)
1247 {
1248  const StripTopology& topol = stripDet->specificTopology();
1249  GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
1250  LocalPoint localHit = plane.toLocal(globalpos);
1251  //track direction
1252  LocalVector locdir=hit.localDirection();
1253  //rotate track in new frame
1254 
1255  GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
1256  LocalVector dir=plane.toLocal(globaldir);
1257  float scale = -localHit.z() / dir.z();
1258 
1259  LocalPoint projectedPos = localHit + scale*dir;
1260 
1261  float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
1262 
1263  LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame
1264 
1265  LocalVector localStripDir( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
1266 
1267  return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
1268 }
1269 
1270 template<unsigned int D>
1273  typedef typename AlgebraicROOTObject<D>::Vector VecD;
1274  typedef typename AlgebraicROOTObject<D,D>::SymMatrix SMatDD;
1275  VecD r = asSVector<D>(rhit->parameters()) - me.measuredParameters<D>(*rhit);
1276 
1277  SMatDD R = asSMatrix<D>(rhit->parametersError()) + me.measuredError<D>(*rhit);
1278  R.Invert();
1279  return ROOT::Math::Similarity(r,R) ;
1280 }
1281 
#define LogDebug(id)
std::string out
Definition: TestTrackHits.h:87
dbl * delta
Definition: mlp_gen.cc:36
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:106
TH1F * hPixSimHitVecSize
T getParameter(std::string const &) const
std::map< std::string, TH1F * > hPullGP_X_ts_mono
TH1F * probXnoshare
int i
Definition: DBlmapReader.cc:9
TH2F * hPrjClsize_vs_Chi2
TH1F * hSt2ClusterSize
std::map< std::string, TH1F * > hPullGM_X_ts_mono
std::map< std::string, TH1F * > hPullGM_Y_ts_stereo
ROOT::Math::SMatrix< double, D1, D1, ROOT::Math::MatRepSym< double, D1 > > SymMatrix
TrackerHitAssociator::Config trackerHitAssociatorConfig_
Definition: TestTrackHits.h:80
std::map< std::string, TH1F * > hPullGP_Y_ts_mono
virtual float stripAngle(float strip) const =0
TH1F * hTotChi2GoodHit
tuple pp
Definition: createTree.py:15
const_iterator end() const
last iterator over the map (read only)
LocalVector momentumAtEntry() const
The momentum of the track that produced the hit, at entry point.
Definition: PSimHit.h:47
std::map< std::string, TH1F * > hPullGP_Z_rs_mono
std::map< std::string, TH1F * > hPullGP_Z_rs_stereo
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
TH2F * hSt2Clsize_vs_Chi2
std::map< std::string, TH1F * > hPullGP_Y_rs_stereo
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::map< std::string, TH1F * > hPullGP_Z_ts_mono
AlgebraicVector measuredParameters(const TrackingRecHit &)
TH1F * probYnoshare
const_iterator find(const key_type &k) const
find element with specified reference key
std::map< std::string, TH1F * > hPullGP_Y_rs
std::map< std::string, TH1F * > hPullGP_Z_tr_mono
edm::Handle< std::vector< Trajectory > > trajCollectionHandle
Definition: TestTrackHits.h:94
const CartesianTrajectoryError cartesianError() const
T y() const
Definition: PV3DBase.h:63
std::map< std::string, TH1F * > hChi2Increment
std::map< std::string, TH1F * > hChi2BadHit
TH2F * hSt1Clsize_vs_Chi2
TH2F * hGoodHit_vs_Chi2
GlobalPoint globalPosition() const
std::map< std::string, TH1F * > hChi2SharedHit
std::map< std::string, TH1F * > hChi2Increment_mono
std::map< std::string, TH1F * > hPullGP_Y_ts
edm::Handle< reco::TrackToTrackingParticleAssociator > trackAssociator
Definition: TestTrackHits.h:93
std::map< std::string, TH1F * > hPullGP_X_rs
std::map< std::string, TH1F * > hPullGP_X_tr_stereo
TH1F * hSimHitVecSize
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:40
std::map< std::string, TH1F * > hPullGP_Y_tr
std::map< std::string, TH1F * > hChi2DeltaHit
TH2F * hPixClsize_vs_Chi2
virtual float strip(const LocalPoint &) const =0
std::stringstream title
TH1F * energyLossRatio
TH1F * probYshared
edm::ESHandle< Propagator > thePropagator
Definition: TestTrackHits.h:90
int iEvent
Definition: GenABIO.cc:230
edm::Handle< edm::View< reco::Track > > trackCollectionHandle
Definition: TestTrackHits.h:95
std::map< std::string, TH1F * > hPullGM_Z_ts_stereo
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:646
std::map< std::string, TH1F * > hPullGP_Z_ts
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
edm::ESHandle< MagneticField > theMF
Definition: TestTrackHits.h:89
Local3DPoint localPosition() const
Definition: PSimHit.h:44
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:544
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
std::map< std::string, TH1F * > hPullGP_X_ts
TH2F * hProcess_vs_Chi2
std::map< std::string, TH1F * > hPullGP_Z_rs
TestTrackHits(const edm::ParameterSet &)
T sqrt(T t)
Definition: SSEVec.h:18
edm::ESHandle< TrajectoryStateUpdator > theUpdator
Definition: TestTrackHits.h:92
TH1F * hClusterSize
std::map< std::string, TH1F * > hChi2Increment_stereo
AlgebraicSymMatrix measuredError(const TrackingRecHit &)
double pt() const
track transverse momentum
Definition: TrackBase.h:616
std::map< std::string, TH1F * > hChi2GoodHit
T z() const
Definition: PV3DBase.h:64
std::string srcName
Definition: TestTrackHits.h:84
std::map< std::string, TH1F * > hPullGM_Z_ts
TH1F * hTotChi2SharedHit
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::map< std::string, TH1F * > hPullGP_Z_tr
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
TH1F * goodbadmerged
TransientTrackingRecHit::ConstRecHitPointer CTTRHp
Definition: CkfDebugger.h:39
std::map< std::string, TH1F * > hPullGP_X_ts_stereo
float probabilityY() const
Definition: SiPixelRecHit.h:90
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:815
virtual void endJob()
virtual void beginRun(edm::Run &run, const edm::EventSetup &)
TH1F * hPrjSimHitVecSize
std::string builderName
Definition: TestTrackHits.h:83
std::string tpName
Definition: TestTrackHits.h:85
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:418
TH1F * hSt2SimHitVecSize
#define LogTrace(id)
const AlgebraicSymMatrix66 & matrix() const
edm::ESHandle< TransientTrackingRecHitBuilder > theBuilder
Definition: TestTrackHits.h:91
LocalVector localDirection() const
Obsolete. Same as momentumAtEntry().unit(), for backward compatibility.
Definition: PSimHit.h:52
std::vector< SimTrack >::const_iterator g4t_iterator
edm::Handle< TrackingParticleCollection > trackingParticleCollectionHandle
Definition: TestTrackHits.h:97
std::map< std::string, TH1F * > hPullGP_Y_rs_mono
Definition: DetId.h:18
TrajectoryStateOnSurface TSOS
Definition: TestHits.cc:19
ROOT::Math::SVector< double, D1 > Vector
TH1F * hTotChi2NSharedHit
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:150
std::map< std::string, TH1F * > hPullGM_Y_ts_mono
TH2F * hClsize_vs_Chi2
const GlobalError position() const
Position error submatrix.
TH1F * hTotChi2Increment
std::map< std::string, TH1F * > hPullGP_X_rs_mono
virtual void analyze(const edm::Event &, const edm::EventSetup &)
std::map< std::string, TH1F * > hPullGP_Y_tr_stereo
std::map< std::string, TH1F * > hPullGP_Z_tr_stereo
const T & get() const
Definition: EventSetup.h:56
T const * product() const
Definition: ESHandle.h:86
TH1F * hSt1SimHitVecSize
std::map< std::string, TH1F * > hPullGP_X_tr
TH1F * hSt1ClusterSize
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:53
std::map< std::string, TH1F * > hPullGP_X_tr_mono
unsigned short processType() const
Definition: PSimHit.h:118
edm::ESHandle< TrackerGeometry > theG
Definition: TestTrackHits.h:88
TH1F * probXshared
std::map< std::string, TH1F * > hChi2NSharedHit
std::vector< SimHitIdpr > associateHitId(const TrackingRecHit &thit) const
GlobalVector globalMomentum() const
std::pair< LocalPoint, LocalVector > projectHit(const PSimHit &, const StripGeomDetUnit *, const BoundPlane &)
std::map< std::string, TH1F * > hPullGM_X_ts
std::map< std::string, TH1F * > hPullGP_Z_ts_stereo
std::map< std::string, TH1F * > hPullGP_X_rs_stereo
std::map< std::string, TH1F * > hPullGP_Y_ts_stereo
edm::Handle< TrajTrackAssociationCollection > trajTrackAssociationCollectionHandle
Definition: TestTrackHits.h:96
std::vector< PSimHit > associateHit(const TrackingRecHit &thit) const
TH1F * hPrjClusterSize
TH1F * hTotChi2BadHit
dbl *** dir
Definition: mlp_gen.cc:35
std::map< std::string, TH1F * > hPullGM_Z_ts_mono
TH1F * hPixClusterSize
float probabilityX() const
Definition: SiPixelRecHit.h:87
double computeChi2Increment(MeasurementExtractor, TransientTrackingRecHit::ConstRecHitPointer)
std::map< std::string, TH2F * > hChi2IncrementVsEta
Detector det() const
get the detector field from this detid
Definition: DetId.h:35
T x() const
Definition: PV3DBase.h:62
const PositionType & position() const
TH1F * hTotChi2DeltaHit
std::string propagatorName
Definition: TestTrackHits.h:82
Definition: Run.h:43
std::map< std::string, TH1F * > hPullGM_X_ts_stereo
std::string updatorName
Definition: TestTrackHits.h:86
std::map< std::string, TH1F * > hPullGM_Y_ts
std::map< std::string, TH1F * > hPullGP_Y_tr_mono
Our base class.
Definition: SiPixelRecHit.h:23