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.
8 #include <TDirectory.h>
9 
13 
16 using namespace std;
17 using namespace edm;
18 
20  conf_(iConfig) {
21  LogTrace("TestTrackHits") << conf_;
22  propagatorName = conf_.getParameter<std::string>("Propagator");
23  builderName = conf_.getParameter<std::string>("TTRHBuilder");
25  tpName = conf_.getParameter<std::string>("tpname");
28 // ParameterSet cuts = conf_.getParameter<ParameterSet>("RecoTracksCuts");
29 // selectRecoTracks = RecoTrackSelector(cuts.getParameter<double>("ptMin"),
30 // cuts.getParameter<double>("minRapidity"),
31 // cuts.getParameter<double>("maxRapidity"),
32 // cuts.getParameter<double>("tip"),
33 // cuts.getParameter<double>("lip"),
34 // cuts.getParameter<int>("minHit"),
35 // cuts.getParameter<double>("maxChi2"));
36 }
37 
39 
41 {
42  iSetup.get<TrackerDigiGeometryRecord>().get(theG);
43  iSetup.get<IdealMagneticFieldRecord>().get(theMF);
47  iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits",trackAssociator);
48 
49  file = new TFile(out.c_str(),"recreate");
50  for (int i=0; i!=6; i++)
51  for (int j=0; j!=9; j++){
52  if (i==0 && j>2) break;
53  if (i==1 && j>1) break;
54  if (i==2 && j>3) break;
55  if (i==3 && j>2) break;
56  if (i==4 && j>5) break;
57  if (i==5 && j>8) break;
58 
59  title.str("");
60  title << "Chi2Increment_" << i+1 << "-" << j+1 ;
61  hChi2Increment[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
62  title.str("");
63  title << "Chi2IncrementVsEta_" << i+1 << "-" << j+1 ;
64  hChi2IncrementVsEta[title.str()] = new TH2F(title.str().c_str(),title.str().c_str(),50,-2.5,2.5,1000,0,100);
65  title.str("");
66  title << "Chi2GoodHit_" << i+1 << "-" << j+1 ;
67  hChi2GoodHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
68  title.str("");
69  title << "Chi2BadHit_" << i+1 << "-" << j+1 ;
70  hChi2BadHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
71  title.str("");
72  title << "Chi2DeltaHit_" << i+1 << "-" << j+1 ;
73  hChi2DeltaHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
74  title.str("");
75  title << "Chi2NSharedHit_" << i+1 << "-" << j+1 ;
76  hChi2NSharedHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
77  title.str("");
78  title << "Chi2SharedHit_" << i+1 << "-" << j+1 ;
79  hChi2SharedHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
80 
81  title.str("");
82  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
83  hPullGP_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
84  title.str("");
85  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
86  hPullGP_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
87  title.str("");
88  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
89  hPullGP_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
90 
91  title.str("");
92  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
93  hPullGM_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
94  title.str("");
95  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
96  hPullGM_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
97  title.str("");
98  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
99  hPullGM_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
100 
101  title.str("");
102  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
103  hPullGP_X_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
104  title.str("");
105  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
106  hPullGP_Y_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
107  title.str("");
108  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
109  hPullGP_Z_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
110 
111  title.str("");
112  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
113  hPullGP_X_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
114  title.str("");
115  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
116  hPullGP_Y_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
117  title.str("");
118  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
119  hPullGP_Z_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
120 
121  if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
122  //mono
123  title.str("");
124  title << "Chi2Increment_mono_" << i+1 << "-" << j+1 ;
125  hChi2Increment_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
126 
127  title.str("");
128  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
129  hPullGP_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
130  title.str("");
131  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
132  hPullGP_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
133  title.str("");
134  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
135  hPullGP_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
136 
137  title.str("");
138  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
139  hPullGM_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
140  title.str("");
141  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
142  hPullGM_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
143  title.str("");
144  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
145  hPullGM_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
146 
147  title.str("");
148  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
149  hPullGP_X_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
150  title.str("");
151  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
152  hPullGP_Y_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
153  title.str("");
154  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
155  hPullGP_Z_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
156 
157  title.str("");
158  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
159  hPullGP_X_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
160  title.str("");
161  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
162  hPullGP_Y_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
163  title.str("");
164  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
165  hPullGP_Z_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
166 
167  //stereo
168  title.str("");
169  title << "Chi2Increment_stereo_" << i+1 << "-" << j+1 ;
170  hChi2Increment_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
171 
172  title.str("");
173  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
174  hPullGP_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
175  title.str("");
176  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
177  hPullGP_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
178  title.str("");
179  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
180  hPullGP_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
181 
182  title.str("");
183  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
184  hPullGM_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
185  title.str("");
186  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
187  hPullGM_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
188  title.str("");
189  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
190  hPullGM_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
191 
192  title.str("");
193  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
194  hPullGP_X_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
195  title.str("");
196  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
197  hPullGP_Y_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
198  title.str("");
199  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
200  hPullGP_Z_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
201 
202  title.str("");
203  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
204  hPullGP_X_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
205  title.str("");
206  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
207  hPullGP_Y_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
208  title.str("");
209  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
210  hPullGP_Z_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
211  }
212  }
213  hTotChi2Increment = new TH1F("TotChi2Increment","TotChi2Increment",1000,0,100);
214  hTotChi2GoodHit = new TH1F("TotChi2GoodHit","TotChi2GoodHit",1000,0,100);
215  hTotChi2BadHit = new TH1F("TotChi2BadHit","TotChi2BadHit",1000,0,100);
216  hTotChi2DeltaHit = new TH1F("TotChi2DeltaHit","TotChi2DeltaHit",1000,0,100);
217  hTotChi2NSharedHit = new TH1F("TotChi2NSharedHit","TotChi2NSharedHit",1000,0,100);
218  hTotChi2SharedHit = new TH1F("TotChi2SharedHit","TotChi2SharedHit",1000,0,100);
219  hProcess_vs_Chi2 = new TH2F("Process_vs_Chi2","Process_vs_Chi2",1000,0,100,17,-0.5,16.5);
220  hClsize_vs_Chi2 = new TH2F("Clsize_vs_Chi2","Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
221  hPixClsize_vs_Chi2= new TH2F("PixClsize_vs_Chi2","PixClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
222  hPrjClsize_vs_Chi2= new TH2F("PrjClsize_vs_Chi2","PrjClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
223  hSt1Clsize_vs_Chi2= new TH2F("St1Clsize_vs_Chi2","St1Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
224  hSt2Clsize_vs_Chi2= new TH2F("St2Clsize_vs_Chi2","St2Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
225  hGoodHit_vs_Chi2 = new TH2F("GoodHit_vs_Chi2","GoodHit_vs_Chi2",10000,0,1000,2,-0.5,1.5);
226  hClusterSize = new TH1F("ClusterSize","ClusterSize",40,-0.5,39.5);
227  hPixClusterSize = new TH1F("PixClusterSize","PixClusterSize",40,-0.5,39.5);
228  hPrjClusterSize = new TH1F("PrjClusterSize","PrjClusterSize",40,-0.5,39.5);
229  hSt1ClusterSize = new TH1F("St1ClusterSize","St1ClusterSize",40,-0.5,39.5);
230  hSt2ClusterSize = new TH1F("St2ClusterSize","St2ClusterSize",40,-0.5,39.5);
231  hSimHitVecSize = new TH1F("hSimHitVecSize","hSimHitVecSize",40,-0.5,39.5);
232  hPixSimHitVecSize = new TH1F("PixSimHitVecSize","PixSimHitVecSize",40,-0.5,39.5);
233  hPrjSimHitVecSize = new TH1F("PrjSimHitVecSize","PrjSimHitVecSize",40,-0.5,39.5);
234  hSt1SimHitVecSize = new TH1F("St1SimHitVecSize","St1SimHitVecSize",40,-0.5,39.5);
235  hSt2SimHitVecSize = new TH1F("St2SimHitVecSize","St2SimHitVecSize",40,-0.5,39.5);
236  goodbadmerged = new TH1F("goodbadmerged","goodbadmerged",5,0.5,5.5);
237  energyLossRatio = new TH1F("energyLossRatio","energyLossRatio",100,0,1);
238  mergedPull = new TH1F("mergedPull","mergedPull",200,0,20);
239  probXgood = new TH1F("probXgood","probXgood",110,0,1.1);
240  probXbad = new TH1F("probXbad","probXbad",110,0,1.1);
241  probXdelta = new TH1F("probXdelta","probXdelta",110,0,1.1);
242  probXshared = new TH1F("probXshared","probXshared",110,0,1.1);
243  probXnoshare = new TH1F("probXnoshare","probXnoshare",110,0,1.1);
244  probYgood = new TH1F("probYgood","probYgood",110,0,1.1);
245  probYbad = new TH1F("probYbad","probYbad",110,0,1.1);
246  probYdelta = new TH1F("probYdelta","probYdelta",110,0,1.1);
247  probYshared = new TH1F("probYshared","probYshared",110,0,1.1);
248  probYnoshare = new TH1F("probYnoshare","probYnoshare",110,0,1.1);
249 }
250 
252 {
253  //Retrieve tracker topology from geometry
255  iSetup.get<IdealGeometryRecord>().get(tTopo);
256 
257 
258  LogDebug("TestTrackHits") << "new event" ;
264  iEvent.getByLabel("offlineBeamSpot",beamSpot);
265 
266  LogTrace("TestTrackHits") << "Tr collection size=" << trackCollectionHandle->size();
267  LogTrace("TestTrackHits") << "TP collection size=" << trackingParticleCollectionHandle->size();
268 
269  hitAssociator = new TrackerHitAssociator(iEvent);
270 
271  reco::RecoToSimCollection recSimColl=trackAssociator->associateRecoToSim(trackCollectionHandle,
273  &iEvent,&iSetup);
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  delete hitAssociator;
942 }
943 
945  //file->Write();
946  TDirectory * chi2d = file->mkdir("Chi2Increment");
947 
948  TDirectory * gp_ts = file->mkdir("GP_TSOS-SimHit");
949  TDirectory * gm_ts = file->mkdir("GM_TSOS-SimHit");
950  TDirectory * gp_tr = file->mkdir("GP_TSOS-RecHit");
951  TDirectory * gp_rs = file->mkdir("GP_RecHit-SimHit");
952 
953  TDirectory * gp_tsx = gp_ts->mkdir("X");
954  TDirectory * gp_tsy = gp_ts->mkdir("Y");
955  TDirectory * gp_tsz = gp_ts->mkdir("Z");
956  TDirectory * gm_tsx = gm_ts->mkdir("X");
957  TDirectory * gm_tsy = gm_ts->mkdir("Y");
958  TDirectory * gm_tsz = gm_ts->mkdir("Z");
959  TDirectory * gp_trx = gp_tr->mkdir("X");
960  TDirectory * gp_try = gp_tr->mkdir("Y");
961  TDirectory * gp_trz = gp_tr->mkdir("Z");
962  TDirectory * gp_rsx = gp_rs->mkdir("X");
963  TDirectory * gp_rsy = gp_rs->mkdir("Y");
964  TDirectory * gp_rsz = gp_rs->mkdir("Z");
965 
966  TDirectory * gp_tsx_mono = gp_ts->mkdir("MONOX");
967  TDirectory * gp_tsy_mono = gp_ts->mkdir("MONOY");
968  TDirectory * gp_tsz_mono = gp_ts->mkdir("MONOZ");
969  TDirectory * gm_tsx_mono = gm_ts->mkdir("MONOX");
970  TDirectory * gm_tsy_mono = gm_ts->mkdir("MONOY");
971  TDirectory * gm_tsz_mono = gm_ts->mkdir("MONOZ");
972  TDirectory * gp_trx_mono = gp_tr->mkdir("MONOX");
973  TDirectory * gp_try_mono = gp_tr->mkdir("MONOY");
974  TDirectory * gp_trz_mono = gp_tr->mkdir("MONOZ");
975  TDirectory * gp_rsx_mono = gp_rs->mkdir("MONOX");
976  TDirectory * gp_rsy_mono = gp_rs->mkdir("MONOY");
977  TDirectory * gp_rsz_mono = gp_rs->mkdir("MONOZ");
978 
979  TDirectory * gp_tsx_stereo = gp_ts->mkdir("STEREOX");
980  TDirectory * gp_tsy_stereo = gp_ts->mkdir("STEREOY");
981  TDirectory * gp_tsz_stereo = gp_ts->mkdir("STEREOZ");
982  TDirectory * gm_tsx_stereo = gm_ts->mkdir("STEREOX");
983  TDirectory * gm_tsy_stereo = gm_ts->mkdir("STEREOY");
984  TDirectory * gm_tsz_stereo = gm_ts->mkdir("STEREOZ");
985  TDirectory * gp_trx_stereo = gp_tr->mkdir("STEREOX");
986  TDirectory * gp_try_stereo = gp_tr->mkdir("STEREOY");
987  TDirectory * gp_trz_stereo = gp_tr->mkdir("STEREOZ");
988  TDirectory * gp_rsx_stereo = gp_rs->mkdir("STEREOX");
989  TDirectory * gp_rsy_stereo = gp_rs->mkdir("STEREOY");
990  TDirectory * gp_rsz_stereo = gp_rs->mkdir("STEREOZ");
991 
992  chi2d->cd();
993  hTotChi2Increment->Write();
994  hTotChi2GoodHit->Write();
995  hTotChi2BadHit->Write();
996  hTotChi2DeltaHit->Write();
997  hTotChi2NSharedHit->Write();
998  hTotChi2SharedHit->Write();
999  hProcess_vs_Chi2->Write();
1000  hClsize_vs_Chi2->Write();
1001  hPixClsize_vs_Chi2->Write();
1002  hPrjClsize_vs_Chi2->Write();
1003  hSt1Clsize_vs_Chi2->Write();
1004  hSt2Clsize_vs_Chi2->Write();
1005  hGoodHit_vs_Chi2->Write();
1006  hClusterSize->Write();
1007  hPixClusterSize->Write();
1008  hPrjClusterSize->Write();
1009  hSt1ClusterSize->Write();
1010  hSt2ClusterSize->Write();
1011  hSimHitVecSize->Write();
1012  hPixSimHitVecSize->Write();
1013  hPrjSimHitVecSize->Write();
1014  hSt1SimHitVecSize->Write();
1015  hSt2SimHitVecSize->Write();
1016  goodbadmerged->Write();
1017  energyLossRatio->Write();
1018  mergedPull->Write();
1019  probXgood->Write();
1020  probXbad->Write();
1021  probXdelta->Write();
1022  probXshared->Write();
1023  probXnoshare->Write();
1024  probYgood->Write();
1025  probYbad->Write();
1026  probYdelta->Write();
1027  probYshared->Write();
1028  probYnoshare->Write();
1029  for (int i=0; i!=6; i++)
1030  for (int j=0; j!=9; j++){
1031  if (i==0 && j>2) break;
1032  if (i==1 && j>1) break;
1033  if (i==2 && j>3) break;
1034  if (i==3 && j>2) break;
1035  if (i==4 && j>5) break;
1036  if (i==5 && j>8) break;
1037  chi2d->cd();
1038  title.str("");
1039  title << "Chi2Increment_" << i+1 << "-" << j+1 ;
1040  hChi2Increment[title.str()]->Write();
1041  title.str("");
1042  title << "Chi2IncrementVsEta_" << i+1 << "-" << j+1 ;
1043  hChi2IncrementVsEta[title.str()]->Write();
1044  title.str("");
1045  title << "Chi2GoodHit_" << i+1 << "-" << j+1 ;
1046  hChi2GoodHit[title.str()]->Write();
1047  title.str("");
1048  title << "Chi2BadHit_" << i+1 << "-" << j+1 ;
1049  hChi2BadHit[title.str()]->Write();
1050  title.str("");
1051  title << "Chi2DeltaHit_" << i+1 << "-" << j+1 ;
1052  hChi2DeltaHit[title.str()]->Write();
1053  title.str("");
1054  title << "Chi2NSharedHit_" << i+1 << "-" << j+1 ;
1055  hChi2NSharedHit[title.str()]->Write();
1056  title.str("");
1057  title << "Chi2SharedHit_" << i+1 << "-" << j+1 ;
1058  hChi2SharedHit[title.str()]->Write();
1059 
1060  gp_ts->cd();
1061  gp_tsx->cd();
1062  title.str("");
1063  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
1064  hPullGP_X_ts[title.str()]->Write();
1065  gp_tsy->cd();
1066  title.str("");
1067  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
1068  hPullGP_Y_ts[title.str()]->Write();
1069  gp_tsz->cd();
1070  title.str("");
1071  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
1072  hPullGP_Z_ts[title.str()]->Write();
1073 
1074  gm_ts->cd();
1075  gm_tsx->cd();
1076  title.str("");
1077  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
1078  hPullGM_X_ts[title.str()]->Write();
1079  gm_tsy->cd();
1080  title.str("");
1081  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
1082  hPullGM_Y_ts[title.str()]->Write();
1083  gm_tsz->cd();
1084  title.str("");
1085  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
1086  hPullGM_Z_ts[title.str()]->Write();
1087 
1088  gp_tr->cd();
1089  gp_trx->cd();
1090  title.str("");
1091  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
1092  hPullGP_X_tr[title.str()]->Write();
1093  gp_try->cd();
1094  title.str("");
1095  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
1096  hPullGP_Y_tr[title.str()]->Write();
1097  gp_trz->cd();
1098  title.str("");
1099  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
1100  hPullGP_Z_tr[title.str()]->Write();
1101 
1102  gp_rs->cd();
1103  gp_rsx->cd();
1104  title.str("");
1105  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
1106  hPullGP_X_rs[title.str()]->Write();
1107  gp_rsy->cd();
1108  title.str("");
1109  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
1110  hPullGP_Y_rs[title.str()]->Write();
1111  gp_rsz->cd();
1112  title.str("");
1113  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
1114  hPullGP_Z_rs[title.str()]->Write();
1115 
1116  if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
1117  chi2d->cd();
1118  title.str("");
1119  title << "Chi2Increment_mono_" << i+1 << "-" << j+1 ;
1120  hChi2Increment_mono[title.str()]->Write();
1121  title.str("");
1122  title << "Chi2Increment_stereo_" << i+1 << "-" << j+1 ;
1123  hChi2Increment_stereo[title.str()]->Write();
1124  //mono
1125  gp_ts->cd();
1126  gp_tsx_mono->cd();
1127  title.str("");
1128  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
1129  hPullGP_X_ts_mono[title.str()]->Write();
1130  gp_tsy_mono->cd();
1131  title.str("");
1132  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
1133  hPullGP_Y_ts_mono[title.str()]->Write();
1134  gp_tsz_mono->cd();
1135  title.str("");
1136  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
1137  hPullGP_Z_ts_mono[title.str()]->Write();
1138 
1139  gm_ts->cd();
1140  gm_tsx_mono->cd();
1141  title.str("");
1142  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
1143  hPullGM_X_ts_mono[title.str()]->Write();
1144  gm_tsy_mono->cd();
1145  title.str("");
1146  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
1147  hPullGM_Y_ts_mono[title.str()]->Write();
1148  gm_tsz_mono->cd();
1149  title.str("");
1150  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
1151  hPullGM_Z_ts_mono[title.str()]->Write();
1152 
1153  gp_tr->cd();
1154  gp_trx_mono->cd();
1155  title.str("");
1156  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
1157  hPullGP_X_tr_mono[title.str()]->Write();
1158  gp_try_mono->cd();
1159  title.str("");
1160  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
1161  hPullGP_Y_tr_mono[title.str()]->Write();
1162  gp_trz_mono->cd();
1163  title.str("");
1164  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
1165  hPullGP_Z_tr_mono[title.str()]->Write();
1166 
1167  gp_rs->cd();
1168  gp_rsx_mono->cd();
1169  title.str("");
1170  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
1171  hPullGP_X_rs_mono[title.str()]->Write();
1172  gp_rsy_mono->cd();
1173  title.str("");
1174  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
1175  hPullGP_Y_rs_mono[title.str()]->Write();
1176  gp_rsz_mono->cd();
1177  title.str("");
1178  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
1179  hPullGP_Z_rs_mono[title.str()]->Write();
1180 
1181  //stereo
1182  gp_ts->cd();
1183  gp_tsx_stereo->cd();
1184  title.str("");
1185  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
1186  hPullGP_X_ts_stereo[title.str()]->Write();
1187  gp_tsy_stereo->cd();
1188  title.str("");
1189  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
1190  hPullGP_Y_ts_stereo[title.str()]->Write();
1191  gp_tsz_stereo->cd();
1192  title.str("");
1193  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
1194  hPullGP_Z_ts_stereo[title.str()]->Write();
1195 
1196  gm_ts->cd();
1197  gm_tsx_stereo->cd();
1198  title.str("");
1199  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
1200  hPullGM_X_ts_stereo[title.str()]->Write();
1201  gm_tsy_stereo->cd();
1202  title.str("");
1203  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
1204  hPullGM_Y_ts_stereo[title.str()]->Write();
1205  gm_tsz_stereo->cd();
1206  title.str("");
1207  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
1208  hPullGM_Z_ts_stereo[title.str()]->Write();
1209 
1210  gp_tr->cd();
1211  gp_trx_stereo->cd();
1212  title.str("");
1213  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
1214  hPullGP_X_tr_stereo[title.str()]->Write();
1215  gp_try_stereo->cd();
1216  title.str("");
1217  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
1218  hPullGP_Y_tr_stereo[title.str()]->Write();
1219  gp_trz_stereo->cd();
1220  title.str("");
1221  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
1222  hPullGP_Z_tr_stereo[title.str()]->Write();
1223 
1224  gp_rs->cd();
1225  gp_rsx_stereo->cd();
1226  title.str("");
1227  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
1228  hPullGP_X_rs_stereo[title.str()]->Write();
1229  gp_rsy_stereo->cd();
1230  title.str("");
1231  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
1232  hPullGP_Y_rs_stereo[title.str()]->Write();
1233  gp_rsz_stereo->cd();
1234  title.str("");
1235  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
1236  hPullGP_Z_rs_stereo[title.str()]->Write();
1237  }
1238  }
1239 
1240  file->Close();
1241 }
1242 
1243 
1244 //needed by to do the residual for matched hits
1245 //taken from SiStripTrackingRecHitsValid.cc
1246 std::pair<LocalPoint,LocalVector>
1247 TestTrackHits::projectHit( const PSimHit& hit, const StripGeomDetUnit* stripDet, const BoundPlane& plane)
1248 {
1249  const StripTopology& topol = stripDet->specificTopology();
1250  GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
1251  LocalPoint localHit = plane.toLocal(globalpos);
1252  //track direction
1253  LocalVector locdir=hit.localDirection();
1254  //rotate track in new frame
1255 
1256  GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
1257  LocalVector dir=plane.toLocal(globaldir);
1258  float scale = -localHit.z() / dir.z();
1259 
1260  LocalPoint projectedPos = localHit + scale*dir;
1261 
1262  float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
1263 
1264  LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame
1265 
1266  LocalVector localStripDir( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
1267 
1268  return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
1269 }
1270 
1271 template<unsigned int D>
1274  typedef typename AlgebraicROOTObject<D>::Vector VecD;
1275  typedef typename AlgebraicROOTObject<D,D>::SymMatrix SMatDD;
1276  VecD r = asSVector<D>(rhit->parameters()) - me.measuredParameters<D>(*rhit);
1277 
1278  SMatDD R = asSMatrix<D>(rhit->parametersError()) + me.measuredError<D>(*rhit);
1279  R.Invert();
1280  return ROOT::Math::Similarity(r,R) ;
1281 }
1282 
#define LogDebug(id)
std::string out
Definition: TestTrackHits.h:88
dbl * delta
Definition: mlp_gen.cc:36
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:114
TH1F * hPixSimHitVecSize
T getParameter(std::string const &) const
std::map< std::string, TH1F * > hPullGP_X_ts_mono
TH1F * probXnoshare
virtual FreeTrajectoryState propagate(const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const final
Definition: Propagator.h:119
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
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
TrackerHitAssociator * hitAssociator
Definition: TestTrackHits.h:81
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:95
const CartesianTrajectoryError cartesianError() const
T y() const
Definition: PV3DBase.h:63
const edm::ParameterSet conf_
Definition: TestTrackHits.h:80
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
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:91
int iEvent
Definition: GenABIO.cc:230
edm::Handle< edm::View< reco::Track > > trackCollectionHandle
Definition: TestTrackHits.h:96
std::map< std::string, TH1F * > hPullGM_Z_ts_stereo
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:699
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:90
Local3DPoint localPosition() const
Definition: PSimHit.h:44
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:597
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:48
edm::ESHandle< TrajectoryStateUpdator > theUpdator
Definition: TestTrackHits.h:93
TH1F * hClusterSize
std::map< std::string, TH1F * > hChi2Increment_stereo
AlgebraicSymMatrix measuredError(const TrackingRecHit &)
double pt() const
track transverse momentum
Definition: TrackBase.h:669
std::map< std::string, TH1F * > hChi2GoodHit
T z() const
Definition: PV3DBase.h:64
std::string srcName
Definition: TestTrackHits.h:85
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:868
virtual void endJob()
virtual void beginRun(edm::Run &run, const edm::EventSetup &)
TH1F * hPrjSimHitVecSize
std::string builderName
Definition: TestTrackHits.h:84
bool first
Definition: L1TdeRCT.cc:75
std::string tpName
Definition: TestTrackHits.h:86
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:402
TH1F * hSt2SimHitVecSize
#define LogTrace(id)
const AlgebraicSymMatrix66 & matrix() const
edm::ESHandle< TransientTrackingRecHitBuilder > theBuilder
Definition: TestTrackHits.h:92
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:98
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:55
edm::ESHandle< TrackAssociatorBase > trackAssociator
Definition: TestTrackHits.h:94
T const * product() const
Definition: ESHandle.h:86
TH1F * hSt1SimHitVecSize
std::map< std::string, TH1F * > hPullGP_X_tr
TH1F * hSt1ClusterSize
std::map< std::string, TH1F * > hPullGP_X_tr_mono
unsigned short processType() const
Definition: PSimHit.h:118
edm::ESHandle< TrackerGeometry > theG
Definition: TestTrackHits.h:89
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:97
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:83
Definition: Run.h:41
std::map< std::string, TH1F * > hPullGM_X_ts_stereo
std::string updatorName
Definition: TestTrackHits.h:87
std::map< std::string, TH1F * > hPullGM_Y_ts
std::map< std::string, TH1F * > hPullGP_Y_tr_mono
Our base class.
Definition: SiPixelRecHit.h:23