CMS 3D CMS Logo

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