CMS 3D CMS Logo

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