18 #include "CbmRichUtil.h"
23 #include "FairMCPoint.h"
28 #include "TClonesArray.h"
33 #include <boost/assign/list_of.hpp>
38 using boost::assign::list_of;
47 : FairTask(
"LitTrackingQA", 1)
51 , fMCTrackCreator(NULL)
63 , fUseConsecutivePointsInSts(true)
80 , fMvdHitMatches(NULL)
84 , fRichProjections(NULL)
85 , fRichRingMatches(NULL)
123 for (Int_t
i = 0;
i < trackVariants.size();
i++) {
125 make_pair(trackVariants[
i], multimap<pair<Int_t, Int_t>, Int_t>()));
136 fHM->
H1(
"hen_EventNo_TrackingQa")->Fill(0.5);
137 Int_t eventNum =
fHM->
H1(
"hen_EventNo_TrackingQa")->GetEntries() - 1;
138 std::cout <<
"CbmLitTrackingQa::Exec: event=" << eventNum << std::endl;
150 TDirectory* oldir = gDirectory;
151 TFile* outFile = FairRootManager::Instance()->GetOutFile();
152 if (outFile != NULL) {
156 gDirectory->cd(oldir->GetPath());
166 FairRootManager* ioman = FairRootManager::Instance();
167 if (NULL == ioman) { Fatal(
"Init",
"CbmRootManager is not instantiated"); }
172 if (NULL ==
fMCTracks) { Fatal(
"Init",
"No MCTrack array!"); }
174 fGlobalTracks = (TClonesArray*) ioman->GetObject(
"GlobalTrack");
175 if (NULL ==
fGlobalTracks) { Fatal(
"Init",
"No GlobalTrack array!"); }
179 if (NULL ==
fMvdPoints) { Fatal(
"Init",
": No MvdPoint array!"); }
181 if (NULL ==
fMvdHitMatches) { Fatal(
"Init",
": No MvdHitMatch array!"); }
185 fStsTracks = (TClonesArray*) ioman->GetObject(
"StsTrack");
186 if (NULL ==
fStsTracks) { Fatal(
"Init",
": No StsTrack array!"); }
187 fStsMatches = (TClonesArray*) ioman->GetObject(
"StsTrackMatch");
188 if (NULL ==
fStsMatches) { Fatal(
"Init",
": No StsTrackMatch array!"); }
192 fRichRings = (TClonesArray*) ioman->GetObject(
"RichRing");
193 if (NULL ==
fRichRings) { Fatal(
"Init",
"No RichRing array!"); }
195 if (NULL ==
fRichProjections) { Fatal(
"Init",
"No RichProjection array!"); }
201 fMuchMatches = (TClonesArray*) ioman->GetObject(
"MuchTrackMatch");
202 if (NULL ==
fMuchMatches) { Fatal(
"Init",
"No MuchTrackMatch array!"); }
206 fTrdMatches = (TClonesArray*) ioman->GetObject(
"TrdTrackMatch");
207 if (NULL ==
fTrdMatches) { Fatal(
"Init",
"No TrdTrackMatch array!"); }
212 if (NULL ==
fTofPoints) { Fatal(
"Init",
"No TofPoint array!"); }
213 fTofMatches = (TClonesArray*) ioman->GetObject(
"TofHitMatch");
214 if (NULL ==
fTofMatches) { Fatal(
"Init",
"No TofHitMatch array!"); }
219 vector<string> tmp = list_of(
"All")(
"Primary")(
"Secondary")(
"Reference")(
221 "PionMinus")(
"KaonPlus")(
"KaonMinus");
226 vector<string> tmp = list_of(
"All")(
"AllReference")(
"Electron")(
227 "ElectronReference")(
"Pion")(
"PionReference");
233 vector<string> tmp = list_of(
"Electron");
240 vector<string> tmp = list_of(
"Electron");
247 vector<string> tmp = list_of(
"All")(
"TrueMatch")(
"WrongMatch");
305 const string& parameter,
306 const string& xTitle,
311 assert(opt ==
"track" || opt ==
"ring" || opt ==
"track_pid"
312 || opt ==
"ring_pid");
313 vector<string> types = list_of(
"Acc")(
"Rec")(
"Eff");
321 for (Int_t iCat = 0; iCat < cat.size(); iCat++) {
322 for (Int_t iType = 0; iType < 3; iType++) {
323 string yTitle = (types[iType] ==
"Eff") ?
"Efficiency [%]" :
"Counter";
325 name +
"_" + cat[iCat] +
"_" + types[iType] +
"_" + parameter;
326 string histTitle = histName +
";" + xTitle +
";" + yTitle;
329 new TH1F(histName.c_str(), histTitle.c_str(), nofBins, minBin, maxBin));
335 const string& parameter,
336 const string& xTitle,
337 const string& yTitle,
345 assert(opt ==
"track" || opt ==
"ring" || opt ==
"track_pid"
346 || opt ==
"ring_pid");
347 vector<string> types = list_of(
"Acc")(
"Rec")(
"Eff");
355 for (Int_t iCat = 0; iCat < cat.size(); iCat++) {
356 for (Int_t iType = 0; iType < 3; iType++) {
357 string zTitle = (types[iType] ==
"Eff") ?
"Efficiency [%]" :
"Counter";
359 name +
"_" + cat[iCat] +
"_" + types[iType] +
"_" + parameter;
360 string histTitle = histName +
";" + xTitle +
";" + yTitle +
";" + zTitle;
362 new TH2F(histName.c_str(),
375 const string& parameter,
376 const string& xTitle,
380 vector<string> types = list_of(
"RecPions")(
"Rec")(
"PionSup");
382 for (Int_t iType = 0; iType < 3; iType++) {
384 (types[iType] ==
"PionSup") ?
"Pion suppression" :
"Counter";
386 + types[iType] +
"_" + parameter;
387 string histTitle = histName +
";" + xTitle +
";" + yTitle;
390 new TH1F(histName.c_str(), histTitle.c_str(), nofBins, minBin, maxBin));
396 const string& xTitle,
397 const string& yTitle,
401 TH1F*
h =
new TH1F(name.c_str(),
402 string(name +
";" + xTitle +
";" + yTitle).c_str(),
410 const string& xTitle,
411 const string& yTitle,
412 const string& zTitle,
420 new TH2F(name.c_str(),
421 (name +
";" + xTitle +
";" + yTitle +
";" + zTitle).c_str(),
432 string type[] = {
"All",
"True",
"Fake",
"TrueOverAll",
"FakeOverAll"};
433 Double_t
min[] = {-0.5, -0.5, -0.5, -0.1, -0.1};
434 Double_t
max[] = {99.5, 99.5, 99.5, 1.1, 1.1};
435 Int_t bins[] = {100, 100, 100, 12, 12};
436 for (Int_t
i = 0;
i < 5;
i++) {
437 string xTitle = (
i == 3 ||
i == 4) ?
"Ratio" :
"Number of hits";
438 string histName =
"hth_" + detName +
"_TrackHits_" + type[
i];
444 const vector<string>& detectors) {
445 vector<string> histos;
446 Int_t nofDetectors = detectors.size();
447 for (Int_t iDet = 0; iDet < nofDetectors; iDet++) {
449 for (Int_t
i = 0;
i <= iDet;
i++) {
450 histEff += detectors[
i];
452 string histNorm = histEff;
453 histos.push_back(
"hte_" + histEff +
"_" + histNorm);
454 for (Int_t
i = iDet + 1;
i < nofDetectors;
i++) {
455 histNorm += detectors[
i];
456 histos.push_back(
"hte_" + histEff +
"_" + histNorm);
464 vector<string> detectors;
472 vector<string> names2;
483 names.insert(names1.begin(), names1.end());
484 names.insert(names2.begin(), names2.end());
485 vector<string> nameVector(names.begin(), names.end());
490 set<string> trackVariants;
492 vector<string> detectors;
498 for (Int_t
i = 0;
i < detectors.size();
i++) {
499 name += detectors[
i];
500 trackVariants.insert(name);
511 for (Int_t
i = 0;
i < detectors.size();
i++) {
512 name += detectors[
i];
513 trackVariants.insert(name);
516 vector<string> trackVariantsVector(trackVariants.begin(),
517 trackVariants.end());
519 trackVariantsVector.push_back(
"Rich");
521 return trackVariantsVector;
525 set<string> variants;
527 vector<string> detectors;
531 for (Int_t
i = 0;
i < detectors.size();
i++) {
532 name += detectors[
i];
533 variants.insert(name);
543 for (Int_t
i = 0;
i < detectors.size();
i++) {
544 name += detectors[
i];
545 variants.insert(name);
548 vector<string> variantsVector(variants.begin(), variants.end());
549 return variantsVector;
553 set<string> trackVariants;
555 vector<string> detectors;
561 for (Int_t
i = 0;
i < detectors.size();
i++) {
562 name += detectors[
i];
563 if (detectors[
i] == detName)
break;
572 Double_t minNofPoints = 0.;
573 Double_t maxNofPoints = 100.;
574 Int_t nofBinsPoints = 100;
597 string histName =
"hte_Much_" + norm;
638 string histName =
"hte_Trd_" + norm;
679 string histName =
"hte_Tof_" + norm;
744 "hte_Rich_Rich",
"RingXc",
"X [cm]", 60, -120, 120,
"ring");
746 "hte_Rich_Rich",
"RingYc",
"|Y| [cm]", 25, 100, 200,
"ring");
774 for (Int_t iHist = 0; iHist < histoNames.size(); iHist++) {
775 string name = histoNames[iHist];
776 string opt = (name.find(
"Rich") == string::npos) ?
"track" :
"ring";
797 if (name.find(
"Rich") != string::npos) {
818 for (Int_t iHist = 0; iHist < psVariants.size(); iHist++) {
819 string name =
"hps_" + psVariants[iHist];
882 CreateH2(
"hng_NofGhosts_Rich_RingXcYc",
892 CreateH1(
"hng_NofGhosts_RichStsMatching_Nh",
898 CreateH1(
"hng_NofGhosts_RichElId_Nh",
904 CreateH1(
"hng_NofGhosts_StsRichMatching_Nh",
921 Int_t nofBinsC = 100000;
922 Double_t maxXC = 100000.;
923 CreateH1(
"hno_NofObjects_GlobalTracks",
930 CreateH1(
"hno_NofObjects_StsTracks",
937 CreateH1(
"hno_NofObjects_TrdTracks",
944 CreateH1(
"hno_NofObjects_MuchTracks",
951 CreateH1(
"hno_NofObjects_RichRings",
957 CreateH1(
"hno_NofObjects_RichProjections",
958 "Projections per event",
966 CreateH1(
"hen_EventNo_TrackingQa",
"",
"", 1, 0, 1.);
974 map<string, multimap<pair<Int_t, Int_t>, Int_t>>::iterator it;
976 multimap<pair<Int_t, Int_t>, Int_t>& mcRecoMap = (*it).second;
986 for (Int_t iTrack = 0; iTrack < nofGlobalTracks; iTrack++) {
1004 Double_t rtDistance = 10;
1017 fHM->
H1(
"hng_NofGhosts_Sts_Nh")->Fill(nofHits);
1025 fHM->
H1(
"hng_NofGhosts_StsRichMatching_Nh")->Fill(nofHits);
1041 fHM->
H1(
"hng_NofGhosts_Trd_Nh")->Fill(nofHits);
1051 Int_t nofHits = muchTrackMatch->
GetNofHits();
1056 fHM->
H1(
"hng_NofGhosts_Much_Nh")->Fill(nofHits);
1070 fHM->
H1(
"hng_NofGhosts_Rich_Nh")->Fill(nofHits);
1076 fHM->
H1(
"hng_NofGhosts_Rich_RingXcYc")
1080 if (rtDistance < 1.)
1081 fHM->
H1(
"hng_NofGhosts_RichStsMatching_Nh")->Fill(nofHits);
1083 Double_t momentumMc = 0.;
1084 if (stsTrackMatch != NULL) {
1087 if (mcTrack != NULL) momentumMc = mcTrack->
GetP();
1091 iTrack, momentumMc)) {
1092 fHM->
H1(
"hng_NofGhosts_RichElId_Nh")->Fill(nofHits);
1099 pair<Int_t, Int_t> stsMCId = {-1, -1}, trdMCId = {-1, -1},
1100 muchMCId = {-1, -1}, richMCId = {-1, -1};
1101 list<pair<Int_t, Int_t>> tofMCIds;
1117 const vector<CbmLink>& tofMCLinks = tofMatch->
GetLinks();
1119 for (vector<CbmLink>::const_iterator tofMCIt = tofMCLinks.begin();
1120 tofMCIt != tofMCLinks.end();
1122 const FairMCPoint* tofPoint =
1125 if (tofPoint != NULL)
1127 pair<Int_t, Int_t>(tofMCIt->GetEntry(), tofPoint->GetTrackID()));
1135 map<string, multimap<pair<Int_t, Int_t>, Int_t>>::iterator it;
1137 string name = (*it).first;
1138 multimap<pair<Int_t, Int_t>, Int_t>& mcRecoMap = (*it).second;
1139 Bool_t sts = (name.find(
"Sts") != string::npos)
1140 ? isStsOk && stsMCId.second != -1
1142 Bool_t trd = (name.find(
"Trd") != string::npos)
1143 ? isTrdOk && stsMCId == trdMCId
1145 Bool_t much = (name.find(
"Much") != string::npos)
1146 ? isMuchOk && stsMCId == muchMCId
1148 Bool_t tof = (name.find(
"Tof") != string::npos)
1150 && find(tofMCIds.begin(), tofMCIds.end(), stsMCId)
1153 Bool_t rich = (name.find(
"Rich") != string::npos)
1154 ? isRichOk && stsMCId == richMCId
1157 if (sts && trd && much && tof && rich) {
1158 pair<pair<Int_t, Int_t>, Int_t> tmp = make_pair(stsMCId, iTrack);
1159 mcRecoMap.insert(tmp);
1167 Int_t nofRings =
fRichRings->GetEntriesFast();
1168 for (Int_t iRing = 0; iRing < nofRings; iRing++) {
1172 pair<Int_t, Int_t> richMCId = {
1176 if (isRichOk && -1 != richMCId.second) {
1177 pair<pair<Int_t, Int_t>, Int_t> tmp = make_pair(richMCId, iRing);
1188 fHM->
H1(
"hth_Mvd_TrackHits_All")->Fill(nofHits);
1194 Int_t nofTrueHits = 0, nofFakeHits = 0;
1195 for (Int_t iHit = 0; iHit < nofHits; iHit++) {
1199 if (NULL == hitMatch)
continue;
1202 const FairMCPoint* point =
1203 static_cast<const FairMCPoint*
>(
fMvdPoints->
Get(0, eventId, pointId));
1204 if (NULL == point)
continue;
1205 Int_t mcTrackId = point->GetTrackID();
1206 if (mcTrackId == stsMcTrackId) {
1212 fHM->
H1(
"hth_Mvd_TrackHits_True")->Fill(nofTrueHits);
1213 fHM->
H1(
"hth_Mvd_TrackHits_Fake")->Fill(nofFakeHits);
1215 fHM->
H1(
"hth_Mvd_TrackHits_TrueOverAll")
1216 ->Fill(Double_t(nofTrueHits) / Double_t(nofHits));
1217 fHM->
H1(
"hth_Mvd_TrackHits_FakeOverAll")
1218 ->Fill(Double_t(nofFakeHits) / Double_t(nofHits));
1232 assert(detName !=
"");
1233 string histName =
"hth_" + detName +
"_TrackHits";
1237 fHM->
H1(histName +
"_TrueOverAll")
1239 fHM->
H1(histName +
"_FakeOverAll")
1244 vector<TH1*> effHistos =
fHM->
H1Vector(
"(hte|hpe)_.*_Eff_.*");
1245 Int_t nofEffHistos = effHistos.size();
1248 Int_t nofExistsMcTracks = 0;
1249 for (Int_t iMCTrack = 0; iMCTrack < nofMcTracks; ++iMCTrack) {
1258 Int_t nofPointsSts =
1260 Int_t nofPointsTrd =
1262 Int_t nofPointsMuch =
1274 Bool_t stsConsecutive =
1293 Double_t angle = std::abs(mom.Theta() * 180 /
TMath::Pi());
1294 Double_t mcP = mcTrack->
GetP();
1296 Double_t mcPt = mcTrack->
GetPt();
1299 map<string, vector<Double_t>> parMap;
1301 vector<Double_t> tmp = list_of(mcP);
1304 vector<Double_t> tmp1 = list_of(mcY);
1307 vector<Double_t> tmp2 = list_of(mcPt);
1308 parMap[
"pt"] = tmp2;
1310 vector<Double_t> tmp3 = list_of(angle);
1311 parMap[
"Angle"] = tmp3;
1313 vector<Double_t> tmp4 = list_of(0);
1318 vector<Double_t> tmp5 = list_of(boa);
1319 parMap[
"BoA"] = tmp5;
1321 vector<Double_t> tmp5X = list_of(ringX);
1322 parMap[
"RingXc"] = tmp5X;
1324 vector<Double_t> tmp5Y = list_of(std::abs(ringY));
1325 parMap[
"RingYc"] = tmp5Y;
1327 vector<Double_t> tmp6 = list_of(ringX)(ringY);
1328 parMap[
"RingXcYc"] = tmp6;
1330 vector<Double_t> tmp7 = list_of(nofHitsRich);
1331 parMap[
"RingNh"] = tmp7;
1333 vector<Double_t> tmp8 = list_of(mcY)(mcPt);
1335 parMap[
"YPt"] = tmp8;
1337 for (Int_t iHist = 0; iHist < nofEffHistos; iHist++) {
1338 TH1* hist = effHistos[iHist];
1339 string histName = hist->GetName();
1341 string effName =
split[1];
1342 string normName =
split[2];
1343 string catName =
split[3];
1344 string histTypeName =
split[0];
1345 string parName =
split[5];
1346 assert(parMap.count(parName) != 0);
1348 vector<Double_t> par = list_of(0);
1349 if (parName ==
"Np") {
1350 vector<Double_t> tmp = list_of(
1353 : (effName ==
"Trd")
1355 : (effName ==
"Much") ? nofPointsMuch
1356 : (effName ==
"Tof") ? nofPointsTof : 0);
1359 par = parMap[parName];
1362 Bool_t sts = (normName.find(
"Sts") != string::npos) ? isStsOk :
true;
1363 Bool_t trd = (normName.find(
"Trd") != string::npos) ? isTrdOk :
true;
1364 Bool_t much = (normName.find(
"Much") != string::npos) ? isMuchOk :
true;
1365 Bool_t tof = (normName.find(
"Tof") != string::npos) ? isTofOk :
true;
1366 Bool_t rich = (normName.find(
"Rich") != string::npos) ? isRichOk :
true;
1368 if (effName ==
"Trd" || effName ==
"Much" || effName ==
"Tof") {
1371 fMcToRecoMap[prevRecName].find(make_pair(iEvent, iMCTrack))
1373 Bool_t accOk = isPrevRec && sts && trd && much && tof && rich;
1385 Bool_t accOk = sts && trd && much && tof && rich;
1387 if (histName.find(
"Rich") == string::npos) {
1415 const multimap<pair<Int_t, Int_t>, Int_t>& mcMap,
1416 const string& histName,
1417 const string& histTypeName,
1418 const string& effName,
1419 const string& catName,
1420 const vector<Double_t>& par) {
1425 Bool_t accOk =
function(
fMCTracks, mcEventId, mcId);
1427 (histTypeName ==
"hte")
1428 ? (mcMap.find(pair<Int_t, Int_t>(mcEventId, mcId)) != mcMap.end()
1430 : (
ElectronId(mcEventId, mcId, mcMap, effName) && accOk);
1431 Int_t nofParams = par.size();
1432 assert(nofParams < 3 && nofParams > 0);
1433 if (nofParams == 1) {
1434 if (accOk) {
fHM->
H1(accHistName)->Fill(par[0]); }
1435 if (recOk) {
fHM->
H1(recHistName)->Fill(par[0]); }
1436 }
else if (nofParams == 2) {
1437 if (accOk) {
fHM->
H1(accHistName)->Fill(par[0], par[1]); }
1438 if (recOk) {
fHM->
H1(recHistName)->Fill(par[0], par[1]); }
1445 const multimap<pair<Int_t, Int_t>, Int_t>& mcMap,
1446 const string& histName,
1447 const string& histTypeName,
1448 const string& effName,
1449 const string& catName,
1450 const vector<Double_t>& par) {
1451 Int_t nofHitsInRing =
1457 Bool_t accOk =
function(
fMCTracks, mcEventId, mcId, nofHitsInRing);
1459 (histTypeName ==
"hte")
1460 ? (mcMap.find(make_pair(mcEventId, mcId)) != mcMap.end() && accOk)
1461 : (
ElectronId(mcEventId, mcId, mcMap, effName) && accOk);
1462 Int_t nofParams = par.size();
1463 assert(nofParams < 3 && nofParams > 0);
1464 if (nofParams == 1) {
1465 if (accOk) {
fHM->
H1(accHistName)->Fill(par[0]); }
1466 if (recOk) {
fHM->
H1(recHistName)->Fill(par[0]); }
1467 }
else if (nofParams == 2) {
1468 if (accOk) {
fHM->
H1(accHistName)->Fill(par[0], par[1]); }
1469 if (recOk) {
fHM->
H1(recHistName)->Fill(par[0], par[1]); }
1476 const multimap<pair<Int_t, Int_t>, Int_t>& mcMap,
1477 const string& effName) {
1478 multimap<pair<Int_t, Int_t>, Int_t>::const_iterator it =
1479 mcMap.find(pair<Int_t, Int_t>(mcEventId, mcId));
1480 if (it == mcMap.end())
return false;
1481 Int_t globalTrackIndex = (*it).second;
1486 Bool_t isRichElectron =
1491 Bool_t isTrdElectron =
1496 Bool_t isTofElectron =
1501 return isRichElectron && isTrdElectron && isTofElectron;
1506 vector<TH1*> histos =
fHM->
H1Vector(
"hps_.*_RecPions_.*");
1507 Int_t nofHistos = histos.size();
1510 for (Int_t iGT = 0; iGT < nofGlobalTracks; iGT++) {
1514 if (stsIndex < 0)
continue;
1517 const FairTrackParam* richProjection =
1519 if (richProjection == NULL || richProjection->GetX() == 0
1520 || richProjection->GetY() == 0)
1530 if (mcIdSts < 0)
continue;
1535 Double_t p = mom.Mag();
1537 if (std::abs(pdgSts) == 211) {
1538 Bool_t isRichElectron =
1542 Bool_t isTrdElectron =
1546 Bool_t isTofElectron =
1551 for (Int_t iHist = 0; iHist < nofHistos; iHist++) {
1552 histos[iHist]->Fill(p);
1553 string name = histos[iHist]->GetName();
1555 string effName =
split[1];
1556 string category =
split[2];
1559 ((effName.find(
"Rich") != string::npos) ? isRichElectron :
true)
1560 && ((effName.find(
"Trd") != string::npos) ? isTrdElectron :
true)
1561 && ((effName.find(
"Tof") != string::npos) ? isTofElectron :
true);
1578 fHM->
H1(
"hno_NofObjects_GlobalTracks")
1581 fHM->
H1(
"hno_NofObjects_StsTracks")->Fill(
fStsTracks->GetEntriesFast());
1584 fHM->
H1(
"hno_NofObjects_RichRings")->Fill(
fRichRings->GetEntriesFast());
1585 Int_t projCount = 0;
1587 for (Int_t iProj = 0; iProj < nProj; iProj++) {
1589 if (NULL == proj || proj->GetX() == 0. || proj->GetY() == 0)
continue;
1592 fHM->
H1(
"hno_NofObjects_RichProjections")->Fill(projCount);