CbmRoot
PairAnalysisCutQA.cxx
Go to the documentation of this file.
1 // //
3 // Authors:
4 // Julian Book <Julian.Book@cern.ch>
5 //
6 /*
7  Allow to monitor how many tracks, pairs, events pass the selection criterion
8  in any of the cuts added to the corresponding filters. Further it automatically
9  calculates the MC matching efficiency seperately for each detector and base PDG
10  particle after each cut.
11 
12  All you need to add to your config is the following:
13 
14  PairAnalysis::SetCutQA();
15 
16 
17 */
18 // //
20 
21 #include "PairAnalysisCutQA.h"
22 
23 #include <TCollection.h>
24 #include <TList.h>
25 #include <TVectorD.h>
26 
27 #include "CbmMCTrack.h"
28 
29 #include "AnalysisCuts.h"
30 #include "PairAnalysisCutGroup.h"
31 
32 #include "PairAnalysisEvent.h"
33 #include "PairAnalysisPairKF.h"
34 #include "PairAnalysisPairLV.h"
35 #include "PairAnalysisTrack.h"
36 
37 #include "PairAnalysisHelper.h"
38 
39 
41 
42 
44  : PairAnalysisCutQA("QAcuts", "QAcuts") {
45  //
46  // Default constructor
47  //
48 }
49 
50 //_____________________________________________________________________
51 PairAnalysisCutQA::PairAnalysisCutQA(const char* name, const char* title)
52  : TNamed(name, title), fQAHistList() {
53  //
54  // Named Constructor
55  //
56  for (Int_t itype = 0; itype < kNtypes; itype++) {
57  fNCuts[itype] = 1;
58  for (Int_t i = 0; i < 20; i++) {
59  fCutNames[i][itype] = "";
60  }
61  }
62  fTypeKeys[kTrack] = "Track";
63  fTypeKeys[kTrack2] = "FinalTrack";
64  fTypeKeys[kTrackMC] = "MCTrack";
65  fTypeKeys[kPair] = "Pair";
66  fTypeKeys[kPrePair] = "PrePair";
67  fTypeKeys[kEvent] = "Event";
68  fQAHistList.SetOwner(kFALSE);
69 }
70 
71 //_____________________________________________________________________
73  //
74  //Default Destructor
75  //
76  fQAHistList.Clear();
77 }
78 
79 //_____________________________________________________________________
81 
82  fQAHistList.SetName(Form("%s", GetName()));
83 
84  THashList* table = new THashList;
85  table->SetOwner(kTRUE);
86  table->SetName("Event");
87  fQAHistList.Add(table);
88 
89  table = new THashList;
90  table->SetOwner(kTRUE);
91  table->SetName("Track");
92  fQAHistList.Add(table);
93 
94  table = new THashList;
95  table->SetOwner(kTRUE);
96  table->SetName("Pair");
97  fQAHistList.Add(table);
98 
99 
100  TH1I* fCutQA = 0x0; // qa histogram for counts
101  TH2I* fPdgCutQA = 0x0; // qa histogram for PDG codes
102  TProfile2D* fEffCutQA = 0x0; // qa histogram for matching efficicy
103 
104  const TVectorD* binsPdg = PairAnalysisHelper::MakeLinBinning(5, 0, 5);
105  const TVectorD* binsDet = PairAnalysisHelper::MakeLinBinning(6, 0, 6);
106  // loop over all types
107  for (Int_t itype = 0; itype < kNtypes; itype++) {
108  // printf("\n type: %d\n",itype);
109  TString logic = "passed";
110  if (itype == kPrePair) logic = "rejected";
111 
112  const TVectorD* binsX =
114  // create histogram based on added cuts
115  fCutQA = new TH1I(
116  fTypeKeys[itype],
117  Form(
118  "%sQA;cuts;# %s %ss", fTypeKeys[itype], logic.Data(), fTypeKeys[itype]),
119  fNCuts[itype],
120  binsX->GetMatrixArray());
121 
122  if (itype == kTrack || itype == kTrack2) {
123  fPdgCutQA = new TH2I(Form("%sPDG", fTypeKeys[itype]),
124  Form("%sPDG;cuts;PDG code;# %s %ss",
125  fTypeKeys[itype],
126  logic.Data(),
127  fTypeKeys[itype]),
128  fNCuts[itype],
129  binsX->GetMatrixArray(),
130  binsPdg->GetNrows() - 1,
131  binsPdg->GetMatrixArray());
132 
133  fEffCutQA =
134  new TProfile2D(Form("%sMatchEff", fTypeKeys[itype]),
135  Form("%sMatchEff;cuts;detector;<#epsilon_{match}^{MC}>",
136  fTypeKeys[itype]),
137  fNCuts[itype],
138  binsX->GetMatrixArray(),
139  binsDet->GetNrows() - 1,
140  binsDet->GetMatrixArray());
141  } else {
142  fPdgCutQA = 0x0;
143  fEffCutQA = 0x0;
144  }
145 
146  // delete surplus vector
147  delete binsX;
148 
149  // Set labels to histograms
150  fCutNames[0][itype] = "no cuts";
151  if (fNCuts[kPrePair] > 1)
152  fCutNames[0][kTrack2] = "pair prefilter";
153  else
154  fCutNames[0][kTrack2] = "1st track filter";
155  // loop over all cuts
156  for (Int_t i = 0; i < fNCuts[itype]; i++) {
157  fCutQA->GetXaxis()->SetBinLabel(i + 1, fCutNames[i][itype]);
158  if (fPdgCutQA)
159  fPdgCutQA->GetXaxis()->SetBinLabel(i + 1, fCutNames[i][itype]);
160  if (fEffCutQA)
161  fEffCutQA->GetXaxis()->SetBinLabel(i + 1, fCutNames[i][itype]);
162  // printf(" itype:%s %d -> cut:%s \n",fTypeKeys[itype],itype,fCutNames[i][itype]);
163  }
164 
165  // pdg label
166  if (fPdgCutQA) {
167  TString pdglbl = "";
168  for (Int_t i = 0; i < binsPdg->GetNrows() - 1; i++) {
169  switch (i + 1) {
170  case 1: pdglbl = "electron"; break; // electron
171  case 2: pdglbl = "muon"; break; // muon
172  case 3: pdglbl = "pion"; break; // pion
173  case 4: pdglbl = "kaon"; break; // kaon
174  case 5: pdglbl = "proton"; break; // proton
175  }
176  fPdgCutQA->GetYaxis()->SetBinLabel(i + 1, pdglbl.Data());
177  }
178  }
179 
180  // detector label
181  if (fEffCutQA) {
182  TString detlbl = "";
183  for (Int_t i = 0; i < binsDet->GetNrows() - 1; i++) {
184  switch (i + 1) {
185  case 1:
187  break;
188  case 2:
190  break;
191  case 3:
193  break;
194  case 4:
196  break;
197  case 5:
199  break;
200  case 6:
202  break;
203  }
204  fEffCutQA->GetYaxis()->SetBinLabel(i + 1, detlbl.Data());
205  }
206  }
207 
208  // add to output list
209  switch (itype) {
210  case kEvent:
211  static_cast<THashList*>(fQAHistList.FindObject("Event"))
212  ->AddLast(fCutQA);
213  break;
214  case kTrack:
215  case kTrack2:
216  case kTrackMC:
217  static_cast<THashList*>(fQAHistList.FindObject("Track"))
218  ->AddLast(fCutQA);
219  if (fPdgCutQA)
220  static_cast<THashList*>(fQAHistList.FindObject("Track"))
221  ->AddLast(fPdgCutQA);
222  if (fEffCutQA)
223  static_cast<THashList*>(fQAHistList.FindObject("Track"))
224  ->AddLast(fEffCutQA);
225  break;
226  case kPair:
227  case kPrePair:
228  static_cast<THashList*>(fQAHistList.FindObject("Pair"))
229  ->AddLast(fCutQA);
230  break;
231  }
232  }
233 
234  // delete surplus
235  delete binsPdg;
236  delete binsDet;
237 }
238 
239 //_____________________________________________________________________
241  //
242  // add track filter cuts to the qa histogram
243  //
244 
245 
246  TIter listIterator(filter->GetCuts());
247  while (AnalysisCuts* thisCut = (AnalysisCuts*) listIterator()) {
248  Bool_t addCut = kTRUE;
249 
250  // add new cut class to the list
251  if (addCut) {
252  fCutNames[fNCuts[kTrack]][kTrack] = thisCut->GetTitle();
253  // printf("add cut %s to %d \n",thisCut->GetTitle(),fNCuts[kTrack]);
254  fNCuts[kTrack]++;
255  }
256 
257  } // pair filter loop
258 }
259 
260 //_____________________________________________________________________
262  //
263  // add MC track filter cuts to the qa histogram
264  //
265 
266 
267  TIter listIterator(filter->GetCuts());
268  while (AnalysisCuts* thisCut = (AnalysisCuts*) listIterator()) {
269  Bool_t addCut = kTRUE;
270 
271  // add new cut class to the list
272  if (addCut) {
273  fCutNames[fNCuts[kTrackMC]][kTrackMC] = thisCut->GetTitle();
274  // printf("add cut %s to %d \n",thisCut->GetTitle(),fNCuts[kTrack]);
275  fNCuts[kTrackMC]++;
276  }
277 
278  } // pair filter loop
279 }
280 
281 //_____________________________________________________________________
283  //
284  // add track filter cuts to the qa histogram
285  //
286  if (!filter) return;
287 
288  TIter listIterator(filter->GetCuts());
289  while (AnalysisCuts* thisCut = (AnalysisCuts*) listIterator()) {
290  Bool_t addCut = kTRUE;
291 
292  // add new cut class to the list
293  if (addCut) {
294  fCutNames[fNCuts[kTrack2]][kTrack2] = thisCut->GetTitle();
295  // printf("add cut %s to %d \n",thisCut->GetTitle(),fNCuts[kTrack]);
296  fNCuts[kTrack2]++;
297  }
298 
299  } // pair filter loop
300 }
301 
302 
303 //_____________________________________________________________________
305  //
306  // add track filter cuts to the qa histogram
307  //
308 
309  TIter listIterator(pairFilter->GetCuts());
310  while (AnalysisCuts* thisCut = (AnalysisCuts*) listIterator()) {
311  Bool_t addCut = kTRUE;
312 
313  // add new cut class to the list
314  if (addCut) {
315  fCutNames[fNCuts[kPair]][kPair] = thisCut->GetTitle();
316  // printf("add cut %s to %d \n",thisCut->GetTitle(),fNCuts[kPair]);
317  fNCuts[kPair]++;
318  }
319 
320  } // trk filter loop
321 }
322 
323 //_____________________________________________________________________
325  //
326  // add track filter cuts to the qa histogram
327  //
328  if (!pairFilter) return;
329 
330  TIter listIterator(pairFilter->GetCuts());
331  while (AnalysisCuts* thisCut = (AnalysisCuts*) listIterator()) {
332  Bool_t addCut = kTRUE;
333 
334  // add new cut class to the list
335  if (addCut) {
336  fCutNames[fNCuts[kPrePair]][kPrePair] = thisCut->GetTitle();
337  // printf("add cut %s to %d \n",thisCut->GetTitle(),fNCuts[kPair]);
338  fNCuts[kPrePair]++;
339  }
340 
341  } // trk filter loop
342 }
343 
344 //_____________________________________________________________________
346  //
347  // add track filter cuts to the qa histogram
348  //
349 
350 
351  TIter listIterator(eventFilter->GetCuts());
352  while (AnalysisCuts* thisCut = (AnalysisCuts*) listIterator()) {
353  Bool_t addCut = kTRUE;
354 
355  // add new cut class to the list
356  if (addCut) {
357  fCutNames[fNCuts[kEvent]][kEvent] = thisCut->GetTitle();
358  // printf("add cut %s to %d \n",thisCut->GetTitle(),fNCuts[kEvent]);
359  fNCuts[kEvent]++;
360  }
361 
362  } // trk filter loop
363 }
364 
365 //_____________________________________________________________________
366 void PairAnalysisCutQA::Fill(UInt_t mask, TObject* obj, UInt_t addIdx) {
367  //
368  // fill the corresponding step in the qa histogram
369  //
370 
371  UInt_t idx = GetObjIndex(obj) + addIdx;
372 
373  // pdg to pdg label
374  Int_t pdg = 0;
375  TString pdglbl = "";
376  if (idx == kTrack || idx == kTrack2) {
377  pdg = (Int_t)(static_cast<PairAnalysisTrack*>(obj)->PdgCode());
378  switch (TMath::Abs(pdg)) {
379  case 11: pdglbl = "electron"; break; // electron
380  case 13: pdglbl = "muon"; break; // muon
381  case 211: pdglbl = "pion"; break; // pion
382  case 321: pdglbl = "kaon"; break; // kaon
383  case 2212: pdglbl = "proton"; break; // proton
384  }
385  }
386 
387  // find histolist
388  THashList* histos = 0x0;
389  switch (idx) {
390  case kEvent:
391  histos = static_cast<THashList*>(fQAHistList.FindObject("Event"));
392  break;
393  case kTrack:
394  case kTrack2:
395  case kTrackMC:
396  histos = static_cast<THashList*>(fQAHistList.FindObject("Track"));
397  break;
398  case kPair:
399  case kPrePair:
400  histos = static_cast<THashList*>(fQAHistList.FindObject("Pair"));
401  break;
402  }
403 
404  // loop over cutmask and check decision
405  Int_t cutstep = 1;
406  for (Int_t iCut = 0; iCut < fNCuts[idx] - 1; ++iCut) {
407  // UInt_t cutMask=1<<iCut; // for each cut
408  UInt_t cutMask = (1 << (iCut + 1)) - 1; // increasing cut match
409 
410  // passed
411  if ((mask & cutMask) == cutMask) {
412  static_cast<TH1I*>(histos->FindObject(fTypeKeys[idx]))->Fill(cutstep);
413  if (pdg)
414  static_cast<TH2I*>(histos->FindObject(Form("%sPDG", fTypeKeys[idx])))
415  ->Fill(cutstep, pdglbl.Data(), 1.);
416 
417  // fill detector dependent
418  if (idx == kTrack || idx == kTrack2) {
419  TProfile2D* detQA = static_cast<TProfile2D*>(
420  histos->FindObject(Form("%sMatchEff", fTypeKeys[idx])));
421  PairAnalysisTrack* t = static_cast<PairAnalysisTrack*>(obj);
422  detQA->Fill(cutstep,
424  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kMvd))));
425  detQA->Fill(cutstep,
427  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kSts))));
428  detQA->Fill(cutstep,
430  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kRich))));
431  detQA->Fill(cutstep,
433  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kTrd))));
434  detQA->Fill(cutstep,
436  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kTof))));
437  detQA->Fill(cutstep,
439  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kMuch))));
440  }
441 
442  ++cutstep;
443  }
444  }
445 }
446 
447 //_____________________________________________________________________
448 void PairAnalysisCutQA::FillAll(TObject* obj, UInt_t addIdx) {
449  //
450  // fill the corresponding step in the qa histogram
451  //
452 
453  UInt_t idx = GetObjIndex(obj) + addIdx;
454 
455  // pdg to pdg label
456  Int_t pdg = 0;
457  TString pdglbl = "";
458  if (idx == kTrack || idx == kTrack2) {
459  pdg = (Int_t)(static_cast<PairAnalysisTrack*>(obj)->PdgCode());
460  switch (TMath::Abs(pdg)) {
461  case 11: pdglbl = "electron"; break; // electron
462  case 13: pdglbl = "muon"; break; // muon
463  case 211: pdglbl = "pion"; break; // pion
464  case 321: pdglbl = "kaon"; break; // kaon
465  case 2212: pdglbl = "proton"; break; // proton
466  }
467  }
468 
469  // find histolist
470  THashList* histos = 0x0;
471  switch (idx) {
472  case kEvent:
473  histos = static_cast<THashList*>(fQAHistList.FindObject("Event"));
474  break;
475  case kTrack:
476  case kTrack2:
477  case kTrackMC:
478  histos = static_cast<THashList*>(fQAHistList.FindObject("Track"));
479  break;
480  case kPair:
481  case kPrePair:
482  histos = static_cast<THashList*>(fQAHistList.FindObject("Pair"));
483  break;
484  }
485 
486  // fill
487  static_cast<TH1I*>(histos->FindObject(fTypeKeys[idx]))->Fill(0);
488  if (pdg)
489  static_cast<TH2I*>(histos->FindObject(Form("%sPDG", fTypeKeys[idx])))
490  ->Fill(0., pdglbl.Data(), 1.);
491 
492  // fill detector dependent
493  if (idx == kTrack || idx == kTrack2) {
494  TProfile2D* detQA = static_cast<TProfile2D*>(
495  histos->FindObject(Form("%sMatchEff", fTypeKeys[idx])));
496  PairAnalysisTrack* t = static_cast<PairAnalysisTrack*>(obj);
497  detQA->Fill(0.,
499  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kMvd))));
500  detQA->Fill(0.,
502  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kSts))));
503  detQA->Fill(0.,
505  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kRich))));
506  detQA->Fill(0.,
508  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kTrd))));
509  detQA->Fill(0.,
511  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kTof))));
512  detQA->Fill(0.,
514  t->TestBit(BIT(14 + ToIntegralType(ECbmModuleId::kMuch))));
515  }
516 }
517 
518 //______________________________________________________________________
519 UInt_t PairAnalysisCutQA::GetObjIndex(TObject* obj) {
520  //
521  // return the corresponding idex
522  //
523  // printf("INFO: object type is a %s \n", obj->IsA()->GetName());
524  if (obj->IsA() == CbmMCTrack::Class())
525  return kTrackMC;
526  else if (obj->IsA() == PairAnalysisTrack::Class())
527  return kTrack;
528  else if (obj->IsA() == PairAnalysisPairLV::Class())
529  return kPair;
530  else if (obj->IsA() == PairAnalysisPairKF::Class())
531  return kPair;
532  else if (obj->IsA() == PairAnalysisEvent::Class())
533  return kEvent;
534  else
535  fprintf(stderr,
536  "ERROR: object type not supported, please let the author know "
537  "about it\n"); //, obj->IsA()->GetName());
538  return -1;
539 }
PairAnalysisCutQA::kTrack
@ kTrack
Definition: PairAnalysisCutQA.h:27
PairAnalysisCutQA::kTrackMC
@ kTrackMC
Definition: PairAnalysisCutQA.h:27
PairAnalysisCutQA::PairAnalysisCutQA
PairAnalysisCutQA()
PairAnalysisHelper.h
PairAnalysisHelper::GetDetName
TString GetDetName(ECbmModuleId det)
Definition: PairAnalysisHelper.cxx:468
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
PairAnalysisCutQA::fCutNames
const char * fCutNames[20][kNtypes]
Definition: PairAnalysisCutQA.h:51
ECbmModuleId::kMvd
@ kMvd
Micro-Vertex Detector.
PairAnalysisTrack
Definition: PairAnalysisTrack.h:37
PairAnalysisHelper::MakeLinBinning
TVectorD * MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
Definition: PairAnalysisHelper.cxx:70
ECbmModuleId::kTof
@ kTof
Time-of-flight Detector.
PairAnalysisCutQA::fQAHistList
THashList fQAHistList
Definition: PairAnalysisCutQA.h:49
PairAnalysisCutQA::FillAll
void FillAll(TObject *obj, UInt_t addIdx=0)
Definition: PairAnalysisCutQA.cxx:448
PairAnalysisCutQA::AddPairFilter
void AddPairFilter(AnalysisFilter *pairFilter)
Definition: PairAnalysisCutQA.cxx:304
PairAnalysisCutQA::kNtypes
@ kNtypes
Definition: PairAnalysisCutQA.h:27
PairAnalysisCutQA::AddTrackFilterMC
void AddTrackFilterMC(AnalysisFilter *trkFilterMC)
Definition: PairAnalysisCutQA.cxx:261
PairAnalysisCutQA::GetObjIndex
UInt_t GetObjIndex(TObject *obj)
Definition: PairAnalysisCutQA.cxx:519
PairAnalysisPairLV.h
PairAnalysisCutQA::AddPrePairFilter
void AddPrePairFilter(AnalysisFilter *pairFilter)
Definition: PairAnalysisCutQA.cxx:324
AnalysisFilter::GetCuts
TList * GetCuts() const
Definition: AnalysisFilter.h:28
PairAnalysisCutQA::kPair
@ kPair
Definition: PairAnalysisCutQA.h:27
ECbmModuleId::kRich
@ kRich
Ring-Imaging Cherenkov Detector.
PairAnalysisCutGroup.h
PairAnalysisCutQA::kEvent
@ kEvent
Definition: PairAnalysisCutQA.h:27
PairAnalysisCutQA
Definition: PairAnalysisCutQA.h:24
PairAnalysisTrack.h
ClassImp
ClassImp(PairAnalysisCutQA) PairAnalysisCutQA
Definition: PairAnalysisCutQA.cxx:40
ECbmModuleId::kTrd
@ kTrd
Transition Radiation Detector.
PairAnalysisCutQA::kTrack2
@ kTrack2
Definition: PairAnalysisCutQA.h:27
CbmMCTrack.h
PairAnalysisCutQA::kPrePair
@ kPrePair
Definition: PairAnalysisCutQA.h:27
AnalysisFilter
Definition: AnalysisFilter.h:15
PairAnalysisCutQA::AddTrackFilter2
void AddTrackFilter2(AnalysisFilter *trkFilter2)
Definition: PairAnalysisCutQA.cxx:282
PairAnalysisCutQA::AddEventFilter
void AddEventFilter(AnalysisFilter *eventFilter)
Definition: PairAnalysisCutQA.cxx:345
ToIntegralType
constexpr auto ToIntegralType(T enumerator) -> typename std::underlying_type< T >::type
Definition: CbmDefs.h:24
PairAnalysisTrack::PdgCode
Int_t PdgCode() const
Definition: PairAnalysisTrack.h:108
ECbmModuleId::kMuch
@ kMuch
Muon detection system.
PairAnalysisCutQA::~PairAnalysisCutQA
virtual ~PairAnalysisCutQA()
Definition: PairAnalysisCutQA.cxx:72
ECbmModuleId::kSts
@ kSts
Silicon Tracking System.
PairAnalysisCutQA.h
PairAnalysisCutQA::fTypeKeys
const char * fTypeKeys[kNtypes]
Definition: PairAnalysisCutQA.h:52
PairAnalysisCutQA::AddTrackFilter
void AddTrackFilter(AnalysisFilter *trkFilter)
Definition: PairAnalysisCutQA.cxx:240
PairAnalysisEvent.h
AnalysisCuts
Definition: AnalysisCuts.h:12
AnalysisCuts.h
PairAnalysisCutQA::Fill
void Fill(UInt_t mask, TObject *obj, UInt_t addIdx=0)
Definition: PairAnalysisCutQA.cxx:366
PairAnalysisCutQA::fNCuts
Int_t fNCuts[kNtypes]
Definition: PairAnalysisCutQA.h:50
PairAnalysisCutQA::Init
void Init()
Definition: PairAnalysisCutQA.cxx:80
PairAnalysisPairKF.h