diff --git a/core/qa/CbmQaTask.cxx b/core/qa/CbmQaTask.cxx
index d01349387c578253215e3005bf5f9895e59ff478..4bc40d55986ceb2c1ff3f7e4a91a0b77b4b04268 100644
--- a/core/qa/CbmQaTask.cxx
+++ b/core/qa/CbmQaTask.cxx
@@ -37,7 +37,6 @@ CbmQaTask::CbmQaTask(const char* name, const char* prefix, int verbose, bool isM
void CbmQaTask::Exec(Option_t* /*option*/)
{
fNofEvents.SetVal(fNofEvents.GetVal() + 1);
- LOG_IF(info, fVerbose > 1) << fName << ": event " << fNofEvents.GetVal();
this->InitTimeSlice();
this->FillHistograms();
}
diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C
index a8d95846eae557a9a50406d6de46bd9887735e93..310a3e3b6261db1ac503c7bf7240416858155cd3 100644
--- a/macro/mcbm/mcbm_qa.C
+++ b/macro/mcbm/mcbm_qa.C
@@ -36,7 +36,7 @@
#endif
void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test",
- TString setupName = "mcbm_beam_2020_03")
+ TString setupName = "mcbm_beam_2020_03", Bool_t bUseMC = kTRUE)
{
// ========================================================================
@@ -48,7 +48,6 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test",
// ------------------------------------------------------------------------
// ----- Environment --------------------------------------------------
- bool bUseMC = true; // flag: true - MC information is used
int verbose = 6; // verbose level
TString myName = "mcbm_qa"; // this macro's name for screen output
TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
@@ -164,9 +163,16 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test",
// ------------------------------------------------------------------------
// ----- MCDataManager (legacy mode) -----------------------------------
- CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 1);
- mcManager->AddFile(traFile);
- run->AddTask(mcManager);
+ if (bUseMC) {
+ std::cout << "-I- " << myName << ": Adding MC manager and MC to reco matching tasks\n";
+
+ auto* mcManager = new CbmMCDataManager("MCDataManager", 1);
+ mcManager->AddFile(traFile);
+ run->AddTask(mcManager);
+
+ auto* matchRecoToMC = new CbmMatchRecoToMC();
+ run->AddTask(matchRecoToMC);
+ }
// ------------------------------------------------------------------------
// ----- MUCH QA ---------------------------------
@@ -179,6 +185,7 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test",
}
// ------------------------------------------------------------------------
+
// ----- CA tracking QA ---------------------------------------------------
// Tracking detector interface initialization
run->AddTask(new CbmTrackingDetectorInterfaceInit());
diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx
index d636d9bbf4c7ef7ec8e51e2f913659eb49a8ea3c..cfaea58318bbff75ccdb57c56147c67c784a54c6 100644
--- a/reco/L1/CbmCaMCModule.cxx
+++ b/reco/L1/CbmCaMCModule.cxx
@@ -470,7 +470,6 @@ void CbmCaMCModule::MatchRecoAndMCTracks()
trkRe.SetMaxPurity(maxPurity);
} // loop over hit map: end
assert(trkRe.GetNofMCTracks() < 2);
-
} // loop over reconstructed tracks: end
}
diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h
index 6097328d6e5cc3b8afea9cd1f18a3ef8d9b3466f..9423ce3d123d2f4c4b222835fa864eb30f542ab6 100644
--- a/reco/L1/L1LinkDef.h
+++ b/reco/L1/L1LinkDef.h
@@ -34,6 +34,7 @@
#pragma link C++ class CbmCaInputQaSts + ;
#pragma link C++ class CbmCaInputQaTrd + ;
#pragma link C++ class CbmCaInputQaTof + ;
+#pragma link C++ enum cbm::ca::ETrackType;
#pragma link C++ class cbm::ca::OutputQa + ;
#pragma link C++ class ca::tools::WindowFinder + ;
diff --git a/reco/L1/catools/CaToolsMCData.cxx b/reco/L1/catools/CaToolsMCData.cxx
index d4b8cdf36c0e203c8baca9ce98570b0fd2ccbc4c..33213add261b31944358880775aeff42d3671bb3 100644
--- a/reco/L1/catools/CaToolsMCData.cxx
+++ b/reco/L1/catools/CaToolsMCData.cxx
@@ -94,7 +94,6 @@ void MCData::Clear()
//
void MCData::InitTrackInfo(const L1Vector<CbmL1HitDebugInfo>& vHits)
{
- LOG(info) << "\033[1;32m!!!! FLAG a\033[0m";
for (auto& aTrk : fvTracks) {
// Assign hits to tracks
aTrk.ClearHitIndexes();
@@ -109,7 +108,6 @@ void MCData::InitTrackInfo(const L1Vector<CbmL1HitDebugInfo>& vHits)
aTrk.InitPointsInfo(fvPoints);
aTrk.InitHitsInfo(vHits);
}
- LOG(info) << "\033[1;32m!!!! FLAG b\033[0m";
}
// ---------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx
index 549a52f2951f3fbc18f4c60483eaedbcc9e8d7ac..af42719b282bd283fc7439ba8e4557612913fd5a 100644
--- a/reco/L1/qa/CbmCaOutputQa.cxx
+++ b/reco/L1/qa/CbmCaOutputQa.cxx
@@ -15,11 +15,41 @@
#include "L1InitManager.h"
using ca::tools::Debugger;
+using ::ca::tools::MCTrack;
using cbm::ca::OutputQa;
// ---------------------------------------------------------------------------------------------------------------------
//
-OutputQa::OutputQa(int verbose, bool isMCUsed) : CbmQaTask("CbmCaOutputQa", "caout", verbose, isMCUsed) {}
+OutputQa::OutputQa(int verbose, bool isMCUsed) : CbmQaTask("CbmCaOutputQa", "caout", verbose, isMCUsed)
+{
+ // Turn on default track classes
+ AddTrackType(ETrackType::kAll);
+ AddTrackType(ETrackType::kGhost);
+ AddTrackType(ETrackType::kPrim);
+ AddTrackType(ETrackType::kSec);
+ AddTrackType(ETrackType::kPrimPIP);
+ AddTrackType(ETrackType::kPrimPIM);
+ AddTrackType(ETrackType::kSecPIP);
+ AddTrackType(ETrackType::kSecPIM);
+ AddTrackType(ETrackType::kPrimMUP);
+ AddTrackType(ETrackType::kPrimMUM);
+ AddTrackType(ETrackType::kSecMUP);
+ AddTrackType(ETrackType::kSecMUM);
+
+ //AddTrackType(ETrackType::kAllE);
+ //AddTrackType(ETrackType::kAllMU);
+ //AddTrackType(ETrackType::kAllPI);
+ //AddTrackType(ETrackType::kAllK);
+ //AddTrackType(ETrackType::kAllPPBAR);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void OutputQa::AddTrackType(ETrackType type, bool flag)
+{
+ LOG_IF(fatal, type >= kEND || type < 0) << fName << ": attempt to add unsupported track type " << type;
+ fvbTrackTypeOn[type] = flag;
+}
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -36,68 +66,169 @@ void OutputQa::FillHistograms()
// ** Fill distributions for reconstructed tracks (including ghost ones) **
// ************************************************************************
- for (size_t iTrkReco = 0; iTrkReco < fvRecoTracks.size(); ++iTrkReco) {
- const auto& recoTrk = fvRecoTracks[iTrkReco];
+ // TODO: (experimental) more general approach for filling track types utilizing cuts on MC and reco tracks
+ if constexpr (kEXPTRACKFILL) {
+ // Fill histograms for each defined track type
+ for (int iType = 0; iType < ETrackType::kEND; ++iType) {
+ if (fvbTrackTypeOn[iType]) {
+ fvpTrackHistograms[iType]->FillRecoTracks();
+ if (fvpTrackHistograms[iType]->IsMCUsed()) { fvpTrackHistograms[iType]->FillMCTracks(); }
+ } // if track type ID filled
+ } // track type ID
+ }
+ else {
+ for (size_t iTrkReco = 0; iTrkReco < fvRecoTracks.size(); ++iTrkReco) {
+ const auto& recoTrk = fvRecoTracks[iTrkReco];
+
+ // Reject tracks, which do not contain hits
+ if (recoTrk.GetNofHits() < 1) { continue; }
- // Reject tracks, which do not contain hits
- if (recoTrk.GetNofHits() < 1) { continue; }
+ if (fvbTrackTypeOn[ETrackType::kAll]) { fvpTrackHistograms[ETrackType::kAll]->FillRecoTrack(iTrkReco); }
- fvpTrackHistograms[ETrackType::kAll]->FillRecoTrack(iTrkReco);
+ if (IsMCUsed()) {
+ if (fvbTrackTypeOn[ETrackType::kGhost] && recoTrk.IsGhost()) {
+ fvpTrackHistograms[ETrackType::kGhost]->FillRecoTrack(iTrkReco);
+ }
+ int iTrkMC = recoTrk.GetMatchedMCTrackIndex();
+ if (iTrkMC > -1) {
+ const auto& mcTrk = fMCData.GetTrack(iTrkMC);
+ int pdg = mcTrk.GetPdgCode();
+ bool isPrimary = mcTrk.IsPrimary();
+
+ // Cut tracks, which did not leave hits in tracker
+ if (mcTrk.GetNofHits() == 0) { continue; }
+
+ if (isPrimary) {
+ if (fvbTrackTypeOn[ETrackType::kPrim]) { fvpTrackHistograms[ETrackType::kPrim]->FillRecoTrack(iTrkReco); }
+
+ if (fvbTrackTypeOn[ETrackType::kPrimPI] && std::abs(pdg) == 211) {
+ fvpTrackHistograms[ETrackType::kPrimPI]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kPrimPIP] && pdg == +211) {
+ fvpTrackHistograms[ETrackType::kPrimPIP]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kPrimPIM] && pdg == -211) {
+ fvpTrackHistograms[ETrackType::kPrimPIM]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kPrimMU] && std::abs(pdg) == 13) {
+ fvpTrackHistograms[ETrackType::kPrimMU]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kPrimMUP] && pdg == +13) {
+ fvpTrackHistograms[ETrackType::kPrimMUP]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kPrimMUM] && pdg == -13) {
+ fvpTrackHistograms[ETrackType::kPrimMUM]->FillRecoTrack(iTrkReco);
+ }
+ }
+ else {
+ if (fvbTrackTypeOn[ETrackType::kSec]) { fvpTrackHistograms[ETrackType::kSec]->FillRecoTrack(iTrkReco); }
+
+ if (fvbTrackTypeOn[ETrackType::kSecPI] && std::abs(pdg) == 211) {
+ fvpTrackHistograms[ETrackType::kSecPI]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kSecPIP] && pdg == +211) {
+ fvpTrackHistograms[ETrackType::kSecPIP]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kSecPIM] && pdg == -211) {
+ fvpTrackHistograms[ETrackType::kSecPIM]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kSecMU] && std::abs(pdg) == 13) {
+ fvpTrackHistograms[ETrackType::kSecMU]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kSecMUP] && pdg == +13) {
+ fvpTrackHistograms[ETrackType::kSecMUP]->FillRecoTrack(iTrkReco);
+ }
+
+ if (fvbTrackTypeOn[ETrackType::kSecMUM] && pdg == -13) {
+ fvpTrackHistograms[ETrackType::kSecMUM]->FillRecoTrack(iTrkReco);
+ }
+ }
+ }
+ }
+ } // loop over recoTrk: end
+
+ // *************************************
+ // ** Fill distributions of MC-tracks **
+ // *************************************
if (IsMCUsed()) {
- if (recoTrk.IsGhost()) { fvpTrackHistograms[ETrackType::kGhost]->FillRecoTrack(iTrkReco); }
- int iTrkMC = recoTrk.GetMatchedMCTrackIndex();
- if (iTrkMC > -1) {
+
+
+ for (int iTrkMC = 0; iTrkMC < fMCData.GetNofTracks(); ++iTrkMC) {
const auto& mcTrk = fMCData.GetTrack(iTrkMC);
- int pdg = mcTrk.GetPdgCode();
- bool isPrimary = mcTrk.IsPrimary();
+
+ // Cut tracks, which did not leave hits in tracker
+ if (mcTrk.GetNofHits() == 0) { continue; }
+
+ // Fill different track categories
+ if (fvbTrackTypeOn[ETrackType::kAll]) { fvpTrackHistograms[ETrackType::kAll]->FillMCTrack(iTrkMC); }
+
+ int pdg = mcTrk.GetPdgCode();
+ bool isPrimary = mcTrk.IsPrimary();
if (isPrimary) {
- fvpTrackHistograms[ETrackType::kPrim]->FillRecoTrack(iTrkReco);
+ if (fvbTrackTypeOn[ETrackType::kPrim]) { fvpTrackHistograms[ETrackType::kPrim]->FillMCTrack(iTrkMC); }
- if (pdg == +211) { fvpTrackHistograms[ETrackType::kPrimPIP]->FillRecoTrack(iTrkReco); }
- if (pdg == -211) { fvpTrackHistograms[ETrackType::kPrimPIM]->FillRecoTrack(iTrkReco); }
- }
- else {
- fvpTrackHistograms[ETrackType::kSec]->FillRecoTrack(iTrkReco);
+ if (fvbTrackTypeOn[ETrackType::kPrimPI] && std::abs(pdg) == 211) {
+ fvpTrackHistograms[ETrackType::kPrimPI]->FillMCTrack(iTrkMC);
+ }
- if (pdg == +211) { fvpTrackHistograms[ETrackType::kSecPIP]->FillRecoTrack(iTrkReco); }
- if (pdg == -211) { fvpTrackHistograms[ETrackType::kSecPIM]->FillRecoTrack(iTrkReco); }
- }
- }
- }
+ if (fvbTrackTypeOn[ETrackType::kPrimPIP] && pdg == +211) {
+ fvpTrackHistograms[ETrackType::kPrimPIP]->FillMCTrack(iTrkMC);
+ }
+ if (fvbTrackTypeOn[ETrackType::kPrimPIM] && pdg == -211) {
+ fvpTrackHistograms[ETrackType::kPrimPIM]->FillMCTrack(iTrkMC);
+ }
- } // loop over recoTrk: end
+ if (fvbTrackTypeOn[ETrackType::kPrimMU] && std::abs(pdg) == 13) {
+ fvpTrackHistograms[ETrackType::kPrimMU]->FillMCTrack(iTrkMC);
+ }
+ if (fvbTrackTypeOn[ETrackType::kPrimMUP] && pdg == +13) {
+ fvpTrackHistograms[ETrackType::kPrimMUP]->FillMCTrack(iTrkMC);
+ }
- // *************************************
- // ** Fill distributions of MC-tracks **
- // *************************************
- if (IsMCUsed()) {
- for (int iTrkMC = 0; iTrkMC < fMCData.GetNofTracks(); ++iTrkMC) {
- const auto& mcTrk = fMCData.GetTrack(iTrkMC);
+ if (fvbTrackTypeOn[ETrackType::kPrimMUM] && pdg == -13) {
+ fvpTrackHistograms[ETrackType::kPrimMUM]->FillMCTrack(iTrkMC);
+ }
+ }
+ else {
+ if (fvbTrackTypeOn[ETrackType::kSec]) { fvpTrackHistograms[ETrackType::kSec]->FillMCTrack(iTrkMC); }
- // Cut tracks, which did not leave hits in tracker
- if (mcTrk.GetNofHits() == 0) { continue; }
+ if (fvbTrackTypeOn[ETrackType::kSecPIP] && std::abs(pdg) == 211) {
+ fvpTrackHistograms[ETrackType::kSecPIP]->FillMCTrack(iTrkMC);
+ }
- // Fill different track categories
- fvpTrackHistograms[ETrackType::kAll]->FillMCTrack(iTrkMC);
+ if (fvbTrackTypeOn[ETrackType::kSecPIP] && pdg == +211) {
+ fvpTrackHistograms[ETrackType::kPrimPIP]->FillMCTrack(iTrkMC);
+ }
- int pdg = mcTrk.GetPdgCode();
- bool isPrimary = mcTrk.IsPrimary();
+ if (fvbTrackTypeOn[ETrackType::kSecPIM] && pdg == -211) {
+ fvpTrackHistograms[ETrackType::kPrimPIM]->FillMCTrack(iTrkMC);
+ }
- if (isPrimary) {
- fvpTrackHistograms[ETrackType::kPrim]->FillMCTrack(iTrkMC);
+ if (fvbTrackTypeOn[ETrackType::kSecMU] && std::abs(pdg) == 13) {
+ fvpTrackHistograms[ETrackType::kSecMU]->FillMCTrack(iTrkMC);
+ }
- if (pdg == +211) { fvpTrackHistograms[ETrackType::kPrimPIP]->FillMCTrack(iTrkMC); }
- if (pdg == -211) { fvpTrackHistograms[ETrackType::kPrimPIM]->FillMCTrack(iTrkMC); }
- }
- else {
- fvpTrackHistograms[ETrackType::kSec]->FillMCTrack(iTrkMC);
+ if (fvbTrackTypeOn[ETrackType::kSecMUP] && pdg == +13) {
+ fvpTrackHistograms[ETrackType::kSecMUP]->FillMCTrack(iTrkMC);
+ }
- if (pdg == +211) { fvpTrackHistograms[ETrackType::kSecPIP]->FillMCTrack(iTrkMC); }
- if (pdg == -211) { fvpTrackHistograms[ETrackType::kSecPIM]->FillMCTrack(iTrkMC); }
+ if (fvbTrackTypeOn[ETrackType::kSecMUM] && pdg == -13) {
+ fvpTrackHistograms[ETrackType::kSecMUM]->FillMCTrack(iTrkMC);
+ }
+ }
}
}
}
@@ -183,24 +314,229 @@ InitStatus OutputQa::InitDataBranches()
InitStatus OutputQa::InitHistograms()
{
auto RegisterTrackQa = [&](const char* typeName, ETrackType type, bool bSuppressMC = false) {
+ if (!fvbTrackTypeOn[type]) { return; }
bool bUseMC = IsMCUsed() && !bSuppressMC;
fvpTrackHistograms[type] = std::make_unique<TrackTypeQa>(typeName, fsPrefix.Data(), bUseMC, fpFolderRoot);
fvpTrackHistograms[type]->RegisterParameters(fpParameters);
fvpTrackHistograms[type]->RegisterRecoHits(fvHits);
fvpTrackHistograms[type]->RegisterRecoTracks(fvRecoTracks);
fvpTrackHistograms[type]->RegisterMCData(fMCData);
+
+ if constexpr (kEXPTRACKFILL) {
+ // Define track cuts
+ switch (type) {
+ //
+ // General
+ case ETrackType::kAll: break;
+ case ETrackType::kGhost:
+ fvpTrackHistograms[type]->SetRecoTrackCut([](const CbmL1Track& t) -> bool { return t.IsGhost(); });
+ break;
+ case kPrim:
+ fvpTrackHistograms[type]->SetMCTrackCut([](const MCTrack& t) -> bool { return t.IsPrimary(); });
+ break;
+ case kSec:
+ fvpTrackHistograms[type]->SetMCTrackCut([](const MCTrack& t) -> bool { return !t.IsPrimary(); });
+ break;
+ //
+ // Electrons
+ case kAllE:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return std::abs(t.GetPdgCode()) == 11; });
+ break;
+ case kPrimE:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && std::abs(t.GetPdgCode()) == 11; });
+ break;
+ case kPrimEP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == -11; });
+ break;
+ case kPrimEM:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == 11; });
+ break;
+ case kSecE:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && std::abs(t.GetPdgCode()) == 11; });
+ break;
+ case kSecEP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == -11; });
+ break;
+ case kSecEM:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == 11; });
+ break;
+ //
+ // Muons
+ case kAllMU:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return std::abs(t.GetPdgCode()) == 13; });
+ break;
+ case kPrimMU:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && std::abs(t.GetPdgCode()) == 13; });
+ break;
+ case kPrimMUP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == -13; });
+ break;
+ case kPrimMUM:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == 13; });
+ break;
+ case kSecMU:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && std::abs(t.GetPdgCode()) == 13; });
+ break;
+ case kSecMUP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == -13; });
+ break;
+ case kSecMUM:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == 13; });
+ break;
+ //
+ // Pions
+ case kAllPI:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return std::abs(t.GetPdgCode()) == 211; });
+ break;
+ case kPrimPI:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && std::abs(t.GetPdgCode()) == 211; });
+ break;
+ case kPrimPIP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == 211; });
+ break;
+ case kPrimPIM:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == -211; });
+ break;
+ case kSecPI:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && std::abs(t.GetPdgCode()) == 211; });
+ break;
+ case kSecPIP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == 211; });
+ break;
+ case kSecPIM:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == -211; });
+ break;
+ //
+ // Kaons
+ case kAllK:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return std::abs(t.GetPdgCode()) == 211; });
+ break;
+ case kPrimK:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && std::abs(t.GetPdgCode()) == 321; });
+ break;
+ case kPrimKP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == 321; });
+ break;
+ case kPrimKM:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == -321; });
+ break;
+ case kSecK:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && std::abs(t.GetPdgCode()) == 321; });
+ break;
+ case kSecKP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == 321; });
+ break;
+ case kSecKM:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == -321; });
+ break;
+ //
+ // Protons
+ case kAllPPBAR:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return std::abs(t.GetPdgCode()) == 2212; });
+ break;
+ case kPrimPPBAR:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && std::abs(t.GetPdgCode()) == 2212; });
+ break;
+ case kPrimP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == 2212; });
+ break;
+ case kPrimPBAR:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return t.IsPrimary() && t.GetPdgCode() == -2212; });
+ break;
+ case kSecPPBAR:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && std::abs(t.GetPdgCode()) == 2212; });
+ break;
+ case kSecP:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == 2212; });
+ break;
+ case kSecPBAR:
+ fvpTrackHistograms[type]->SetMCTrackCut(
+ [](const MCTrack& t) -> bool { return !t.IsPrimary() && t.GetPdgCode() == -2212; });
+ break;
+ }
+ }
+
fvpTrackHistograms[type]->Init();
};
- RegisterTrackQa("ghost", ETrackType::kGhost, /*suppress MC*/ true);
+ for (int i = 0; i < ETrackType::kEND; ++i) {
+ LOG(info) << i << ' ' << fvpTrackHistograms[i].get() << ' ' << fvbTrackTypeOn[i];
+ }
+
RegisterTrackQa("all", ETrackType::kAll);
if (IsMCUsed()) {
+ RegisterTrackQa("ghost", ETrackType::kGhost, /*suppress MC*/ true);
RegisterTrackQa("prim", ETrackType::kPrim);
RegisterTrackQa("sec", ETrackType::kSec);
+ RegisterTrackQa("all_pi", ETrackType::kAllPI);
+ RegisterTrackQa("prim_pi", ETrackType::kPrimPI);
RegisterTrackQa("prim_pip", ETrackType::kPrimPIP);
RegisterTrackQa("prim_pim", ETrackType::kPrimPIM);
+ RegisterTrackQa("sec_pi", ETrackType::kSecPI);
RegisterTrackQa("sec_pip", ETrackType::kSecPIP);
RegisterTrackQa("sec_pim", ETrackType::kSecPIM);
+ RegisterTrackQa("all_e", ETrackType::kAllE);
+ RegisterTrackQa("prim_e", ETrackType::kPrimE);
+ RegisterTrackQa("prim_ep", ETrackType::kPrimEP);
+ RegisterTrackQa("prim_em", ETrackType::kPrimEM);
+ RegisterTrackQa("sec_e", ETrackType::kSecE);
+ RegisterTrackQa("sec_ep", ETrackType::kSecEP);
+ RegisterTrackQa("sec_em", ETrackType::kSecEM);
+ RegisterTrackQa("all_mu", ETrackType::kAllMU);
+ RegisterTrackQa("prim_mu", ETrackType::kPrimMU);
+ RegisterTrackQa("prim_mup", ETrackType::kPrimMUP);
+ RegisterTrackQa("prim_mum", ETrackType::kPrimMUM);
+ RegisterTrackQa("sec_mu", ETrackType::kSecMU);
+ RegisterTrackQa("sec_mup", ETrackType::kSecMUP);
+ RegisterTrackQa("sec_mum", ETrackType::kSecMUM);
+ RegisterTrackQa("all_k", ETrackType::kAllK);
+ RegisterTrackQa("prim_k", ETrackType::kPrimK);
+ RegisterTrackQa("prim_kp", ETrackType::kPrimKP);
+ RegisterTrackQa("prim_km", ETrackType::kPrimKM);
+ RegisterTrackQa("sec_k", ETrackType::kSecK);
+ RegisterTrackQa("sec_kp", ETrackType::kSecKP);
+ RegisterTrackQa("sec_km", ETrackType::kSecKM);
+ RegisterTrackQa("all_ppbar", ETrackType::kAllPPBAR);
+ RegisterTrackQa("prim_ppbar", ETrackType::kPrimPPBAR);
+ RegisterTrackQa("prim_p", ETrackType::kPrimP);
+ RegisterTrackQa("prim_pbar", ETrackType::kPrimPBAR);
+ RegisterTrackQa("sec_ppbar", ETrackType::kSecPPBAR);
+ RegisterTrackQa("sec_p", ETrackType::kSecP);
+ RegisterTrackQa("sec_pbar", ETrackType::kSecPBAR);
}
return kSUCCESS;
@@ -230,6 +566,7 @@ InitStatus OutputQa::InitTimeSlice()
LOG_IF(info, fVerbose > 1) << fName << ": Data sample consists of " << nHits << " hits, " << nRecoTracks
<< " reco tracks, " << nMCTracks << " MC tracks, " << nMCPoints << " MC points";
+
return kSUCCESS;
}
diff --git a/reco/L1/qa/CbmCaOutputQa.h b/reco/L1/qa/CbmCaOutputQa.h
index cd7ebc203692728682545da907c23fd2b3843adf..2cac1f21bc4c6d637c0dc8a918714d852c74eccb 100644
--- a/reco/L1/qa/CbmCaOutputQa.h
+++ b/reco/L1/qa/CbmCaOutputQa.h
@@ -25,62 +25,92 @@
namespace cbm::ca
{
+ /// @brief Enumeration of track category
+ enum ETrackType
+ {
+ kAll, ///< all tracks
+ kGhost, ///< ghost tracks (no MC is used)
+ kPrim, ///< primary
+ kSec, ///< secondary
+ kAllE, ///< all e-/e+
+ kPrimE, ///< primary e-/e+
+ kPrimEP, ///< primary e+
+ kPrimEM, ///< primary e-
+ kSecE, ///< secondary e-/e+
+ kSecEP, ///< secondary e+
+ kSecEM, ///< secondary e-
+ kAllMU, ///< all mu+/mu-
+ kPrimMU, ///< primary mu+/mu-
+ kPrimMUP, ///< primary mu+
+ kPrimMUM, ///< primary mu-
+ kSecMU, ///< secondary mu+/mu-
+ kSecMUP, ///< secondary mu+
+ kSecMUM, ///< secondary mu-
+ kAllPI, ///< all pi+/pi-
+ kPrimPI, ///< primary pi+/pi-
+ kPrimPIP, ///< primary pi+
+ kPrimPIM, ///< primary pi-
+ kSecPI, ///< secondary pi+/pi-
+ kSecPIP, ///< secondary pi+
+ kSecPIM, ///< secondary pi-
+ kAllK, ///< all K+/K-
+ kPrimK, ///< primary K+/K-
+ kPrimKP, ///< primary K+
+ kPrimKM, ///< primary K-
+ kSecK, ///< secondary K+/K-
+ kSecKP, ///< secondary K+
+ kSecKM, ///< secondary K-
+ kAllPPBAR, ///< all p/p-bar
+ kPrimPPBAR, ///< primary p/p-bar
+ kPrimP, ///< primary p
+ kPrimPBAR, ///< primary p-bar
+ kSecPPBAR, ///< secondary p/p-bar
+ kSecP, ///< secondary p
+ kSecPBAR, ///< secondary p-bar
+ // kPrimD,
+ // kPrimDBAR,
+ // kSecD,
+ // kSecDBAR,
+ // kPrimT,
+ // kPrimTBAR,
+ // kSecT,
+ // kSecTBAR,
+ // kPrim3HE,
+ // kPrim3HEBAR,
+ // kSec3HE,
+ // kSec3HEBAR,
+ // kPrim4HE,
+ // kPrim4HEBAR,
+ // kSec4HE,
+ // kSec4HEBAR,
+ kEND
+ };
+
/// @brief QA-task for CA tracking output results
///
class OutputQa : public CbmQaTask {
- public:
- /// @brief Enumeration of track category
- enum ETrackType
- {
- kAll, ///< all tracks
- kGhost, ///< ghost tracks (no MC is used)
- kPrim, ///< primary tracks
- kSec, ///< secondary tracks
- // kPrimEP, ///< primary electron tracks
- // kPrimEM, ///< primary positron tracks
- // kSecEP, ///< secondary electron tracks
- // kSecEM, ///< secondary positron tracks
- // kPrimMUP, ///< primary .....
- // kPrimMUM,
- // kSecMUP,
- // kSecMUM,
- kPrimPIP,
- kPrimPIM,
- kSecPIP,
- kSecPIM,
- // kPrimKP,
- // kPrimKM,
- // kSecKP,
- // kSecKM,
- // kPrimP,
- // kPrimPBAR,
- // kSecP,
- // kSecPBAR,
- // kPrimD,
- // kPrimDBAR,
- // kSecD,
- // kSecDBAR,
- // kPrimT,
- // kPrimTBAR,
- // kSecT,
- // kSecTBAR,
- // kPrim3HE,
- // kPrim3HEBAR,
- // kSec3HE,
- // kSec3HEBAR,
- // kPrim4HE,
- // kPrim4HEBAR,
- // kSec4HE,
- // kSec4HEBAR,
- kEND
- };
-
+ // WIP: Temporary flag to switch between two different approaches of filling track type histograms
+ //
+ // 1) Experimental approach (kEXPTRACKFILL = true) utilizes unified filling based on defined reco and MC cuts.
+ // 2) Feature is to be studied more precisely (descrepancy in primary/secondary track with a standard approach)
+ // 3) Experimental approach runs in ~10% slower, then the standard
+ static constexpr bool kEXPTRACKFILL = false;
+ public:
/// @brief Constructor from parameters
/// @param verbose Verbosity level
/// @param isMCUsed Flag, if MC information is available for this task
OutputQa(int verbose, bool isMCUsed);
+ /// @brief Adds track type
+ /// @param type Track type
+ /// @param flag Flag: true/false
+ ///
+ /// Adds a track type for building distributions. By default, only all, primary, secondary and ghost track
+ /// distributions are booked.
+ void AddTrackType(ETrackType type, bool flag = true);
+
+
/// @brief Enables debugger
/// @param filename Name of output ROOT file
///
@@ -149,6 +179,11 @@ namespace cbm::ca
/// @brief De-initializes histograms
void DeInit() {}
+ /// @brief Fills histograms for a given track types
+ /// @param type Track type
+ /// @param
+ void FillTrackTypeHistograms();
+
private:
// Flags for detector subsystems being used
bool fbUseMvd = false; ///< is MVD used
@@ -179,8 +214,7 @@ namespace cbm::ca
/// Histograms of different track types
std::array<std::unique_ptr<TrackTypeQa>, ETrackType::kEND> fvpTrackHistograms;
-
- std::array<std::tuple<Style_t, Marker_t, Color_t>, ETrackType::kEND> fvpTrackHistogramStyles;
+ std::array<bool, ETrackType::kEND> fvbTrackTypeOn = {0}; ///< Track type is on
};
} // namespace cbm::ca
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.cxx b/reco/L1/qa/CbmCaTrackTypeQa.cxx
index eed783fb873e9fbf047ce6912d99920d5f5f9c44..77a329d086ea4634bb07e57baf92a1d53199f387 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.cxx
+++ b/reco/L1/qa/CbmCaTrackTypeQa.cxx
@@ -179,6 +179,33 @@ void TrackTypeQa::FillRecoTrack(int iTrkReco, double weight)
}
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackTypeQa::FillRecoTracks()
+{
+ for (int iTrkReco = 0; iTrkReco < (int) fpvRecoTracks->size(); ++iTrkReco) {
+ const auto& recoTrk = (*fpvRecoTracks)[iTrkReco];
+
+ // Reject tracks, which do not contain hits
+ if (recoTrk.GetNofHits() < 1) { continue; }
+
+ // Reject tracks, which do not pass the applied cuts
+ if (!fRecoTrackCut(recoTrk)) { continue; }
+
+ // Reject tracks, which have MC tracks not passing the applied cuts
+ if (fbUseMC) {
+ int iTrkMC = recoTrk.GetMatchedMCTrackIndex();
+ if (iTrkMC > -1) {
+ const auto& mcTrk = fpMCData->GetTrack(iTrkMC);
+ if (!fMCTrackCut(mcTrk)) { continue; }
+ }
+ }
+
+ // Fill histograms
+ this->FillRecoTrack(iTrkReco);
+ }
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
void TrackTypeQa::FillMCTrack(int iTrkMC, double weight)
@@ -210,3 +237,25 @@ void TrackTypeQa::FillMCTrack(int iTrkMC, double weight)
fph_eff_pMC_yMC->Fill(mcTrack.GetRapidity(), mcTrack.GetP(), bReco);
fph_eff_thetaMC_phiMC->Fill(mcTrack.GetPhi(), mcTrack.GetPhi(), bReco);
}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackTypeQa::FillMCTracks()
+{
+ assert(fbUseMC);
+ for (int iTrkMC = 0; iTrkMC < fpMCData->GetNofTracks(); ++iTrkMC) {
+ const auto& mcTrk = fpMCData->GetTrack(iTrkMC);
+
+ // Cut tracks, which did not leave hits in tracker
+ if (mcTrk.GetNofHits() == 0) { continue; }
+
+ // Cut tracks, which do not pass the applied cut
+ if (!fMCTrackCut(mcTrk)) { continue; }
+
+ // TODO: Investigate, how to apply reconstructed track cuts here
+ // NOTE: Reconstructed track cut analogs must be applied to mc tracks
+
+ // Fill histograms
+ this->FillMCTrack(iTrkMC);
+ }
+}
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.h b/reco/L1/qa/CbmCaTrackTypeQa.h
index 13acc148490ab9147d28031f757609c450ebd12f..73d284fcd74c8d0def7c4fff747495589f5f175c 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.h
+++ b/reco/L1/qa/CbmCaTrackTypeQa.h
@@ -23,7 +23,8 @@
namespace ca::tools
{
class MCData;
-}
+ class MCTrack;
+} // namespace ca::tools
class CbmL1Track;
class CbmL1HitDebugInfo;
@@ -69,18 +70,31 @@ namespace cbm::ca
/// @brief Initializes histograms
void Init();
+ /// @brief Flag on MC usage by the track type
+ /// @return true Histograms relied on MC info are filled
+ /// @return false Histograms relied on MC info are not filled
+ bool IsMCUsed() const { return fbUseMC; }
+
/// @brief Fills histograms with various track information
/// @note To be called in the loop over reconstructed tracks full sample
/// @param iTrkReco Index of reconstructed track
/// @param weight Weight
void FillRecoTrack(int iTrkReco, double weight = 1);
+ /// @brief Fills histograms with recon tracks
+ /// @note To be called once per event/TS
+ void FillRecoTracks();
+
/// @brief Fills histograms with mc track information
/// @note To be called in the loop over MC tracks full sample
/// @param iTrkMC Index of MC track
/// @param weight Weight
void FillMCTrack(int iTrkMC, double weight = 1);
+ /// @brief Fills histograms with recon tracks
+ /// @note To be called once per event/TS
+ void FillMCTracks();
+
/// @brief Registers the reconstructed tracks container
/// @param vTracks Reference to reconstructed tracks container
void RegisterRecoTracks(L1Vector<CbmL1Track>& vTracks) { fpvRecoTracks = &vTracks; }
@@ -101,6 +115,18 @@ namespace cbm::ca
/// @param title Title of track type
void SetTitle(const char* title) { fsTitle = title; }
+ /// @brief Sets selection cuts on MC tracks
+ /// @param cut Functor <bool(const MCTrack&)> defining cut on MC track
+ ///
+ /// When function returns false,
+ void SetMCTrackCut(std::function<bool(const ::ca::tools::MCTrack&)> cut) { fMCTrackCut = cut; }
+
+ /// @brief Sets selection cuts on reconstructed tracks
+ /// @param cut Functor <bool(const CbmL1Track&)> defining cut on reconstructed track
+ ///
+ /// When function returns false,
+ void SetRecoTrackCut(std::function<bool(const CbmL1Track&)> cut) { fRecoTrackCut = cut; }
+
// ************************
// ** List of histograms **
// ************************
@@ -177,6 +203,15 @@ namespace cbm::ca
::ca::tools::MCData* fpMCData = nullptr; ///< Pointer to MC data object
std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to parameters object
+ // ** Cuts on tracks for a given track class **
+
+ /// Cut function on MC tracks
+ std::function<bool(const ::ca::tools::MCTrack&)> fMCTrackCut = [](const ::ca::tools::MCTrack&) { return true; };
+
+ /// Cut function on reconstructed tracks
+ std::function<bool(const CbmL1Track&)> fRecoTrackCut = [](const CbmL1Track&) { return true; };
+
+
// **************************
// ** Histogram properties **
// **************************