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