CMS 3D CMS Logo

TestSmoothHits.cc
Go to the documentation of this file.
1 #include "TestSmoothHits.h"
2 
11 #include <TDirectory.h>
13 
17 
20 using namespace std;
21 using namespace edm;
22 
24  trackerHitAssociatorConfig_(consumesCollector()){
25  LogTrace("TestSmoothHits") << iConfig<< std::endl;
26  propagatorName = iConfig.getParameter<std::string>("Propagator");
27  builderName = iConfig.getParameter<std::string>("TTRHBuilder");
28  srcName = iConfig.getParameter<std::string>("src");
29  fname = iConfig.getParameter<std::string>("Fitter");
30  sname = iConfig.getParameter<std::string>("Smoother");
31  mineta = iConfig.getParameter<double>("mineta");
32  maxeta = iConfig.getParameter<double>("maxeta");
33 }
34 
36 
38 {
39 
40  iSetup.get<TrackerDigiGeometryRecord>().get(theG);
41  iSetup.get<IdealMagneticFieldRecord>().get(theMF);
44  iSetup.get<TrajectoryFitter::Record>().get(fname, fit);
45  iSetup.get<TrajectoryFitter::Record>().get(sname, smooth);
46 
47  file = new TFile("testSmoothHits.root","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  title.str("");
57  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
58  hPullGP_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
59  title.str("");
60  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
61  hPullGP_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
62  title.str("");
63  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
64  hPullGP_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
65  title.str("");
66  title << "Chi2Increment_" << i+1 << "-" << j+1;
67  hChi2Increment[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
68 
69  title.str("");
70  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
71  hPullGM_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
72  title.str("");
73  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
74  hPullGM_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
75  title.str("");
76  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
77  hPullGM_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
78 
79  title.str("");
80  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
81  hPullGP_X_tr[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 << "_tr";
84  hPullGP_Y_tr[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 << "_tr";
87  hPullGP_Z_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
88 
89  title.str("");
90  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
91  hPullGP_X_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
92  title.str("");
93  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
94  hPullGP_Y_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
95  title.str("");
96  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
97  hPullGP_Z_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
98 
99  if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
100  //mono
101  title.str("");
102  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
103  hPullGP_X_ts_mono[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 << "_ts_mono";
106  hPullGP_Y_ts_mono[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 << "_ts_mono";
109  hPullGP_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
110 
111  title.str("");
112  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
113  hPullGM_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
114  title.str("");
115  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
116  hPullGM_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
117  title.str("");
118  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
119  hPullGM_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
120 
121  title.str("");
122  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
123  hPullGP_X_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
124  title.str("");
125  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
126  hPullGP_Y_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
127  title.str("");
128  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
129  hPullGP_Z_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
130 
131  title.str("");
132  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
133  hPullGP_X_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
134  title.str("");
135  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
136  hPullGP_Y_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
137  title.str("");
138  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
139  hPullGP_Z_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
140 
141  //stereo
142  title.str("");
143  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
144  hPullGP_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
145  title.str("");
146  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
147  hPullGP_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
148  title.str("");
149  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
150  hPullGP_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
151 
152  title.str("");
153  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
154  hPullGM_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
155  title.str("");
156  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
157  hPullGM_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
158  title.str("");
159  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
160  hPullGM_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
161 
162  title.str("");
163  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
164  hPullGP_X_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
165  title.str("");
166  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
167  hPullGP_Y_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
168  title.str("");
169  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
170  hPullGP_Z_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
171 
172  title.str("");
173  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
174  hPullGP_X_rs_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 << "_rs_stereo";
177  hPullGP_Y_rs_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 << "_rs_stereo";
180  hPullGP_Z_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
181  }
182  }
183  hTotChi2Increment = new TH1F("TotChi2Increment","TotChi2Increment",1000,0,100);
184  hChi2_vs_Process = new TH2F("Chi2_vs_Process","Chi2_vs_Process",1000,0,100,17,-0.5,16.5);
185  hChi2_vs_clsize = new TH2F("Chi2_vs_clsize","Chi2_vs_clsize",1000,0,100,17,-0.5,16.5);
186 }
187 
188 
190 {
191  //Retrieve tracker topology from geometry
193  iSetup.get<TrackerTopologyRcd>().get(tTopo);
194 
195 
196  LogTrace("TestSmoothHits") << "new event" << std::endl;
197 
200 
202 
203  for (TrackCandidateCollection::const_iterator i=theTCCollection->begin(); i!=theTCCollection->end();i++){
204 
205  LogTrace("TestSmoothHits") << "new candidate" << std::endl;
206 
207  const TrackCandidate * theTC = &(*i);
209  const TrackCandidate::range& recHitVec=theTC->recHits();
210 
211  //convert PTrajectoryStateOnDet to TrajectoryStateOnSurface
212 
213 
214  DetId detId(state.detId());
215  TrajectoryStateOnSurface theTSOS=
217 
218  if (theTSOS.globalMomentum().eta()>maxeta || theTSOS.globalMomentum().eta()<mineta) continue;
219 
220  //convert the TrackingRecHit vector to a TransientTrackingRecHit vector
222 
224  i!=recHitVec.second; i++){
225  hits.push_back(theBuilder->build(&(*i) ));
226  }
227 
228  //call the fitter
229  std::vector<Trajectory> fitted = fit->fit(theTC->seed(), hits, theTSOS);
230  //call the smoother
231  std::vector<Trajectory> result;
232  for(std::vector<Trajectory>::iterator it = fitted.begin(); it != fitted.end(); it++) {
233  std::vector<Trajectory> smoothed = smooth->trajectories(*it);
234  result.insert(result.end(), smoothed.begin(), smoothed.end());
235  }
236  if (result.empty()) continue;
237  std::vector<TrajectoryMeasurement> vtm = result[0].measurements();
238 
239  TSOS lastState = theTSOS;
240  for (std::vector<TrajectoryMeasurement>::iterator tm=vtm.begin(); tm!=vtm.end();tm++){
241 
243  if ((rhit)->isValid()==0&&rhit->det()!=nullptr) continue;
244  LogTrace("TestSmoothHits") << "new hit" ;
245 
246  int subdetId = rhit->det()->geographicalId().subdetId();
247  DetId id = rhit->det()->geographicalId();
248  int layerId = tTopo->layer(id);
249  LogTrace("TestSmoothHits") << "subdetId=" << subdetId << " layerId=" << layerId ;
250 
251  double delta = 99999;
252  LocalPoint rhitLPv = rhit->localPosition();
253 
254  std::vector<PSimHit> assSimHits = hitAssociator.associateHit(*(rhit->hit()));
255  if (assSimHits.empty()) continue;
256  PSimHit shit;
257  for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
258  if ((m->localPosition()-rhitLPv).mag()<delta) {
259  shit=*m;
260  delta = (m->localPosition()-rhitLPv).mag();
261  }
262  }
263 
264  TSOS currentState = combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
265  TSOS updatedState = tm->updatedState();
266 
267  //plot chi2 increment
268  double chi2increment = tm->estimate();
269  LogTrace("TestSmoothHits") << "tm->estimate()=" << tm->estimate();
270  title.str("");
271  title << "Chi2Increment_" << subdetId << "-" << layerId;
272  hChi2Increment[title.str()]->Fill( chi2increment );
273  hTotChi2Increment->Fill( chi2increment );
274  hChi2_vs_Process->Fill( chi2increment, shit.processType() );
275  if (dynamic_cast<const SiPixelRecHit*>(rhit->hit()))
276  hChi2_vs_clsize->Fill( chi2increment, ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size() );
277  if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit()))
278  hChi2_vs_clsize->Fill( chi2increment, ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() );
279 
280  //test hits
281  const Surface * surf = &( (rhit)->det()->surface() );
282  LocalVector shitLMom;
283  LocalPoint shitLPos;
284  if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
285  double rechitmatchedx = rhit->localPosition().x();
286  double rechitmatchedy = rhit->localPosition().y();
287  double mindist = 999999;
288  float distx, disty;
289  std::pair<LocalPoint,LocalVector> closestPair;
290  const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) ((const GluedGeomDet *)(rhit)->det())->stereoDet();
291  const BoundPlane& plane = (rhit)->det()->surface();
292  for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++) {
293  //project simhit;
294  std::pair<LocalPoint,LocalVector> hitPair = projectHit((*m),stripDet,plane);
295  distx = fabs(rechitmatchedx - hitPair.first.x());
296  disty = fabs(rechitmatchedy - hitPair.first.y());
297  double dist = distx*distx+disty*disty;
298  if(sqrt(dist)<mindist){
299  mindist = dist;
300  closestPair = hitPair;
301  }
302  }
303  shitLPos = closestPair.first;
304  shitLMom = closestPair.second;
305  } else {
306  shitLPos = shit.localPosition();
307  shitLMom = shit.momentumAtEntry();
308  }
309  GlobalVector shitGMom = surf->toGlobal(shitLMom);
310  GlobalPoint shitGPos = surf->toGlobal(shitLPos);
311 // if (dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())) {
312 // double rechitmatchedx = rhit->localPosition().x();
313 // double rechitmatchedy = rhit->localPosition().y();
314 // double mindist = 999999;
315 // double distx, disty;
316 // std::pair<LocalPoint,LocalVector> closestPair;
317 // const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) ((const GluedGeomDet *)(rhit)->det())->stereoDet();
318 // const BoundPlane& plane = (rhit)->det()->surface();
319 // for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
320 // const PSimHit& hit = *m;
321 // const StripTopology& topol = stripDet->specificTopology();
322 // GlobalPoint globalpos = stripDet->surface().toGlobal(hit.localPosition());
323 // LocalPoint localHit = plane.toLocal(globalpos);
324 // //track direction
325 // LocalVector locdir=hit.localDirection();
326 // //rotate track in new frame
327 // GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
328 // LocalVector dir=plane.toLocal(globaldir);
329 // float scale = -localHit.z() / dir.z();
330 // LocalPoint projectedPos = localHit + scale*dir;
331 // float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
332 // LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame
333 // LocalVector localStripDir = LocalVector( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
334 // std::pair<LocalPoint,LocalVector> hitPair( projectedPos, localStripDir);
335 // distx = fabs(rechitmatchedx - hitPair.first.x());
336 // disty = fabs(rechitmatchedy - hitPair.first.y());
337 // double dist = distx*distx+disty*disty;
338 // if(sqrt(dist)<mindist){
339 // mindist = dist;
340 // closestPair = hitPair;
341 // }
342 // }
343 // shitLPos = closestPair.first;
344 // shitLMom = closestPair.second;
345 // } else {
346 // shitLPos = shit.localPosition();
347 // shitLMom = shit.momentumAtEntry();
348 // }
349 // GlobalVector shitGMom = surf->toGlobal(shitLMom);
350 // GlobalPoint shitGPos = surf->toGlobal(shitLPos);
351 
352  GlobalVector tsosGMom = currentState.globalMomentum();
353  GlobalError tsosGMEr(currentState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
354  GlobalPoint tsosGPos = currentState.globalPosition();
355  GlobalError tsosGPEr = currentState.cartesianError().position();
356 
357  GlobalPoint rhitGPos = (rhit)->globalPosition();
358  GlobalError rhitGPEr = (rhit)->globalPositionError();
359 
360  double pullGPX_rs = (rhitGPos.x()-shitGPos.x())/sqrt(rhitGPEr.cxx());
361  double pullGPY_rs = (rhitGPos.y()-shitGPos.y())/sqrt(rhitGPEr.cyy());
362  double pullGPZ_rs = (rhitGPos.z()-shitGPos.z())/sqrt(rhitGPEr.czz());
363  //double pullGPX_rs = (rhitGPos.x()-shitGPos.x());
364  //double pullGPY_rs = (rhitGPos.y()-shitGPos.y());
365  //double pullGPZ_rs = (rhitGPos.z()-shitGPos.z());
366 
367  LogTrace("TestSmoothHits") << "rs" << std::endl;
368 
369  title.str("");
370  title << "PullGP_X_" << subdetId << "-" << layerId << "_rs";
371  hPullGP_X_rs[title.str()]->Fill( pullGPX_rs );
372  title.str("");
373  title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs";
374  hPullGP_Y_rs[title.str()]->Fill( pullGPY_rs );
375  title.str("");
376  title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs";
377  hPullGP_Z_rs[title.str()]->Fill( pullGPZ_rs );
378 
379  double pullGPX_tr = (tsosGPos.x()-rhitGPos.x())/sqrt(tsosGPEr.cxx()+rhitGPEr.cxx());
380  double pullGPY_tr = (tsosGPos.y()-rhitGPos.y())/sqrt(tsosGPEr.cyy()+rhitGPEr.cyy());
381  double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z())/sqrt(tsosGPEr.czz()+rhitGPEr.czz());
382  //double pullGPX_tr = (tsosGPos.x()-rhitGPos.x());
383  //double pullGPY_tr = (tsosGPos.y()-rhitGPos.y());
384  //double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z());
385 
386  LogTrace("TestSmoothHits") << "tr" << std::endl;
387 
388  title.str("");
389  title << "PullGP_X_" << subdetId << "-" << layerId << "_tr";
390  hPullGP_X_tr[title.str()]->Fill( pullGPX_tr );
391  title.str("");
392  title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr";
393  hPullGP_Y_tr[title.str()]->Fill( pullGPY_tr );
394  title.str("");
395  title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr";
396  hPullGP_Z_tr[title.str()]->Fill( pullGPZ_tr );
397 
398  double pullGPX_ts = (tsosGPos.x()-shitGPos.x())/sqrt(tsosGPEr.cxx());
399  double pullGPY_ts = (tsosGPos.y()-shitGPos.y())/sqrt(tsosGPEr.cyy());
400  double pullGPZ_ts = (tsosGPos.z()-shitGPos.z())/sqrt(tsosGPEr.czz());
401  //double pullGPX_ts = (tsosGPos.x()-shitGPos.x());
402  //double pullGPY_ts = (tsosGPos.y()-shitGPos.y());
403  //double pullGPZ_ts = (tsosGPos.z()-shitGPos.z());
404 
405  LogTrace("TestSmoothHits") << "ts1" << std::endl;
406 
407  title.str("");
408  title << "PullGP_X_" << subdetId << "-" << layerId << "_ts";
409  hPullGP_X_ts[title.str()]->Fill( pullGPX_ts );
410  title.str("");
411  title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts";
412  hPullGP_Y_ts[title.str()]->Fill( pullGPY_ts );
413  title.str("");
414  title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts";
415  hPullGP_Z_ts[title.str()]->Fill( pullGPZ_ts );
416 
417  double pullGMX_ts = (tsosGMom.x()-shitGMom.x())/sqrt(tsosGMEr.cxx());
418  double pullGMY_ts = (tsosGMom.y()-shitGMom.y())/sqrt(tsosGMEr.cyy());
419  double pullGMZ_ts = (tsosGMom.z()-shitGMom.z())/sqrt(tsosGMEr.czz());
420  //double pullGMX_ts = (tsosGMom.x()-shitGMom.x());
421  //double pullGMY_ts = (tsosGMom.y()-shitGMom.y());
422  //double pullGMZ_ts = (tsosGMom.z()-shitGMom.z());
423 
424  LogTrace("TestSmoothHits") << "ts2" << std::endl;
425 
426  title.str("");
427  title << "PullGM_X_" << subdetId << "-" << layerId << "_ts";
428  hPullGM_X_ts[title.str()]->Fill( pullGMX_ts );
429  title.str("");
430  title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts";
431  hPullGM_Y_ts[title.str()]->Fill( pullGMY_ts );
432  title.str("");
433  title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts";
434  hPullGM_Z_ts[title.str()]->Fill( pullGMZ_ts );
435 
436  if (dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())) {
437  //mono
438  LogTrace("TestSmoothHits") << "MONO HIT" << std::endl;
439  auto m = dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->monoHit();
440  CTTRHp tMonoHit =
441  theBuilder->build(&m);
442  if (tMonoHit==nullptr) continue;
443  vector<PSimHit> assMonoSimHits = hitAssociator.associateHit(*tMonoHit->hit());
444  if (assMonoSimHits.empty()) continue;
445  const PSimHit sMonoHit = *(assSimHits.begin());
446  const Surface * monoSurf = &( tMonoHit->det()->surface() );
447  if (monoSurf==nullptr) continue;
448  TSOS monoState = thePropagator->propagate(lastState,*monoSurf);
449  if (monoState.isValid()==0) continue;
450 
451  LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
452  GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
453  LocalPoint monoShitLPos = sMonoHit.localPosition();
454  GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
455 
456  GlobalVector monoTsosGMom = monoState.globalMomentum();
457  GlobalError monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
458  GlobalPoint monoTsosGPos = monoState.globalPosition();
459  GlobalError monoTsosGPEr = monoState.cartesianError().position();
460 
461  GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
462  GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
463 
464  double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x())/sqrt(monoRhitGPEr.cxx());
465  double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y())/sqrt(monoRhitGPEr.cyy());
466  double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z())/sqrt(monoRhitGPEr.czz());
467  //double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x());
468  //double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y());
469  //double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z());
470 
471  title.str("");
472  title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
473  hPullGP_X_rs_mono[title.str()]->Fill( pullGPX_rs_mono );
474  title.str("");
475  title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
476  hPullGP_Y_rs_mono[title.str()]->Fill( pullGPY_rs_mono );
477  title.str("");
478  title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
479  hPullGP_Z_rs_mono[title.str()]->Fill( pullGPZ_rs_mono );
480 
481  double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x())/sqrt(monoTsosGPEr.cxx()+monoRhitGPEr.cxx());
482  double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y())/sqrt(monoTsosGPEr.cyy()+monoRhitGPEr.cyy());
483  double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z())/sqrt(monoTsosGPEr.czz()+monoRhitGPEr.czz());
484  //double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x());
485  //double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y());
486  //double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z());
487 
488  title.str("");
489  title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
490  hPullGP_X_tr_mono[title.str()]->Fill( pullGPX_tr_mono );
491  title.str("");
492  title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
493  hPullGP_Y_tr_mono[title.str()]->Fill( pullGPY_tr_mono );
494  title.str("");
495  title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
496  hPullGP_Z_tr_mono[title.str()]->Fill( pullGPZ_tr_mono );
497 
498  double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x())/sqrt(monoTsosGPEr.cxx());
499  double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y())/sqrt(monoTsosGPEr.cyy());
500  double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z())/sqrt(monoTsosGPEr.czz());
501  //double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x());
502  //double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y());
503  //double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z());
504 
505  title.str("");
506  title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
507  hPullGP_X_ts_mono[title.str()]->Fill( pullGPX_ts_mono );
508  title.str("");
509  title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
510  hPullGP_Y_ts_mono[title.str()]->Fill( pullGPY_ts_mono );
511  title.str("");
512  title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
513  hPullGP_Z_ts_mono[title.str()]->Fill( pullGPZ_ts_mono );
514 
515  double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x())/sqrt(monoTsosGMEr.cxx());
516  double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y())/sqrt(monoTsosGMEr.cyy());
517  double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z())/sqrt(monoTsosGMEr.czz());
518  //double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x());
519  //double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y());
520  //double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z());
521 
522  title.str("");
523  title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
524  hPullGM_X_ts_mono[title.str()]->Fill( pullGMX_ts_mono );
525  title.str("");
526  title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
527  hPullGM_Y_ts_mono[title.str()]->Fill( pullGMY_ts_mono );
528  title.str("");
529  title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
530  hPullGM_Z_ts_mono[title.str()]->Fill( pullGMZ_ts_mono );
531 
532  //stereo
533  LogTrace("TestSmoothHits") << "STEREO HIT" << std::endl;
534  auto s = dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->stereoHit();
535  CTTRHp tStereoHit =
536  theBuilder->build(&s);
537  if (tStereoHit==nullptr) continue;
538  vector<PSimHit> assStereoSimHits = hitAssociator.associateHit(*tStereoHit->hit());
539  if (assStereoSimHits.empty()) continue;
540  const PSimHit sStereoHit = *(assSimHits.begin());
541  const Surface * stereoSurf = &( tStereoHit->det()->surface() );
542  if (stereoSurf==nullptr) continue;
543  TSOS stereoState = thePropagator->propagate(lastState,*stereoSurf);
544  if (stereoState.isValid()==0) continue;
545 
546  LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
547  GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
548  LocalPoint stereoShitLPos = sStereoHit.localPosition();
549  GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
550 
551  GlobalVector stereoTsosGMom = stereoState.globalMomentum();
552  GlobalError stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
553  GlobalPoint stereoTsosGPos = stereoState.globalPosition();
554  GlobalError stereoTsosGPEr = stereoState.cartesianError().position();
555 
556  GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
557  GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
558 
559  double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/sqrt(stereoRhitGPEr.cxx());
560  double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/sqrt(stereoRhitGPEr.cyy());
561  double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/sqrt(stereoRhitGPEr.czz());
562  //double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x());
563  //double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y());
564  //double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z());
565 
566  title.str("");
567  title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
568  hPullGP_X_rs_stereo[title.str()]->Fill( pullGPX_rs_stereo );
569  title.str("");
570  title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
571  hPullGP_Y_rs_stereo[title.str()]->Fill( pullGPY_rs_stereo );
572  title.str("");
573  title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
574  hPullGP_Z_rs_stereo[title.str()]->Fill( pullGPZ_rs_stereo );
575 
576  double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x())/sqrt(stereoTsosGPEr.cxx()+stereoRhitGPEr.cxx());
577  double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y())/sqrt(stereoTsosGPEr.cyy()+stereoRhitGPEr.cyy());
578  double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z())/sqrt(stereoTsosGPEr.czz()+stereoRhitGPEr.czz());
579  //double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x());
580  //double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y());
581  //double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z());
582 
583  title.str("");
584  title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
585  hPullGP_X_tr_stereo[title.str()]->Fill( pullGPX_tr_stereo );
586  title.str("");
587  title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
588  hPullGP_Y_tr_stereo[title.str()]->Fill( pullGPY_tr_stereo );
589  title.str("");
590  title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
591  hPullGP_Z_tr_stereo[title.str()]->Fill( pullGPZ_tr_stereo );
592 
593  double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x())/sqrt(stereoTsosGPEr.cxx());
594  double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y())/sqrt(stereoTsosGPEr.cyy());
595  double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z())/sqrt(stereoTsosGPEr.czz());
596  //double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x());
597  //double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y());
598  //double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z());
599 
600  title.str("");
601  title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
602  hPullGP_X_ts_stereo[title.str()]->Fill( pullGPX_ts_stereo );
603  title.str("");
604  title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
605  hPullGP_Y_ts_stereo[title.str()]->Fill( pullGPY_ts_stereo );
606  title.str("");
607  title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
608  hPullGP_Z_ts_stereo[title.str()]->Fill( pullGPZ_ts_stereo );
609 
610  double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x())/sqrt(stereoTsosGMEr.cxx());
611  double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y())/sqrt(stereoTsosGMEr.cyy());
612  double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z())/sqrt(stereoTsosGMEr.czz());
613  //double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x());
614  //double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y());
615  //double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z());
616 
617  title.str("");
618  title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
619  hPullGM_X_ts_stereo[title.str()]->Fill( pullGMX_ts_stereo );
620  title.str("");
621  title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
622  hPullGM_Y_ts_stereo[title.str()]->Fill( pullGMY_ts_stereo );
623  title.str("");
624  title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
625  hPullGM_Z_ts_stereo[title.str()]->Fill( pullGMZ_ts_stereo );
626  }
627  lastState = updatedState;
628  //#endif
629  }
630  }
631  LogTrace("TestSmoothHits") << "end of event" << std::endl;
632 }
633 
635  //file->Write();
636  TDirectory * chi2i = file->mkdir("Chi2_Increment");
637 
638  TDirectory * gp_ts = file->mkdir("GP_TSOS-SimHit");
639  TDirectory * gm_ts = file->mkdir("GM_TSOS-SimHit");
640  TDirectory * gp_tr = file->mkdir("GP_TSOS-RecHit");
641  TDirectory * gp_rs = file->mkdir("GP_RecHit-SimHit");
642 
643  TDirectory * gp_tsx = gp_ts->mkdir("X");
644  TDirectory * gp_tsy = gp_ts->mkdir("Y");
645  TDirectory * gp_tsz = gp_ts->mkdir("Z");
646  TDirectory * gm_tsx = gm_ts->mkdir("X");
647  TDirectory * gm_tsy = gm_ts->mkdir("Y");
648  TDirectory * gm_tsz = gm_ts->mkdir("Z");
649  TDirectory * gp_trx = gp_tr->mkdir("X");
650  TDirectory * gp_try = gp_tr->mkdir("Y");
651  TDirectory * gp_trz = gp_tr->mkdir("Z");
652  TDirectory * gp_rsx = gp_rs->mkdir("X");
653  TDirectory * gp_rsy = gp_rs->mkdir("Y");
654  TDirectory * gp_rsz = gp_rs->mkdir("Z");
655 
656  TDirectory * gp_tsx_mono = gp_ts->mkdir("MONOX");
657  TDirectory * gp_tsy_mono = gp_ts->mkdir("MONOY");
658  TDirectory * gp_tsz_mono = gp_ts->mkdir("MONOZ");
659  TDirectory * gm_tsx_mono = gm_ts->mkdir("MONOX");
660  TDirectory * gm_tsy_mono = gm_ts->mkdir("MONOY");
661  TDirectory * gm_tsz_mono = gm_ts->mkdir("MONOZ");
662  TDirectory * gp_trx_mono = gp_tr->mkdir("MONOX");
663  TDirectory * gp_try_mono = gp_tr->mkdir("MONOY");
664  TDirectory * gp_trz_mono = gp_tr->mkdir("MONOZ");
665  TDirectory * gp_rsx_mono = gp_rs->mkdir("MONOX");
666  TDirectory * gp_rsy_mono = gp_rs->mkdir("MONOY");
667  TDirectory * gp_rsz_mono = gp_rs->mkdir("MONOZ");
668 
669  TDirectory * gp_tsx_stereo = gp_ts->mkdir("STEREOX");
670  TDirectory * gp_tsy_stereo = gp_ts->mkdir("STEREOY");
671  TDirectory * gp_tsz_stereo = gp_ts->mkdir("STEREOZ");
672  TDirectory * gm_tsx_stereo = gm_ts->mkdir("STEREOX");
673  TDirectory * gm_tsy_stereo = gm_ts->mkdir("STEREOY");
674  TDirectory * gm_tsz_stereo = gm_ts->mkdir("STEREOZ");
675  TDirectory * gp_trx_stereo = gp_tr->mkdir("STEREOX");
676  TDirectory * gp_try_stereo = gp_tr->mkdir("STEREOY");
677  TDirectory * gp_trz_stereo = gp_tr->mkdir("STEREOZ");
678  TDirectory * gp_rsx_stereo = gp_rs->mkdir("STEREOX");
679  TDirectory * gp_rsy_stereo = gp_rs->mkdir("STEREOY");
680  TDirectory * gp_rsz_stereo = gp_rs->mkdir("STEREOZ");
681 
682  chi2i->cd();
683  hTotChi2Increment->Write();
684  hChi2_vs_Process->Write();
685  hChi2_vs_clsize->Write();
686  for (int i=0; i!=6; i++)
687  for (int j=0; j!=9; j++){
688  if (i==0 && j>2) break;
689  if (i==1 && j>1) break;
690  if (i==2 && j>3) break;
691  if (i==3 && j>2) break;
692  if (i==4 && j>5) break;
693  if (i==5 && j>8) break;
694  chi2i->cd();
695  title.str("");
696  title << "Chi2Increment_" << i+1 << "-" << j+1;
697  hChi2Increment[title.str()]->Write();
698 
699  gp_ts->cd();
700  gp_tsx->cd();
701  title.str("");
702  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
703  hPullGP_X_ts[title.str()]->Write();
704  gp_tsy->cd();
705  title.str("");
706  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
707  hPullGP_Y_ts[title.str()]->Write();
708  gp_tsz->cd();
709  title.str("");
710  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
711  hPullGP_Z_ts[title.str()]->Write();
712 
713  gm_ts->cd();
714  gm_tsx->cd();
715  title.str("");
716  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
717  hPullGM_X_ts[title.str()]->Write();
718  gm_tsy->cd();
719  title.str("");
720  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
721  hPullGM_Y_ts[title.str()]->Write();
722  gm_tsz->cd();
723  title.str("");
724  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
725  hPullGM_Z_ts[title.str()]->Write();
726 
727  gp_tr->cd();
728  gp_trx->cd();
729  title.str("");
730  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
731  hPullGP_X_tr[title.str()]->Write();
732  gp_try->cd();
733  title.str("");
734  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
735  hPullGP_Y_tr[title.str()]->Write();
736  gp_trz->cd();
737  title.str("");
738  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
739  hPullGP_Z_tr[title.str()]->Write();
740 
741  gp_rs->cd();
742  gp_rsx->cd();
743  title.str("");
744  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
745  hPullGP_X_rs[title.str()]->Write();
746  gp_rsy->cd();
747  title.str("");
748  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
749  hPullGP_Y_rs[title.str()]->Write();
750  gp_rsz->cd();
751  title.str("");
752  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
753  hPullGP_Z_rs[title.str()]->Write();
754 
755  if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
756  //mono
757  gp_ts->cd();
758  gp_tsx_mono->cd();
759  title.str("");
760  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
761  hPullGP_X_ts_mono[title.str()]->Write();
762  gp_tsy_mono->cd();
763  title.str("");
764  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
765  hPullGP_Y_ts_mono[title.str()]->Write();
766  gp_tsz_mono->cd();
767  title.str("");
768  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
769  hPullGP_Z_ts_mono[title.str()]->Write();
770 
771  gm_ts->cd();
772  gm_tsx_mono->cd();
773  title.str("");
774  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
775  hPullGM_X_ts_mono[title.str()]->Write();
776  gm_tsy_mono->cd();
777  title.str("");
778  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
779  hPullGM_Y_ts_mono[title.str()]->Write();
780  gm_tsz_mono->cd();
781  title.str("");
782  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
783  hPullGM_Z_ts_mono[title.str()]->Write();
784 
785  gp_tr->cd();
786  gp_trx_mono->cd();
787  title.str("");
788  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
789  hPullGP_X_tr_mono[title.str()]->Write();
790  gp_try_mono->cd();
791  title.str("");
792  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
793  hPullGP_Y_tr_mono[title.str()]->Write();
794  gp_trz_mono->cd();
795  title.str("");
796  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
797  hPullGP_Z_tr_mono[title.str()]->Write();
798 
799  gp_rs->cd();
800  gp_rsx_mono->cd();
801  title.str("");
802  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
803  hPullGP_X_rs_mono[title.str()]->Write();
804  gp_rsy_mono->cd();
805  title.str("");
806  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
807  hPullGP_Y_rs_mono[title.str()]->Write();
808  gp_rsz_mono->cd();
809  title.str("");
810  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
811  hPullGP_Z_rs_mono[title.str()]->Write();
812 
813  //stereo
814  gp_ts->cd();
815  gp_tsx_stereo->cd();
816  title.str("");
817  title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
818  hPullGP_X_ts_stereo[title.str()]->Write();
819  gp_tsy_stereo->cd();
820  title.str("");
821  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
822  hPullGP_Y_ts_stereo[title.str()]->Write();
823  gp_tsz_stereo->cd();
824  title.str("");
825  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
826  hPullGP_Z_ts_stereo[title.str()]->Write();
827 
828  gm_ts->cd();
829  gm_tsx_stereo->cd();
830  title.str("");
831  title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
832  hPullGM_X_ts_stereo[title.str()]->Write();
833  gm_tsy_stereo->cd();
834  title.str("");
835  title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
836  hPullGM_Y_ts_stereo[title.str()]->Write();
837  gm_tsz_stereo->cd();
838  title.str("");
839  title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
840  hPullGM_Z_ts_stereo[title.str()]->Write();
841 
842  gp_tr->cd();
843  gp_trx_stereo->cd();
844  title.str("");
845  title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
846  hPullGP_X_tr_stereo[title.str()]->Write();
847  gp_try_stereo->cd();
848  title.str("");
849  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
850  hPullGP_Y_tr_stereo[title.str()]->Write();
851  gp_trz_stereo->cd();
852  title.str("");
853  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
854  hPullGP_Z_tr_stereo[title.str()]->Write();
855 
856  gp_rs->cd();
857  gp_rsx_stereo->cd();
858  title.str("");
859  title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
860  hPullGP_X_rs_stereo[title.str()]->Write();
861  gp_rsy_stereo->cd();
862  title.str("");
863  title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
864  hPullGP_Y_rs_stereo[title.str()]->Write();
865  gp_rsz_stereo->cd();
866  title.str("");
867  title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
868  hPullGP_Z_rs_stereo[title.str()]->Write();
869  }
870  }
871 
872  file->Close();
873 }
874 
875 //needed by to do the residual for matched hits
876 //taken from SiStripTrackingRecHitsValid.cc
877 std::pair<LocalPoint,LocalVector>
878 TestSmoothHits::projectHit( const PSimHit& hit, const StripGeomDetUnit* stripDet, const BoundPlane& plane)
879 {
880  const StripTopology& topol = stripDet->specificTopology();
881  GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
882  LocalPoint localHit = plane.toLocal(globalpos);
883  //track direction
884  LocalVector locdir=hit.localDirection();
885  //rotate track in new frame
886 
887  GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
888  LocalVector dir=plane.toLocal(globaldir);
889  float scale = -localHit.z() / dir.z();
890 
891  LocalPoint projectedPos = localHit + scale*dir;
892 
893  float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
894 
895  LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame
896 
897  LocalVector localStripDir( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
898 
899  return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
900 }
901 
904 
std::map< std::string, TH1F * > hPullGP_Z_rs_stereo
dbl * delta
Definition: mlp_gen.cc:36
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:106
T getParameter(std::string const &) const
std::string srcName
std::map< std::string, TH1F * > hPullGP_X_ts_mono
std::map< std::string, TH1F * > hPullGM_Y_ts_mono
edm::Handle< TrackCandidateCollection > theTCCollection
std::stringstream title
range recHits() const
TestSmoothHits(const edm::ParameterSet &)
LocalVector momentumAtEntry() const
The momentum of the track that produced the hit, at entry point.
Definition: PSimHit.h:47
std::map< std::string, TH1F * > hChi2Increment
std::map< std::string, TH1F * > hPullGP_X_rs_mono
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::string propagatorName
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
TrajectorySeed const & seed() const
std::map< std::string, TH1F * > hPullGM_Z_ts_stereo
const CartesianTrajectoryError cartesianError() const
T y() const
Definition: PV3DBase.h:63
virtual float strip(const LocalPoint &) const =0
std::vector< ConstRecHitPointer > RecHitContainer
GlobalPoint globalPosition() const
std::map< std::string, TH1F * > hPullGP_Y_rs
std::map< std::string, TH1F * > hPullGM_X_ts_mono
std::map< std::string, TH1F * > hPullGP_Z_tr
std::map< std::string, TH1F * > hPullGP_X_ts
std::pair< const_iterator, const_iterator > range
std::string fname
std::map< std::string, TH1F * > hPullGP_Y_tr_stereo
std::map< std::string, TH1F * > hPullGP_Z_ts
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:42
std::map< std::string, TH1F * > hPullGP_X_rs_stereo
edm::ESHandle< TrajectorySmoother > smooth
void endJob() override
virtual TrajectoryContainer trajectories(const Trajectory &traj) const
std::map< std::string, TH1F * > hPullGP_Y_rs_mono
std::map< std::string, TH1F * > hPullGM_Y_ts_stereo
PTrajectoryStateOnDet const & trajectoryStateOnDet() const
void beginRun(edm::Run const &run, const edm::EventSetup &) override
TH2F * hChi2_vs_Process
std::map< std::string, TH1F * > hPullGP_Z_rs_mono
edm::ESHandle< MagneticField > theMF
TH2F * hChi2_vs_clsize
TrackerHitAssociator::Config trackerHitAssociatorConfig_
int iEvent
Definition: GenABIO.cc:230
std::map< std::string, TH1F * > hPullGP_Y_ts_mono
std::map< std::string, TH1F * > hPullGM_Y_ts
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
std::map< std::string, TH1F * > hPullGP_Z_tr_mono
Local3DPoint localPosition() const
Definition: PSimHit.h:44
virtual float stripAngle(float strip) const =0
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
std::map< std::string, TH1F * > hPullGM_Z_ts
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::map< std::string, TH1F * > hPullGP_X_tr
TransientTrackingRecHit::ConstRecHitPointer CTTRHp
Definition: CkfDebugger.h:39
unsigned int detId() const
std::map< std::string, TH1F * > hPullGP_Y_rs_stereo
std::string sname
edm::ESHandle< TransientTrackingRecHitBuilder > theBuilder
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:475
#define LogTrace(id)
std::map< std::string, TH1F * > hPullGP_Y_tr_mono
virtual RecHitPointer build(const TrackingRecHit *p) const =0
build a tracking rechit from an existing rechit
const AlgebraicSymMatrix66 & matrix() const
std::map< std::string, TH1F * > hPullGP_X_ts_stereo
std::map< std::string, TH1F * > hPullGP_Y_ts
LocalVector localDirection() const
Obsolete. Same as momentumAtEntry().unit(), for backward compatibility.
Definition: PSimHit.h:52
std::map< std::string, TH1F * > hPullGP_Z_ts_stereo
Definition: DetId.h:18
std::map< std::string, TH1F * > hPullGM_Z_ts_mono
TrajectoryStateOnSurface transientState(const PTrajectoryStateOnDet &ts, const Surface *surface, const MagneticField *field)
~TestSmoothHits() override
const GlobalError position() const
Position error submatrix.
std::map< std::string, TH1F * > hPullGP_X_rs
TrajectoryStateOnSurface TSOS
TH1F * hTotChi2Increment
edm::ESHandle< Propagator > thePropagator
const T & get() const
Definition: EventSetup.h:59
std::map< std::string, TH1F * > hPullGP_X_tr_stereo
void analyze(const edm::Event &, const edm::EventSetup &) override
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:53
unsigned short processType() const
Definition: PSimHit.h:118
std::string builderName
unsigned int layer(const DetId &id) const
edm::ESHandle< TrackerGeometry > theG
std::map< std::string, TH1F * > hPullGP_X_tr_mono
std::pair< LocalPoint, LocalVector > projectHit(const PSimHit &, const StripGeomDetUnit *, const BoundPlane &)
HLT enums.
GlobalVector globalMomentum() const
const TrackerGeomDet * idToDet(DetId) const override
std::map< std::string, TH1F * > hPullGP_Z_tr_stereo
std::map< std::string, TH1F * > hPullGP_Z_rs
std::vector< PSimHit > associateHit(const TrackingRecHit &thit) const
dbl *** dir
Definition: mlp_gen.cc:35
std::map< std::string, TH1F * > hPullGP_Y_ts_stereo
std::map< std::string, TH1F * > hPullGP_Z_ts_mono
Detector det() const
get the detector field from this detid
Definition: DetId.h:36
std::map< std::string, TH1F * > hPullGM_X_ts
T x() const
Definition: PV3DBase.h:62
std::map< std::string, TH1F * > hPullGM_X_ts_stereo
TransientTrackingRecHit::ConstRecHitPointer CTTRHp
T const * product() const
Definition: ESHandle.h:86
std::map< std::string, TH1F * > hPullGP_Y_tr
Definition: Run.h:43
Our base class.
Definition: SiPixelRecHit.h:23