10 #include "FairRunAna.h"
15 #include "KFMCTrack.h"
17 #include "KFParticleTopoReconstructor.h"
18 #include "KFTopoPerformance.h"
22 #include "TClonesArray.h"
23 #include "TDatabasePDG.h"
24 #include "TDirectory.h"
42 #include "CbmTofHit.h"
54 : FairTask(name, iVerbose)
55 , fStsTrackBranchName(
"StsTrack")
56 , fGlobalTrackBranchName(
"GlobalTrack")
57 , fRichBranchName(
"RichRing")
58 , fTrdBranchName(
"TrdTrack")
59 , fTrdHitBranchName(
"TrdHit")
60 , fTofBranchName(
"TofHit")
61 , fMuchTrackBranchName(
"MuchTrack")
62 , fMCTracksBranchName(
"MCTrack")
63 , fStsTrackMatchBranchName(
"StsTrackMatch")
64 , fRichRingMatchBranchName(
"RichRingMatch")
65 , fTrdTrackMatchBranchName(
"TrdTrackMatch")
66 , fTofHitMatchBranchName(
"TofHitMatch")
67 , fMuchTrackMatchBranchName(
"MuchTrackMatch")
69 , fGlobalTrackArray(0)
76 , fStsTrackMatchArray(0)
77 , fRichRingMatchArray(0)
78 , fTrdTrackMatchArray(0)
79 , fTofHitMatchArray(0)
80 , fMuchTrackMatchArray(0)
81 , fOutFileName(outFileName)
86 TFile* curFile = gFile;
87 TDirectory* curDirectory = gDirectory;
99 gDirectory->mkdir(
"STS");
100 gDirectory->cd(
"STS");
102 TString histoName[
NStsHisto] = {
"NHits",
"chi2/NDF",
"prob"};
103 TString axisName[
NStsHisto] = {
"N hits",
"#chi^{2}/NDF",
"prob"};
105 float xMin[
NStsHisto] = {-0.5, 0.f, 0.f};
106 float xMax[
NStsHisto] = {15.5, 20.f, 1.f};
108 TString subdirs[8] = {
109 "Tracks",
"e",
"mu",
"pi",
"K",
"p",
"fragments",
"ghost"};
111 for (
int iDir = 0; iDir < 8; iDir++) {
112 gDirectory->mkdir(subdirs[iDir].Data());
113 gDirectory->cd(subdirs[iDir].Data());
115 gDirectory->mkdir(
"TrackFitQA");
116 gDirectory->cd(
"TrackFitQA");
119 TString pull =
"pull";
121 TString parName[5] = {
"X",
"Y",
"Tx",
"Ty",
"QP"};
124 float xMaxFit[5] = {0.05, 0.045, 0.01, 0.01, 0.1};
126 for (
int iH = 0; iH < 5; iH++) {
127 hStsFitHisto[iDir][iH] =
new TH1F((res + parName[iH]).Data(),
128 (res + parName[iH]).Data(),
132 hStsFitHisto[iDir][iH + 5] =
new TH1F((pull + parName[iH]).Data(),
133 (pull + parName[iH]).Data(),
139 gDirectory->cd(
"..");
142 hStsHisto[iDir][iH] =
new TH1F(histoName[iH].Data(),
143 histoName[iH].Data(),
147 hStsHisto[iDir][iH]->GetXaxis()->SetTitle(axisName[iH].Data());
150 gDirectory->cd(
"..");
153 gDirectory->cd(
"..");
155 gDirectory->mkdir(
"MuCh");
156 gDirectory->cd(
"MuCh");
159 "NHits",
"FirstStation",
"LastStation",
"chi2/NDF",
"prob"};
161 "N hits",
"First Station",
"Last Station",
"#chi^{2}/NDF",
"prob"};
162 int nBins[
NMuchHisto] = {16, 16, 16, 100, 100};
163 float xMin[
NMuchHisto] = {-0.5f, -0.5f, -0.5f, 0.f, 0.f};
164 float xMax[
NMuchHisto] = {15.5f, 15.5f, 15.5f, 20.f, 1.f};
166 TString subdirs[3] = {
"mu",
"Background",
"Ghost"};
168 for (
int iDir = 0; iDir < 3; iDir++) {
169 gDirectory->mkdir(subdirs[iDir].Data());
170 gDirectory->cd(subdirs[iDir].Data());
172 hMuchHisto[iDir][iH] =
new TH1F(histoName[iH].Data(),
173 histoName[iH].Data(),
177 hMuchHisto[iDir][iH]->GetXaxis()->SetTitle(axisName[iH].Data());
179 gDirectory->cd(
"..");
182 gDirectory->cd(
"..");
184 gDirectory->mkdir(
"RICH");
185 gDirectory->cd(
"RICH");
187 TString subdirs[10] = {
"AllTracks",
199 for (
int iDir = 0; iDir < 10; iDir++) {
200 gDirectory->mkdir(subdirs[iDir]);
201 gDirectory->cd(subdirs[iDir]);
204 histoName2D[iH], histoName2D[iH], 1000, 0, 15., 1000, 0, 10.);
209 gDirectory->cd(
"..");
212 gDirectory->cd(
"..");
214 gDirectory->mkdir(
"TRD");
215 gDirectory->cd(
"TRD");
217 TString histoName[
NTrdHisto] = {
"Wkn",
"ANN"};
218 TString axisName[
NTrdHisto] = {
"Wkn",
"ANN"};
223 TString subdirs[14] = {
"AllTracks",
238 for (
int iDir = 0; iDir < 14; iDir++) {
239 gDirectory->mkdir(subdirs[iDir].Data());
240 gDirectory->cd(subdirs[iDir].Data());
242 hTrdHisto[iDir][iH] =
new TH1F(histoName[iH].Data(),
243 histoName[iH].Data(),
247 hTrdHisto[iDir][iH]->GetXaxis()->SetTitle(axisName[iH].Data());
250 new TH2F(
"dE/dx",
"dE/dx", 1500, 0., 15., 1000, 0., 1000.);
251 hTrdHisto2D[iDir][0]->GetYaxis()->SetTitle(
"dE/dx [keV/(cm)]");
252 hTrdHisto2D[iDir][0]->GetXaxis()->SetTitle(
"p [GeV/c]");
253 gDirectory->cd(
"..");
256 gDirectory->cd(
"..");
258 gDirectory->mkdir(
"TOF");
259 gDirectory->cd(
"TOF");
261 TString histoName2D[2] = {
"M2P",
"M2dEdX"};
262 TString xAxisName[
NTofHisto2D] = {
"p [GeV/c]",
"dE/dx [keV/(cm)]"};
263 TString yAxisName[
NTofHisto2D] = {
"m^{2} [GeV^{2}/c^{4}]",
264 "m^{2} [GeV^{2}/c^{4}]"};
272 TString subdirs[14] = {
"AllTracks",
287 for (
int iDir = 0; iDir < 14; iDir++) {
288 gDirectory->mkdir(subdirs[iDir].Data());
289 gDirectory->cd(subdirs[iDir].Data());
292 hTofHisto2D[iDir][iH] =
new TH2F(histoName2D[iH].Data(),
293 histoName2D[iH].Data(),
300 hTofHisto2D[iDir][iH]->GetXaxis()->SetTitle(xAxisName[iH]);
301 hTofHisto2D[iDir][iH]->GetYaxis()->SetTitle(yAxisName[iH]);
304 gDirectory->cd(
"..");
307 gDirectory->cd(
"..");
310 gDirectory = curDirectory;
327 FairRootManager* ioman = FairRootManager::Instance();
330 Warning(
"CbmKFTrackQA::Init",
"RootManager not instantiated!");
337 Warning(
"CbmKFTrackQA::Init",
"track-array not found!");
344 Warning(
"CbmKFTrackQA::Init",
"global track array not found!");
349 Warning(
"CbmKFTrackQA::Init",
"TOF hit-array not found!");
354 Warning(
"CbmKFTrackQA::Init",
"TRD track-array not found!");
358 Warning(
"CbmKFTrackQA::Init",
"TRD hit-array not found!");
362 Warning(
"CbmKFTrackQA::Init",
"Rich ring array not found!");
366 Warning(
"CbmKFTrackQA::Init",
"mc track array not found!");
374 Warning(
"CbmKFTrackQA::Init",
"track match array not found!");
382 Warning(
"CbmKFTrackQA::Init",
"RichRing match array not found!");
387 Warning(
"CbmKFTrackQA::Init",
"TofHit match array not found!");
393 Warning(
"CbmKFTrackQA::Init",
"TrdTrack match array not found!");
399 Warning(
"CbmKFTrackQA::Init",
"Much track match array not found!");
404 Warning(
"CbmKFTrackQA::Init",
"Much track-array not found!");
410 if (mcManager == 0) {
411 Warning(
"CbmKFTrackQA::Init",
"mc manager not found!");
417 Warning(
"CbmKFTrackQA::Init",
"tof points not found!");
428 for (Int_t iMC = 0; iMC <
nMCTracks; iMC++) {
442 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg)))
443 q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0;
444 else if (pdg == 1000010020)
446 else if (pdg == -1000010020)
448 else if (pdg == 1000010030)
450 else if (pdg == -1000010030)
452 else if (pdg == 1000020030)
454 else if (pdg == -1000020030)
456 else if (pdg == 1000020040)
458 else if (pdg == -1000020040)
462 Double_t p = cbmMCTrack->
GetP();
471 vector<int> trackMatch(ntrackMatches, -1);
473 for (
int iTr = 0; iTr < ntrackMatches; iTr++) {
477 Float_t bestWeight = 0.f;
478 Float_t totalWeight = 0.f;
479 Int_t mcTrackId = -1;
480 for (
int iLink = 0; iLink < stsTrackMatch->
GetNofLinks(); iLink++) {
487 if (bestWeight / totalWeight < 0.7)
continue;
488 if (mcTrackId >=
nMCTracks || mcTrackId < 0) {
489 std::cout <<
"Sts Matching is wrong! StsTackId = " << mcTrackId
490 <<
" N mc tracks = " <<
nMCTracks << std::endl;
494 mcTracks[mcTrackId].SetReconstructed();
495 trackMatch[iTr] = mcTrackId;
501 for (
int iTr = 0; iTr <
fStsTrackArray->GetEntriesFast(); iTr++) {
503 vRTracks[iTr] = *stsTrack;
505 if (trackMatch[iTr] > -1) pdg[iTr] =
mcTracks[trackMatch[iTr]].PDG();
511 vector<float> vChiToPrimVtx;
512 vector<L1FieldRegion> vField;
513 fitter.
Fit(vRTracks, pdg);
514 fitter.
GetChiToVertex(vRTracks, vField, vChiToPrimVtx, kfVertex, 3000000);
516 for (
unsigned int iTr = 0; iTr < vRTracks.size(); iTr++) {
517 if (trackMatch[iTr] < 0)
continue;
519 const KFMCTrack& mcTrack =
mcTracks[trackMatch[iTr]];
520 if (mcTrack.MotherId() > -1)
continue;
523 const FairTrackParam* parameters = vRTracks[iTr].GetParamFirst();
525 Double_t recoParam[5] = {parameters->GetX(),
529 parameters->GetQp()};
530 Double_t recoError[5] = {parameters->GetCovariance(0, 0),
531 parameters->GetCovariance(1, 1),
532 parameters->GetCovariance(2, 2),
533 parameters->GetCovariance(3, 3),
534 parameters->GetCovariance(4, 4)};
535 Double_t mcParam[5] = {mcTrack.X(),
537 mcTrack.Px() / mcTrack.Pz(),
538 mcTrack.Py() / mcTrack.Pz(),
544 for (
int iParam = 0; iParam < 5; iParam++) {
545 Double_t residual = recoParam[iParam] - mcParam[iParam];
547 Double_t pReco =
fabs(1. / recoParam[iParam]);
548 Double_t pMC =
fabs(1. / mcParam[iParam]);
557 if (recoError[iParam] >= 0.) {
558 Double_t pull = residual /
sqrt(recoError[iParam]);
568 vector<int> trackMuchMatch;
571 trackMuchMatch.resize(nMuchTrackMatches, -1);
573 for (
int iTr = 0; iTr < nMuchTrackMatches; iTr++) {
577 Float_t bestWeight = 0.f;
578 Float_t totalWeight = 0.f;
579 Int_t mcTrackId = -1;
580 for (
int iLink = 0; iLink < muchTrackMatch->
GetNofLinks(); iLink++) {
587 if (bestWeight / totalWeight < 0.7)
continue;
588 if (mcTrackId >=
nMCTracks || mcTrackId < 0) {
589 std::cout <<
"Much Matching is wrong! MuchTackId = " << mcTrackId
590 <<
" N mc tracks = " <<
nMCTracks << std::endl;
594 trackMuchMatch[iTr] = mcTrackId;
599 Warning(
"KF Track QA",
"No GlobalTrack array!");
608 int stsTrackMCIndex = trackMatch[stsTrackIndex];
609 Double_t* stsHistoData =
new Double_t[
NStsHisto];
615 if (stsTrackMCIndex > -1) {
619 hStsHisto[iDir][iH]->Fill(stsHistoData[iH]);
620 hStsHisto[0][iH]->Fill(stsHistoData[iH]);
626 delete[] stsHistoData;
630 if (muchIndex > -1) {
634 int muchTrackMCIndex = trackMuchMatch[muchIndex];
636 Double_t* muchHistoData =
new Double_t[
NMuchHisto];
647 if (stsTrackMCIndex < 0 || stsTrackMCIndex != muchTrackMCIndex)
651 if (TMath::Abs(
mcTracks[stsTrackMCIndex].PDG()) == 13)
658 delete[] muchHistoData;
664 stsPar->Momentum(mom);
666 Double_t p = mom.Mag();
667 Double_t pt = mom.Perp();
668 Double_t pz =
sqrt(p * p - pt * pt);
672 if (richIndex > -1) {
675 int richTrackMCIndex = -1;
680 float bestWeight = 0.f;
681 float totalWeight = 0.f;
682 int bestMCTrackId = -1;
683 for (
int iLink = 0; iLink < richRingMatch->
GetNofLinks();
691 if (bestWeight / totalWeight >= 0.7)
692 richTrackMCIndex = bestMCTrackId;
697 Double_t axisA = richRing->
GetAaxis();
698 Double_t axisB = richRing->
GetBaxis();
704 int iTrackCategory = -1;
705 if (stsTrackMCIndex < 0)
707 else if (richTrackMCIndex < 0)
709 else if (stsTrackMCIndex != richTrackMCIndex)
714 if (iTrackCategory < 10) {
726 vector<int> trackTrdMatch(ntrackTrdMatches, -1);
736 Int_t trdTrackMCIndex = -1;
739 for (
int iTr = 0; iTr < ntrackTrdMatches; iTr++) {
741 Float_t bestWeight = 0.f;
742 Float_t totalWeight = 0.f;
743 for (
int iLink = 0; iLink < trdTrackMatch->
GetNofLinks();
751 if (bestWeight / totalWeight < 0.7)
continue;
752 if (trdTrackMCIndex >=
nMCTracks || trdTrackMCIndex < 0) {
753 std::cout <<
"Trd Matching is wrong! TrdTackId = "
754 << trdTrackMCIndex <<
" N mc tracks = " <<
nMCTracks
759 mcTracks[trdTrackMCIndex].SetReconstructed();
760 trackTrdMatch[iTr] = trdTrackMCIndex;
763 Double_t* trdHistoData =
new Double_t[
NTrdHisto];
768 for (Int_t iTRD = 0; iTRD < trdTrack->
GetNofHits(); iTRD++) {
777 int iTrackCategory = -1;
778 if (stsTrackMCIndex < 0)
780 else if (trdTrackMCIndex < 0)
782 else if (stsTrackMCIndex != trdTrackMCIndex)
788 hTrdHisto[0][iH]->Fill(trdHistoData[iH]);
789 hTrdHisto[iTrackCategory][iH]->Fill(trdHistoData[iH]);
791 dedx = 1e6 * (pz / p) * eloss;
793 hTrdHisto2D[iTrackCategory][0]->Fill(mom.Mag(), dedx);
795 delete[] trdHistoData;
805 if (tofIndex > -1 && stsIndex > -1) {
812 Double_t time = tofHit->
GetTime();
813 Double_t q = stsPar->GetQp() > 0 ? 1. : -1.;
816 * (1. / ((l / time / 29.9792458) * (l / time / 29.9792458)) - 1.);
821 if (tofHit && tofHitMatch && stsMatch) {
826 Int_t tofMCTrackId = tofPoint->GetTrackID();
828 int iHitCategory = -1;
829 if (stsTrackMCIndex < 0)
831 else if (tofMCTrackId < 0)
833 else if (stsTrackMCIndex != tofMCTrackId)
839 hTofHisto2D[iHitCategory][1]->Fill(eloss * 1e6, m2);
848 TDirectory* curr = gDirectory;
849 TFile* currentFile = gFile;
865 if (!obj->IsFolder())
868 TDirectory* cur = gDirectory;
869 TFile* currentFile = gFile;
871 TDirectory* sub = cur->GetDirectory(obj->GetName());
873 TList* listSub = (
static_cast<TDirectory*
>(obj))->GetList();
875 while (TObject* obj1 = it())
884 map<int, int>::iterator it;
893 if (TMath::Abs(getZ - 145) <= 2.0)
return 1;
894 if (TMath::Abs(getZ - 155) <= 2.0)
return 2;
895 if (TMath::Abs(getZ - 165) <= 2.0)
return 3;
896 if (TMath::Abs(getZ - 195) <= 2.0)
return 4;
897 if (TMath::Abs(getZ - 205) <= 2.0)
return 5;
898 if (TMath::Abs(getZ - 215) <= 2.0)
return 6;
899 if (TMath::Abs(getZ - 245) <= 2.0)
return 7;
900 if (TMath::Abs(getZ - 255) <= 2.0)
return 8;
901 if (TMath::Abs(getZ - 265) <= 2.0)
return 9;
902 if (TMath::Abs(getZ - 305) <= 2.0)
return 10;
903 if (TMath::Abs(getZ - 315) <= 2.0)
return 11;
904 if (TMath::Abs(getZ - 325) <= 2.0)
return 12;
905 if (TMath::Abs(getZ - 370) <= 2.0)
return 13;
906 if (TMath::Abs(getZ - 380) <= 2.0)
return 14;
907 if (TMath::Abs(getZ - 390) <= 2.0)
return 15;
908 if (TMath::Abs(getZ - 500) <= 2.0)
return 16;
909 if (TMath::Abs(getZ - 510) <= 2.0)
return 17;
910 if (TMath::Abs(getZ - 520) <= 2.0)
return 18;