BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnPidQa.cxx
Go to the documentation of this file.
1/**
2 * \file BmnPidQa.cxx
3 * \author Andrey Lebedev <andrey.lebedev@gsi.de> - original author for CBM experiment
4 * \author Sergey Merts <sergey.merts@gmail.com> - modifications for BMN experiment
5 * \date 2007-2021
6 */
7
8#include "BmnPidQa.h"
9
10#include "BmnAcceptanceFunction.h"
11#include "BmnDchHit.h"
12#include "BmnEnums.h"
13#include "BmnGemStripHit.h"
14#include "BmnGemTrack.h"
15#include "BmnGlobalTrack.h"
16#include "BmnHistManager.h"
17#include "BmnMCPoint.h"
18#include "BmnMatch.h"
19#include "BmnMath.h"
20#include "BmnPidQaReport.h"
21#include "BmnSiliconHit.h"
22#include "BmnTrackMatch.h"
23#include "BmnUtils.h"
24#include "BmnVertex.h"
25#include "CbmBaseHit.h"
26#include "CbmGlobalTrack.h"
27#include "CbmMCTrack.h"
28#include "CbmStsTrack.h"
29#include "CbmTofHit.h"
30#include "FairMCEventHeader.h"
31#include "FairMCPoint.h"
32#include "FairRunAna.h"
33#include "TClonesArray.h"
34#include "TF1.h"
35#include "TFitResult.h"
36#include "TH1.h"
37#include "TH2F.h"
38
39#include <fstream>
40#include <iostream>
41
42using namespace std;
43using namespace TMath;
44using lit::FindAndReplace;
45using lit::Split;
46
47BmnPidQa::BmnPidQa(TString name, TString storageName)
48 : FairTask("BmnPidQA", 1)
49 , fOutName(name)
50 , fStorageName(storageName)
51 , fHM(NULL)
52 , fOutputDir("./")
53 , fPrimes(kFALSE)
54 , fPRangeMin(0.)
55 , fPRangeMax(6.)
56 , fPRangeBins(500)
57 , fBetaRangeMin(0.)
58 , fBetaRangeMax(1.1)
59 , fBetaRangeBins(500)
60 , fVelocRangeMin(1.5e+8)
61 , fVelocRangeMax(3.2e+8)
62 , fVelocRangeBins(500)
63 , fTimeRangeMin(0.)
64 , fTimeRangeMax(50)
65 , fTimeRangeBins(500)
66 , fMassRangeMin(0.)
67 , fMassRangeMax(1)
68 , fMassRangeBins(100)
69 , fNHitsRangeMin(0)
70 , fNHitsRangeMax(15)
71 , fNHitsRangeBins(500)
72 , fMCTracks(NULL)
73 , fGlobalTracks(NULL)
74{
75 db = TDatabasePDG::Instance();
76 // Add several types of particles to the database.
77 // p, pi, K, e are already in it
78 db->AddParticle("D", "D", 1.876123928, true, 0, 3, "Core", 1000010020, 1000010020, 1000010020);
79 db->AddParticle("T", "T", 2.809432115, true, 0, 3, "Core", 1000010030, 1000010030, 1000010030);
80 db->AddParticle("He3", "He3", 2.809413523, true, 0, 6, "Core", 1000020030, 1000020030, 1000020030);
81 db->AddParticle("He4", "He4", 3.728401326, true, 0, 6, "Core", 1000020040, 1000020040, 1000020040);
82
83 for (PidParticles iSort = (PidParticles)0; iSort != EndPidEnum; iSort = (PidParticles)(iSort + 1)) {
84 Int_t pdg = EnumToPdg(iSort);
85 TParticlePDG* iParticle = db->GetParticle(pdg);
86 fParticles.push_back(iParticle);
87 Double_t mass = iParticle->Mass();
88 if (fMassTable.count(mass) == 0) {
89 fMassTable.insert(pair<Double_t, string>(mass, iParticle->GetName()));
90 fPidStatistics400.insert(pair<string, vector<Int_t>>(iParticle->GetName(), vector<Int_t>(4, 0)));
91 }
92 }
93 fPidStatistics700 = fPidStatistics400;
94}
95
96BmnPidQa::~BmnPidQa()
97{
98 if (fHM)
99 delete fHM;
100}
101
102InitStatus BmnPidQa::Init()
103{
104 fHM = new BmnHistManager();
105 CreateHistograms();
106 ReadDataBranches();
107 return kSUCCESS;
108}
109
110void BmnPidQa::Exec(Option_t* opt)
111{
112 ProcessGlobal();
113}
114
115void BmnPidQa::Finish()
116{
117
118 fHM->WriteToFile();
119 BmnSimulationReport* report = new BmnPidQaReport(fOutName, fStorageName, fMassTable);
120 report->SetOnlyPrimes(fPrimes);
121 report->Create(fHM, fOutputDir);
122 delete report;
123
124 MassTablePrint();
125 cout << endl;
126 cout << "Name" << '\t' << "Total" << '\t' << "True" << '\t' << "False" << endl;
127 cout << "TOF400 stat" << endl;
128 PidStatisticsPrint400();
129 cout << endl;
130 cout << "TOF700 stat" << endl;
131 PidStatisticsPrint700();
132}
133
134void BmnPidQa::ReadDataBranches()
135{
136
137 FairRootManager* ioman = FairRootManager::Instance();
138 cout << "ReadDataBranches()" << endl;
139 if (NULL == ioman)
140 Fatal("Init", "BmnRootManager is not instantiated");
141
142 fMCTracks = (TClonesArray*)ioman->GetObject("MCTrack");
143 if (NULL == fMCTracks)
144 Fatal("Init", "No MCTrack array!");
145
146 fGlobalTracks = (TClonesArray*)ioman->GetObject("BmnGlobalTrack");
147 fGlobalTrackMatches = (TClonesArray*)ioman->GetObject("BmnGlobalTrackMatch");
148
149 fTof400Hits = (TClonesArray*)ioman->GetObject("BmnTof400Hit");
150 fTof700Hits = (TClonesArray*)ioman->GetObject("BmnTof700Hit");
151}
152
153void BmnPidQa::CreateH1(const string& name,
154 const string& xTitle,
155 const string& yTitle,
156 Int_t nofBins,
157 Double_t minBin,
158 Double_t maxBin)
159{
160 TH1F* h = new TH1F(name.c_str(), string(name + ";" + xTitle + ";" + yTitle).c_str(), nofBins, minBin, maxBin);
161 fHM->Add(name, h);
162}
163
164void BmnPidQa::CreateH2(const string& name,
165 const string& xTitle,
166 const string& yTitle,
167 const string& zTitle,
168 Int_t nofBinsX,
169 Double_t minBinX,
170 Double_t maxBinX,
171 Int_t nofBinsY,
172 Double_t minBinY,
173 Double_t maxBinY)
174{
175 TH2F* h = new TH2F(name.c_str(), (name + ";" + xTitle + ";" + yTitle + ";" + zTitle).c_str(), nofBinsX, minBinX,
176 maxBinX, nofBinsY, minBinY, maxBinY);
177 fHM->Add(name, h);
178}
179
180void BmnPidQa::CreateTrackHitsHistogram(const string& detName)
181{
182 string type[] = {"All", "True", "Fake", "TrueOverAll", "FakeOverAll"};
183 Double_t min[] = {0., 0., 0., 0., 0.};
184 Double_t max[] = {20, 20, 20, 1., 1.};
185 Int_t bins[] = {20, 20, 20, 20, 20};
186 for (Int_t i = 0; i < 5; i++) {
187 string xTitle = (i == 3 || i == 4) ? "Ratio" : "Number of hits";
188 string histName = "hth_" + detName + "_TrackHits_" + type[i];
189 CreateH1(histName.c_str(), xTitle, "Yeild", bins[i], min[i], max[i]);
190 }
191}
192
193void BmnPidQa::CreateHistograms()
194{
195
196 cout << "CreateHistograms()" << endl;
197 CreateH2("Banana-plot TOF-400", "TOF-400 P_{rec}/q, GeV/c/e", "Beta, c", "N", fPRangeBins * 3, fPRangeMin,
198 fPRangeMax, fBetaRangeBins * 2, fBetaRangeMin, fBetaRangeMax);
199 CreateH2("Banana-plot TOF-700", "TOF-700 P_{rec}/q, GeV/c/e", "Beta, c", "N", fPRangeBins * 3, fPRangeMin,
200 fPRangeMax, fBetaRangeBins * 2, fBetaRangeMin, fBetaRangeMax);
201
202 for (auto iter = fMassTable.begin(); iter != fMassTable.end(); ++iter) {
203 string nameOfParticle = (*iter).second;
204 cout << "Particle name is " << nameOfParticle << endl;
205
206 // Momentum histograms
207
208 // Velocity hist
209 CreateH2("Total velocity from P TOF-400 for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Velocity, m/s",
210 "N", fPRangeBins * 3, fPRangeMin, fPRangeMax, fVelocRangeBins, fVelocRangeMin, fVelocRangeMax);
211
212 CreateH2("Total velocity from P TOF-700 for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Velocity, m/s",
213 "N", fPRangeBins * 3, fPRangeMin, fPRangeMax, fVelocRangeBins, fVelocRangeMin, fVelocRangeMax);
214
215 // Time hist
216
217 CreateH2("Total time from P TOF-400 for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Time, ns", "N",
218 fPRangeBins * 3, fPRangeMin, fPRangeMax, fTimeRangeBins, fTimeRangeMin, fTimeRangeMax);
219
220 CreateH2("Total time from P TOF-700 for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Time, ns", "N",
221 fPRangeBins * 3, fPRangeMin, fPRangeMax, fTimeRangeBins, fTimeRangeMin, fTimeRangeMax);
222
223 //
224 CreateH1("TOF-400 total_vs_P for " + nameOfParticle, "True total, P_{rec}/q, GeV/c/e", "N", fPRangeBins,
225 fPRangeMin, fPRangeMax);
226 CreateH1("TOF-700 total_vs_P for " + nameOfParticle, "True total, P_{rec}/q, GeV/c", "N", fPRangeBins,
227 fPRangeMin, fPRangeMax);
228 // Histogram stores rate of true-reconstructions
229 CreateH1("TOF-400 true_vs_P for " + nameOfParticle, "Pid hits, P_{rec}/q, GeV/c/e", "N", fPRangeBins,
230 fPRangeMin, fPRangeMax);
231 CreateH1("TOF-700 true_vs_P for " + nameOfParticle, "Pid hits, P_{rec}/q, GeV/c/e", "N", fPRangeBins,
232 fPRangeMin, fPRangeMax);
233 // Histogram stores rate of wrong-reconstructions
234 CreateH1("TOF-400 false_vs_P for " + nameOfParticle, "Pid miss, P_{rec}/q, GeV/c/e", "N", fPRangeBins,
235 fPRangeMin, fPRangeMax);
236 CreateH1("TOF-700 false_vs_P for " + nameOfParticle, "Pid miss, P_{rec}/q, GeV/c/e", "N", fPRangeBins,
237 fPRangeMin, fPRangeMax);
238 // True + False hist for contamination of PID
239 CreateH1("TOF-400 true-false_vs_P for " + nameOfParticle, "Contamination, P_{rec}/q, GeV/c/e", "tof400, %",
240 fPRangeBins, fPRangeMin, fPRangeMax);
241 CreateH1("TOF-700 true-false_vs_P for " + nameOfParticle, "Contamination, P_{rec}/q, GeV/c/e", "tof700, %",
242 fPRangeBins, fPRangeMin, fPRangeMax);
243
244 // Rigidity histograms
245 //
246 // Velocity true/false hist
247 CreateH2("True velocity from P TOF-400 for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Velocity, m/s",
248 "N", fPRangeBins * 3, fPRangeMin, fPRangeMax, fVelocRangeBins, fVelocRangeMin, fVelocRangeMax);
249
250 CreateH2("True velocity from P TOF-700 for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Velocity, m/s",
251 "N", fPRangeBins * 3, fPRangeMin, fPRangeMax, fVelocRangeBins, fVelocRangeMin, fVelocRangeMax);
252
253 CreateH2("False velocity from P TOF-400 for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Velocity, m/s",
254 "N", fPRangeBins * 3, fPRangeMin, fPRangeMax, fVelocRangeBins, fVelocRangeMin, fVelocRangeMax);
255
256 CreateH2("False velocity from P TOF-700 for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Velocity, m/s",
257 "N", fPRangeBins * 3, fPRangeMin, fPRangeMax, fVelocRangeBins, fVelocRangeMin, fVelocRangeMax);
258
259 // Time true/false hist
260
261 CreateH2("True time from P TOF-400 for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Time, ns", "N",
262 fPRangeBins * 3, fPRangeMin, fPRangeMax, fTimeRangeBins, fTimeRangeMin, fTimeRangeMax);
263
264 CreateH2("True time from P TOF-700 for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Time, ns", "N",
265 fPRangeBins * 3, fPRangeMin, fPRangeMax, fTimeRangeBins, fTimeRangeMin, fTimeRangeMax);
266
267 CreateH2("False time from P TOF-400 for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Time, ns", "N",
268 fPRangeBins * 3, fPRangeMin, fPRangeMax, fTimeRangeBins, fTimeRangeMin, fTimeRangeMax);
269
270 CreateH2("False time from P TOF-700 for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Time, ns", "N",
271 fPRangeBins * 3, fPRangeMin, fPRangeMax, fTimeRangeBins, fTimeRangeMin, fTimeRangeMax);
272
273 // Mass^2 hist
274
275 fMassRangeMax = (*iter).first * (*iter).first * 3;
276
277 CreateH2("Total mass^2 from P TOF-400 for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Mass, GeV^2", "N",
278 fPRangeBins * 3, fPRangeMin, fPRangeMax, fMassRangeBins, fMassRangeMin, fMassRangeMax);
279
280 CreateH2("Total mass^2 from P TOF-700 for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Mass, GeV^2", "N",
281 fPRangeBins * 3, fPRangeMin, fPRangeMax, fMassRangeBins, fMassRangeMin, fMassRangeMax);
282
283 CreateH2("True mass^2 from P TOF-400 for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Mass, GeV^2", "N",
284 fPRangeBins * 3, fPRangeMin, fPRangeMax, fMassRangeBins, fMassRangeMin, fMassRangeMax);
285
286 CreateH2("True mass^2 from P TOF-700 for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Mass, GeV^2", "N",
287 fPRangeBins * 3, fPRangeMin, fPRangeMax, fMassRangeBins, fMassRangeMin, fMassRangeMax);
288
289 CreateH2("False mass^2 from P TOF-400 for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Mass, GeV^2", "N",
290 fPRangeBins * 3, fPRangeMin, fPRangeMax, fMassRangeBins, fMassRangeMin, fMassRangeMax);
291
292 CreateH2("False mass^2 from P TOF-700 for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Mass, GeV^2", "N",
293 fPRangeBins * 3, fPRangeMin, fPRangeMax, fMassRangeBins, fMassRangeMin, fMassRangeMax);
294
295 CreateH2("TOF-400 total rigidity-momentum for " + nameOfParticle, "P_{rec}/q, GeV/c/e", "Beta, c", "N_{total}",
296 fPRangeBins * 3, fPRangeMin, fPRangeMax, fBetaRangeBins * 2, fBetaRangeMin, fBetaRangeMax);
297 CreateH2("TOF-700 total rigidity-momentum for " + nameOfParticle, "P_{rec}/q, GeV/c/e", "Beta, c", "N_{total}",
298 fPRangeBins * 3, fPRangeMin, fPRangeMax, fBetaRangeBins * 2, fBetaRangeMin, fBetaRangeMax);
299 CreateH2("TOF-400 true rigidity-momentum for " + nameOfParticle, "P_{rec}/q, GeV/c/e", "Beta, c", "N_{true}",
300 fPRangeBins * 3, fPRangeMin, fPRangeMax, fBetaRangeBins * 2, fBetaRangeMin, fBetaRangeMax);
301 CreateH2("TOF-700 true rigidity-momentum for " + nameOfParticle, "P_{rec}/q, GeV/c/e", "Beta, c", "N_{true}",
302 fPRangeBins * 3, fPRangeMin, fPRangeMax, fBetaRangeBins * 2, fBetaRangeMin, fBetaRangeMax);
303 CreateH2("TOF-400 false rigidity-momentum for " + nameOfParticle, "P_{rec}/q, GeV/c/e", "Beta, c", "N_{false}",
304 fPRangeBins * 3, fPRangeMin, fPRangeMax, fBetaRangeBins * 2, fBetaRangeMin, fBetaRangeMax);
305 CreateH2("TOF-700 false rigidity-momentum for " + nameOfParticle, "P_{rec}/q, GeV/c/e", "Beta, c", "N_{false}",
306 fPRangeBins * 3, fPRangeMin, fPRangeMax, fBetaRangeBins * 2, fBetaRangeMin, fBetaRangeMax);
307
308 /*
309 //Species undivided in time
310 CreateH2("Undivided TOF-400 rigidity-momentum for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Beta, c",
311 "N", fPRangeBins*3, fPRangeMin, fPRangeMax, fBetaRangeBins*2, fBetaRangeMin, fBetaRangeMax);
312
313 CreateH2("Undivided TOF-400 time from P for " + nameOfParticle, "TOF-400 P_{rec}/q, GeV/c/e", "Time, ns", "N",
314 fPRangeBins*3, fPRangeMin, fPRangeMax, fTimeRangeBins, fTimeRangeMin, fTimeRangeMax);
315
316 CreateH2("Undivided TOF-700 rigidity-momentum for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Beta, c",
317 "N", fPRangeBins*3, fPRangeMin, fPRangeMax, fBetaRangeBins*2, fBetaRangeMin, fBetaRangeMax);
318
319 CreateH2("Undivided TOF-700 time from P for " + nameOfParticle, "TOF-700 P_{rec}/q, GeV/c/e", "Time, ns", "N",
320 fPRangeBins*3, fPRangeMin, fPRangeMax, fTimeRangeBins, fTimeRangeMin, fTimeRangeMax);
321
322 */
323 // Number of hits histograms
324 //
325 CreateH1("TOF-400 total_vs_NOfHits for " + nameOfParticle, "True total, N_{of hits}", "N", fNHitsRangeBins,
326 fNHitsRangeMin, fNHitsRangeMax);
327 CreateH1("TOF-700 total_vs_NOfHits for " + nameOfParticle, "True total, N_{of hits}", "N", fNHitsRangeBins,
328 fNHitsRangeMin, fNHitsRangeMax);
329 // Histogram stores rate of true-reconstructions
330 CreateH1("TOF-400 true_vs_NOfHits for " + nameOfParticle, "Pid hits, N_{of hits}", "N", fNHitsRangeBins,
331 fNHitsRangeMin, fNHitsRangeMax);
332 CreateH1("TOF-700 true_vs_NOfHits for " + nameOfParticle, "Pid hits, N_{of hits}", "N", fNHitsRangeBins,
333 fNHitsRangeMin, fNHitsRangeMax);
334 // Histogram stores rate of wrong-reconstructions
335 CreateH1("TOF-400 false_vs_NOfHits for " + nameOfParticle, "Pid miss, N_{of hits}", "N", fNHitsRangeBins,
336 fNHitsRangeMin, fNHitsRangeMax);
337 CreateH1("TOF-700 false_vs_NOfHits for " + nameOfParticle, "Pid miss, N_{of hits}", "N", fNHitsRangeBins,
338 fNHitsRangeMin, fNHitsRangeMax);
339 // True + False hist for contamination of PID
340 CreateH1("TOF-400 true-false_vs_NOfHits for " + nameOfParticle, "Contamination, N_{of hits}", "tof400, %",
341 fNHitsRangeBins, fNHitsRangeMin, fNHitsRangeMax);
342 CreateH1("TOF-700 true-false_vs_NOfHits for " + nameOfParticle, "Contamination, N_{of hits}", "tof700, %",
343 fNHitsRangeBins, fNHitsRangeMin, fNHitsRangeMax);
344 }
345}
346
347void BmnPidQa::ProcessGlobal()
348{
349 for (Int_t iTrack = 0; iTrack < fGlobalTracks->GetEntriesFast(); iTrack++) {
350 // Reciving of reco tracks
351 BmnGlobalTrack* glTrack = (BmnGlobalTrack*)(fGlobalTracks->At(iTrack));
352 if (!glTrack)
353 continue;
354
355 // Banana-plots prepearing
356 if (glTrack->GetBeta(1) > -999.0)
357 fHM->H2("Banana-plot TOF-400")->Fill(glTrack->GetP(), glTrack->GetBeta(1));
358
359 if (glTrack->GetBeta(2) > -999.0)
360 fHM->H2("Banana-plot TOF-700")->Fill(glTrack->GetP(), glTrack->GetBeta(2));
361
362 // Reciving of trackmatches
363 BmnTrackMatch* globTrackMatch = (BmnTrackMatch*)(fGlobalTrackMatches->At(iTrack));
364 if (!globTrackMatch)
365 continue;
366 if (globTrackMatch->GetNofLinks() == 0)
367 continue;
368
369 // Getting of most possible MCTrack index corresponding to the reco track
370 Int_t globMCId = globTrackMatch->GetMatchedLink().GetIndex();
371
372 // Getting of MCTrack by match index
373 CbmMCTrack* mcTrack = (CbmMCTrack*)(fMCTracks->At(globMCId));
374 if (!mcTrack)
375 continue;
376 if (!(mcTrack->GetMotherId() == -1) && GetOnlyPrimes()) {
377 continue;
378 }
379
380 Double_t time(0);
381 Double_t length(0);
382 Double_t velocity(0);
383 Int_t mcPdg = mcTrack->GetPdgCode();
384
385 // TOF400:
386 // Comparation of particle masses in MC and Reco tracks
387 if ((glTrack->GetBeta(1) > -999.0)) {
388
389 Int_t recoPdg1 = EnumToPdg(glTrack->GetParticleTof400());
390 Int_t indexTOF400 = glTrack->GetTof1HitIndex();
391 Double_t mass = glTrack->GetMass2(1);
392
393 if (indexTOF400 != -1) {
394 BmnHit* hit = (BmnHit*)fTof400Hits->At(indexTOF400);
395 length = (hit->GetLength());
396 time = (hit->GetTimeStamp());
397 velocity = length * 1e+7 / time; // length in cm, time in ns => m/s
398 }
399 // Calcuation of corresponding histogram momentum bin
400
401 Double_t p = glTrack->GetP();
402 Int_t nOfHits = glTrack->GetNHits();
403 Double_t beta = glTrack->GetBeta(1);
404 TParticlePDG* recoParticle = GetParticleExtend(recoPdg1);
405 TParticlePDG* mcParticle = GetParticleExtend(mcPdg);
406 if ((mcParticle == 0) || (recoParticle == 0))
407 continue;
408 Double_t recoMass = recoParticle->Mass();
409 Double_t mcMass = mcParticle->Mass();
410 string recoName = GetParticleName(recoMass);
411 string mcName = GetParticleName(mcMass);
412
413 // Hist for undivided in time TOF400
414 /*
415 if((recoPdg1==0)&&((fMassTable.count(mcMass))!=0)){
416 ++fPidStatistics400[mcName][3]; // ++ undivided;
417 fHM->H2("Undivided TOF-400 rigidity-momentum for " + mcName)->Fill(p, beta);
418
419 fHM->H2("Undivided TOF-400 time from P for " + mcName)->Fill(p, time);
420 }*/
421
422 if (abs(mcMass - recoMass) < 0.000001) {
423 if (fMassTable.count(mcMass) != 0) {
424 ++fPidStatistics400[mcName][1]; // +true pid;
425 fHM->H1("TOF-400 true_vs_P for " + mcName)->Fill(p);
426 fHM->H1("TOF-400 true_vs_NOfHits for " + mcName)->Fill(nOfHits);
427 fHM->H2("TOF-400 true rigidity-momentum for " + mcName)->Fill(p, beta);
428 if (indexTOF400 != -1) {
429 fHM->H2("True velocity from P TOF-400 for " + mcName)->Fill(p, velocity);
430 fHM->H2("True time from P TOF-400 for " + mcName)->Fill(p, time);
431 fHM->H2("True mass^2 from P TOF-400 for " + mcName)->Fill(p, mass);
432 }
433 }
434 } else {
435 if (fMassTable.count(recoMass) != 0) {
436 ++fPidStatistics400[recoName][2]; // +false pid;
437 fHM->H1("TOF-400 false_vs_P for " + recoName)->Fill(p);
438 fHM->H1("TOF-400 false_vs_NOfHits for " + recoName)->Fill(nOfHits);
439 fHM->H2("TOF-400 false rigidity-momentum for " + recoName)->Fill(p, beta);
440 if (indexTOF400 != -1) {
441 fHM->H2("False velocity from P TOF-400 for " + recoName)->Fill(p, velocity);
442 fHM->H2("False time from P TOF-400 for " + recoName)->Fill(p, time);
443 fHM->H2("False mass^2 from P TOF-400 for " + recoName)->Fill(p, mass);
444 }
445 }
446 }
447 if (fMassTable.count(mcMass) != 0) {
448 ++fPidStatistics400[mcName][0]; // +total true particle number (mc);
449 fHM->H1("TOF-400 total_vs_P for " + mcName)->Fill(p);
450 fHM->H1("TOF-400 total_vs_NOfHits for " + mcName)->Fill(nOfHits);
451 fHM->H2("TOF-400 total rigidity-momentum for " + mcName)->Fill(p, beta);
452 if (indexTOF400 != -1) {
453 fHM->H2("Total velocity from P TOF-400 for " + mcName)->Fill(p, velocity);
454 fHM->H2("Total time from P TOF-400 for " + mcName)->Fill(p, time);
455 fHM->H2("Total mass^2 from P TOF-400 for " + mcName)->Fill(p, mass);
456 }
457 }
458 }
459
460 // TOF700:
461
462 if ((glTrack->GetBeta(2) > -999.0)) {
463
464 Int_t recoPdg2 = EnumToPdg(glTrack->GetParticleTof700());
465 Double_t mass = glTrack->GetMass2(2);
466 Int_t indexTOF700 = glTrack->GetTof2HitIndex();
467
468 if (indexTOF700 != -1) {
469 BmnHit* hit = (BmnHit*)fTof700Hits->At(indexTOF700);
470 length = (hit->GetLength());
471 time = (hit->GetTimeStamp());
472 velocity = length * 1e+7 / time; // length in cm, time in ns => m/s
473 }
474
475 // Calcuation of corresponding historgram momentum bin
476 Double_t p = glTrack->GetP();
477 Int_t nOfHits = glTrack->GetNHits();
478 Double_t beta = glTrack->GetBeta(2);
479 TParticlePDG* recoParticle = GetParticleExtend(recoPdg2);
480 TParticlePDG* mcParticle = GetParticleExtend(mcPdg);
481 if ((mcParticle == 0) || (recoParticle == 0))
482 continue;
483 Double_t recoMass = recoParticle->Mass();
484 Double_t mcMass = mcParticle->Mass();
485 string recoName = GetParticleName(recoMass);
486 string mcName = GetParticleName(mcMass);
487
488 // Hist for undivided in time TOF700
489 /* if((recoPdg2==0)&&((fMassTable.count(mcMass))!=0)) {
490 ++fPidStatistics700[mcName][3]; // ++ undivided;
491 fHM->H2("Undivided TOF-700 rigidity-momentum for " + mcName)->Fill(p, beta);
492 fHM->H2("Undivided TOF-700 time from P for " + mcName)->Fill(p, time);
493 }*/
494
495 if (abs(mcMass - recoMass) < 0.000001) {
496 if (fMassTable.count(mcMass) != 0) {
497 ++fPidStatistics700[mcName][1]; // +true pid;
498 fHM->H1("TOF-700 true_vs_P for " + mcName)->Fill(p);
499 fHM->H1("TOF-700 true_vs_NOfHits for " + mcName)->Fill(nOfHits);
500 fHM->H2("TOF-700 true rigidity-momentum for " + mcName)->Fill(p, beta);
501 if (indexTOF700 != -1) {
502 fHM->H2("True velocity from P TOF-700 for " + mcName)->Fill(p, velocity);
503 fHM->H2("True time from P TOF-700 for " + mcName)->Fill(p, time);
504 fHM->H2("True mass^2 from P TOF-700 for " + mcName)->Fill(p, mass);
505 }
506 }
507 } else {
508 if (fMassTable.count(recoMass) != 0) {
509 ++fPidStatistics700[recoName][2]; // +false pid;
510 fHM->H1("TOF-700 false_vs_P for " + recoName)->Fill(p);
511 fHM->H1("TOF-700 false_vs_NOfHits for " + recoName)->Fill(nOfHits);
512 fHM->H2("TOF-700 false rigidity-momentum for " + recoName)->Fill(p, beta);
513 if (indexTOF700 != -1) {
514 fHM->H2("False velocity from P TOF-700 for " + recoName)->Fill(p, velocity);
515 fHM->H2("False time from P TOF-700 for " + recoName)->Fill(p, time);
516 fHM->H2("False mass^2 from P TOF-700 for " + recoName)->Fill(p, mass);
517 }
518 }
519 }
520 if (fMassTable.count(mcMass) != 0) {
521 ++fPidStatistics700[mcName][0]; // +total true particle number (mc);
522 fHM->H1("TOF-700 total_vs_P for " + mcName)->Fill(p);
523 fHM->H1("TOF-700 total_vs_NOfHits for " + mcName)->Fill(nOfHits);
524 fHM->H2("TOF-700 total rigidity-momentum for " + (mcName))->Fill(p, beta);
525 if (indexTOF700 != -1) {
526 fHM->H2("Total velocity from P TOF-700 for " + mcName)->Fill(p, velocity);
527 fHM->H2("Total time from P TOF-700 for " + mcName)->Fill(p, time);
528 fHM->H2("Total mass^2 from P TOF-700 for " + mcName)->Fill(p, mass);
529 }
530 }
531 }
532 }
533}
534
535Int_t BmnPidQa::EnumToPdg(PidParticles part)
536{
537 switch (part) {
538 case PidProton:
539 return 2212;
540 case PidPion:
541 return 211;
542 case PidKaon:
543 return 321;
544 case PidElectron:
545 return 11;
546 case PidDeuteron:
547 return 1000010020;
548 case PidTriton:
549 return 1000010030;
550 case PidHelium3:
551 return 1000020030;
552 case PidHelium4:
553 return 1000020040;
554 default:
555 cout << "No such particle in PidParticles\n";
556 return -1;
557 }
558}
559
560TParticlePDG* BmnPidQa::GetParticleExtend(Int_t pdgCode)
561{
562 for (auto sort = fParticles.begin(); sort != fParticles.end(); ++sort) {
563 if ((*sort)->PdgCode() == pdgCode)
564 return *sort;
565 }
566
567 return db->GetParticle(pdgCode);
568}
569
570string BmnPidQa::GetParticleName(Double_t mass)
571{
572 for (auto iter = fMassTable.begin(); iter != fMassTable.end(); ++iter) {
573 Double_t tableMass = iter->first;
574 if (abs(tableMass - mass) < 0.000001)
575 return iter->second;
576 }
577 return "NULL";
578}
579
580void BmnPidQa::MassTablePrint()
581{
582 for (auto iter = fMassTable.begin(); iter != fMassTable.end(); ++iter) {
583 cout << iter->first << ": " << iter->second << " " << endl;
584 }
585}
586void BmnPidQa::PidStatisticsPrint400()
587{
588 for (auto iter = fPidStatistics400.begin(); iter != fPidStatistics400.end(); ++iter)
589 cout << iter->first << " " << iter->second[0] << " " << iter->second[1] << " " << iter->second[2] << endl;
590}
591void BmnPidQa::PidStatisticsPrint700()
592{
593 for (auto iter = fPidStatistics700.begin(); iter != fPidStatistics700.end(); ++iter)
594 cout << iter->first << " " << iter->second[0] << " " << iter->second[1] << " " << iter->second[2] << endl;
595}
BmnPidQaReport.h
Create report for tracking QA.
BmnPidQa.h
FairTask for pid performance calculation.
i
int i
Definition P4_F32vec4.h:22
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:30
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:31
BmnAcceptanceFunction.h
Global function to define the track acceptance. Used in QA.
BmnMCPoint.h
Monte-Carlo point.
Int_t
int Int_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
Option_t
const char Option_t
Fatal
void Fatal(const char *location, const char *msgfmt,...)
BmnGlobalTrack
BmnGlobalTrack.h.
Definition BmnGlobalTrack.h:14
BmnGlobalTrack::GetTof2HitIndex
Int_t GetTof2HitIndex() const
Definition BmnGlobalTrack.h:37
BmnGlobalTrack::GetTof1HitIndex
Int_t GetTof1HitIndex() const
Definition BmnGlobalTrack.h:36
BmnGlobalTrack::GetMass2
Float_t GetMass2(Int_t tofID) const
Definition BmnGlobalTrack.h:43
BmnHistManager
Histogram manager.
Definition BmnHistManager.h:42
BmnHit::GetLength
Double_t GetLength()
Definition BmnHit.h:33
BmnMatch::GetMatchedLink
const BmnLink & GetMatchedLink() const
Definition BmnMatch.h:40
BmnMatch::GetNofLinks
Int_t GetNofLinks() const
Definition BmnMatch.h:41
BmnPidQaReport
Create report for pid QA.
Definition BmnPidQaReport.h:29
BmnPidQa::GetParticleName
string GetParticleName(Double_t mass)
Definition BmnPidQa.cxx:570
BmnPidQa::BmnPidQa
BmnPidQa()
Constructor.
Definition BmnPidQa.h:35
BmnPidQa::~BmnPidQa
virtual ~BmnPidQa()
Destructor.
Definition BmnPidQa.cxx:96
BmnPidQa::fNHitsRangeMin
Double_t fNHitsRangeMin
Definition BmnPidQa.h:138
BmnPidQa::ReadDataBranches
void ReadDataBranches()
Read data branches from input data files.
Definition BmnPidQa.cxx:134
BmnPidQa::fTof400Hits
TClonesArray * fTof400Hits
Definition BmnPidQa.h:164
BmnPidQa::fBetaRangeMin
Double_t fBetaRangeMin
Definition BmnPidQa.h:122
BmnPidQa::fOutName
TString fOutName
Definition BmnPidQa.h:110
BmnPidQa::fParticles
vector< TParticlePDG * > fParticles
Definition BmnPidQa.h:144
BmnPidQa::CreateTrackHitsHistogram
void CreateTrackHitsHistogram(const string &detName)
Definition BmnPidQa.cxx:180
BmnPidQa::fOutputDir
string fOutputDir
Definition BmnPidQa.h:114
BmnPidQa::fTimeRangeMin
Double_t fTimeRangeMin
Definition BmnPidQa.h:130
BmnPidQa::fPRangeBins
Int_t fPRangeBins
Definition BmnPidQa.h:120
BmnPidQa::CreateH2
void CreateH2(const string &name, const string &xTitle, const string &yTitle, const string &zTitle, Int_t nofBinsX, Double_t minBinX, Double_t maxBinX, Int_t nofBinsY, Double_t minBinY, Double_t maxBinY)
Definition BmnPidQa.cxx:164
BmnPidQa::fPrimes
Bool_t fPrimes
Definition BmnPidQa.h:116
BmnPidQa::fMassRangeMax
Double_t fMassRangeMax
Definition BmnPidQa.h:135
BmnPidQa::EnumToPdg
Int_t EnumToPdg(PidParticles part)
Definition BmnPidQa.cxx:535
BmnPidQa::fNHitsRangeMax
Double_t fNHitsRangeMax
Definition BmnPidQa.h:139
BmnPidQa::CreateH1
void CreateH1(const string &name, const string &xTitle, const string &yTitle, Int_t nofBins, Double_t minBin, Double_t maxBin)
Definition BmnPidQa.cxx:153
BmnPidQa::Finish
virtual void Finish()
Derived from FairTask.
Definition BmnPidQa.cxx:115
BmnPidQa::fNHitsRangeBins
Int_t fNHitsRangeBins
Definition BmnPidQa.h:140
BmnPidQa::ProcessGlobal
void ProcessGlobal()
Definition BmnPidQa.cxx:347
BmnPidQa::Init
virtual InitStatus Init()
Derived from FairTask.
Definition BmnPidQa.cxx:102
BmnPidQa::fStorageName
TString fStorageName
Definition BmnPidQa.h:111
BmnPidQa::fBetaRangeBins
Int_t fBetaRangeBins
Definition BmnPidQa.h:124
BmnPidQa::fPidStatistics700
unordered_map< string, vector< Int_t > > fPidStatistics700
Definition BmnPidQa.h:156
BmnPidQa::MassTablePrint
void MassTablePrint()
Definition BmnPidQa.cxx:580
BmnPidQa::fVelocRangeMin
Double_t fVelocRangeMin
Definition BmnPidQa.h:126
BmnPidQa::PidStatisticsPrint700
void PidStatisticsPrint700()
Definition BmnPidQa.cxx:591
BmnPidQa::fGlobalTrackMatches
TClonesArray * fGlobalTrackMatches
Definition BmnPidQa.h:163
BmnPidQa::fMCTracks
TClonesArray * fMCTracks
Definition BmnPidQa.h:161
BmnPidQa::PidStatisticsPrint400
void PidStatisticsPrint400()
Definition BmnPidQa.cxx:586
BmnPidQa::fPRangeMin
Double_t fPRangeMin
Definition BmnPidQa.h:118
BmnPidQa::GetOnlyPrimes
Bool_t GetOnlyPrimes() const
Definition BmnPidQa.h:75
BmnPidQa::fVelocRangeBins
Int_t fVelocRangeBins
Definition BmnPidQa.h:128
BmnPidQa::GetParticleExtend
TParticlePDG * GetParticleExtend(Int_t pdgCode)
Definition BmnPidQa.cxx:560
BmnPidQa::CreateHistograms
void CreateHistograms()
Definition BmnPidQa.cxx:193
BmnPidQa::fHM
BmnHistManager * fHM
Definition BmnPidQa.h:113
BmnPidQa::fMassRangeBins
Int_t fMassRangeBins
Definition BmnPidQa.h:136
BmnPidQa::db
TDatabasePDG * db
Definition BmnPidQa.h:142
BmnPidQa::fGlobalTracks
TClonesArray * fGlobalTracks
Definition BmnPidQa.h:162
BmnPidQa::Exec
virtual void Exec(Option_t *opt)
Derived from FairTask.
Definition BmnPidQa.cxx:110
BmnPidQa::fBetaRangeMax
Double_t fBetaRangeMax
Definition BmnPidQa.h:123
BmnPidQa::fPRangeMax
Double_t fPRangeMax
Definition BmnPidQa.h:119
BmnPidQa::fVelocRangeMax
Double_t fVelocRangeMax
Definition BmnPidQa.h:127
BmnPidQa::fTimeRangeMax
Double_t fTimeRangeMax
Definition BmnPidQa.h:131
BmnPidQa::fMassTable
unordered_map< Double_t, string > fMassTable
Definition BmnPidQa.h:146
BmnPidQa::fMassRangeMin
Double_t fMassRangeMin
Definition BmnPidQa.h:134
BmnPidQa::fPidStatistics400
unordered_map< string, vector< Int_t > > fPidStatistics400
Definition BmnPidQa.h:156
BmnPidQa::fTimeRangeBins
Int_t fTimeRangeBins
Definition BmnPidQa.h:132
BmnPidQa::fTof700Hits
TClonesArray * fTof700Hits
Definition BmnPidQa.h:165
BmnSimulationReport
Base class for simulation reports.
Definition BmnSimulationReport.h:34
BmnSimulationReport::Create
void Create(BmnHistManager *histManager, const string &outputDir)
Main function which creates report data.
Definition BmnSimulationReport.cxx:35
BmnSimulationReport::SetOnlyPrimes
void SetOnlyPrimes(const Bool_t prime)
Definition BmnSimulationReport.h:104
BmnTrack::GetNHits
Int_t GetNHits() const
Accessors.
Definition BmnTrack.h:49
BmnTrack::GetP
Double_t GetP()
Definition BmnTrack.h:67
CbmMCTrack
CbmMCTrack.h.
Definition CbmMCTrack.h:26
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition CbmMCTrack.h:57
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Accessors.
Definition CbmMCTrack.h:56
TDatabasePDG::Instance
static TDatabasePDG * Instance()
TNamed::GetName
const char * GetName() const override
TParticlePDG::Mass
Double_t Mass() const
beta
double beta(double x, double y)
abs
RVec< PromoteType< T > > abs(const RVec< T > &v)
lit::FindAndReplace
string FindAndReplace(const string &name, const string &oldSubstr, const string &newSubstr)
Definition BmnUtils.cxx:19
lit::Split
vector< string > Split(const string &name, char delimiter)
Definition BmnUtils.cxx:27
std
Definition json.hpp:4446