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