BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
CbmKFParticlesFinderQA.cxx
Go to the documentation of this file.
1/*
2 *====================================================================
3 *
4 * CBM KF Track Quality
5 *
6 * Authors: M.Zyzak
7 *
8 * e-mail :
9 *
10 *====================================================================
11 *
12 * KF Fit performance
13 *
14 *====================================================================
15 */
16#include "CbmKFParticlesFinderQA.h"
17
18#include "CbmKF.h"
19#include "CbmKFParticlesFinder.h"
20/*#include "CbmKFParticle.h"*/
21#include "CbmKFVertex.h"
22#include "CbmL1PFFitter.h"
23#include "CbmMCTrack.h"
24#include "CbmMvdHitMatch.h"
25#include "CbmMvdPoint.h"
26#include "CbmStsCluster.h"
27#include "CbmStsDigi.h"
28#include "CbmStsHit.h"
29#include "CbmStsPoint.h"
30#include "CbmStsTrack.h"
31#include "CbmTrackMatch.h"
32#include "CbmVertex.h"
33#include "FairRunAna.h"
34#include "KFParticle.h"
35#include "KFParticleMatch.h"
36#include "KFParticleSIMD.h"
37#include "L1Field.h"
38#include "TAxis.h"
39#include "TClonesArray.h"
40#include "TDirectory.h"
41#include "TFile.h"
42#include "TH1F.h"
43#include "TH2F.h"
44
45#include <algorithm>
46#include <cmath>
47#include <fstream>
48#include <iostream>
49#include <map>
50using namespace std;
51using std::vector;
52
53CbmKFParticlesFinderQA::CbmKFParticlesFinderQA(CbmKFParticlesFinder* pf,
54 Int_t iVerbose,
55 int findParticlesMode,
56 int perf,
57 const char* name,
58 const char* title)
59 : FairTask(name, iVerbose)
60 , fFindParticlesMode(findParticlesMode)
61 , fPerformance(perf)
62 , fMinNStations(4)
63 , fMinRecoMom(0.1)
64 , fPF(pf)
65 , fMCTrackPoints()
66 , vMCParticles()
67 , MCtoRParticleId()
68 , RtoMCParticleId()
69 , vIsBkgWithSamePDG()
70 , flistStsPts(0)
71 , flistMvdPts(0)
72 , flistMCTracks(0)
73 , flistStsTracksMatch(0)
74 , flistStsTracks(0)
75 , flistStsHits(0)
76 , flistMvdHits(0)
77 , flistMvdHitMatches(0)
78 , flistStsClusters(0)
79 , flistStsDigi(0)
80 , fPrimVtx(0)
81 , fRecParticles(0)
82 , fMCParticles(0)
83 , fMatchParticles(0)
84 , fSaveParticles(0)
85 , fSaveMCParticles(0)
86 , outfileName("CbmParticlesFinderQA.root")
87 , fEfffileName("Efficiency.txt")
88 , histodir(0)
89 , vStsHitMatch()
90 , vStsPointMatch()
91 , vMvdPointMatch()
92 , vMCTrackMatch()
93 , fParteff()
94 , fNEvents(0)
95 , hFitDaughtersQA()
96 , hFitQA()
97 , hPartParam()
98 , hPartParamBG()
99 , hPartParamGhost()
100 ,
101 // hPartParamCorrBG(),
102 hPartParamSignal()
103 , hPartParamQA()
104 , hPartParam2D()
105 , hPartParam2DBG()
106 , hPartParam2DGhost()
107 , hPartParam2DSignal()
108 , hPVFitQa()
109 , fMotherPdgToIndex()
110 , hMotherPdg()
111 , hTrackParameters()
112{
113 TString partNames[nHistoMotherPdg] = {"e+", "e-", "#mu+", "#mu-", "#pi+", "#pi-",
114 "K+", "K-", "p+", "p-", "K0s", "#Lambda",
115 "#Lambda_bar", "#Xi+", "#Xi-", "#Omega+", "#Omega-"};
116 int partPdg[nHistoMotherPdg] = {-11, 11, -13, 13, 211, -211, 321, -321, 2212,
117 -2212, 310, 3122, -3122, 3312, -3312, 3334, -3334};
118 for (int iMP = 0; iMP < nHistoMotherPdg; iMP++)
119 fMotherPdgToIndex[partPdg[iMP]] = iMP;
120
121 TDirectory* currentDir = gDirectory;
122 gDirectory->mkdir("KFParticlesFinder");
123 gDirectory->cd("KFParticlesFinder");
124 histodir = gDirectory;
125 gDirectory->mkdir("Particles");
126 gDirectory->cd("Particles");
127 {
128 for (int iPart = 0; iPart < fParteff.nParticles; ++iPart) {
129 gDirectory->mkdir(fParteff.partName[iPart].Data());
130 gDirectory->cd(fParteff.partName[iPart].Data());
131 {
132 TString res = "res";
133 TString pull = "pull";
134
135 gDirectory->mkdir("DaughtersQA");
136 gDirectory->cd("DaughtersQA");
137 {
138 TString parName[nFitQA / 2] = {"X", "Y", "Z", "Px", "Py", "Pz", "E", "M"};
139 int nBins = 100;
140 float xMax[nFitQA / 2] = {0.15, 0.15, 0.03, 0.01, 0.01, 0.06, 0.06, 0.01};
141
142 for (int iH = 0; iH < nFitQA / 2; iH++) {
143 hFitDaughtersQA[iPart][iH] = new TH1F((res + parName[iH]).Data(), (res + parName[iH]).Data(),
144 nBins, -xMax[iH], xMax[iH]);
145 hFitDaughtersQA[iPart][iH + 8] =
146 new TH1F((pull + parName[iH]).Data(), (pull + parName[iH]).Data(), nBins, -6, 6);
147 }
148 }
149 gDirectory->cd(".."); // particle directory
150
151 gDirectory->mkdir("FitQA");
152 gDirectory->cd("FitQA");
153 {
154 TString parName[nFitQA / 2] = {"X", "Y", "Z", "Px", "Py", "Pz", "E", "M"};
155 int nBins = 50;
156 float xMax[nFitQA / 2] = {0.15, 0.15, 1.2, 0.02, 0.02, 0.15, 0.15, 0.006};
157
158 for (int iH = 0; iH < nFitQA / 2; iH++) {
159 hFitQA[iPart][iH] = new TH1F((res + parName[iH]).Data(), (res + parName[iH]).Data(), nBins,
160 -xMax[iH], xMax[iH]);
161 hFitQA[iPart][iH + 8] =
162 new TH1F((pull + parName[iH]).Data(), (pull + parName[iH]).Data(), nBins, -6, 6);
163 }
164 }
165 gDirectory->cd(".."); // particle directory
166
167 gDirectory->mkdir("Parameters");
168 gDirectory->cd("Parameters");
169 {
170 TString parName[nHistoPartParam] = {"M", "p", "p_{t}", "y", "DL",
171 "c#tau", "chi2ndf", "chi2ndfTopo", "prob", "#theta",
172 "phi", "z", "l/dl"};
173 TString parAxisName[nHistoPartParam] = {"m [GeV/c^{2}]",
174 "p [GeV/c]",
175 "p_{t} [GeV/c]",
176 "y",
177 "Decay length [cm]",
178 "Life time c#tau [cm]",
179 "chi2/ndf",
180 "chi2/ndf topo",
181 "prob",
182 "#theta [rad]",
183 "phi [rad]",
184 "z [cm]",
185 "l/dl"};
186 int nBins[nHistoPartParam] = {1000, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100};
187 float xMin[nHistoPartParam] = {
188 fParteff.partMHistoMin[iPart], 0., 0., 0., -5., 0., 0., 0., 0., -2., -2., -5., -1.};
189 float xMax[nHistoPartParam] = {
190 fParteff.partMHistoMax[iPart], 10., 3., 6., 55., 30., 20., 20., 1., 2., 2., 55., 35.};
191
192 for (int iH = 0; iH < nHistoPartParam; iH++) {
193 hPartParam[iPart][iH] =
194 new TH1F(parName[iH].Data(), parName[iH].Data(), nBins[iH], xMin[iH], xMax[iH]);
195 hPartParam[iPart][iH]->GetXaxis()->SetTitle(parAxisName[iH].Data());
196 }
197
198 hPartParam2D[iPart][0] =
199 new TH2F("y-p_{t}", "y-p_{t}", nBins[3], xMin[3], xMax[3], nBins[2], xMin[2], xMax[2]);
200 hPartParam2D[iPart][0]->GetXaxis()->SetTitle("y");
201 hPartParam2D[iPart][0]->GetYaxis()->SetTitle("p_{t} [GeV/c]");
202
203 gDirectory->mkdir("Signal");
204 gDirectory->cd("Signal");
205 {
206 for (int iH = 0; iH < nHistoPartParam; iH++) {
207 hPartParamSignal[iPart][iH] =
208 new TH1F((parName[iH]).Data(), (parName[iH]).Data(), nBins[iH], xMin[iH], xMax[iH]);
209 hPartParamSignal[iPart][iH]->GetXaxis()->SetTitle(parAxisName[iH].Data());
210 }
211
212 hPartParam2DSignal[iPart][0] =
213 new TH2F("y-p_{t}", "y-p_{t}", nBins[3], xMin[3], xMax[3], nBins[2], xMin[2], xMax[2]);
214 hPartParam2DSignal[iPart][0]->GetXaxis()->SetTitle("y");
215 hPartParam2DSignal[iPart][0]->GetYaxis()->SetTitle("p_{t} [GeV/c]");
216
217 gDirectory->mkdir("QA");
218 gDirectory->cd("QA");
219 {
220 int nBinsQA = 50;
221 float xMaxQA[nHistoPartParamQA] = {0.01, 0.001, 0.001};
222 for (int iH = 0; iH < nHistoPartParamQA; iH++) {
223 hPartParamQA[iPart][iH] =
224 new TH1F((res + parName[iH]).Data(), (res + parName[iH]).Data(), nBinsQA,
225 -xMaxQA[iH], xMaxQA[iH]);
226 hPartParamQA[iPart][iH + nHistoPartParamQA] =
227 new TH1F((pull + parName[iH]).Data(), (pull + parName[iH]).Data(), nBinsQA, -6, 6);
228 }
229 }
230 gDirectory->cd(".."); // particle directory / Parameters / Signal
231 }
232 gDirectory->cd(".."); // particle directory / Parameters
233 gDirectory->mkdir("Background");
234 gDirectory->cd("Background");
235 {
236 for (int iH = 0; iH < nHistoPartParam; iH++) {
237 hPartParamBG[iPart][iH] =
238 new TH1F((parName[iH]).Data(), (parName[iH]).Data(), nBins[iH], xMin[iH], xMax[iH]);
239 hPartParamBG[iPart][iH]->GetXaxis()->SetTitle(parAxisName[iH].Data());
240 }
241
242 hPartParam2DBG[iPart][0] =
243 new TH2F("y-p_{t}", "y-p_{t}", nBins[3], xMin[3], xMax[3], nBins[2], xMin[2], xMax[2]);
244 hPartParam2DBG[iPart][0]->GetXaxis()->SetTitle("y");
245 hPartParam2DBG[iPart][0]->GetYaxis()->SetTitle("p_{t} [GeV/c]");
246 }
247 // gDirectory->cd(".."); // particle directory
248 // gDirectory->mkdir("CorrBG");
249 // gDirectory->cd("CorrBG");
250 // {
251 // for(int iH=0; iH<nHistoPartParam; iH++)
252 // {
253 // hPartParamCorrBG[iPart][iH] = new
254 // TH1F((parName[iH]).Data(),(parName[iH]).Data(),
255 // nBins[iH],xMin[iH],xMax[iH]);
256 // hPartParamCorrBG[iPart][iH]->GetXaxis()->SetTitle(parAxisName[iH].Data());
257 // }
258 //
259 // hPartParam2DCorrBG[iPart][0] = new TH2F("y-p_{t}","y-p_{t}",
260 // nBins[3],xMin[3],xMax[3],
261 // nBins[2],xMin[2],xMax[2]);
262 // hPartParam2DCorrBG[iPart][0]->GetXaxis()->SetTitle("y");
263 // hPartParam2DCorrBG[iPart][0]->GetYaxis()->SetTitle("p_{t} [GeV/c]");
264 // }
265 gDirectory->cd(".."); // particle directory
266 gDirectory->mkdir("Ghost");
267 gDirectory->cd("Ghost");
268 {
269 for (int iH = 0; iH < nHistoPartParam; iH++) {
270 hPartParamGhost[iPart][iH] =
271 new TH1F((parName[iH]).Data(), (parName[iH]).Data(), nBins[iH], xMin[iH], xMax[iH]);
272 hPartParamGhost[iPart][iH]->GetXaxis()->SetTitle(parAxisName[iH].Data());
273 }
274
275 hPartParam2DGhost[iPart][0] =
276 new TH2F("y-p_{t}", "y-p_{t}", nBins[3], xMin[3], xMax[3], nBins[2], xMin[2], xMax[2]);
277 hPartParam2DGhost[iPart][0]->GetXaxis()->SetTitle("y");
278 hPartParam2DGhost[iPart][0]->GetYaxis()->SetTitle("p_{t} [GeV/c]");
279 }
280 gDirectory->cd(".."); // particle directory
281 }
282 gDirectory->cd(".."); // particle directory
283 }
284 gDirectory->cd(".."); // Particles
285 }
286 }
287 gDirectory->cd(".."); // L1
288 gDirectory->mkdir("PrimaryVertexQA");
289 gDirectory->cd("PrimaryVertexQA");
290 {
291 struct
292 {
293 string name;
294 string title;
295 Int_t n;
296 Double_t l, r;
297 } Table[nHistosPV] = {{"PVResX", "x_{rec}-x_{mc}, cm", 100, -0.006f, 0.006f},
298 {"PVResY", "y_{rec}-y_{mc}, cm", 100, -0.006f, 0.006f},
299 {"PVResZ", "z_{rec}-z_{mc}, cm", 100, -0.06f, 0.06f},
300 {"PVPullX", "Pull X", 100, -6.f, 6.f},
301 {"PVPullY", "Pull Y", 100, -6.f, 6.f},
302 {"PVPullZ", "Pull Z", 100, -6.f, 6.f}};
303 for (int iHPV = 0; iHPV < nHistosPV; ++iHPV) {
304 hPVFitQa[iHPV] = new TH1F(Table[iHPV].name.data(), Table[iHPV].title.data(), Table[iHPV].n, Table[iHPV].l,
305 Table[iHPV].r);
306 }
307 }
308 gDirectory->cd(".."); // L1
309 gDirectory->mkdir("MotherPDG");
310 gDirectory->cd("MotherPDG");
311 {
312 for (int iMP = 0; iMP < nHistoMotherPdg; iMP++)
313 hMotherPdg[iMP] = new TH1F(partNames[iMP].Data(), partNames[iMP].Data(), 10000, -5000, 5000);
314 }
315 gDirectory->cd(".."); // particle directory
316
317 gDirectory->mkdir("TrackParameters");
318 gDirectory->cd("TrackParameters");
319 {
320 TString chi2Name = "Chi2Prim";
321 for (int iPart = 0; iPart < CbmKFPartEfficiencies::nParticles; iPart++) {
322 TString chi2NamePart = "Chi2Prim";
323 chi2NamePart += "_";
324 chi2NamePart += fParteff.partName[iPart].Data();
325 hTrackParameters[iPart] = new TH1F(chi2NamePart.Data(), chi2NamePart.Data(), 1000, 0, 100);
326 }
327 hTrackParameters[CbmKFPartEfficiencies::nParticles] = new TH1F("Chi2Prim", "Chi2Prim", 1000, 0, 100);
328 }
329 gDirectory->cd(".."); // particle directory
330
331 gDirectory = currentDir;
332}
333
334CbmKFParticlesFinderQA::~CbmKFParticlesFinderQA()
335{
336 if (fSaveParticles)
337 fRecParticles->Delete();
338 if (fSaveMCParticles) {
339 fMCParticles->Delete();
340 fMatchParticles->Delete();
341 }
342}
343
344InitStatus CbmKFParticlesFinderQA::ReInit()
345{
346 return Init();
347}
348
349InitStatus CbmKFParticlesFinderQA::Init()
350{
351 FairRootManager* fManger = FairRootManager::Instance();
352
353 flistMCTracks = dynamic_cast<TClonesArray*>(fManger->GetObject("MCTrack"));
354 flistStsPts = dynamic_cast<TClonesArray*>(fManger->GetObject("StsPoint"));
355 flistStsTracks = dynamic_cast<TClonesArray*>(fManger->GetObject("StsTrack"));
356 flistStsTracksMatch = dynamic_cast<TClonesArray*>(fManger->GetObject("StsTrackMatch"));
357 flistStsHits = dynamic_cast<TClonesArray*>(fManger->GetObject("StsHit"));
358 flistStsClusters = dynamic_cast<TClonesArray*>(fManger->GetObject("StsCluster"));
359 flistStsDigi = dynamic_cast<TClonesArray*>(fManger->GetObject("StsDigi"));
360 fPrimVtx = (CbmVertex*)fManger->GetObject("PrimaryVertex");
361
362 if (CbmKF::Instance()->GetNMvdStations() == 0) {
363 flistMvdPts = 0;
364 flistMvdHits = 0;
365 flistMvdHitMatches = 0;
366 } else {
367 flistMvdPts = dynamic_cast<TClonesArray*>(fManger->GetObject("MvdPoint"));
368 flistMvdHits = dynamic_cast<TClonesArray*>(fManger->GetObject("MvdHit"));
369 flistMvdHitMatches = dynamic_cast<TClonesArray*>(fManger->GetObject("MvdHitMatch"));
370 }
371
372 if (fSaveParticles) {
373 // create and register TClonesArray with output reco particles
374 fRecParticles = new TClonesArray("KFParticle", 100);
375 fManger->Register("RecoParticles", "KFParticle", fRecParticles, kTRUE);
376 }
377
378 if (fSaveMCParticles) {
379 // create and register TClonesArray with output MC particles
380 fMCParticles = new TClonesArray("KFMCParticle", 100);
381 fManger->Register("KFMCParticles", "KFParticle", fMCParticles, kTRUE);
382
383 // create and register TClonesArray with matching between reco and MC particles
384 fMatchParticles = new TClonesArray("KFParticleMatch", 100);
385 fManger->Register("KFParticleMatch", "KFParticle", fMatchParticles, kTRUE);
386 }
387
388 return kSUCCESS;
389}
390
391void CbmKFParticlesFinderQA::Exec(Option_t* option)
392{
393 if (fSaveParticles)
394 fRecParticles->Delete();
395 if (fSaveMCParticles) {
396 fMCParticles->Delete();
397 fMatchParticles->Delete();
398 }
399
400 vStsHitMatch.clear();
401 vStsPointMatch.clear();
402 vMvdPointMatch.clear();
403 vMCTrackMatch.clear();
404 vStsHitMatch.resize(flistStsHits->GetEntriesFast());
405 vStsPointMatch.resize(flistStsPts->GetEntriesFast());
406 if (flistMvdPts)
407 vMvdPointMatch.resize(flistMvdPts->GetEntriesFast());
408 vMCTrackMatch.resize(flistMCTracks->GetEntriesFast());
409 for (unsigned int iP = 0; iP < vStsPointMatch.size(); iP++)
410 vStsPointMatch[iP] = -1;
411 for (unsigned int iP = 0; iP < vMvdPointMatch.size(); iP++)
412 vMvdPointMatch[iP] = -1;
413 for (unsigned int iTr = 0; iTr < vMCTrackMatch.size(); iTr++)
414 vMCTrackMatch[iTr] = -1;
415
416 if (flistMvdHits) {
417 for (int iH = 0; iH < flistMvdHits->GetEntriesFast(); iH++) {
418 CbmMvdHitMatch* match = static_cast<CbmMvdHitMatch*>(flistMvdHitMatches->At(iH));
419 if (match) {
420 int iMC = match->GetPointId();
421 if (iMC < 0)
422 continue;
423 vMvdPointMatch[iMC] = iH;
424 }
425 }
426 }
427
428 StsHitMatch();
429
430 for (int iTr = 0; iTr < flistStsTracksMatch->GetEntriesFast(); iTr++) {
431 CbmTrackMatch* StsTrackMatch = (CbmTrackMatch*)flistStsTracksMatch->At(iTr);
432 if (StsTrackMatch->GetNofMCTracks() == 0)
433 continue;
434 const int mcTrackId = StsTrackMatch->GetMCTrackId();
435 vMCTrackMatch[mcTrackId] = iTr;
436 }
437
438 fMCTrackPoints.clear();
439 fMCTrackPoints.resize(flistMCTracks->GetEntriesFast());
440 if (CbmKF::Instance()->GetNMvdStations() > 0)
441 for (Int_t iMvd = 0; iMvd < flistMvdPts->GetEntriesFast(); iMvd++) {
442 CbmMvdPoint* MvdPoint = (CbmMvdPoint*)flistMvdPts->At(iMvd);
443 if (!MvdPoint)
444 continue;
445 fMCTrackPoints[MvdPoint->GetTrackID()].MvdArray.push_back(MvdPoint);
446 int iHit = vMvdPointMatch[iMvd];
447 if (iHit < 0)
448 continue;
449 CbmMvdHit* MvdHit = (CbmMvdHit*)flistStsHits->At(iHit);
450 if (!MvdHit)
451 continue;
452 fMCTrackPoints[MvdPoint->GetTrackID()].MvdHitsArray.push_back(MvdHit);
453 }
454 for (Int_t iSts = 0; iSts < flistStsPts->GetEntriesFast(); iSts++) {
455 CbmStsPoint* StsPoint = (CbmStsPoint*)flistStsPts->At(iSts);
456 fMCTrackPoints[StsPoint->GetTrackID()].StsArray.push_back(StsPoint);
457 int iHit = vStsPointMatch[iSts];
458 if (iHit < 0)
459 continue;
460 CbmStsHit* StsHit = (CbmStsHit*)flistStsHits->At(iHit);
461 fMCTrackPoints[StsPoint->GetTrackID()].StsHitsArray.push_back(StsHit);
462 }
463
464 GetMCParticles();
465 FindReconstructableMCParticles();
466 MatchParticles();
467 PartEffPerformance();
468 PartHistoPerformance();
469
470 // save particles to file
471 if (fSaveParticles) {
472 for (unsigned int iP = 0; iP < fPF->GetParticles().size(); iP++) {
473 new ((*fRecParticles)[iP]) KFParticle(fPF->GetParticles()[iP]);
474 }
475 }
476
477 if (fSaveMCParticles) {
478 for (unsigned int iP = 0; iP < fPF->GetParticles().size(); iP++) {
479 new ((*fMatchParticles)[iP]) KFParticleMatch();
480 KFParticleMatch* p = (KFParticleMatch*)(fMatchParticles->At(iP));
481
482 Short_t matchType = 0;
483 int iMCPart = -1;
484 if (!RtoMCParticleId[iP].IsMatchedWithPdg()) // background
485 {
486 if (RtoMCParticleId[iP].IsMatched()) {
487 iMCPart = RtoMCParticleId[iP].GetBestMatchWithPdg();
488 matchType = 1;
489 }
490 } else {
491 iMCPart = RtoMCParticleId[iP].GetBestMatchWithPdg();
492 matchType = 2;
493 }
494
495 p->SetMatch(iMCPart);
496 p->SetMatchType(matchType);
497 }
498
499 for (unsigned int iP = 0; iP < vMCParticles.size(); iP++) {
500 new ((*fMCParticles)[iP]) KFMCParticle(vMCParticles[iP]);
501 // KFMCParticle *p = (KFMCParticle*)( fMCParticles->At(iP) );
502 // *p = vMCParticles[iP];
503 }
504 }
505}
506
507void CbmKFParticlesFinderQA::Finish()
508{
509 if (!(outfileName == "")) {
510 TDirectory* curr = gDirectory;
511 TFile* currentFile = gFile;
512 // Open output file and write histograms
513 TFile* outfile;
514
515 outfile = new TFile(outfileName.Data(), "RECREATE");
516 outfile->cd();
517 WriteHistos(histodir);
518 outfile->Close();
519 outfile->Delete();
520 gFile = currentFile;
521 gDirectory = curr;
522 } else {
523 WriteHistosCurFile(histodir);
524 }
525 std::fstream eff(fEfffileName.Data(), fstream::out);
526 eff << fParteff;
527 eff.close();
528}
529
530void CbmKFParticlesFinderQA::WriteHistos(TObject* obj)
531{
532 if (!obj->IsFolder())
533 obj->Write();
534 else {
535 TDirectory* cur = gDirectory;
536 TDirectory* sub = cur->mkdir(obj->GetName());
537 sub->cd();
538 TList* listSub = (static_cast<TDirectory*>(obj))->GetList();
539 TIter it(listSub);
540 while (TObject* obj1 = it())
541 WriteHistos(obj1);
542 cur->cd();
543 }
544}
545
546void CbmKFParticlesFinderQA::WriteHistosCurFile(TObject* obj)
547{
548 if (!obj->IsFolder())
549 obj->Write();
550 else {
551 TDirectory* cur = gDirectory;
552 TDirectory* sub = cur->GetDirectory(obj->GetName());
553 sub->cd();
554 TList* listSub = (static_cast<TDirectory*>(obj))->GetList();
555 TIter it(listSub);
556 while (TObject* obj1 = it())
557 WriteHistosCurFile(obj1);
558 cur->cd();
559 }
560}
561
562void CbmKFParticlesFinderQA::StsHitMatch()
563{
564 const bool useLinks = 1; // 0 - use HitMatch, one_to_one; 1 - use FairLinks, many_to_many. Set 0 to switch to old
565 // definition of efficiency.
566 // TODO: fix trunk problem with links. Set useLinks = 1
567
568 for (int iH = 0; iH < flistStsHits->GetEntriesFast(); iH++) {
569 CbmStsHit* stsHit = dynamic_cast<CbmStsHit*>(flistStsHits->At(iH));
570 vStsHitMatch[iH] = -1;
571 if (useLinks) {
572 if (flistStsClusters) {
573 // if ( NLinks != 2 ) cout << "HitMatch: Error. Hit wasn't matched with 2 clusters." << endl;
574 // 1-st cluster
575 vector<int> stsPointIds; // save here all mc-points matched with first cluster
576 vector<int> stsPointIds_hit;
577 int iL = 0;
578 FairLink link = stsHit->GetLink(iL);
579 CbmStsCluster* stsCluster = dynamic_cast<CbmStsCluster*>(flistStsClusters->At(link.GetIndex()));
580 int NLinks2 = stsCluster->GetNLinks();
581 for (int iL2 = 0; iL2 < NLinks2; iL2++) {
582 FairLink link2 = stsCluster->GetLink(iL2);
583 CbmStsDigi* stsDigi = dynamic_cast<CbmStsDigi*>(flistStsDigi->At(link2.GetIndex()));
584 const int NLinks3 = stsDigi->GetNLinks();
585 for (int iL3 = 0; iL3 < NLinks3; iL3++) {
586 FairLink link3 = stsDigi->GetLink(iL3);
587 int stsPointId = link3.GetIndex();
588 stsPointIds.push_back(stsPointId);
589 } // for mcPoint
590 } // for digi
591 // 2-nd cluster
592 iL = 1;
593 link = stsHit->GetLink(iL);
594 stsCluster = dynamic_cast<CbmStsCluster*>(flistStsClusters->At(link.GetIndex()));
595 NLinks2 = stsCluster->GetNLinks();
596 for (int iL2 = 0; iL2 < NLinks2; iL2++) {
597 FairLink link2 = stsCluster->GetLink(iL2);
598 CbmStsDigi* stsDigi = dynamic_cast<CbmStsDigi*>(flistStsDigi->At(link2.GetIndex()));
599 const int NLinks3 = stsDigi->GetNLinks();
600 for (int iL3 = 0; iL3 < NLinks3; iL3++) {
601 FairLink link3 = stsDigi->GetLink(iL3);
602 int stsPointId = link3.GetIndex();
603
604 if (!find(&(stsPointIds[0]), &(stsPointIds[stsPointIds.size()]), stsPointId))
605 continue; // check if first cluster matched with same mc-point
606 // CbmStsPoint *pt = dynamic_cast<CbmStsPoint*>( listStsPts->At(stsPointId) );
607 // if(!pt) continue;
608 // CbmMCTrack *MCTrack = dynamic_cast<CbmMCTrack*>( listMCTracks->At( pt->GetTrackID()
609 // ) );
610 // if ( !MCTrack ) continue;
611 // if ( abs(MCTrack->GetPdgCode()) >= 10000 ) continue;
612 bool save = 1;
613 for (unsigned int iP = 0; iP < stsPointIds_hit.size(); iP++)
614 if (vStsHitMatch[iH] == stsPointIds_hit[iP])
615 save = 0;
616 stsPointIds_hit.push_back(stsPointId);
617 vStsHitMatch[iH] = stsPointId;
618 vStsPointMatch[stsPointId] = iH;
619 } // for mcPoint
620 } // for digi
621 } // if clusters
622 } // if useLinks
623 } // for hits
624}
625
626void CbmKFParticlesFinderQA::GetMCParticles()
627{
628 vMCParticles.clear();
629 // convert
630 for (int i = 0; i < flistMCTracks->GetEntriesFast(); i++) {
631 CbmMCTrack& mtra = *(static_cast<CbmMCTrack*>(flistMCTracks->At(i)));
632 KFMCParticle part;
633 part.SetMCTrackID(i);
634 part.SetMotherId(mtra.GetMotherId());
635 part.SetPDG(mtra.GetPdgCode());
636 vMCParticles.push_back(part);
637 }
638 // find relations
639 const unsigned int nMCParticles = vMCParticles.size();
640 for (unsigned int iP = 0; iP < nMCParticles; iP++) {
641 KFMCParticle& part = vMCParticles[iP];
642 for (unsigned int iP2 = 0; iP2 < nMCParticles; iP2++) {
643 KFMCParticle& part2 = vMCParticles[iP2];
644
645 if (part.GetMotherId() == part2.GetMCTrackID()) {
646 part2.AddDaughter(iP);
647 }
648 }
649 }
650}
651
652void CbmKFParticlesFinderQA::FindReconstructableMCParticles()
653{
654 const unsigned int nMCParticles = vMCParticles.size();
655
656 for (unsigned int iP = 0; iP < nMCParticles; iP++) {
657 KFMCParticle& part = vMCParticles[iP];
658 CheckMCParticleIsReconstructable(part);
659 }
660}
661
662void CbmKFParticlesFinderQA::CheckMCParticleIsReconstructable(KFMCParticle& part)
663{
664 if (part.IsReconstructable())
665 return;
666
667 // tracks
668 if (/*part.NDaughters() == 0*/ part.GetPDG() == -211 || part.GetPDG() == 211 || part.GetPDG() == 2212
669 || part.GetPDG() == -2212 || part.GetPDG() == 321 || part.GetPDG() == -321 || part.GetPDG() == 11
670 || part.GetPDG() == -11 || part.GetPDG() == 13 || part.GetPDG() == -13)
671 { // TODO other particles
672
673 switch (fFindParticlesMode) {
674 case 1:
675 part.SetAsReconstructable();
676 break;
677 case 2: {
678 int iMCPart = part.GetMCTrackID();
679 CbmMCTrack* mcTr = (static_cast<CbmMCTrack*>(flistMCTracks->At(iMCPart)));
680 if (fMCTrackPoints[iMCPart].IsReconstructable(mcTr, fMinNStations, fPerformance, fMinRecoMom))
681 part.SetAsReconstructable();
682 } break;
683 case 3: {
684 int iMCPart = part.GetMCTrackID();
685 if (vMCTrackMatch[iMCPart] >= 0)
686 part.SetAsReconstructable();
687 } break;
688 default:
689 part.SetAsReconstructable();
690 };
691 }
692 // mother particles
693 else
694 {
695 for (int iPart = 0; iPart < fParteff.nParticles; iPart++) {
696 if (part.GetPDG() == fParteff.partPDG[iPart]) {
697 const unsigned int nDaughters = fParteff.partDaughterPdg[iPart].size();
698 if (part.GetDaughterIds().size() != nDaughters)
699 return;
700 int pdg[nDaughters];
701
702 for (unsigned int iD = 0; iD < nDaughters; iD++)
703 pdg[iD] = vMCParticles[part.GetDaughterIds()[iD]].GetPDG();
704
705 bool isDaughterFound[nDaughters];
706 for (unsigned int iDMC = 0; iDMC < nDaughters; iDMC++)
707 isDaughterFound[iDMC] = 0;
708
709 bool isReco = 1;
710 for (unsigned int iDMC = 0; iDMC < nDaughters; iDMC++)
711 for (unsigned int iD = 0; iD < nDaughters; iD++)
712 if (pdg[iD] == fParteff.partDaughterPdg[iPart][iDMC])
713 isDaughterFound[iDMC] = 1;
714
715 for (unsigned int iDMC = 0; iDMC < nDaughters; iDMC++)
716 isReco = isReco && isDaughterFound[iDMC];
717
718 if (!isReco)
719 return;
720 }
721 }
722
723 const vector<int>& dIds = part.GetDaughterIds();
724 const unsigned int nD = dIds.size();
725 bool reco = 1;
726 for (unsigned int iD = 0; iD < nD && reco; iD++) {
727 KFMCParticle& dp = vMCParticles[dIds[iD]];
728 CheckMCParticleIsReconstructable(dp);
729 reco &= dp.IsReconstructable();
730 }
731 if (reco)
732 part.SetAsReconstructable();
733 }
734}
735
736/// Procedure for match Reconstructed and MC Particles. Should be called before Performances
737void CbmKFParticlesFinderQA::MatchParticles()
738{
739 // get all reco particles ( temp )
740 MCtoRParticleId.clear();
741 RtoMCParticleId.clear();
742 MCtoRParticleId.resize(vMCParticles.size());
743 RtoMCParticleId.resize(fPF->GetParticles().size());
744 vIsBkgWithSamePDG.clear();
745 vIsBkgWithSamePDG.resize(fPF->GetParticles().size());
746
747 // match tracks ( particles which are direct copy of tracks )
748 for (unsigned int iRP = 0; iRP < fPF->GetParticles().size(); iRP++) {
749 // CbmKFParticle &rPart = fPF->GetParticles()[iRP];
750 KFParticle& rPart = fPF->GetParticles()[iRP];
751
752 if (rPart.NDaughters() != 1)
753 continue;
754
755 const int rTrackId = rPart.DaughterIds()[0];
756 CbmTrackMatch* StsTrackMatch = (CbmTrackMatch*)flistStsTracksMatch->At(rTrackId);
757 if (StsTrackMatch->GetNofMCTracks() == 0)
758 continue;
759 const int mcTrackId = StsTrackMatch->GetMCTrackId();
760
761 for (unsigned int iMP = 0; iMP < vMCParticles.size(); iMP++) {
762 KFMCParticle& mPart = vMCParticles[iMP];
763 if (mPart.GetMCTrackID() == mcTrackId) { // match is found
764 if (mPart.GetPDG() == rPart.GetPDG()) {
765 MCtoRParticleId[iMP].ids.push_back(iRP);
766 RtoMCParticleId[iRP].ids.push_back(iMP);
767 } else {
768 MCtoRParticleId[iMP].idsMI.push_back(iRP);
769 RtoMCParticleId[iRP].idsMI.push_back(iMP);
770 }
771 }
772 }
773 }
774
775 // match created mother particles
776 for (unsigned int iRP = 0; iRP < fPF->GetParticles().size(); iRP++) {
777 // CbmKFParticle &rPart = fPF->GetParticles()[iRP];
778 KFParticle& rPart = fPF->GetParticles()[iRP];
779 const unsigned int NRDaughters = rPart.NDaughters();
780 if (NRDaughters < 2)
781 continue;
782
783 unsigned int iD = 0;
784 int mmId = -2; // mother MC id
785 {
786 const int rdId = rPart.DaughterIds()[iD];
787 if (!RtoMCParticleId[rdId].IsMatched())
788 continue;
789 const int mdId = RtoMCParticleId[rdId].GetBestMatch();
790 mmId = vMCParticles[mdId].GetMotherId();
791 }
792 iD++;
793 for (; iD < NRDaughters; iD++) {
794 const int rdId = rPart.DaughterIds()[iD];
795 if (!RtoMCParticleId[rdId].IsMatched())
796 break;
797 const int mdId = RtoMCParticleId[rdId].GetBestMatch();
798 if (vMCParticles[mdId].GetMotherId() != mmId)
799 break;
800 }
801 if (iD == NRDaughters && mmId != -1) { // match is found and it is not primary vertex
802 KFMCParticle& mmPart = vMCParticles[mmId];
803 if (mmPart.GetPDG() == rPart.GetPDG() && mmPart.NDaughters() == rPart.NDaughters()) {
804 MCtoRParticleId[mmId].ids.push_back(iRP);
805 RtoMCParticleId[iRP].ids.push_back(mmId);
806 } else {
807 MCtoRParticleId[mmId].idsMI.push_back(iRP);
808 RtoMCParticleId[iRP].idsMI.push_back(mmId);
809 }
810 }
811 }
812}
813
814void CbmKFParticlesFinderQA::PartEffPerformance()
815{
816 CbmKFPartEfficiencies partEff; // efficiencies for current event
817
818 const int NRP = fPF->GetParticles().size();
819 for (int iP = 0; iP < NRP; ++iP) {
820 // const CbmKFParticle &part = fPF->GetParticles()[iP];
821 const KFParticle& part = fPF->GetParticles()[iP];
822 const int pdg = part.GetPDG();
823
824 const bool isBG = RtoMCParticleId[iP].idsMI.size() != 0;
825 const bool isGhost = !RtoMCParticleId[iP].IsMatched();
826
827 for (int iPart = 0; iPart < fParteff.nParticles; iPart++)
828 if (pdg == fParteff.partPDG[iPart])
829 partEff.IncReco(isGhost, isBG, fParteff.partName[iPart].Data());
830 }
831
832 const int NMP = vMCParticles.size();
833 for (int iP = 0; iP < NMP; ++iP) {
834 const KFMCParticle& part = vMCParticles[iP];
835 if (!part.IsReconstructable())
836 continue;
837 const int pdg = part.GetPDG();
838 const int mId = part.GetMotherId();
839
840 const bool isReco = MCtoRParticleId[iP].ids.size() != 0;
841 const int nClones = MCtoRParticleId[iP].ids.size() - 1;
842
843 for (int iPart = 0; iPart < fParteff.nParticles; iPart++) {
844 if (pdg == fParteff.partPDG[iPart]) {
845 partEff.Inc(isReco, nClones, fParteff.partName[iPart].Data());
846 if (mId == -1)
847 partEff.Inc(isReco, nClones, (fParteff.partName[iPart] + "_prim").Data());
848 else
849 partEff.Inc(isReco, nClones, (fParteff.partName[iPart] + "_sec").Data());
850 }
851 }
852 }
853
854 fNEvents++;
855
856 fParteff += partEff;
857
858 partEff.CalcEff();
859 fParteff.CalcEff();
860
861 // cout.precision(3);
862 if (fVerbose) {
863 if (fNEvents % 100 == 0) {
864 cout << " ---- KF Particle finder --- " << endl;
865 // cout << "L1 STAT : " << fNEvents << " EVENT " << endl << endl;
866 // partEff.PrintEff();
867 // cout << endl;
868 cout << "ACCUMULATED STAT : " << fNEvents << " EVENTS " << endl << endl;
869 fParteff.PrintEff();
870
871 cout << endl;
872 // cout<<"CA Track Finder: " << L1_CATIME/L1_fNEvents << " s/ev" << endl << endl;
873 }
874 }
875}
876
877void CbmKFParticlesFinderQA::PartHistoPerformance()
878{
879 CbmKFVertex vtx;
880 if (fPrimVtx)
881 vtx = CbmKFVertex(*fPrimVtx);
882
883 for (unsigned int iP = 0; iP < fPF->GetParticles().size(); iP++) {
884 int iParticle = fParteff.GetParticleIndex(fPF->GetParticles()[iP].GetPDG());
885 if (iParticle < 0)
886 continue;
887
888 Double_t M, ErrM;
889 Double_t dL, ErrdL; // decay length
890 Double_t cT, ErrcT; // c*tau
891 Double_t P, ErrP;
892 Double_t Pt;
893 Double_t Rapidity;
894 Double_t Theta;
895 Double_t Phi;
896 Double_t Z;
897 // CbmKFParticle TempPart = fPF->GetParticles()[iP];
898 KFParticle TempPart = fPF->GetParticles()[iP];
899 TempPart.GetMass(M, ErrM);
900
901 TempPart.GetMomentum(P, ErrP);
902 Pt = TempPart.GetPt();
903 Rapidity = TempPart.GetRapidity();
904 TempPart.GetDecayLength(dL, ErrdL);
905 TempPart.GetLifeTime(cT, ErrcT);
906 Double_t chi2 = TempPart.GetChi2();
907 Int_t ndf = TempPart.GetNDF();
908 Double_t prob = TMath::Prob(chi2, ndf);
909 Theta = TempPart.GetTheta();
910 Phi = TempPart.GetPhi();
911 Z = TempPart.GetZ();
912
913 KFParticleSIMD tempSIMDPart(TempPart);
914 fvec l, dl;
915 KFParticleSIMD pv(vtx);
916 tempSIMDPart.GetDistanceToVertexLine(pv, l, dl);
917 dL = sqrt(TempPart.X() * TempPart.X() + TempPart.Y() * TempPart.Y() + TempPart.Z() * TempPart.Z());
918
919 tempSIMDPart.SetProductionVertex(pv);
920 fvec chi2topo = tempSIMDPart.Chi2() / tempSIMDPart.GetNDF();
921
922 hPartParam[iParticle][0]->Fill(M);
923 hPartParam[iParticle][1]->Fill(P);
924 hPartParam[iParticle][2]->Fill(Pt);
925 hPartParam[iParticle][3]->Fill(Rapidity);
926 hPartParam[iParticle][4]->Fill(dL);
927 hPartParam[iParticle][5]->Fill(cT);
928 hPartParam[iParticle][6]->Fill(chi2 / ndf);
929 hPartParam[iParticle][7]->Fill(chi2topo[0]);
930 hPartParam[iParticle][8]->Fill(prob);
931 hPartParam[iParticle][9]->Fill(Theta);
932 hPartParam[iParticle][10]->Fill(Phi);
933 hPartParam[iParticle][11]->Fill(Z);
934 hPartParam[iParticle][12]->Fill(l[0] / dl[0]);
935
936 hPartParam2D[iParticle][0]->Fill(Rapidity, Pt, 1);
937
938 if (!RtoMCParticleId[iP].IsMatchedWithPdg()) // background
939 {
940 if (!RtoMCParticleId[iP].IsMatched()) {
941 hPartParamGhost[iParticle][0]->Fill(M);
942 hPartParamGhost[iParticle][1]->Fill(P);
943 hPartParamGhost[iParticle][2]->Fill(Pt);
944 hPartParamGhost[iParticle][3]->Fill(Rapidity);
945 hPartParamGhost[iParticle][4]->Fill(dL);
946 hPartParamGhost[iParticle][5]->Fill(cT);
947 hPartParamGhost[iParticle][6]->Fill(chi2 / ndf);
948 hPartParamGhost[iParticle][7]->Fill(chi2topo[0]);
949 hPartParamGhost[iParticle][8]->Fill(prob);
950 hPartParamGhost[iParticle][9]->Fill(Theta);
951 hPartParamGhost[iParticle][10]->Fill(Phi);
952 hPartParamGhost[iParticle][11]->Fill(Z);
953 hPartParamGhost[iParticle][12]->Fill(l[0] / dl[0]);
954
955 hPartParam2DGhost[iParticle][0]->Fill(Rapidity, Pt, 1);
956
957 // if(vIsBkgWithSamePDG[iP])
958 // {
959 // hPartParamCorrBG[iParticle][ 0]->Fill(M);
960 // hPartParamCorrBG[iParticle][ 1]->Fill(P);
961 // hPartParamCorrBG[iParticle][ 2]->Fill(Pt);
962 // hPartParamCorrBG[iParticle][ 3]->Fill(Rapidity);
963 // hPartParamCorrBG[iParticle][ 4]->Fill(dL);
964 // hPartParamCorrBG[iParticle][ 5]->Fill(cT);
965 // hPartParamCorrBG[iParticle][ 6]->Fill(chi2/ndf);
966 // hPartParamCorrBG[iParticle][ 7]->Fill(prob);
967 // hPartParamCorrBG[iParticle][ 8]->Fill(Theta);
968 // hPartParamCorrBG[iParticle][ 9]->Fill(Phi);
969 // hPartParamCorrBG[iParticle][10]->Fill(Z);
970 //
971 // hPartParam2DCorrBG[iParticle][0]->Fill(Rapidity,Pt,1);
972 // }
973 } else {
974 hPartParamBG[iParticle][0]->Fill(M);
975 hPartParamBG[iParticle][1]->Fill(P);
976 hPartParamBG[iParticle][2]->Fill(Pt);
977 hPartParamBG[iParticle][3]->Fill(Rapidity);
978 hPartParamBG[iParticle][4]->Fill(dL);
979 hPartParamBG[iParticle][5]->Fill(cT);
980 hPartParamBG[iParticle][6]->Fill(chi2 / ndf);
981 hPartParamBG[iParticle][7]->Fill(chi2topo[0]);
982 hPartParamBG[iParticle][8]->Fill(prob);
983 hPartParamBG[iParticle][9]->Fill(Theta);
984 hPartParamBG[iParticle][10]->Fill(Phi);
985 hPartParamBG[iParticle][11]->Fill(Z);
986 hPartParamBG[iParticle][12]->Fill(l[0] / dl[0]);
987
988 hPartParam2DBG[iParticle][0]->Fill(Rapidity, Pt, 1);
989 }
990 continue;
991 }
992 hPartParamSignal[iParticle][0]->Fill(M);
993 hPartParamSignal[iParticle][1]->Fill(P);
994 hPartParamSignal[iParticle][2]->Fill(Pt);
995 hPartParamSignal[iParticle][3]->Fill(Rapidity);
996 hPartParamSignal[iParticle][4]->Fill(dL);
997 hPartParamSignal[iParticle][5]->Fill(cT);
998 hPartParamSignal[iParticle][6]->Fill(chi2 / ndf);
999 hPartParamSignal[iParticle][7]->Fill(chi2topo[0]);
1000 hPartParamSignal[iParticle][8]->Fill(prob);
1001 hPartParamSignal[iParticle][9]->Fill(Theta);
1002 hPartParamSignal[iParticle][10]->Fill(Phi);
1003 hPartParamSignal[iParticle][11]->Fill(Z);
1004 hPartParamSignal[iParticle][12]->Fill(l[0] / dl[0]);
1005
1006 hPartParam2DSignal[iParticle][0]->Fill(Rapidity, Pt, 1);
1007
1008 int iMCPart = RtoMCParticleId[iP].GetBestMatchWithPdg();
1009 KFMCParticle& mcPart = vMCParticles[iMCPart];
1010 // Fit quality of the mother particle
1011 {
1012 int iMCTrack = mcPart.GetMCTrackID();
1013 CbmMCTrack& mcTrack = *(static_cast<CbmMCTrack*>(flistMCTracks->At(iMCTrack)));
1014 int mcDaughterId = mcPart.GetDaughterIds()[0];
1015 CbmMCTrack& mcDaughter = *(static_cast<CbmMCTrack*>(flistMCTracks->At(mcDaughterId)));
1016
1017 Double_t decayVtx[3] = {mcDaughter.GetStartX(), mcDaughter.GetStartY(), mcDaughter.GetStartZ()};
1018 Double_t recParam[8] = {0};
1019 Double_t errParam[8] = {0};
1020
1021 for (int iPar = 0; iPar < 3; iPar++) {
1022 recParam[iPar] = TempPart.GetParameter(iPar);
1023 Double_t error = TempPart.GetCovariance(iPar, iPar);
1024 if (error < 0.) {
1025 error = 1.e20;
1026 }
1027 errParam[iPar] = TMath::Sqrt(error);
1028 }
1029 TempPart.TransportToPoint(decayVtx);
1030 for (int iPar = 3; iPar < 7; iPar++) {
1031 recParam[iPar] = TempPart.GetParameter(iPar);
1032 Double_t error = TempPart.GetCovariance(iPar, iPar);
1033 if (error < 0.) {
1034 error = 1.e20;
1035 }
1036 errParam[iPar] = TMath::Sqrt(error);
1037 }
1038
1039 Double_t Emc = sqrt(mcTrack.GetP() * mcTrack.GetP() + mcTrack.GetMass() * mcTrack.GetMass());
1040 Double_t res[8] = {0}, pull[8] = {0},
1041 mcParam[8] = {decayVtx[0], decayVtx[1], decayVtx[2], mcTrack.GetPx(),
1042 mcTrack.GetPy(), mcTrack.GetPz(), Emc, mcTrack.GetMass()};
1043 for (int iPar = 0; iPar < 7; iPar++) {
1044 res[iPar] = recParam[iPar] - mcParam[iPar];
1045 if (fabs(errParam[iPar]) > 1.e-20)
1046 pull[iPar] = res[iPar] / errParam[iPar];
1047 }
1048 res[7] = M - mcParam[7];
1049 if (fabs(ErrM) > 1.e-20)
1050 pull[7] = res[7] / ErrM;
1051
1052 for (int iPar = 0; iPar < 8; iPar++) {
1053 hFitQA[iParticle][iPar]->Fill(res[iPar]);
1054 hFitQA[iParticle][iPar + 8]->Fill(pull[iPar]);
1055 }
1056
1057 // Double_t mcT = mcDaughter.GetStartT() - mcTrack.GetStartT();
1058 }
1059 // Fit quality of daughters
1060 for (int iD = 0; iD < mcPart.NDaughters(); ++iD) {
1061 int mcDaughterId = mcPart.GetDaughterIds()[iD];
1062 // if(!MCtoRParticleId[mcDaughterId].IsMatchedWithPdg()) continue;
1063 if (!MCtoRParticleId[mcDaughterId].IsMatched())
1064 continue;
1065 CbmMCTrack& mcTrack = *(static_cast<CbmMCTrack*>(flistMCTracks->At(mcDaughterId)));
1066 // int recDaughterId = MCtoRParticleId[mcDaughterId].GetBestMatchWithPdg();
1067 int recDaughterId = MCtoRParticleId[mcDaughterId].GetBestMatch();
1068 // CbmKFParticle Daughter = fPF->GetParticles()[recDaughterId];
1069 KFParticle Daughter = fPF->GetParticles()[recDaughterId];
1070 Daughter.GetMass(M, ErrM);
1071
1072 Double_t decayVtx[3] = {mcTrack.GetStartX(), mcTrack.GetStartY(), mcTrack.GetStartZ()};
1073 Daughter.TransportToPoint(decayVtx);
1074
1075 Double_t Emc = sqrt(mcTrack.GetP() * mcTrack.GetP() + mcTrack.GetMass() * mcTrack.GetMass());
1076 Double_t res[8] = {0}, pull[8] = {0}, mcParam[8] = {mcTrack.GetStartX(),
1077 mcTrack.GetStartY(),
1078 mcTrack.GetStartZ(),
1079 mcTrack.GetPx(),
1080 mcTrack.GetPy(),
1081 mcTrack.GetPz(),
1082 Emc,
1083 mcTrack.GetMass()};
1084 for (int iPar = 0; iPar < 7; iPar++) {
1085 Double_t error = Daughter.GetCovariance(iPar, iPar);
1086 if (error < 0.) {
1087 error = 1.e20;
1088 }
1089 error = TMath::Sqrt(error);
1090 res[iPar] = Daughter.GetParameter(iPar) - mcParam[iPar];
1091 if (fabs(error) > 1.e-20)
1092 pull[iPar] = res[iPar] / error;
1093 }
1094 res[7] = M - mcParam[7];
1095 if (fabs(ErrM) > 1.e-20)
1096 pull[7] = res[7] / ErrM;
1097
1098 for (int iPar = 0; iPar < 8; iPar++) {
1099 hFitDaughtersQA[iParticle][iPar]->Fill(res[iPar]);
1100 hFitDaughtersQA[iParticle][iPar + 8]->Fill(pull[iPar]);
1101 }
1102 }
1103 }
1104 // Find MC parameters of the primary vertex
1105 float mcPVx[3] = {0.f};
1106 for (int iMC = 0; iMC < flistMCTracks->GetEntriesFast(); ++iMC) {
1107 CbmMCTrack* mcTrack = (CbmMCTrack*)flistMCTracks->At(iMC);
1108 if (mcTrack->GetMotherId() < 0) {
1109 mcPVx[0] = mcTrack->GetStartX();
1110 mcPVx[1] = mcTrack->GetStartY();
1111 mcPVx[2] = mcTrack->GetStartZ();
1112 break;
1113 }
1114 }
1115
1116 float dRPVr[3] = {vtx.GetRefX() - mcPVx[0], vtx.GetRefY() - mcPVx[1], vtx.GetRefZ() - mcPVx[2]};
1117 float dRPVp[3] = {dRPVr[0] / sqrt(vtx.GetCovMatrix()[0]), dRPVr[1] / sqrt(vtx.GetCovMatrix()[2]),
1118 dRPVr[2] / sqrt(vtx.GetCovMatrix()[5])};
1119 for (unsigned int iHPV = 0; iHPV < 3; ++iHPV)
1120 hPVFitQa[iHPV]->Fill(dRPVr[iHPV]);
1121 for (unsigned int iHPV = 3; iHPV < 6; ++iHPV)
1122 hPVFitQa[iHPV]->Fill(dRPVp[iHPV - 3]);
1123
1124 // fill mother pdg histo
1125 for (int i = 0; i < flistMCTracks->GetEntriesFast(); i++) {
1126 CbmMCTrack& mtra = *(static_cast<CbmMCTrack*>(flistMCTracks->At(i)));
1127 int motherId = mtra.GetMotherId();
1128 int pdg = mtra.GetPdgCode();
1129 int histoIndex = GetMotherPdgIndex(pdg);
1130 if (histoIndex != -1) {
1131 if (motherId > -1) {
1132 CbmMCTrack& mother = *(static_cast<CbmMCTrack*>(flistMCTracks->At(mtra.GetMotherId())));
1133 int motherPdg = mother.GetPdgCode();
1134 hMotherPdg[histoIndex]->Fill(motherPdg);
1135 } else
1136 hMotherPdg[histoIndex]->Fill(-1);
1137 }
1138 }
1139
1140 // fill histo with chi2_prim
1141 for (unsigned int iP = 0; iP < fPF->GetParticles().size(); iP++) {
1142 const KFParticle& rPart = fPF->GetParticles()[iP];
1143 const unsigned int NRDaughters = rPart.NDaughters();
1144 if (NRDaughters > 1)
1145 break;
1146 if (RtoMCParticleId[iP].GetBestMatch() < 0)
1147 continue;
1148 KFMCParticle& mPart = vMCParticles[RtoMCParticleId[iP].GetBestMatch()];
1149
1150 if (mPart.GetMotherId() < 0) {
1151 hTrackParameters[CbmKFPartEfficiencies::nParticles]->Fill(fPF->GetChiPrim()[iP]);
1152 continue;
1153 }
1154
1155 KFMCParticle& mMotherPart = vMCParticles[mPart.GetMotherId()];
1156
1157 int iParticle = fParteff.GetParticleIndex(mMotherPart.GetPDG());
1158 float chiPrim = fPF->GetChiPrim()[iP];
1159 if (iParticle > -1 && iParticle < CbmKFPartEfficiencies::nParticles)
1160 hTrackParameters[iParticle]->Fill(chiPrim);
1161 }
1162} // void CbmKFParticlesFinderQA::ParticlesEfficienciesPerformance()
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
fvec
F32vec4 fvec
Definition P4_F32vec4.h:231
Int_t
int Int_t
Short_t
short Short_t
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
Option_t
const char Option_t
gDirectory
#define gDirectory
gFile
#define gFile
CbmKFPartEfficiencies::Inc
void Inc(bool isReco, int nClones, TString name)
Definition CbmKFPartEfficiencies.h:444
CbmKFPartEfficiencies::IncReco
void IncReco(bool isGhost, bool isBg, TString name)
Definition CbmKFPartEfficiencies.h:455
CbmKFParticlesFinderQA::hPartParam
TH1F * hPartParam[CbmKFPartEfficiencies::nParticles][nHistoPartParam]
Definition CbmKFParticlesFinderQA.h:127
CbmKFParticlesFinderQA::RtoMCParticleId
vector< KFMatchParticles > RtoMCParticleId
Definition CbmKFParticlesFinderQA.h:81
CbmKFParticlesFinderQA::vMvdPointMatch
std::vector< int > vMvdPointMatch
Definition CbmKFParticlesFinderQA.h:113
CbmKFParticlesFinderQA::hPartParam2DSignal
TH2F * hPartParam2DSignal[CbmKFPartEfficiencies::nParticles][nHistoPartParam2D]
Definition CbmKFParticlesFinderQA.h:141
CbmKFParticlesFinderQA::fMotherPdgToIndex
std::map< int, int > fMotherPdgToIndex
Definition CbmKFParticlesFinderQA.h:147
CbmKFParticlesFinderQA::hPartParam2DGhost
TH2F * hPartParam2DGhost[CbmKFPartEfficiencies::nParticles][nHistoPartParam2D]
Definition CbmKFParticlesFinderQA.h:139
CbmKFParticlesFinderQA::fParteff
CbmKFPartEfficiencies fParteff
Definition CbmKFParticlesFinderQA.h:117
CbmKFParticlesFinderQA::hPartParamBG
TH1F * hPartParamBG[CbmKFPartEfficiencies::nParticles][nHistoPartParam]
Definition CbmKFParticlesFinderQA.h:128
CbmKFParticlesFinderQA::hPartParam2D
TH2F * hPartParam2D[CbmKFPartEfficiencies::nParticles][nHistoPartParam2D]
Definition CbmKFParticlesFinderQA.h:137
CbmKFParticlesFinderQA::fMCTrackPoints
vector< CbmKFTrErrMCPoints > fMCTrackPoints
Definition CbmKFParticlesFinderQA.h:78
CbmKFParticlesFinderQA::vMCParticles
vector< KFMCParticle > vMCParticles
Definition CbmKFParticlesFinderQA.h:79
CbmKFParticlesFinderQA::hFitQA
TH1F * hFitQA[CbmKFPartEfficiencies::nParticles][nFitQA]
Definition CbmKFParticlesFinderQA.h:123
CbmKFParticlesFinderQA::hMotherPdg
TH1F * hMotherPdg[nHistoMotherPdg]
Definition CbmKFParticlesFinderQA.h:148
CbmKFParticlesFinderQA::WriteHistosCurFile
void WriteHistosCurFile(TObject *obj)
Definition CbmKFParticlesFinderQA.cxx:546
CbmKFParticlesFinderQA::vStsPointMatch
std::vector< int > vStsPointMatch
Definition CbmKFParticlesFinderQA.h:112
CbmKFParticlesFinderQA::MatchParticles
void MatchParticles()
Procedure for match Reconstructed and MC Particles. Should be called before Performances.
Definition CbmKFParticlesFinderQA.cxx:737
CbmKFParticlesFinderQA::CheckMCParticleIsReconstructable
void CheckMCParticleIsReconstructable(KFMCParticle &part)
Definition CbmKFParticlesFinderQA.cxx:662
CbmKFParticlesFinderQA::hPartParamGhost
TH1F * hPartParamGhost[CbmKFPartEfficiencies::nParticles][nHistoPartParam]
Definition CbmKFParticlesFinderQA.h:129
CbmKFParticlesFinderQA::fPF
CbmKFParticlesFinder * fPF
Definition CbmKFParticlesFinderQA.h:76
CbmKFParticlesFinderQA::hPartParamQA
TH1F * hPartParamQA[CbmKFPartEfficiencies::nParticles][nHistoPartParamQA *2]
Definition CbmKFParticlesFinderQA.h:134
CbmKFParticlesFinderQA::MCtoRParticleId
vector< KFMatchParticles > MCtoRParticleId
Definition CbmKFParticlesFinderQA.h:80
CbmKFParticlesFinderQA::hFitDaughtersQA
TH1F * hFitDaughtersQA[CbmKFPartEfficiencies::nParticles][nFitQA]
Definition CbmKFParticlesFinderQA.h:122
CbmKFParticlesFinderQA::hPartParamSignal
TH1F * hPartParamSignal[CbmKFPartEfficiencies::nParticles][nHistoPartParam]
Definition CbmKFParticlesFinderQA.h:131
CbmKFParticlesFinderQA::CbmKFParticlesFinderQA
CbmKFParticlesFinderQA(CbmKFParticlesFinder *pf=0, Int_t iVerbose=1, int findParticlesMode=1, int perf=3, const char *name="CbmKFParticlesFinderQA", const char *title="Cbm KF Particles Finder QA")
Definition CbmKFParticlesFinderQA.cxx:53
CbmKFParticlesFinderQA::hPartParam2DBG
TH2F * hPartParam2DBG[CbmKFPartEfficiencies::nParticles][nHistoPartParam2D]
Definition CbmKFParticlesFinderQA.h:138
CbmKFVertex::GetRefZ
Double_t & GetRefZ()
Definition CbmKFVertex.h:39
CbmKFVertex::GetRefX
Double_t & GetRefX()
I/O parameters for vertex fit routines.
Definition CbmKFVertex.h:37
CbmKFVertex::GetCovMatrix
Double_t * GetCovMatrix()
Definition CbmKFVertex.h:40
CbmKFVertex::GetRefY
Double_t & GetRefY()
Definition CbmKFVertex.h:38
CbmKF::Instance
static CbmKF * Instance()
Definition CbmKF.h:34
CbmMCTrack
CbmMCTrack.h.
Definition CbmMCTrack.h:26
CbmMCTrack::GetMass
Double_t GetMass() const
Definition CbmMCTrack.cxx:127
CbmMCTrack::GetPy
Double_t GetPy() const
Definition CbmMCTrack.h:59
CbmMCTrack::GetPz
Double_t GetPz() const
Definition CbmMCTrack.h:60
CbmMCTrack::GetPx
Double_t GetPx() const
Definition CbmMCTrack.h:58
CbmMCTrack::GetStartZ
Double_t GetStartZ() const
Definition CbmMCTrack.h:63
CbmMCTrack::GetStartY
Double_t GetStartY() const
Definition CbmMCTrack.h:62
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition CbmMCTrack.h:57
CbmMCTrack::GetP
Double_t GetP() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetStartX
Double_t GetStartX() const
Definition CbmMCTrack.h:61
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Accessors.
Definition CbmMCTrack.h:56
CbmMvdHitMatch
CbmMvdHitMatch header file.
Definition CbmMvdHitMatch.h:20
CbmMvdHit
CbmMvdHit.h.
Definition CbmMvdHit.h:26
CbmMvdPoint
CbmMvdPoint.h.
Definition CbmMvdPoint.h:27
CbmStsCluster
CbmStsCluster.h.
Definition CbmStsCluster.h:32
CbmStsDigi
CbmStsDigi.h.
Definition CbmStsDigi.h:32
CbmTrackMatch
CbmStsTrackMatch.h.
Definition CbmTrackMatch.h:16
CbmTrackMatch::GetMCTrackId
Int_t GetMCTrackId() const
Index of matched MC track.
Definition CbmTrackMatch.h:35
CbmTrackMatch::GetNofMCTracks
Int_t GetNofMCTracks() const
Number of MCTracks with common hits.
Definition CbmTrackMatch.h:47
CbmVertex
CbmVertex.h.
Definition CbmVertex.h:21
KFMCParticle::IsReconstructable
bool IsReconstructable() const
Definition KFMCParticle.h:42
KFMCParticle::GetPDG
int GetPDG() const
Definition KFMCParticle.h:41
KFMCParticle::SetMCTrackID
void SetMCTrackID(int id)
Definition KFMCParticle.h:36
KFMCParticle::GetDaughterIds
const std::vector< int > & GetDaughterIds() const
Definition KFMCParticle.h:32
KFMCParticle::NDaughters
int NDaughters() const
Definition KFMCParticle.h:31
KFMCParticle::GetMotherId
int GetMotherId() const
Definition KFMCParticle.h:40
KFMCParticle::SetAsReconstructable
void SetAsReconstructable()
Definition KFMCParticle.h:43
KFMCParticle::SetPDG
void SetPDG(int pdg)
Definition KFMCParticle.h:35
KFMCParticle::SetMotherId
void SetMotherId(int id)
Definition KFMCParticle.h:37
KFMCParticle::GetMCTrackID
int GetMCTrackID() const
Definition KFMCParticle.h:39
KFMCParticle::AddDaughter
void AddDaughter(int i)
Definition KFMCParticle.cxx:29
KFParticleBaseSIMD::GetDistanceToVertexLine
void GetDistanceToVertexLine(const KFParticleBaseSIMD &Vertex, fvec &l, fvec &dl, fvec *isParticleFromVertex=0) const
Definition KFParticleBaseSIMD.cxx:1772
KFParticleBase::DaughterIds
const std::vector< int > & DaughterIds() const
Definition KFParticleBase.h:261
KFParticleBase::GetPDG
int GetPDG() const
Definition KFParticleBase.h:267
KFParticleBase::NDaughters
int NDaughters() const
Definition KFParticleBase.h:260
KFParticleMatch::SetMatchType
void SetMatchType(Short_t i)
Definition KFParticleMatch.h:33
KFParticleMatch::SetMatch
void SetMatch(Int_t i)
Definition KFParticleMatch.h:31
KFParticleSIMD::GetNDF
fvec GetNDF() const
Definition KFParticleSIMD.h:504
KFParticleSIMD::SetProductionVertex
void SetProductionVertex(const KFParticleSIMD &Vtx)
Definition KFParticleSIMD.h:864
KFParticleSIMD::Chi2
const fvec & Chi2() const
Definition KFParticleSIMD.h:149
KFParticle::GetZ
Double_t GetZ() const
Definition KFParticle.h:442
KFParticle::GetPhi
Double_t GetPhi() const
Definition KFParticle.h:529
KFParticle::GetTheta
Double_t GetTheta() const
Definition KFParticle.h:192
KFParticle::Y
const Double_t & Y() const
Definition KFParticle.h:129
KFParticle::GetMass
Double_t GetMass() const
Definition KFParticle.h:547
KFParticle::GetNDF
Int_t GetNDF() const
Definition KFParticle.h:482
KFParticle::GetMomentum
Double_t GetMomentum() const
Definition KFParticle.h:538
KFParticle::GetChi2
Double_t GetChi2() const
Definition KFParticle.h:477
KFParticle::X
const Double_t & X() const
Definition KFParticle.h:128
KFParticle::Z
const Double_t & Z() const
Definition KFParticle.h:130
KFParticle::GetCovariance
Double_t GetCovariance(int i) const
Definition KFParticle.h:492
KFParticle::GetParameter
Double_t GetParameter(int i) const
Definition KFParticle.h:487
KFParticle::GetDecayLength
Double_t GetDecayLength() const
Definition KFParticle.h:556
KFParticle::TransportToPoint
void TransportToPoint(const Double_t xyz[])
Definition KFParticle.h:897
KFParticle::GetPt
Double_t GetPt() const
Definition KFParticle.h:511
KFParticle::GetLifeTime
Double_t GetLifeTime() const
Definition KFParticle.h:574
KFParticle::GetRapidity
Double_t GetRapidity() const
Definition KFParticle.h:191
TDirectoryFile::cd
Bool_t cd() override
TDirectory::GetDirectory
virtual TDirectory * GetDirectory(const char *namecycle, Bool_t printError=false, const char *funcname="GetDirectory")
TDirectory::cd
virtual Bool_t cd()
TDirectory::mkdir
virtual TDirectory * mkdir(const char *name, const char *title="", Bool_t returnExistingDirectory=kFALSE)
TFile::Delete
void Delete(const char *namecycle="") override
TFile::Close
void Close(Option_t *option="") override
TObject::IsFolder
virtual Bool_t IsFolder() const
TObject::GetName
virtual const char * GetName() const
TObject::Write
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
TString::Data
const char * Data() const
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
n
*brief user defined string literal for JSON values *sa std::size_t n
Definition json.hpp:20362
TMath::Prob
Double_t Prob(Double_t chi2, Int_t ndf)
TMath::Sqrt
Double_t Sqrt(Double_t x)
std
Definition json.hpp:4446