CMS 3D CMS Logo

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