StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
QualityPlotter.h
1 #ifndef QUALITY_PLOTTER_H
2 #define QUALITY_PLOTTER_H
3 
4 #include <map>
5 #include <string>
6 
7 #include "TH1.h"
8 #include "TH1F.h"
9 #include "TH2F.h"
10 #include "TH3F.h"
11 #include "TVector3.h"
12 
13 #include "GenFit/FitStatus.h"
14 
15 #include "StFwdTrackMaker/FwdTrackerConfig.h"
16 #include "StFwdTrackMaker/Common.h"
17 #include "StFwdTrackMaker/include/Tracker/FwdHit.h"
18 
19 
20 class QualityPlotter {
21  public:
22  QualityPlotter(FwdTrackerConfig &_cfg) : cfg(_cfg) {
23  }
24 
25  void makeHistograms(size_t maxI) {
26  using namespace std;
27  maxIterations = maxI;
28  hist["FinalEff"] = new TH1F("FinalEff", ";Eff", 100, 0, 2.0);
29  hist["FinalN7Eff"] = new TH1F("FinalN7Eff", ";Eff", 100, 0, 2.0);
30  hist["FinalN6Eff"] = new TH1F("FinalN6Eff", ";Eff", 100, 0, 2.0);
31  hist["FinalN5Eff"] = new TH1F("FinalN5Eff", ";Eff", 100, 0, 2.0);
32  hist["FinalN4Eff"] = new TH1F("FinalN4Eff", ";Eff", 100, 0, 2.0);
33  hist["FinalN7Quality"] = new TH1F("FinalN7Quality", ";Quality", 150, 0, 1.5);
34  hist["AllQuality"] = new TH1F("AllQuality", ";Quality", 150, 0, 1.5);
35  hist["DurationPerEvent"] = new TH1F("DurationPerEvent", ";Duration(ms) per Event", 1e5, 0, 1e5);
36 
37  hist["nHitsOnTrack"] = new TH1F("nHitsOnTrack", ";nHit", 10, 0, 10);
38  hist["nHitsOnTrackMc"] = new TH1F("nHitsOnTrackMc", ";nHit", 10, 0, 10);
39 
40  hist["InvPtRes"] = new TH1F("InvPtRes", ";(p_{T}^{MC} - p_{T}^{RC}) / p_{T}^{MC}", 500, -5, 5);
41  hist["InvPtResVsNHits"] = new TH2F("InvPtResVsNHits", ";(p_{T}^{MC} - p_{T}^{RC}) / p_{T}^{MC}", 10, 0, 10, 500, -5, 5);
42  hist["PtRes"] = new TH1F("PtRes", ";(p_{T}^{MC} - p_{T}^{RC}) / p_{T}^{MC}", 500, -5, 5);
43  hist["PtResVsTrue"] = new TH2F("PtResVsTrue", ";q^{MC} #times p_{T}^{MC};(p_{T}^{MC} - p_{T}^{RC}) / p_{T}^{MC}", 100, -5, 5, 200, -2, 2);
44  hist["InvPtResVsTrue"] = new TH2F("InvPtResVsTrue", ";q^{MC} #times p_{T}^{MC}; #sigma_{p_{T}^{-1}}", 100, -5, 5, 200, -2, 2);
45  hist["DeltaPt"] = new TH1F("DeltaPt", ";p_{T}^{RC} - p_{T}^{MC} (GeV/c)", 500, -5, 5);
46  hist["InvPtResVsEta"] = new TH2F("InvPtResVsEta", ";#eta^{MC}; #sigma_{p_{T}^{-1}}", 300, 2, 5, 200, -2, 2);
47 
48  hist["RecoPtVsMcPt"] = new TH2F("RecoPtVsMcPt", "; p_{T}^{MC}; p_{T}^{RC}", 100, 0, 5, 100, 0, 5);
49  hist["QMatrix"] = new TH2F("QMatrix", ";GEN;RECO", 6, -3, 3, 6, -3, 3);
50  hist["FitPValue"] = new TH1F("FitPValue", ";", 1000, 0, 10);
51  hist["FitChi2"] = new TH1F("FitChi2", ";", 2000, 0, 200);
52  hist["FitChi2Ndf"] = new TH1F("FitChi2Ndf", ";", 1500, 0, 150);
53  hist["FitNFailedHits"] = new TH1F("FitNFailedHits", ";", 10, 0, 10);
54 
55  hist["RightQVsMcPt"] = new TH1F("RightQVsMcPt", ";p_{T}^{GEN}; N", 100, -5, 5);
56  hist["WrongQVsMcPt"] = new TH1F("WrongQVsMcPt", ";p_{T}^{GEN}; N", 100, -5, 5);
57  hist["AllQVsMcPt"] = new TH1F("AllQVsMcPt", ";p_{T}^{GEN}; N", 100, -5, 5);
58 
59  hist["McPt"] = new TH1F("McPt", ";p_{T}^{MC} (GeV/c)", 100, 0, 10);
60  hist["McPt_4hits"] = new TH1F("McPt_4hits", ";p_{T}^{MC} (GeV/c)", 100, 0, 10);
61  hist["McPt_5hits"] = new TH1F("McPt_5hits", ";p_{T}^{MC} (GeV/c)", 100, 0, 10);
62  hist["McPt_6hits"] = new TH1F("McPt_6hits", ";p_{T}^{MC} (GeV/c)", 100, 0, 10);
63  hist["McPt_7hits"] = new TH1F("McPt_7hits", ";p_{T}^{MC} (GeV/c)", 100, 0, 10);
64  hist["McPtFound"] = new TH1F("McPtFound", ";p_{T}^{MC} (GeV/c)", 100, 0, 10);
65  hist["McPtFoundAllQ"] = new TH1F("McPtFoundAllQ", ";p_{T}^{MC} (GeV/c)", 100, 0, 10);
66 
67  hist["McEta"] = new TH1F("McEta", ";#eta^{MC} (GeV/c)", 100, 2, 5);
68  hist["McEta_4hits"] = new TH1F("McEta_4hits", ";#eta^{MC} (GeV/c)", 100, 2, 5);
69  hist["McEta_5hits"] = new TH1F("McEta_5hits", ";#eta^{MC} (GeV/c)", 100, 2, 5);
70  hist["McEta_6hits"] = new TH1F("McEta_6hits", ";#eta^{MC} (GeV/c)", 100, 2, 5);
71  hist["McEta_7hits"] = new TH1F("McEta_7hits", ";#eta^{MC} (GeV/c)", 100, 2, 5);
72  hist["McEtaFound"] = new TH1F("McEtaFound", ";#eta^{MC} (GeV/c)", 100, 2, 5);
73  hist["McEtaFoundAllQ"] = new TH1F("McEtaFoundAllQ", ";#eta^{MC} (GeV/c)", 100, 2, 5);
74 
75  hist["McPhi"] = new TH1F("McPhi", ";#phi^{MC} (GeV/c)", 128, -3.2, 3.2);
76  hist["McPhi_4hits"] = new TH1F("McPhi_4hits", ";#phi^{MC} (GeV/c)", 128, -3.2, 3.2);
77  hist["McPhi_5hits"] = new TH1F("McPhi_5hits", ";#phi^{MC} (GeV/c)", 128, -3.2, 3.2);
78  hist["McPhi_6hits"] = new TH1F("McPhi_6hits", ";#phi^{MC} (GeV/c)", 128, -3.2, 3.2);
79  hist["McPhi_7hits"] = new TH1F("McPhi_7hits", ";#phi^{MC} (GeV/c)", 128, -3.2, 3.2);
80  hist["McPhiFound"] = new TH1F("McPhiFound", ";#phi^{MC} (GeV/c)", 128, -3.2, 3.2);
81  hist["McPhiFoundAllQ"] = new TH1F("McPhiFoundAllQ", ";#phi^{MC} (GeV/c)", 128, -3.2, 3.2);
82  hist["nFailedFits"] = new TH1F("nFailedFits", ";;nFailedFits", 5, 0, 5);
83 
84  // 2D efficiency
85  hist["McPtPhiFound"] = new TH2F("McPtPhiFound", "; p_{T}^{MC} (GeV/c);#phi^{MC} (GeV/c)", 20, 0, 10, 64, -3.2, 3.2);
86  hist["McPtPhiFoundAllQ"] = new TH2F("McPtPhiFoundAllQ", "; p_{T}^{MC} (GeV/c);#phi^{MC} (GeV/c)", 20, 0, 10, 64, -3.2, 3.2);
87  hist["McPtPhi_4hits"] = new TH2F("McPtPhi_4hits", "; p_{T}^{MC} (GeV/c);#phi^{MC} (GeV/c)", 20, 0, 10, 64, -3.2, 3.2);
88 
89  // 3D Efficiency
90  hist["McPtEtaPhiFound"] = new TH3F("McPtEtaPhiFound", "; p_{T}^{MC} (GeV/c); #eta; #phi^{MC} (GeV/c)", 50, 0, 5, 30, 2, 5, 64, -3.2, 3.2);
91  hist["McPtEtaPhiFoundAllQ"] = new TH3F("McPtEtaPhiFoundAllQ", "; p_{T}^{MC} (GeV/c); #eta; #phi^{MC} (GeV/c)", 50, 0, 5, 30, 2, 5, 64, -3.2, 3.2);
92  hist["McPtEtaPhi_4hits"] = new TH3F("McPtEtaPhi_4hits", "; p_{T}^{MC} (GeV/c); #eta; #phi^{MC} (GeV/c)", 50, 0, 5, 30, 2, 5, 64, -3.2, 3.2);
93 
94  // for ( size_t track_len : { 4, 5, 6, 7 } )
95  {
96  size_t track_len = 4;
97  // only Q > 0.9
98  string n = TString::Format("McPtFound%lu", track_len).Data();
99  hist[n] = new TH1F(n.c_str(), ";p_{T}^{MC} (GeV/c)", 100, 0, 10);
100  n = TString::Format("McEtaFound%lu", track_len).Data();
101  hist[n] = new TH1F(n.c_str(), ";#eta^{MC} (GeV/c)", 100, 2, 5);
102  n = TString::Format("McPhiFound%lu", track_len).Data();
103  hist[n] = new TH1F(n.c_str(), ";#phi^{MC} (GeV/c)", 128, -3.2, 3.2);
104 
105  n = TString::Format("McPtPhiFound%lu", track_len).Data();
106  hist[n] = new TH2F(n.c_str(), "; p_{T}^{MC} (GeV/c);#phi^{MC} (GeV/c)", 20, 0, 10, 64, -3.2, 3.2);
107 
108  n = TString::Format("McPtEtaPhiFound%lu", track_len).Data();
109  hist[n] = new TH3F(n.c_str(), "; p_{T}^{MC} (GeV/c); #eta; #phi^{MC} (GeV/c)", 50, 0, 5, 30, 2, 5, 64, -3.2, 3.2);
110 
111  // all Q versions
112  n = TString::Format("McPtFound%luAllQ", track_len).Data();
113  hist[n] = new TH1F(n.c_str(), ";p_{T}^{MC} (GeV/c)", 100, 0, 10);
114  n = TString::Format("McEtaFound%luAllQ", track_len).Data();
115  hist[n] = new TH1F(n.c_str(), ";#eta^{MC} (GeV/c)", 100, 2, 5);
116  n = TString::Format("McPhiFound%luAllQ", track_len).Data();
117  hist[n] = new TH1F(n.c_str(), ";#phi^{MC} (GeV/c)", 128, -3.2, 3.2);
118 
119  n = TString::Format("McPtPhiFound%luAllQ", track_len).Data();
120  hist[n] = new TH2F(n.c_str(), "; p_{T}^{MC} (GeV/c);#phi^{MC} (GeV/c)", 20, 0, 10, 64, -3.2, 3.2);
121 
122  n = TString::Format("McPtEtaPhiFound%luAllQ", track_len).Data();
123  hist[n] = new TH3F(n.c_str(), "; p_{T}^{MC} (GeV/c); #eta; #phi^{MC} (GeV/c)", 50, 0, 5, 30, 2, 5, 64, -3.2, 3.2);
124  }
125 
126  hist["McHitMap"] = new TH1F("McHitMap", ";VID", 15, 0, 15);
127 
128  for (size_t i = 0; i < 15; i++) { // hack to prevent crash...
129  string n = "McHitMapLayer" + to_string(i);
130  hist[n] = new TH2F(n.c_str(), ("Layer " + to_string(i) + ";x;y").c_str(), 200, 100, 100, 200, 100, 100);
131  }
132 
133  for (size_t i = 0; i < maxIterations; i++) {
134  string n = "NTracksAfter" + to_string(i);
135  hist[n] = new TH1F(n.c_str(), ("N Tracks After Iteration" + to_string(i)).c_str(), 500, 0, 500);
136 
137  n = "RunningFractionFoundVsIt" + to_string(i);
138  hist[n] = new TH1F(n.c_str(), (";Found ( All It <= " + to_string(i) + " ) / Total Found").c_str(), 110, 0, 1.1);
139 
140  n = "FractionFoundVsIt" + to_string(i);
141  hist[n] = new TH1F(n.c_str(), (";Found ( It = " + to_string(i) + " ) / Total Found").c_str(), 110, 0, 1.1);
142 
143  n = "DurationIt" + to_string(i);
144  hist[n] = new TH1F(n.c_str(), (";Duration(ms) for Iteration " + to_string(i)).c_str(), 1e5, 0, 1e5);
145  }
146  }
147 
148  void writeHistograms() {
149  for (auto nh : hist) {
150  nh.second->SetDirectory(gDirectory);
151  nh.second->Write();
152  }
153  }
154 
155  TH1 *get(std::string hn) {
156  if (hist.count(hn) > 0)
157  return hist[hn];
158  return nullptr; //careful
159  }
160 
161  void startIteration() {
162  // start the timer
163  itStart = FwdTrackerUtils::nowNanoSecond();
164  }
165  void afterIteration(size_t iteration, std::vector<Seed_t> acceptedTracks) {
166 
167  size_t nTracks = acceptedTracks.size();
168  nTracksAfterIteration.push_back(nTracks); // assume that we call the iterations in order
169 
170  long long itEnd = FwdTrackerUtils::nowNanoSecond();
171  long long duration = (itEnd - itStart) * 1e-6; // milliseconds
172  this->get("DurationIt" + to_string(iteration))->Fill(duration);
173  }
174 
175  void startEvent() {
176  eventStart = FwdTrackerUtils::nowNanoSecond();
177  }
178  void summarizeEvent(std::vector<Seed_t> foundTracks, std::map<int, shared_ptr<McTrack>> &mcTrackMap, std::vector<TVector3> fitMoms, std::vector<genfit::FitStatus> fitStatus) {
179 
180  using namespace std;
181 
182  long long duration = (FwdTrackerUtils::nowNanoSecond() - eventStart) * 1e-6; // milliseconds
183  this->get("DurationPerEvent")->Fill(duration);
184 
185  // make a map of the number of tracks found for each # of hits
186  map<size_t, size_t> tracks_found_by_nHits;
187 
188  for (auto t : foundTracks) {
189  tracks_found_by_nHits[t.size()]++;
190  }
191 
192  for (size_t i = 0; i < 9; i++) {
193  if (tracks_found_by_nHits.count(i) > 0)
194  this->get("nHitsOnTrack")->Fill(i, tracks_found_by_nHits[i]);
195  }
196 
197  // the total number of tracks found (all nHits)
198  size_t nTotal = foundTracks.size();
199 
200  // Now compute fraction of tracks found vs. iterations
201  float runningFrac = 0;
202 
203  for (size_t i = 0; i < maxIterations; i++) {
204  if (nTracksAfterIteration.size() < i + 1) {
205  break;
206  }
207 
208  float frac = (float)nTracksAfterIteration[i] / (float)nTotal;
209  this->get("FractionFoundVsIt" + to_string(i))->Fill(frac);
210  runningFrac += frac;
211  this->get("RunningFractionFoundVsIt" + to_string(i))->Fill(runningFrac);
212  }
213 
214  // fill McInfo
215  for (auto kv : mcTrackMap) {
216 
217  if (kv.second == nullptr)
218  continue;
219 
220  this->get("nHitsOnTrackMc")->Fill(kv.second->mHits.size());
221  this->get("McPt")->Fill(kv.second->mPt);
222  this->get("McEta")->Fill(kv.second->mEta);
223  this->get("McPhi")->Fill(kv.second->mPhi);
224 
225  if (kv.second->mHits.size() >= 4) {
226  this->get("McPt_4hits")->Fill(kv.second->mPt);
227  this->get("McEta_4hits")->Fill(kv.second->mEta);
228  this->get("McPhi_4hits")->Fill(kv.second->mPhi);
229 
230  this->get("McPtPhi_4hits")->Fill(kv.second->mPt, kv.second->mPhi);
231  ((TH3 *)this->get("McPtEtaPhi_4hits"))->Fill(kv.second->mPt, kv.second->mEta, kv.second->mPhi);
232  }
233 
234  if (kv.second->mHits.size() >= 5) {
235  this->get("McPt_5hits")->Fill(kv.second->mPt);
236  this->get("McEta_5hits")->Fill(kv.second->mEta);
237  this->get("McPhi_5hits")->Fill(kv.second->mPhi);
238  }
239 
240  if (kv.second->mHits.size() >= 6) {
241  this->get("McPt_6hits")->Fill(kv.second->mPt);
242  this->get("McEta_6hits")->Fill(kv.second->mEta);
243  this->get("McPhi_6hits")->Fill(kv.second->mPhi);
244  }
245 
246  if (kv.second->mHits.size() >= 7) {
247  this->get("McPt_7hits")->Fill(kv.second->mPt);
248  this->get("McEta_7hits")->Fill(kv.second->mEta);
249  this->get("McPhi_7hits")->Fill(kv.second->mPhi);
250  }
251 
252  for (auto h : kv.second->mHits) {
253  auto fh = static_cast<FwdHit *>(h);
254  this->get("McHitMap")->Fill(abs(fh->_vid));
255  std::string n = "McHitMapLayer" + std::to_string(fh->getLayer());
256  this->get(n)->Fill(fh->getX(), fh->getY());
257  }
258  }
259 
260  float avgQuality = 0;
261  // calculate quality of tracks
262  size_t track_index = 0;
263 
264  for (auto t : foundTracks) {
265 
266  map<int, float> qual_map;
267 
268  for (auto hit : t) {
269  qual_map[static_cast<FwdHit *>(hit)->_tid]++;
270  }
271 
272  for (auto &kv : qual_map) {
273  kv.second = kv.second / t.size();
274  }
275 
276  // now get the quality corresponding to most hits on single track
277  // we need to remeber the key also, so we can lookup a hit which corresponds to maxq - to retrieve track info
278  float quality = 0;
279  int mctid = -1;
280 
281  for (auto kv : qual_map) {
282  if (kv.second >= quality) {
283  quality = kv.second;
284  mctid = kv.first;
285  }
286  }
287 
288  avgQuality += quality;
289  this->get("AllQuality")->Fill(quality);
290 
291  if (mctid > 0 && quality >= 3.0 / 4.0 - 0.001) {
292  this->get("McPtFoundAllQ")->Fill(mcTrackMap[mctid]->mPt);
293  this->get("McEtaFoundAllQ")->Fill(mcTrackMap[mctid]->mEta);
294  this->get("McPhiFoundAllQ")->Fill(mcTrackMap[mctid]->mPhi);
295 
296  // for ( size_t min_track_len : { 4, 5, 6, 7 } )
297  {
298  size_t min_track_len = 4;
299  if (t.size() >= min_track_len) {
300  this->get(TString::Format("McPtPhiFound%luAllQ", min_track_len).Data())->Fill(mcTrackMap[mctid]->mPt, mcTrackMap[mctid]->mPhi);
301  ((TH3 *)this->get(TString::Format("McPtEtaPhiFound%luAllQ", min_track_len).Data()))->Fill(mcTrackMap[mctid]->mPt, mcTrackMap[mctid]->mEta, mcTrackMap[mctid]->mPhi);
302  this->get(TString::Format("McPtFound%luAllQ", min_track_len).Data())->Fill(mcTrackMap[mctid]->mPt);
303  this->get(TString::Format("McEtaFound%luAllQ", min_track_len).Data())->Fill(mcTrackMap[mctid]->mEta);
304  this->get(TString::Format("McPhiFound%luAllQ", min_track_len).Data())->Fill(mcTrackMap[mctid]->mPhi);
305  }
306  }
307 
308  } // quality >= 3/4
309 
310  // if ( mctid > 0 && quality >= 3.0/4.0 ) {
311  if (mctid > 0 && quality >= 4.0 / 4.0) {
312  // for ( size_t min_track_len : { 4, 5, 6, 7 } )
313  {
314  size_t min_track_len = 4;
315  if (t.size() >= min_track_len) {
316  this->get(TString::Format("McPtPhiFound%lu", min_track_len).Data())->Fill(mcTrackMap[mctid]->mPt, mcTrackMap[mctid]->mPhi);
317  ((TH3 *)this->get(TString::Format("McPtEtaPhiFound%lu", min_track_len).Data()))->Fill(mcTrackMap[mctid]->mPt, mcTrackMap[mctid]->mEta, mcTrackMap[mctid]->mPhi);
318  this->get(TString::Format("McPtFound%lu", min_track_len).Data())->Fill(mcTrackMap[mctid]->mPt);
319  this->get(TString::Format("McEtaFound%lu", min_track_len).Data())->Fill(mcTrackMap[mctid]->mEta);
320  this->get(TString::Format("McPhiFound%lu", min_track_len).Data())->Fill(mcTrackMap[mctid]->mPhi);
321  }
322  }
323 
324  this->get("McPtFound")->Fill(mcTrackMap[mctid]->mPt);
325  this->get("McPtPhiFound")->Fill(mcTrackMap[mctid]->mPt, mcTrackMap[mctid]->mPhi);
326  ((TH3 *)this->get("McPtEtaPhiFound"))->Fill(mcTrackMap[mctid]->mPt, mcTrackMap[mctid]->mEta, mcTrackMap[mctid]->mPhi);
327  this->get("McEtaFound")->Fill(mcTrackMap[mctid]->mEta);
328  this->get("McPhiFound")->Fill(mcTrackMap[mctid]->mPhi);
329 
330  float mcpt = mcTrackMap[mctid]->mPt;
331  float mceta = mcTrackMap[mctid]->mEta;
332  int mcq = (int)mcTrackMap[mctid]->mQ;
333  float rcpt = 0;
334 
335  int rcq = 0;
336  float pval = 999;
337  float chi2 = 999;
338  float rchi2 = 999;
339  int nFailedHits = 0;
340 
341  if (track_index < fitStatus.size()) {
342  rcq = (int)fitStatus[track_index].getCharge();
343  pval = fitStatus[track_index].getPVal();
344  chi2 = fitStatus[track_index].getChi2();
345  rchi2 = fitStatus[track_index].getChi2() / fitStatus[track_index].getNdf();
346  nFailedHits = fitStatus[track_index].getNFailedPoints();
347  rcpt = fitMoms[track_index].Pt();
348 
349  if (fitStatus[track_index].isFitConvergedFully() == false) {
350  rcq = -10;
351  pval = 999;
352  chi2 = -10;
353  rchi2 = -10;
354  nFailedHits = 9;
355  }
356  }
357 
358  float dPt = mcpt - rcpt;
359  float dInvPt = (1.0 / mcpt) - (1.0 / rcpt);
360 
361  if (t.size() >= 4 && rcpt > 0.01) {
362  this->get("DeltaPt")->Fill(dPt);
363  this->get("PtRes")->Fill(dPt / mcpt);
364  this->get("InvPtRes")->Fill(dInvPt / (1.0 / mcpt));
365  this->get("InvPtResVsNHits")->Fill(t.size(), dInvPt / (1.0 / mcpt));
366  this->get("PtResVsTrue")->Fill(mcpt * mcq, dPt / mcpt);
367  this->get("InvPtResVsTrue")->Fill(mcpt * mcq, dInvPt / (1.0 / mcpt));
368  this->get("InvPtResVsEta")->Fill(mceta, dInvPt / (1.0 / mcpt));
369  this->get("RecoPtVsMcPt")->Fill(mcpt, rcpt);
370  this->get("FitPValue")->Fill(pval);
371  this->get("FitChi2")->Fill(chi2);
372  this->get("FitChi2Ndf")->Fill(rchi2);
373  this->get("FitNFailedHits")->Fill(nFailedHits);
374 
375  if (abs(rcq) == 1) {
376  this->get("QMatrix")->Fill(mcq, rcq);
377  }
378 
379  if (mcq == rcq)
380  this->get("RightQVsMcPt")->Fill(mcpt * mcq);
381  else if (rcq != -10)
382  this->get("WrongQVsMcPt")->Fill(mcpt * mcq);
383 
384  if (rcq != -10)
385  this->get("AllQVsMcPt")->Fill(mcpt * mcq);
386  }
387 
388  if (rcpt < 0.01) {
389  this->get("nFailedFits")->Fill(1);
390  }
391 
392  } else if (mctid == 0) {
393  // this->get( "McPtFound" )->Fill( 0 );
394  // this->get( "McEtaFound" )->Fill( 0 );
395  // this->get( "McPhiFound" )->Fill( -10 );
396  }
397 
398  // fill the pT versus efficiency for found tracks
399  // this->
400  track_index++;
401  } // found track
402 
403  avgQuality /= (float)nTotal;
404  this->get("FinalN7Quality")->Fill(avgQuality);
405 
406  nTracksAfterIteration.clear();
407  } // summarize event
408 
409  void finish() {
410  hist["NQMatrix"] = (TH2 *)this->get("QMatrix")->Clone("NQMatrix");
411  this->get("NQMatrix")->Scale(1.0 / this->get("NQMatrix")->GetEntries());
412 
413  this->hist["EffVsMcPt"] = (TH1 *)this->get("McPtFound")->Clone("EffVsMcPt");
414  this->get("EffVsMcPt")->Divide(this->get("McPt"));
415 
416  this->hist["EffVsMcEta"] = (TH1 *)this->get("McEtaFound")->Clone("EffVsMcEta");
417  this->get("EffVsMcEta")->Divide(this->get("McEta"));
418 
419  this->hist["EffVsMcPhi"] = (TH1 *)this->get("McPhiFound")->Clone("EffVsMcPhi");
420  this->get("EffVsMcPhi")->Divide(this->get("McPhi"));
421 
422  // for ( size_t i : { 4 } ) {
423  {
424  size_t i = 4;
425  string n = TString::Format("McPt_%luhits", i).Data();
426  this->hist["EffVs" + n] = (TH1 *)this->get(TString::Format("McPtFound%lu", i).Data())->Clone(("EffVs" + n).c_str());
427  this->get("EffVs" + n)->Divide(this->get(n));
428 
429  this->hist["EffVs" + n + "_AllQ"] = (TH1 *)this->get(TString::Format("McPtFound%luAllQ", i).Data())->Clone(("EffVs" + n + "_AllQ").c_str());
430  this->get("EffVs" + n + "_AllQ")->Divide(this->get(n));
431 
432  n = TString::Format("McEta_%luhits", i).Data();
433  this->hist["EffVs" + n] = (TH1 *)this->get(TString::Format("McEtaFound%lu", i).Data())->Clone(("EffVs" + n).c_str());
434  this->get("EffVs" + n)->Divide(this->get(n));
435 
436  this->hist["EffVs" + n + "_AllQ"] = (TH1 *)this->get(TString::Format("McEtaFound%luAllQ", i).Data())->Clone(("EffVs" + n + "_AllQ").c_str());
437  this->get("EffVs" + n + "_AllQ")->Divide(this->get(n));
438 
439  n = TString::Format("McPhi_%luhits", i).Data();
440  this->hist["EffVs" + n] = (TH1 *)this->get(TString::Format("McPhiFound%lu", i).Data())->Clone(("EffVs" + n).c_str());
441  this->get("EffVs" + n)->Divide(this->get(n));
442 
443  this->hist["EffVs" + n + "_AllQ"] = (TH1 *)this->get(TString::Format("McPhiFound%luAllQ", i).Data())->Clone(("EffVs" + n + "_AllQ").c_str());
444  this->get("EffVs" + n + "_AllQ")->Divide(this->get(n));
445 
446  n = TString::Format("McPtPhi_%luhits", i).Data();
447  this->hist["EffVs" + n] = (TH1 *)this->get(TString::Format("McPtPhiFound%lu", i).Data())->Clone(("EffVs" + n).c_str());
448  this->get("EffVs" + n)->Divide(this->get(n));
449 
450  this->hist["EffVs" + n + "_AllQ"] = (TH1 *)this->get(TString::Format("McPtPhiFound%luAllQ", i).Data())->Clone(("EffVs" + n + "_AllQ").c_str());
451  this->get("EffVs" + n + "_AllQ")->Divide(this->get(n));
452 
453  n = TString::Format("McPtEtaPhi_%luhits", i).Data();
454  this->hist["EffVs" + n] = (TH1 *)this->get(TString::Format("McPtEtaPhiFound%lu", i).Data())->Clone(("EffVs" + n).c_str());
455  this->get("EffVs" + n)->Divide(this->get(n));
456 
457  this->hist["EffVs" + n + "_AllQ"] = (TH1 *)this->get(TString::Format("McPtEtaPhiFound%luAllQ", i).Data())->Clone(("EffVs" + n + "_AllQ").c_str());
458  this->get("EffVs" + n + "_AllQ")->Divide(this->get(n));
459  }
460 
461  this->hist["EffVsMcPtAllQ"] = (TH1 *)this->get("McPtFoundAllQ")->Clone("EffVsMcPtAllQ");
462  this->get("EffVsMcPtAllQ")->Divide(this->get("McPt"));
463 
464  this->hist["EffVsMcEtaAllQ"] = (TH1 *)this->get("McEtaFoundAllQ")->Clone("EffVsMcEtaAllQ");
465  this->get("EffVsMcEtaAllQ")->Divide(this->get("McEta"));
466 
467  this->hist["EffVsMcPhiAllQ"] = (TH1 *)this->get("McPhiFoundAllQ")->Clone("EffVsMcPhiAllQ");
468  this->get("EffVsMcPhiAllQ")->Divide(this->get("McPhi"));
469 
470  // make the charge misid hists
471  auto hRQ = this->get("RightQVsMcPt");
472  auto hWQ = this->get("WrongQVsMcPt");
473  TH1 *hAllQSum = (TH1 *)hRQ->Clone("hAllQSum");
474  hAllQSum->Add(hWQ);
475 
476  this->hist["ChargeMisIdVsMcPt"] = (TH1 *)hWQ->Clone("ChargeMisIdVsMcPt");
477  this->hist["ChargeMisIdVsMcPt"]->Divide(hAllQSum);
478  }
479 
480  private:
481  FwdTrackerConfig &cfg;
482  std::map<std::string, TH1 *> hist;
483 
484  vector<size_t> nTracksAfterIteration;
485  size_t maxIterations;
486  long long itStart;
487  long long eventStart;
488 };
489 
490 #endif
FwdHit
Definition: FwdHit.h:68
Generated by   doxygen 1.8.5