diff --git a/core/qa/CbmQaTask.h b/core/qa/CbmQaTask.h index 23d4d885f98680d798980dde9c2520108a5824c9..732d33e5224f85239eaef4ee49e51704c3170e61 100644 --- a/core/qa/CbmQaTask.h +++ b/core/qa/CbmQaTask.h @@ -101,6 +101,11 @@ class CbmQaTask : public FairTask, public CbmQaIO { /// \return false Some of the objects are inconsistent bool CompareQaObjects(); + /// \brief Disables event-by-event execution + /// + /// By default, if the branch + void DisableEventMode(); + /// \brief FairTask: Defines action of the task in the event/TS void Exec(Option_t* /*option*/) override; @@ -155,7 +160,6 @@ class CbmQaTask : public FairTask, public CbmQaIO { /// \brief Sets name of the setup void SetSetupName(const char* setup) { fsSetupName = setup; } - protected: /// \brief De-initialize the task virtual void DeInit() {} diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C index 20c0a2ac75d03c036cbd8eb15f9e63b49adaf39d..2eed8d369af3f7c034dc314546e1ba3281e36916 100644 --- a/macro/run/run_qa.C +++ b/macro/run/run_qa.C @@ -38,11 +38,11 @@ #endif void run_qa(TString data = "test", TString setup = "sis100_electron", Int_t nEvents = -1, TString dataTra2 = "", - TString dataTra3 = "", TString configName = ""); + TString dataTra3 = "", TString configName = "", TString sEvBuildRaw = ""); void run_qa(TString dataTra, TString dataRaw, TString dataReco, TString dataPar, TString dataSink, TString setup = "sis100_electron", Int_t nEvents = -1, TString dataTra2 = "", TString dataTra3 = "", - TString configName = "") + TString configName = "", TString sEvBuildRaw = "") { gROOT->SetBatch(kTRUE); @@ -89,7 +89,7 @@ void run_qa(TString dataTra, TString dataRaw, TString dataReco, TString dataPar, // ------------------------------------------------------------------------ // ----- Some global switches ----------------------------------------- - //Bool_t eventBased = !sEvBuildRaw.IsNull(); + Bool_t bEventBasedReco = !sEvBuildRaw.IsNull(); bool bUseMvd = geo->IsActive(ECbmModuleId::kMvd); bool bUseSts = geo->IsActive(ECbmModuleId::kSts); bool bUseRich = geo->IsActive(ECbmModuleId::kRich); @@ -259,6 +259,7 @@ void run_qa(TString dataTra, TString dataRaw, TString dataReco, TString dataPar, caParFile.ReplaceAll(".root", ".L1Parameters.dat"); auto* pCaOutputQa = new cbm::ca::OutputQa(verbose, bUseMC); + pCaOutputQa->SetProcessFullTs(!bEventBasedReco); pCaOutputQa->SetStsTrackingMode(); pCaOutputQa->ReadParameters(caParFile.Data()); if (configName.Length() != 0) { @@ -335,7 +336,8 @@ void run_qa(TString dataTra, TString dataRaw, TString dataReco, TString dataPar, // ------------------------------------------------------------------------ } -void run_qa(TString data, TString setup, Int_t nEvents, TString dataTra2, TString dataTra3, TString configName) +void run_qa(TString data, TString setup, Int_t nEvents, TString dataTra2, TString dataTra3, TString configName, + TString sEvBuildRaw) { - run_qa(data, data, data, data, data, setup, nEvents, dataTra2, dataTra3, configName); + run_qa(data, data, data, data, data, setup, nEvents, dataTra2, dataTra3, configName, sEvBuildRaw); } diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx index 12007e5b166efdc4b534be4f6875e96a346dec3d..a84a66a92b95a09792334063d57dd78c4190a2ff 100644 --- a/reco/L1/qa/CbmCaOutputQa.cxx +++ b/reco/L1/qa/CbmCaOutputQa.cxx @@ -1,4 +1,4 @@ -/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt +/* Copyright (C) 2023-2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt SPDX-License-Identifier: GPL-3.0-only Authors: Sergei Zharko [committer] */ @@ -80,909 +80,770 @@ OutputQa::OutputQa(int verbose, bool isMCUsed) : CbmQaTask("CbmCaOutputQa", verb // --------------------------------------------------------------------------------------------------------------------- // -void OutputQa::FillHistograms() +void OutputQa::Check() { - fMonitor.IncrementCounter(EMonitorKey::kEvent); - - // ************************************************************************ - // ** Fill distributions for reconstructed tracks (including ghost ones) ** - // ************************************************************************ + // Create summary table + if (IsMCUsed()) { + fpMCModule->Finish(); - for (size_t iTrkReco = 0; iTrkReco < fvRecoTracks.size(); ++iTrkReco) { - const auto& recoTrk = fvRecoTracks[iTrkReco]; + // TODO: Add cuts on entries from fmSummaryTableEntries + std::vector<ETrackType> vTypesToPlot; + int nRows = std::count_if(fmSummaryTableEntries.begin(), fmSummaryTableEntries.end(), + [&](const auto& f) { return fvbTrackTypeOn[f] && fvpTrackHistograms[f]->IsMCUsed(); }) + + 1; - // Reject tracks, which do not contain hits - if (recoTrk.GetNofHits() < 1) { - continue; + CbmQaTable* aTable = MakeQaObject<CbmQaTable>("summary_table", "Tracking summary table", nRows + 1, 9); + int iRow = 0; + std::vector<std::string> colNames = {"Efficiency", "Killed", "Length", "Fakes", "Clones", + "Reco/Evt", "MC/Evt", "Nst(hit)", "Nst(point)"}; + aTable->SetNamesOfCols(colNames); + aTable->SetColWidth(14); + double nEvents = static_cast<double>(GetEventNumber()); + LOG(info) << "Number of events: " << GetEventNumber(); + for (auto trType : fmSummaryTableEntries) { + if (!fvbTrackTypeOn[trType] || !fvpTrackHistograms[trType]->IsMCUsed()) { + continue; + } + aTable->SetRowName(iRow, fvpTrackHistograms[trType]->GetTitle()); + aTable->SetCell(iRow, 0, fvpTrackHistograms[trType]->GetIntegratedEff()); + aTable->SetCell(iRow, 1, fvpTrackHistograms[trType]->GetKilledRate()); + aTable->SetCell(iRow, 2, fvpTrackHistograms[trType]->GetAverageRecoLength()); + aTable->SetCell(iRow, 3, fvpTrackHistograms[trType]->GetAverageFakeLength()); + aTable->SetCell(iRow, 4, fvpTrackHistograms[trType]->GetClonesRate()); + aTable->SetCell(iRow, 5, fvpTrackHistograms[trType]->GetNofRecoTracksMatched() / nEvents); + aTable->SetCell(iRow, 6, fvpTrackHistograms[trType]->GetNofMCTracks() / nEvents); + aTable->SetCell(iRow, 7, fvpTrackHistograms[trType]->GetAverageNofStationsWithHit()); + aTable->SetCell(iRow, 8, fvpTrackHistograms[trType]->GetAverageNofStationsWithPoint()); + ++iRow; + } + double nGhosts = 0.; + if (fvpTrackHistograms[ETrackType::kGhost] && fvpTrackHistograms[ETrackType::kAll]) { + nGhosts = fvpTrackHistograms[ETrackType::kGhost]->fph_reco_p->GetEntries(); + aTable->SetRowName(iRow, "N ghosts"); + aTable->SetCell(iRow, 0, nGhosts); + aTable->SetRowName(iRow + 1, "Ghost rate"); + aTable->SetCell(iRow + 1, 0, nGhosts / fMonitor.GetCounterValue(EMonitorKey::kTrack)); } + LOG(info) << '\n' << aTable->ToString(3); + } - FillRecoTrack(ETrackType::kAll, iTrkReco); + LOG(info) << '\n' << fMonitor.ToString(); +} - if (IsMCUsed()) { - // NOTE: The ghost status of track is now defined by its purity, thus it can still contain MC information - if (recoTrk.IsGhost()) { - FillRecoTrack(ETrackType::kGhost, iTrkReco); - } +// --------------------------------------------------------------------------------------------------------------------- +// +void OutputQa::DrawEvent() +{ + constexpr double kZmin = 0.; + constexpr double kZmax = 300.; + constexpr double kXmin = -100.; + constexpr double kXmax = +100.; + constexpr double kYmin = -100.; + constexpr double kYmax = +100.; + constexpr std::array<Color_t, 11> kColorMC = {205, 209, 213, 217, 225, 208, 213, 216, 219, 224, 227}; + constexpr std::array<Color_t, 3> kColorGhost = {201, 202, 203}; + constexpr int kCanvX = 1920; + constexpr int kCanvY = 1080; + constexpr double kLMargin = 0.05; // Left margin + constexpr double kVEMargin = 0.15; // Vertical margin (top + bottom) + constexpr double kRMarginDispl = 0.4; + constexpr double kVIMargin = 0.0001; + constexpr Marker_t kMarkerPointWHit = 25; + constexpr Marker_t kMarkerPointWOHit = 5; + constexpr Marker_t kMarkerHitWPoint = 24; + constexpr Marker_t kMarkerHitWOPoint = 28; + constexpr double kFontSize = 0.035; + constexpr Width_t kLineWidth = 1; + constexpr Style_t kLineMCTrackReconstructed = 9; + constexpr Style_t kLineMCTrackNotReconstructed = 10; + constexpr Style_t kLineRecoTrackGhost = 2; + constexpr Style_t kLineRecoTrackNotGhost = 1; - int iTrkMC = recoTrk.GetMatchedMCTrackIndex(); - if (iTrkMC > -1) { - const auto& mcTrk = fMCData.GetTrack(iTrkMC); - int pdg = mcTrk.GetPdgCode(); - bool isPrimary = mcTrk.IsPrimary(); + int iEvent = GetEventNumber(); + auto* pCanv = MakeQaObject<CbmQaCanvas>(Form("events/event_%d", iEvent), Form("event_%d", iEvent), kCanvX, kCanvY); + pCanv->Divide(1, 2, kVIMargin, kVIMargin); + pCanv->cd(1); + gPad->SetMargin(kLMargin, kRMarginDispl, kVIMargin, kVEMargin); + gPad->SetGrid(); + auto* pHistX = new TH2F("hFrameX", ";z [cm];x [cm]", 2, kZmin, kZmax, 2, kXmin, kXmax); + pHistX->GetYaxis()->SetTitleOffset(0.6); + pHistX->GetXaxis()->SetLabelSize(kFontSize); + pHistX->GetYaxis()->SetLabelSize(kFontSize); + pHistX->GetXaxis()->SetTitleSize(kFontSize); + pHistX->GetYaxis()->SetTitleSize(kFontSize); + pHistX->SetStats(false); + pHistX->Draw(); + pCanv->cd(2); + gPad->SetMargin(kLMargin, kRMarginDispl, kVEMargin, kVIMargin); + gPad->SetGrid(); + auto* pHistY = new TH2F("hFrameY", ";z [cm];y [cm]", 2, kZmin, kZmax, 2, kYmin, kYmax); + pHistY->GetYaxis()->SetTitleOffset(0.6); + pHistY->GetXaxis()->SetLabelSize(kFontSize); + pHistY->GetYaxis()->SetLabelSize(kFontSize); + pHistY->GetXaxis()->SetTitleSize(kFontSize); + pHistY->GetYaxis()->SetTitleSize(kFontSize); + pHistY->SetStats(false); + pHistY->Draw(); - // Cut tracks, which did not leave hits in tracker - if (mcTrk.GetNofHits() == 0) { - continue; - } + pCanv->cd(1); - if (isPrimary) { - FillRecoTrack(ETrackType::kPrim, iTrkReco); - bool bFast = mcTrk.GetP() > CbmL1Constants::MinFastMom; - bool bLong = mcTrk.GetTotNofStationsWithHit() == fpParameters->GetNstationsActive(); - if (bFast) { - FillRecoTrack(ETrackType::kPrimFast, iTrkReco); - } - if (bLong) { - FillRecoTrack(ETrackType::kPrimLong, iTrkReco); - } - if (bLong && bFast) { - FillRecoTrack(ETrackType::kPrimLongFast, iTrkReco); - } - } - else { - FillRecoTrack(ETrackType::kSec, iTrkReco); - if (mcTrk.GetP() > CbmL1Constants::MinFastMom) { - FillRecoTrack(ETrackType::kSecFast, iTrkReco); - } - } + auto* pHeader = new TLegend(kLMargin, 1 - kVEMargin + 0.01, 0.99, 0.99); + pHeader->SetNColumns(6); + pHeader->SetTextSize(kFontSize); + pHeader->SetMargin(0.1); + pHeader->AddEntry(static_cast<TObject*>(nullptr), Form("event #%d", iEvent), ""); + pHeader->AddEntry(new TMarker(0, 0, kMarkerPointWHit), "point w/ hit", "p"); + pHeader->AddEntry(new TMarker(0, 0, kMarkerPointWOHit), "point w/o hit", "p"); + pHeader->AddEntry(new TMarker(0, 0, kMarkerHitWPoint), "hit w/ point", "p"); + pHeader->AddEntry(new TMarker(0, 0, kMarkerHitWOPoint), "hit w/o point", "p"); + pHeader->Draw("same"); - // Track distributions for different particle species - switch (std::abs(pdg)) { - case 211: // pion - FillRecoTrack(ETrackType::kAllPI, iTrkReco); - if (isPrimary) { - FillRecoTrack(ETrackType::kPrimPI, iTrkReco); - FillRecoTrack((pdg > 0) ? ETrackType::kPrimPIP : ETrackType::kPrimPIM, iTrkReco); - } - else { - FillRecoTrack(ETrackType::kSecPI, iTrkReco); - FillRecoTrack((pdg > 0) ? ETrackType::kSecPIP : ETrackType::kSecPIM, iTrkReco); - } - break; - case 2212: // proton - FillRecoTrack(ETrackType::kAllPPBAR, iTrkReco); - if (isPrimary) { - FillRecoTrack(ETrackType::kPrimPPBAR, iTrkReco); - FillRecoTrack((pdg > 0) ? ETrackType::kPrimP : ETrackType::kPrimPBAR, iTrkReco); - } - else { - FillRecoTrack(ETrackType::kSecPPBAR, iTrkReco); - FillRecoTrack((pdg > 0) ? ETrackType::kSecP : ETrackType::kSecPBAR, iTrkReco); - } - break; - case 321: // kaon - FillRecoTrack(ETrackType::kAllK, iTrkReco); - if (isPrimary) { - FillRecoTrack(ETrackType::kPrimK, iTrkReco); - FillRecoTrack((pdg > 0) ? ETrackType::kPrimKP : ETrackType::kPrimKM, iTrkReco); - } - else { - FillRecoTrack(ETrackType::kSecK, iTrkReco); - FillRecoTrack((pdg > 0) ? ETrackType::kSecKP : ETrackType::kSecKM, iTrkReco); - } - break; - case 11: // electron - FillRecoTrack(ETrackType::kAllE, iTrkReco); - if (isPrimary) { - FillRecoTrack(ETrackType::kPrimE, iTrkReco); - FillRecoTrack((pdg < 0) ? ETrackType::kPrimEP : ETrackType::kPrimEM, iTrkReco); - } - else { - FillRecoTrack(ETrackType::kSecE, iTrkReco); - FillRecoTrack((pdg < 0) ? ETrackType::kSecEP : ETrackType::kSecEM, iTrkReco); - } - break; - case 13: // muon - FillRecoTrack(ETrackType::kAllMU, iTrkReco); - if (isPrimary) { - FillRecoTrack(ETrackType::kPrimMU, iTrkReco); - FillRecoTrack((pdg > 0) ? ETrackType::kPrimMUP : ETrackType::kPrimMUM, iTrkReco); - } - else { - FillRecoTrack(ETrackType::kSecMU, iTrkReco); - FillRecoTrack((pdg > 0) ? ETrackType::kSecMUP : ETrackType::kSecMUM, iTrkReco); - } - break; - } // switch abs(pdg): end - } - } - } // loop over recoTrk: end + auto* pLegendReco = new TLegend(1 - kRMarginDispl, kVIMargin, 0.99, 1 - kVEMargin, "Reco tracks"); + pLegendReco->SetMargin(0.1); + pLegendReco->SetTextSize(kFontSize); + //pLegendReco->SetHeader("Reco Tracks", "C"); + pCanv->cd(2); + auto* pLegendMC = new TLegend(1 - kRMarginDispl, kVEMargin, 0.99, 1 - kVIMargin, "MC tracks"); + pLegendMC->SetMargin(0.1); + pLegendMC->SetTextSize(kFontSize); + //pLegendMC->SetHeader("MC Tracks", "C"); + int iColorRec = 0; + int iColorGhost = 0; - // ************************************* - // ** Fill distributions of MC-tracks ** - // ************************************* + // Draw MC tracks + std::map<int, Color_t> mMCtrkColors; // Trk ID vs. color if (IsMCUsed()) { - for (int iTrkMC = 0; iTrkMC < fMCData.GetNofTracks(); ++iTrkMC) { - const auto& mcTrk = fMCData.GetTrack(iTrkMC); - - // ----- CUTS ON MC TRACKS - // Cut tracks, which did not leave hits in tracker - if (mcTrk.GetNofHits() == 0) { + // Draw MC tracks + int nMCTracks = fMCData.GetNofTracks(); + for (int iTmc = 0; iTmc < nMCTracks; ++iTmc) { + const auto& trk = fMCData.GetTrack(iTmc); + int nPoints = trk.GetNofPoints(); + if (nPoints == 0) { continue; } - - // Cut tracks, which cannot be reconstructed - if (!mcTrk.IsReconstructable()) { - continue; + std::vector<double> trkPointX(nPoints); + std::vector<double> trkPointY(nPoints); + std::vector<double> trkPointZ(nPoints); + for (int iPLoc = 0; iPLoc < nPoints; ++iPLoc) { + const auto& point = fMCData.GetPoint(trk.GetPointIndexes()[iPLoc]); + trkPointX[iPLoc] = point.GetX(); + trkPointY[iPLoc] = point.GetY(); + trkPointZ[iPLoc] = point.GetZ(); } - int pdg = mcTrk.GetPdgCode(); - bool isPrimary = mcTrk.IsPrimary(); + Color_t currColor = 1; + Style_t currStyle = trk.IsReconstructed() ? kLineMCTrackReconstructed : kLineMCTrackNotReconstructed; + currColor = kColorMC[iColorRec]; + iColorRec = (iColorRec + 1) % static_cast<int>(kColorMC.size()); + mMCtrkColors[trk.GetId()] = currColor; + { + pCanv->cd(1); + auto* gr = new TGraph(nPoints, trkPointZ.data(), trkPointX.data()); + gr->SetMarkerStyle(1); + gr->SetMarkerColor(currColor); + gr->SetLineColor(currColor); + gr->SetLineStyle(currStyle); + gr->SetLineWidth(kLineWidth); + gr->Draw("LSAME"); - // Fill different track categories - FillMCTrack(ETrackType::kAll, iTrkMC); - if (isPrimary) { - FillMCTrack(ETrackType::kPrim, iTrkMC); - bool bFast = mcTrk.GetP() > CbmL1Constants::MinFastMom; - bool bLong = mcTrk.GetTotNofStationsWithHit() == fpParameters->GetNstationsActive(); - if (bFast) { - FillMCTrack(ETrackType::kPrimFast, iTrkMC); - } - if (bLong) { - FillMCTrack(ETrackType::kPrimLong, iTrkMC); + std::stringstream msg; + msg << "ID=" << trk.GetId() << ", "; + msg << "PDG=" << trk.GetPdgCode() << ", "; + msg << "p=" << trk.GetP() << " GeV/c, "; + msg << "rec-able=" << trk.IsReconstructable() << ", "; + msg << "rec-ed=" << trk.IsReconstructed() << ", "; + if (trk.GetNofRecoTracks() > 0) { + msg << "reco_IDs=("; + for (int iTr : trk.GetRecoTrackIndexes()) { + msg << iTr << ' '; + } + msg << "), "; } - if (bLong && bFast) { - FillMCTrack(ETrackType::kPrimLongFast, iTrkMC); + if (trk.GetNofTouchTracks() > 0) { + msg << "touch_IDs=("; + for (int iTr : trk.GetTouchTrackIndexes()) { + msg << iTr << ' '; + } + msg << "), "; } + pLegendMC->AddEntry(gr, msg.str().c_str(), "l"); } - else { - FillMCTrack(ETrackType::kSec, iTrkMC); - if (mcTrk.GetP() > CbmL1Constants::MinFastMom) { - FillMCTrack(ETrackType::kSecFast, iTrkMC); - } + { + pCanv->cd(2); + auto* gr = new TGraph(nPoints, trkPointZ.data(), trkPointY.data()); + gr->SetMarkerStyle(1); + gr->SetMarkerColor(currColor); + gr->SetLineColor(currColor); + gr->SetLineStyle(currStyle); + gr->SetLineWidth(kLineWidth); + gr->Draw("LSAME"); } + } - // Track distributions for different particle species - switch (std::abs(pdg)) { - case 211: // pion - FillMCTrack(ETrackType::kAllPI, iTrkMC); - if (isPrimary) { - FillMCTrack(ETrackType::kPrimPI, iTrkMC); - FillMCTrack((pdg > 0) ? ETrackType::kPrimPIP : ETrackType::kPrimPIM, iTrkMC); - } - else { - FillMCTrack(ETrackType::kSecPI, iTrkMC); - FillMCTrack((pdg > 0) ? ETrackType::kSecPIP : ETrackType::kSecPIM, iTrkMC); - } - break; - case 2212: // proton - FillMCTrack(ETrackType::kAllPPBAR, iTrkMC); - if (isPrimary) { - FillMCTrack(ETrackType::kPrimPPBAR, iTrkMC); - FillMCTrack((pdg > 0) ? ETrackType::kPrimP : ETrackType::kPrimPBAR, iTrkMC); - } - else { - FillMCTrack(ETrackType::kSecPPBAR, iTrkMC); - FillMCTrack((pdg > 0) ? ETrackType::kSecP : ETrackType::kSecPBAR, iTrkMC); - } - break; - case 321: // kaon - FillMCTrack(ETrackType::kAllK, iTrkMC); - if (isPrimary) { - FillMCTrack(ETrackType::kPrimK, iTrkMC); - FillMCTrack((pdg > 0) ? ETrackType::kPrimKP : ETrackType::kPrimKM, iTrkMC); - } - else { - FillMCTrack(ETrackType::kSecK, iTrkMC); - FillMCTrack((pdg > 0) ? ETrackType::kSecKP : ETrackType::kSecKM, iTrkMC); - } - break; - case 11: // electron - FillMCTrack(ETrackType::kAllE, iTrkMC); - if (isPrimary) { - FillMCTrack(ETrackType::kPrimE, iTrkMC); - FillMCTrack((pdg < 0) ? ETrackType::kPrimEP : ETrackType::kPrimEM, iTrkMC); - } - else { - FillMCTrack(ETrackType::kSecE, iTrkMC); - FillMCTrack((pdg < 0) ? ETrackType::kSecEP : ETrackType::kSecEM, iTrkMC); - } - break; - case 13: // muon - FillMCTrack(ETrackType::kAllMU, iTrkMC); - if (isPrimary) { - FillMCTrack(ETrackType::kPrimMU, iTrkMC); - FillMCTrack((pdg > 0) ? ETrackType::kPrimMUP : ETrackType::kPrimMUM, iTrkMC); - } - else { - FillMCTrack(ETrackType::kSecMU, iTrkMC); - FillMCTrack((pdg > 0) ? ETrackType::kSecMUP : ETrackType::kSecMUM, iTrkMC); - } - break; - } // switch abs(pdg): end - } // iTrkMC - } // IsMCUsed() -} - -// --------------------------------------------------------------------------------------------------------------------- -// -InitStatus OutputQa::InitCanvases() -{ - /// Set of track types to compare - std::vector<ETrackType> vCmpTypesGeneral = {kAll, kPrim, kSec}; - std::vector<ETrackType> vCmpTypesPrim = {kPrim, kPrimE, kPrimMU, kPrimPI, kPrimK, kPrimPPBAR}; - std::vector<ETrackType> vCmpTypesSec = {kSec, kSecE, kSecMU, kSecPI, kSecK, kSecPPBAR}; - std::vector<ETrackType> vCmpTypesPions = {kAllPI, kPrimPIP, kPrimPIM, kSecPIP, kSecPIM}; - std::vector<ETrackType> vCmpTypesKaons = {kAllK, kPrimKP, kPrimKM, kSecKP, kSecKM}; - std::vector<ETrackType> vCmpTypesProtons = {kAllPPBAR, kPrimP, kPrimPBAR, kSecP, kSecPBAR}; - - /// @brief Function to draw generic canvas of histogram comparison - auto DrawTrackDistributions = [&](CbmQaCanvas* pCanv, std::function<TH1F*(ETrackType)> Hist) { - pCanv->Divide2D(6); - pCanv->cd(1); - gPad->SetLogy(); - DrawSetOf<TH1F>(vCmpTypesGeneral, Hist); - pCanv->cd(2); - gPad->SetLogy(); - DrawSetOf<TH1F>(vCmpTypesPrim, Hist); - pCanv->cd(3); - gPad->SetLogy(); - DrawSetOf<TH1F>(vCmpTypesSec, Hist); - pCanv->cd(4); - gPad->SetLogy(); - DrawSetOf<TH1F>(vCmpTypesPions, Hist); - pCanv->cd(5); - gPad->SetLogy(); - DrawSetOf<TH1F>(vCmpTypesKaons, Hist); - pCanv->cd(6); - gPad->SetLogy(); - DrawSetOf<TH1F>(vCmpTypesProtons, Hist); - }; - - /// @brief Function to draw generic canvas of efficiencies comparison - auto DrawTrackEfficiens = [&](CbmQaCanvas* pCanv, std::function<TProfile*(ETrackType)> Prof) { - pCanv->Divide2D(3); - pCanv->cd(1); - DrawSetOf<TProfile>(vCmpTypesGeneral, Prof); - pCanv->cd(2); - DrawSetOf<TProfile>(vCmpTypesPrim, Prof); - pCanv->cd(3); - DrawSetOf<TProfile>(vCmpTypesSec, Prof); - }; - - - if (IsMCUsed()) { - // ** Reconstructed track distributions ** - // Reconstructed pseudorapidity - auto* pc_reco_eta = - MakeQaObject<CbmQaCanvas>("reco_eta", "Reconstructed track pseudorapidity", kCXSIZEPX * 3, kCYSIZEPX * 2); - DrawTrackDistributions(pc_reco_eta, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_reco_eta; }); - - // MC pseudorapidity - auto* pc_reco_etaMC = - MakeQaObject<CbmQaCanvas>("reco_etaMC", "Reconstructed track MC pseudorapidity", kCXSIZEPX * 3, kCYSIZEPX * 2); - DrawTrackDistributions(pc_reco_etaMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_reco_etaMC; }); - - // MC momentum - auto* pc_reco_pMC = - MakeQaObject<CbmQaCanvas>("reco_pMC", "Reconstructed track MC momentum", kCXSIZEPX * 3, kCYSIZEPX * 2); - DrawTrackDistributions(pc_reco_pMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_reco_pMC; }); - - // MC rapidity - auto* pc_reco_yMC = - MakeQaObject<CbmQaCanvas>("reco_yMC", "Reconstructed track MC rapidity", kCXSIZEPX * 3, kCYSIZEPX * 2); - DrawTrackDistributions(pc_reco_yMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_reco_yMC; }); - - // ** MC track distributions ** - - // MC momentum - auto* pc_mc_pMC = - MakeQaObject<CbmQaCanvas>("mc_pMC", "MC reconstructable track MC momentum", kCXSIZEPX * 3, kCYSIZEPX * 2); - DrawTrackDistributions(pc_mc_pMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_mc_pMC; }); - - // MC rapidity - auto* pc_mc_yMC = - MakeQaObject<CbmQaCanvas>("mc_yMC", "MC reconstructable track MC rapidity", kCXSIZEPX * 3, kCYSIZEPX * 2); - DrawTrackDistributions(pc_mc_yMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_mc_yMC; }); - - // MC rapidity vs. MC momentum - // auto* pc_mc_pMC_yMC = - MakeQaObject<CbmQaCanvas>("mc_pMC_yMC", "MC track MC mom. vs. rapidity ", kCXSIZEPX * 3, kCYSIZEPX * 2); - DrawSetOf<TH2F>(vCmpTypesGeneral, [&](ETrackType t) -> TH2F* { return fvpTrackHistograms[t]->fph_reco_pMC_yMC; }); - - // ** Efficiencies ** - - // MC momentum - auto* pc_eff_pMC = MakeQaObject<CbmQaCanvas>("eff_pMC", "Tracking Eff. vs. MC momentum", kCXSIZEPX * 3, kCYSIZEPX); - DrawTrackEfficiens(pc_eff_pMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_pMC; }); + // Draw MC points + int nPoints = fMCData.GetNofPoints(); + for (int iP = 0; iP < nPoints; ++iP) { + const auto& point = fMCData.GetPoint(iP); + bool bHasHit = point.GetHitIndexes().size() > 0; + Marker_t style = bHasHit ? kMarkerPointWHit : kMarkerPointWOHit; + Color_t color = mMCtrkColors.at(point.GetTrackId()); + { + pCanv->cd(1); + auto* marker = new TMarker(point.GetZ(), point.GetX(), style); + marker->SetMarkerColor(color); + marker->Draw("same"); - auto* pc_eff_yMC = MakeQaObject<CbmQaCanvas>("eff_yMC", "Tracking Eff. vs. MC rapidity", kCXSIZEPX * 3, kCYSIZEPX); - DrawTrackEfficiens(pc_eff_yMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_yMC; }); + auto* pText = new TText(point.GetZ() + 2., point.GetX() + 2., Form("%d", point.GetStationId())); + pText->SetTextColor(color); + pText->SetTextSize(kFontSize); + pText->Draw("same"); + } + { + pCanv->cd(2); + auto* marker = new TMarker(point.GetZ(), point.GetY(), style); + marker->SetMarkerColor(color); + marker->Draw("same"); - auto* pc_eff_thetaMC = - MakeQaObject<CbmQaCanvas>("eff_thetaMC", "Tracking Eff. vs. MC polar angle", kCXSIZEPX * 3, kCYSIZEPX); - DrawTrackEfficiens(pc_eff_thetaMC, - [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_thetaMC; }); + auto* pText = new TText(point.GetZ() + 2., point.GetY() + 2., Form("%d", point.GetStationId())); + pText->SetTextColor(color); + pText->SetTextSize(kFontSize); + pText->Draw("same"); + } + } - auto* pc_eff_phiMC = - MakeQaObject<CbmQaCanvas>("eff_phiMC", "Tracking Eff. vs. MC azimuthal angle", kCXSIZEPX * 3, kCYSIZEPX); - DrawTrackEfficiens(pc_eff_phiMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_phiMC; }); + // Draw reconstructed tracks + int nRecoTracks = fvRecoTracks.size(); + std::vector<char> vbHitUsed(fvHits.size()); + std::vector<Color_t> vRecoTrkColors(fvHits.size()); // Reco hit ID vs. color + if (nRecoTracks > 0) { + for (int iTr = 0; iTr < nRecoTracks; ++iTr) { + const auto& trk = fvRecoTracks[iTr]; + Color_t currColor = 1; + Style_t currStyle = trk.IsGhost() ? kLineRecoTrackGhost : kLineRecoTrackNotGhost; + if (trk.IsGhost()) { + currColor = kColorGhost[iColorGhost]; + iColorGhost = (iColorGhost + 1) % static_cast<int>(kColorGhost.size()); + } + else { + int iTmc = trk.GetMatchedMCTrackIndex(); + currColor = iTmc > -1 ? mMCtrkColors[iTmc] : 1; + } - auto* pc_eff_etaMC = - MakeQaObject<CbmQaCanvas>("eff_etaMC", "Tracking Eff. vs. MC pseudorapidity", kCXSIZEPX * 3, kCYSIZEPX); - DrawTrackEfficiens(pc_eff_etaMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_etaMC; }); + int nHits = trk.GetNofHits(); + std::vector<double> trkHitX(nHits); + std::vector<double> trkHitY(nHits); + std::vector<double> trkHitZ(nHits); + for (int iHLoc = 0; iHLoc < nHits; ++iHLoc) { + int iH = trk.GetHitIndexes()[iHLoc]; + const auto& hit = fvHits[iH]; + vbHitUsed[iH] = true; + trkHitX[iHLoc] = hit.GetX(); + trkHitY[iHLoc] = hit.GetY(); + trkHitZ[iHLoc] = hit.GetZ(); + vRecoTrkColors[iH] = currColor; + } + { + pCanv->cd(1); + auto* gr = new TGraph(nHits, trkHitZ.data(), trkHitX.data()); + gr->SetMarkerStyle(1); + gr->SetMarkerColor(currColor); + gr->SetLineColor(currColor); + gr->SetLineStyle(currStyle); + gr->SetLineWidth(kLineWidth + 2); + gr->Draw("LSAME"); + } + { + pCanv->cd(2); + auto* gr = new TGraph(nHits, trkHitZ.data(), trkHitY.data()); + gr->SetMarkerStyle(1); + gr->SetMarkerColor(currColor); + gr->SetLineColor(currColor); + gr->SetLineStyle(currStyle); + gr->SetLineWidth(kLineWidth + 2); + gr->Draw("LSAME"); + std::stringstream msg; + msg << "ID=" << trk.index << ", "; + msg << "MC_IDs=("; + for (int iTmc : trk.GetMCTrackIndexes()) { + msg << iTmc << ' '; + } + msg << "), "; + msg << "purity=" << trk.GetMaxPurity(); + pLegendReco->AddEntry(gr, msg.str().c_str(), "l"); + } + } + } + // Draw hits + int nHits = fvHits.size(); + if (nHits > 0) { + for (int iH = 0; iH < nHits; ++iH) { + const auto& hit = fvHits[iH]; + bool bFake = hit.GetBestMcPointId() == -1; + bool bUsed = vbHitUsed[iH]; - // ** Pulls and residuals ** - // NOTE: stored in a subdirectory for a given track type and point type - for (int iType = 0; iType < ETrackType::kEND; ++iType) { - if (fvbTrackTypeOn[iType] && fvpTrackHistograms[iType]->IsMCUsed()) { - fvpTrackHistograms[iType]->fpFitQaFirstHit->CreateResidualPlot(); - fvpTrackHistograms[iType]->fpFitQaFirstHit->CreatePullPlot(); - fvpTrackHistograms[iType]->fpFitQaLastHit->CreateResidualPlot(); - fvpTrackHistograms[iType]->fpFitQaLastHit->CreatePullPlot(); + Marker_t style = bFake ? kMarkerHitWOPoint : kMarkerHitWPoint; + Color_t color = bUsed ? vRecoTrkColors[iH] : 1; + { + pCanv->cd(1); + auto* marker = new TMarker(hit.GetZ(), hit.GetX(), style); + marker->SetMarkerColor(color); + marker->Draw("same"); + } + { + pCanv->cd(2); + auto* marker = new TMarker(hit.GetZ(), hit.GetY(), style); + marker->SetMarkerColor(color); + marker->Draw("same"); + } } } } - - return kSUCCESS; + pCanv->cd(1); + pLegendReco->Draw(); + pCanv->cd(2); + pLegendMC->Draw(); } // --------------------------------------------------------------------------------------------------------------------- // -InitStatus OutputQa::InitDataBranches() +void OutputQa::FillHistograms() { + fMonitor.IncrementCounter(EMonitorKey::kEvent); - LOG_IF(fatal, !fpParameters.get()) - << fName << ": CA parameters object is not defined. Please, provide initializer or read parameters from binary " - << "via OutputQa::ReadParameters(filename) from the qa macro"; - - // Turn off detectors that are not used in the reconstruction setup - - fbUseMvd = fbUseMvd && (fpParameters->GetNstationsActive(ca::EDetectorID::kMvd) > 0); - fbUseSts = fbUseSts && (fpParameters->GetNstationsActive(ca::EDetectorID::kSts) > 0); - fbUseMuch = fbUseMuch && (fpParameters->GetNstationsActive(ca::EDetectorID::kMuch) > 0); - fbUseTrd = fbUseTrd && (fpParameters->GetNstationsActive(ca::EDetectorID::kTrd) > 0); - fbUseTof = fbUseTof && (fpParameters->GetNstationsActive(ca::EDetectorID::kTof) > 0); - - // Turn off detectors, which hits are not presented in input tree - - auto* fairManager = FairRootManager::Instance(); - fbUseMvd = fbUseMvd && fairManager->GetObject("MvdHit"); - fbUseSts = fbUseSts && fairManager->GetObject("StsHit"); - fbUseMuch = fbUseMuch && fairManager->GetObject("MuchPixelHit"); - fbUseTrd = fbUseTrd && fairManager->GetObject("TrdHit"); - fbUseTof = fbUseTof && fairManager->GetObject("TofHit"); - - - LOG(info) << fName << ": Detector subsystems used:"; - LOG(info) << "\tMVD: " << (fbUseMvd ? "ON" : "OFF"); - LOG(info) << "\tSTS: " << (fbUseSts ? "ON" : "OFF"); - LOG(info) << "\tMuCh: " << (fbUseMuch ? "ON" : "OFF"); - LOG(info) << "\tTRD: " << (fbUseTrd ? "ON" : "OFF"); - LOG(info) << "\tTOF: " << (fbUseTof ? "ON" : "OFF"); - - LOG(info) << fName << ": Initializing data branches"; - - // Initialize IO data manager - if (!fpDataManager.get()) { - fpDataManager = std::make_shared<ca::DataManager>(); - } - - // Initialize time slice reader instance - fpTSReader->SetTrackingMode(fTrackingMode); - fpTSReader->SetDetector(ca::EDetectorID::kMvd, fbUseMvd); - fpTSReader->SetDetector(ca::EDetectorID::kSts, fbUseSts); - fpTSReader->SetDetector(ca::EDetectorID::kMuch, fbUseMuch); - fpTSReader->SetDetector(ca::EDetectorID::kTrd, fbUseTrd); - fpTSReader->SetDetector(ca::EDetectorID::kTof, fbUseTof); - - fpTSReader->RegisterParameters(fpParameters); - fpTSReader->RegisterTracksContainer(fvRecoTracks); - fpTSReader->RegisterQaHitContainer(fvHits); - fpTSReader->RegisterHitIndexContainer(fvHitIds); - fpTSReader->SetSortQaHits(true); - if (!fpTSReader->InitRun()) { - return kFATAL; - } + // ************************************************************************ + // ** Fill distributions for reconstructed tracks (including ghost ones) ** + // ************************************************************************ - // Initialize MC module - if (IsMCUsed()) { - fpMCModule = std::make_shared<MCModule>(fVerbose, fPerformanceMode); - fpMCModule->SetDetector(ca::EDetectorID::kMvd, fbUseMvd); - fpMCModule->SetDetector(ca::EDetectorID::kSts, fbUseSts); - fpMCModule->SetDetector(ca::EDetectorID::kMuch, fbUseMuch); - fpMCModule->SetDetector(ca::EDetectorID::kTrd, fbUseTrd); - fpMCModule->SetDetector(ca::EDetectorID::kTof, fbUseTof); + for (size_t iTrkReco = 0; iTrkReco < fvRecoTracks.size(); ++iTrkReco) { + const auto& recoTrk = fvRecoTracks[iTrkReco]; - fpMCModule->RegisterMCData(fMCData); - fpMCModule->RegisterRecoTrackContainer(fvRecoTracks); - fpMCModule->RegisterHitIndexContainer(fvHitIds); - fpMCModule->RegisterQaHitContainer(fvHits); - fpMCModule->RegisterParameters(fpParameters); - fpMCModule->RegisterFirstHitIndexes(fpTSReader->GetHitFirstIndexDet()); - if (!fpMCModule->InitRun()) { - return kFATAL; + // Reject tracks, which do not contain hits + if (recoTrk.GetNofHits() < 1) { + continue; } - } - // Initialize monitor - fMonitor.SetCounterName(EMonitorKey::kEvent, "N events"); - fMonitor.SetCounterName(EMonitorKey::kTrack, "N reco tracks"); - fMonitor.SetCounterName(EMonitorKey::kHit, "N hits"); - fMonitor.SetCounterName(EMonitorKey::kMcTrack, "N MC tracks"); - fMonitor.SetCounterName(EMonitorKey::kMcPoint, "N MC points"); - fMonitor.SetRatioKeys({EMonitorKey::kEvent, EMonitorKey::kTrack}); + FillRecoTrack(ETrackType::kAll, iTrkReco); - return kSUCCESS; -} + if (IsMCUsed()) { + // NOTE: The ghost status of track is now defined by its purity, thus it can still contain MC information + if (recoTrk.IsGhost()) { + FillRecoTrack(ETrackType::kGhost, iTrkReco); + } -// --------------------------------------------------------------------------------------------------------------------- -// -InitStatus OutputQa::InitHistograms() -{ - auto RegisterTrackQa = [&](const char* typeName, const char* title, ETrackType type, bool bSuppressMC = false) { - if (!fvbTrackTypeOn[type]) { - return; - } - bool bUseMC = IsMCUsed() && !bSuppressMC; - fvsTrackTypeName[type] = typeName; - fvpTrackHistograms[type] = std::make_unique<TrackTypeQa>(typeName, fsPrefix.Data(), bUseMC, fpvObjList); - fvpTrackHistograms[type]->SetRootFolderName(fsRootFolderName + "/" + typeName); - fvpTrackHistograms[type]->SetTitle(title); - fvpTrackHistograms[type]->RegisterParameters(fpParameters); - fvpTrackHistograms[type]->RegisterRecoHits(fvHits); - fvpTrackHistograms[type]->RegisterRecoTracks(fvRecoTracks); - fvpTrackHistograms[type]->RegisterMCData(fMCData); - fvpTrackHistograms[type]->SetDrawAtt(fvTrackDrawAtts[type].fColor, fvTrackDrawAtts[type].fMarker); - fvpTrackHistograms[type]->Init(); - }; + int iTrkMC = recoTrk.GetMatchedMCTrackIndex(); + if (iTrkMC > -1) { + const auto& mcTrk = fMCData.GetTrack(iTrkMC); + int pdg = mcTrk.GetPdgCode(); + bool isPrimary = mcTrk.IsPrimary(); - //for (int i = 0; i < ETrackType::kEND; ++i) { - //LOG(info) << i << ' ' << fvpTrackHistograms[i].get() << ' ' << fvbTrackTypeOn[i]; - //} + // Cut tracks, which did not leave hits in tracker + if (mcTrk.GetNofHits() == 0) { + continue; + } - // TODO: Replace these parameters into the AddTrackType method!!! - RegisterTrackQa("all", "all", ETrackType::kAll); - if (IsMCUsed()) { - RegisterTrackQa("prim_long_fast", "primary long fast", ETrackType::kPrimLongFast); - RegisterTrackQa("prim_long", "primary long", ETrackType::kPrimLong); - RegisterTrackQa("ghost", "ghost", ETrackType::kGhost, true); - RegisterTrackQa("prim", "primary", ETrackType::kPrim); - RegisterTrackQa("prim_fast", "primary fast", ETrackType::kPrimFast); - RegisterTrackQa("sec", "secondary", ETrackType::kSec); - RegisterTrackQa("sec_fast", "secondary fast", ETrackType::kSecFast); - RegisterTrackQa("all_pi", "all #pi^{#pm}", ETrackType::kAllPI); - RegisterTrackQa("prim_pi", "primary #pi^{#pm}", ETrackType::kPrimPI); - RegisterTrackQa("prim_pip", "primary #pi^{#plus}", ETrackType::kPrimPIP); - RegisterTrackQa("prim_pim", "primary #pi^{#minus}", ETrackType::kPrimPIM); - RegisterTrackQa("sec_pi", "secondary #pi^{#pm}", ETrackType::kSecPI); - RegisterTrackQa("sec_pip", "secondary #pi^{#plus}", ETrackType::kSecPIP); - RegisterTrackQa("sec_pim", "secondary #pi^{#minus}", ETrackType::kSecPIM); - RegisterTrackQa("all_e", "all e^{#pm}", ETrackType::kAllE); - RegisterTrackQa("prim_e", "primary e^{#pm}", ETrackType::kPrimE); - RegisterTrackQa("prim_ep", "primary e^{#plus}", ETrackType::kPrimEP); - RegisterTrackQa("prim_em", "primary e^{#minus}", ETrackType::kPrimEM); - RegisterTrackQa("sec_e", "secondary e^{#pm}", ETrackType::kSecE); - RegisterTrackQa("sec_ep", "secondary e^{#plus}", ETrackType::kSecEP); - RegisterTrackQa("sec_em", "secondary e^{#minus}", ETrackType::kSecEM); - RegisterTrackQa("all_mu", "all #mu^{#pm}", ETrackType::kAllMU); - RegisterTrackQa("prim_mu", "primary #mu^{#pm}", ETrackType::kPrimMU); - RegisterTrackQa("prim_mup", "primary #mu^{#plus}", ETrackType::kPrimMUP); - RegisterTrackQa("prim_mum", "primary #mu^{#minus}", ETrackType::kPrimMUM); - RegisterTrackQa("sec_mu", "secondary #mu^{#pm}", ETrackType::kSecMU); - RegisterTrackQa("sec_mup", "secondary #mu^{#plus}", ETrackType::kSecMUP); - RegisterTrackQa("sec_mum", "secondary #mu^{#minus}", ETrackType::kSecMUM); - RegisterTrackQa("all_k", "all K^{#pm}", ETrackType::kAllK); - RegisterTrackQa("prim_k", "primary K^{#pm}", ETrackType::kPrimK); - RegisterTrackQa("prim_kp", "primary K^{#plus}", ETrackType::kPrimKP); - RegisterTrackQa("prim_km", "primary K^{#minus}", ETrackType::kPrimKM); - RegisterTrackQa("sec_k", "secondary K^{#pm}", ETrackType::kSecK); - RegisterTrackQa("sec_kp", "secondary K^{#plus}", ETrackType::kSecKP); - RegisterTrackQa("sec_km", "secondary K^{#minus}", ETrackType::kSecKM); - RegisterTrackQa("all_ppbar", "all p/#bar{p}", ETrackType::kAllPPBAR); - RegisterTrackQa("prim_ppbar", "primary p/#bar{p}", ETrackType::kPrimPPBAR); - RegisterTrackQa("prim_p", "primary p", ETrackType::kPrimP); - RegisterTrackQa("prim_pbar", "primary #bar{p}", ETrackType::kPrimPBAR); - RegisterTrackQa("sec_ppbar", "secondary p/#bar{p}", ETrackType::kSecPPBAR); - RegisterTrackQa("sec_p", "secondary p", ETrackType::kSecP); - RegisterTrackQa("sec_pbar", "secondary #bar{p}", ETrackType::kSecPBAR); - } + if (isPrimary) { + FillRecoTrack(ETrackType::kPrim, iTrkReco); + bool bFast = mcTrk.GetP() > CbmL1Constants::MinFastMom; + bool bLong = mcTrk.GetTotNofStationsWithHit() == fpParameters->GetNstationsActive(); + if (bFast) { + FillRecoTrack(ETrackType::kPrimFast, iTrkReco); + } + if (bLong) { + FillRecoTrack(ETrackType::kPrimLong, iTrkReco); + } + if (bLong && bFast) { + FillRecoTrack(ETrackType::kPrimLongFast, iTrkReco); + } + } + else { + FillRecoTrack(ETrackType::kSec, iTrkReco); + if (mcTrk.GetP() > CbmL1Constants::MinFastMom) { + FillRecoTrack(ETrackType::kSecFast, iTrkReco); + } + } - // Init default track types for the summary table - if (!fmSummaryTableEntries.size()) { - // clang-format off - fmSummaryTableEntries = { - ETrackType::kPrimLongFast, - ETrackType::kPrimLong, - ETrackType::kPrimFast, - ETrackType::kAll, - ETrackType::kPrim, - ETrackType::kSec - }; - // clang-format on - } + // Track distributions for different particle species + switch (std::abs(pdg)) { + case 211: // pion + FillRecoTrack(ETrackType::kAllPI, iTrkReco); + if (isPrimary) { + FillRecoTrack(ETrackType::kPrimPI, iTrkReco); + FillRecoTrack((pdg > 0) ? ETrackType::kPrimPIP : ETrackType::kPrimPIM, iTrkReco); + } + else { + FillRecoTrack(ETrackType::kSecPI, iTrkReco); + FillRecoTrack((pdg > 0) ? ETrackType::kSecPIP : ETrackType::kSecPIM, iTrkReco); + } + break; + case 2212: // proton + FillRecoTrack(ETrackType::kAllPPBAR, iTrkReco); + if (isPrimary) { + FillRecoTrack(ETrackType::kPrimPPBAR, iTrkReco); + FillRecoTrack((pdg > 0) ? ETrackType::kPrimP : ETrackType::kPrimPBAR, iTrkReco); + } + else { + FillRecoTrack(ETrackType::kSecPPBAR, iTrkReco); + FillRecoTrack((pdg > 0) ? ETrackType::kSecP : ETrackType::kSecPBAR, iTrkReco); + } + break; + case 321: // kaon + FillRecoTrack(ETrackType::kAllK, iTrkReco); + if (isPrimary) { + FillRecoTrack(ETrackType::kPrimK, iTrkReco); + FillRecoTrack((pdg > 0) ? ETrackType::kPrimKP : ETrackType::kPrimKM, iTrkReco); + } + else { + FillRecoTrack(ETrackType::kSecK, iTrkReco); + FillRecoTrack((pdg > 0) ? ETrackType::kSecKP : ETrackType::kSecKM, iTrkReco); + } + break; + case 11: // electron + FillRecoTrack(ETrackType::kAllE, iTrkReco); + if (isPrimary) { + FillRecoTrack(ETrackType::kPrimE, iTrkReco); + FillRecoTrack((pdg < 0) ? ETrackType::kPrimEP : ETrackType::kPrimEM, iTrkReco); + } + else { + FillRecoTrack(ETrackType::kSecE, iTrkReco); + FillRecoTrack((pdg < 0) ? ETrackType::kSecEP : ETrackType::kSecEM, iTrkReco); + } + break; + case 13: // muon + FillRecoTrack(ETrackType::kAllMU, iTrkReco); + if (isPrimary) { + FillRecoTrack(ETrackType::kPrimMU, iTrkReco); + FillRecoTrack((pdg > 0) ? ETrackType::kPrimMUP : ETrackType::kPrimMUM, iTrkReco); + } + else { + FillRecoTrack(ETrackType::kSecMU, iTrkReco); + FillRecoTrack((pdg > 0) ? ETrackType::kSecMUP : ETrackType::kSecMUM, iTrkReco); + } + break; + } // switch abs(pdg): end + } + } + } // loop over recoTrk: end - return kSUCCESS; -} -// --------------------------------------------------------------------------------------------------------------------- -// -void OutputQa::Check() -{ - // Create summary table + // ************************************* + // ** Fill distributions of MC-tracks ** + // ************************************* if (IsMCUsed()) { - fpMCModule->Finish(); + for (int iTrkMC = 0; iTrkMC < fMCData.GetNofTracks(); ++iTrkMC) { + const auto& mcTrk = fMCData.GetTrack(iTrkMC); - // TODO: Add cuts on entries from fmSummaryTableEntries - std::vector<ETrackType> vTypesToPlot; - int nRows = std::count_if(fmSummaryTableEntries.begin(), fmSummaryTableEntries.end(), - [&](const auto& f) { return fvbTrackTypeOn[f] && fvpTrackHistograms[f]->IsMCUsed(); }) - + 1; + // ----- CUTS ON MC TRACKS + // Cut tracks, which did not leave hits in tracker + if (mcTrk.GetNofHits() == 0) { + continue; + } - CbmQaTable* aTable = MakeQaObject<CbmQaTable>("summary_table", "Tracking summary table", nRows + 1, 9); - int iRow = 0; - std::vector<std::string> colNames = {"Efficiency", "Killed", "Length", "Fakes", "Clones", - "Reco/Evt", "MC/Evt", "Nst(hit)", "Nst(point)"}; - aTable->SetNamesOfCols(colNames); - aTable->SetColWidth(14); - double nEvents = static_cast<double>(GetEventNumber()); - LOG(info) << "Number of events: " << GetEventNumber(); - for (auto trType : fmSummaryTableEntries) { - if (!fvbTrackTypeOn[trType] || !fvpTrackHistograms[trType]->IsMCUsed()) { + // Cut tracks, which cannot be reconstructed + if (!mcTrk.IsReconstructable()) { continue; } - aTable->SetRowName(iRow, fvpTrackHistograms[trType]->GetTitle()); - aTable->SetCell(iRow, 0, fvpTrackHistograms[trType]->GetIntegratedEff()); - aTable->SetCell(iRow, 1, fvpTrackHistograms[trType]->GetKilledRate()); - aTable->SetCell(iRow, 2, fvpTrackHistograms[trType]->GetAverageRecoLength()); - aTable->SetCell(iRow, 3, fvpTrackHistograms[trType]->GetAverageFakeLength()); - aTable->SetCell(iRow, 4, fvpTrackHistograms[trType]->GetClonesRate()); - aTable->SetCell(iRow, 5, fvpTrackHistograms[trType]->GetNofRecoTracksMatched() / nEvents); - aTable->SetCell(iRow, 6, fvpTrackHistograms[trType]->GetNofMCTracks() / nEvents); - aTable->SetCell(iRow, 7, fvpTrackHistograms[trType]->GetAverageNofStationsWithHit()); - aTable->SetCell(iRow, 8, fvpTrackHistograms[trType]->GetAverageNofStationsWithPoint()); - ++iRow; - } - double nGhosts = 0.; - if (fvpTrackHistograms[ETrackType::kGhost] && fvpTrackHistograms[ETrackType::kAll]) { - nGhosts = fvpTrackHistograms[ETrackType::kGhost]->fph_reco_p->GetEntries(); - aTable->SetRowName(iRow, "N ghosts"); - aTable->SetCell(iRow, 0, nGhosts); - aTable->SetRowName(iRow + 1, "Ghost rate"); - aTable->SetCell(iRow + 1, 0, nGhosts / fMonitor.GetCounterValue(EMonitorKey::kTrack)); - } - LOG(info) << '\n' << aTable->ToString(3); - } + int pdg = mcTrk.GetPdgCode(); + bool isPrimary = mcTrk.IsPrimary(); - LOG(info) << '\n' << fMonitor.ToString(); + // Fill different track categories + FillMCTrack(ETrackType::kAll, iTrkMC); + if (isPrimary) { + FillMCTrack(ETrackType::kPrim, iTrkMC); + bool bFast = mcTrk.GetP() > CbmL1Constants::MinFastMom; + bool bLong = mcTrk.GetTotNofStationsWithHit() == fpParameters->GetNstationsActive(); + if (bFast) { + FillMCTrack(ETrackType::kPrimFast, iTrkMC); + } + if (bLong) { + FillMCTrack(ETrackType::kPrimLong, iTrkMC); + } + if (bLong && bFast) { + FillMCTrack(ETrackType::kPrimLongFast, iTrkMC); + } + } + else { + FillMCTrack(ETrackType::kSec, iTrkMC); + if (mcTrk.GetP() > CbmL1Constants::MinFastMom) { + FillMCTrack(ETrackType::kSecFast, iTrkMC); + } + } + + // Track distributions for different particle species + switch (std::abs(pdg)) { + case 211: // pion + FillMCTrack(ETrackType::kAllPI, iTrkMC); + if (isPrimary) { + FillMCTrack(ETrackType::kPrimPI, iTrkMC); + FillMCTrack((pdg > 0) ? ETrackType::kPrimPIP : ETrackType::kPrimPIM, iTrkMC); + } + else { + FillMCTrack(ETrackType::kSecPI, iTrkMC); + FillMCTrack((pdg > 0) ? ETrackType::kSecPIP : ETrackType::kSecPIM, iTrkMC); + } + break; + case 2212: // proton + FillMCTrack(ETrackType::kAllPPBAR, iTrkMC); + if (isPrimary) { + FillMCTrack(ETrackType::kPrimPPBAR, iTrkMC); + FillMCTrack((pdg > 0) ? ETrackType::kPrimP : ETrackType::kPrimPBAR, iTrkMC); + } + else { + FillMCTrack(ETrackType::kSecPPBAR, iTrkMC); + FillMCTrack((pdg > 0) ? ETrackType::kSecP : ETrackType::kSecPBAR, iTrkMC); + } + break; + case 321: // kaon + FillMCTrack(ETrackType::kAllK, iTrkMC); + if (isPrimary) { + FillMCTrack(ETrackType::kPrimK, iTrkMC); + FillMCTrack((pdg > 0) ? ETrackType::kPrimKP : ETrackType::kPrimKM, iTrkMC); + } + else { + FillMCTrack(ETrackType::kSecK, iTrkMC); + FillMCTrack((pdg > 0) ? ETrackType::kSecKP : ETrackType::kSecKM, iTrkMC); + } + break; + case 11: // electron + FillMCTrack(ETrackType::kAllE, iTrkMC); + if (isPrimary) { + FillMCTrack(ETrackType::kPrimE, iTrkMC); + FillMCTrack((pdg < 0) ? ETrackType::kPrimEP : ETrackType::kPrimEM, iTrkMC); + } + else { + FillMCTrack(ETrackType::kSecE, iTrkMC); + FillMCTrack((pdg < 0) ? ETrackType::kSecEP : ETrackType::kSecEM, iTrkMC); + } + break; + case 13: // muon + FillMCTrack(ETrackType::kAllMU, iTrkMC); + if (isPrimary) { + FillMCTrack(ETrackType::kPrimMU, iTrkMC); + FillMCTrack((pdg > 0) ? ETrackType::kPrimMUP : ETrackType::kPrimMUM, iTrkMC); + } + else { + FillMCTrack(ETrackType::kSecMU, iTrkMC); + FillMCTrack((pdg > 0) ? ETrackType::kSecMUP : ETrackType::kSecMUM, iTrkMC); + } + break; + } // switch abs(pdg): end + } // iTrkMC + } // IsMCUsed() } // --------------------------------------------------------------------------------------------------------------------- // -void OutputQa::DrawEvent() +InitStatus OutputQa::InitCanvases() { - constexpr double kZmin = 0.; - constexpr double kZmax = 300.; - constexpr double kXmin = -100.; - constexpr double kXmax = +100.; - constexpr double kYmin = -100.; - constexpr double kYmax = +100.; - constexpr std::array<Color_t, 11> kColorMC = {205, 209, 213, 217, 225, 208, 213, 216, 219, 224, 227}; - constexpr std::array<Color_t, 3> kColorGhost = {201, 202, 203}; - constexpr int kCanvX = 1920; - constexpr int kCanvY = 1080; - constexpr double kLMargin = 0.05; // Left margin - constexpr double kVEMargin = 0.15; // Vertical margin (top + bottom) - constexpr double kRMarginDispl = 0.4; - constexpr double kVIMargin = 0.0001; - constexpr Marker_t kMarkerPointWHit = 25; - constexpr Marker_t kMarkerPointWOHit = 5; - constexpr Marker_t kMarkerHitWPoint = 24; - constexpr Marker_t kMarkerHitWOPoint = 28; - constexpr double kFontSize = 0.035; - constexpr Width_t kLineWidth = 1; - constexpr Style_t kLineMCTrackReconstructed = 9; - constexpr Style_t kLineMCTrackNotReconstructed = 10; - constexpr Style_t kLineRecoTrackGhost = 2; - constexpr Style_t kLineRecoTrackNotGhost = 1; + /// Set of track types to compare + std::vector<ETrackType> vCmpTypesGeneral = {kAll, kPrim, kSec}; + std::vector<ETrackType> vCmpTypesPrim = {kPrim, kPrimE, kPrimMU, kPrimPI, kPrimK, kPrimPPBAR}; + std::vector<ETrackType> vCmpTypesSec = {kSec, kSecE, kSecMU, kSecPI, kSecK, kSecPPBAR}; + std::vector<ETrackType> vCmpTypesPions = {kAllPI, kPrimPIP, kPrimPIM, kSecPIP, kSecPIM}; + std::vector<ETrackType> vCmpTypesKaons = {kAllK, kPrimKP, kPrimKM, kSecKP, kSecKM}; + std::vector<ETrackType> vCmpTypesProtons = {kAllPPBAR, kPrimP, kPrimPBAR, kSecP, kSecPBAR}; - int iEvent = GetEventNumber(); - auto* pCanv = MakeQaObject<CbmQaCanvas>(Form("events/event_%d", iEvent), Form("event_%d", iEvent), kCanvX, kCanvY); - pCanv->Divide(1, 2, kVIMargin, kVIMargin); - pCanv->cd(1); - gPad->SetMargin(kLMargin, kRMarginDispl, kVIMargin, kVEMargin); - gPad->SetGrid(); - auto* pHistX = new TH2F("hFrameX", ";z [cm];x [cm]", 2, kZmin, kZmax, 2, kXmin, kXmax); - pHistX->GetYaxis()->SetTitleOffset(0.6); - pHistX->GetXaxis()->SetLabelSize(kFontSize); - pHistX->GetYaxis()->SetLabelSize(kFontSize); - pHistX->GetXaxis()->SetTitleSize(kFontSize); - pHistX->GetYaxis()->SetTitleSize(kFontSize); - pHistX->SetStats(false); - pHistX->Draw(); - pCanv->cd(2); - gPad->SetMargin(kLMargin, kRMarginDispl, kVEMargin, kVIMargin); - gPad->SetGrid(); - auto* pHistY = new TH2F("hFrameY", ";z [cm];y [cm]", 2, kZmin, kZmax, 2, kYmin, kYmax); - pHistY->GetYaxis()->SetTitleOffset(0.6); - pHistY->GetXaxis()->SetLabelSize(kFontSize); - pHistY->GetYaxis()->SetLabelSize(kFontSize); - pHistY->GetXaxis()->SetTitleSize(kFontSize); - pHistY->GetYaxis()->SetTitleSize(kFontSize); - pHistY->SetStats(false); - pHistY->Draw(); + /// @brief Function to draw generic canvas of histogram comparison + auto DrawTrackDistributions = [&](CbmQaCanvas* pCanv, std::function<TH1F*(ETrackType)> Hist) { + pCanv->Divide2D(6); + pCanv->cd(1); + gPad->SetLogy(); + DrawSetOf<TH1F>(vCmpTypesGeneral, Hist); + pCanv->cd(2); + gPad->SetLogy(); + DrawSetOf<TH1F>(vCmpTypesPrim, Hist); + pCanv->cd(3); + gPad->SetLogy(); + DrawSetOf<TH1F>(vCmpTypesSec, Hist); + pCanv->cd(4); + gPad->SetLogy(); + DrawSetOf<TH1F>(vCmpTypesPions, Hist); + pCanv->cd(5); + gPad->SetLogy(); + DrawSetOf<TH1F>(vCmpTypesKaons, Hist); + pCanv->cd(6); + gPad->SetLogy(); + DrawSetOf<TH1F>(vCmpTypesProtons, Hist); + }; + + /// @brief Function to draw generic canvas of efficiencies comparison + auto DrawTrackEfficiens = [&](CbmQaCanvas* pCanv, std::function<TProfile*(ETrackType)> Prof) { + pCanv->Divide2D(3); + pCanv->cd(1); + DrawSetOf<TProfile>(vCmpTypesGeneral, Prof); + pCanv->cd(2); + DrawSetOf<TProfile>(vCmpTypesPrim, Prof); + pCanv->cd(3); + DrawSetOf<TProfile>(vCmpTypesSec, Prof); + }; - pCanv->cd(1); - auto* pHeader = new TLegend(kLMargin, 1 - kVEMargin + 0.01, 0.99, 0.99); - pHeader->SetNColumns(6); - pHeader->SetTextSize(kFontSize); - pHeader->SetMargin(0.1); - pHeader->AddEntry(static_cast<TObject*>(nullptr), Form("event #%d", iEvent), ""); - pHeader->AddEntry(new TMarker(0, 0, kMarkerPointWHit), "point w/ hit", "p"); - pHeader->AddEntry(new TMarker(0, 0, kMarkerPointWOHit), "point w/o hit", "p"); - pHeader->AddEntry(new TMarker(0, 0, kMarkerHitWPoint), "hit w/ point", "p"); - pHeader->AddEntry(new TMarker(0, 0, kMarkerHitWOPoint), "hit w/o point", "p"); - pHeader->Draw("same"); + if (IsMCUsed()) { + // ** Reconstructed track distributions ** + // Reconstructed pseudorapidity + auto* pc_reco_eta = + MakeQaObject<CbmQaCanvas>("reco_eta", "Reconstructed track pseudorapidity", kCXSIZEPX * 3, kCYSIZEPX * 2); + DrawTrackDistributions(pc_reco_eta, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_reco_eta; }); - auto* pLegendReco = new TLegend(1 - kRMarginDispl, kVIMargin, 0.99, 1 - kVEMargin, "Reco tracks"); - pLegendReco->SetMargin(0.1); - pLegendReco->SetTextSize(kFontSize); - //pLegendReco->SetHeader("Reco Tracks", "C"); - pCanv->cd(2); - auto* pLegendMC = new TLegend(1 - kRMarginDispl, kVEMargin, 0.99, 1 - kVIMargin, "MC tracks"); - pLegendMC->SetMargin(0.1); - pLegendMC->SetTextSize(kFontSize); - //pLegendMC->SetHeader("MC Tracks", "C"); + // MC pseudorapidity + auto* pc_reco_etaMC = + MakeQaObject<CbmQaCanvas>("reco_etaMC", "Reconstructed track MC pseudorapidity", kCXSIZEPX * 3, kCYSIZEPX * 2); + DrawTrackDistributions(pc_reco_etaMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_reco_etaMC; }); - int iColorRec = 0; - int iColorGhost = 0; + // MC momentum + auto* pc_reco_pMC = + MakeQaObject<CbmQaCanvas>("reco_pMC", "Reconstructed track MC momentum", kCXSIZEPX * 3, kCYSIZEPX * 2); + DrawTrackDistributions(pc_reco_pMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_reco_pMC; }); - // Draw MC tracks - std::map<int, Color_t> mMCtrkColors; // Trk ID vs. color - if (IsMCUsed()) { - // Draw MC tracks - int nMCTracks = fMCData.GetNofTracks(); - for (int iTmc = 0; iTmc < nMCTracks; ++iTmc) { - const auto& trk = fMCData.GetTrack(iTmc); - int nPoints = trk.GetNofPoints(); - if (nPoints == 0) { - continue; - } - std::vector<double> trkPointX(nPoints); - std::vector<double> trkPointY(nPoints); - std::vector<double> trkPointZ(nPoints); - for (int iPLoc = 0; iPLoc < nPoints; ++iPLoc) { - const auto& point = fMCData.GetPoint(trk.GetPointIndexes()[iPLoc]); - trkPointX[iPLoc] = point.GetX(); - trkPointY[iPLoc] = point.GetY(); - trkPointZ[iPLoc] = point.GetZ(); - } - Color_t currColor = 1; - Style_t currStyle = trk.IsReconstructed() ? kLineMCTrackReconstructed : kLineMCTrackNotReconstructed; - currColor = kColorMC[iColorRec]; - iColorRec = (iColorRec + 1) % static_cast<int>(kColorMC.size()); - mMCtrkColors[trk.GetId()] = currColor; - { - pCanv->cd(1); - auto* gr = new TGraph(nPoints, trkPointZ.data(), trkPointX.data()); - gr->SetMarkerStyle(1); - gr->SetMarkerColor(currColor); - gr->SetLineColor(currColor); - gr->SetLineStyle(currStyle); - gr->SetLineWidth(kLineWidth); - gr->Draw("LSAME"); + // MC rapidity + auto* pc_reco_yMC = + MakeQaObject<CbmQaCanvas>("reco_yMC", "Reconstructed track MC rapidity", kCXSIZEPX * 3, kCYSIZEPX * 2); + DrawTrackDistributions(pc_reco_yMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_reco_yMC; }); - std::stringstream msg; - msg << "ID=" << trk.GetId() << ", "; - msg << "PDG=" << trk.GetPdgCode() << ", "; - msg << "p=" << trk.GetP() << " GeV/c, "; - msg << "rec-able=" << trk.IsReconstructable() << ", "; - msg << "rec-ed=" << trk.IsReconstructed() << ", "; - if (trk.GetNofRecoTracks() > 0) { - msg << "reco_IDs=("; - for (int iTr : trk.GetRecoTrackIndexes()) { - msg << iTr << ' '; - } - msg << "), "; - } - if (trk.GetNofTouchTracks() > 0) { - msg << "touch_IDs=("; - for (int iTr : trk.GetTouchTrackIndexes()) { - msg << iTr << ' '; - } - msg << "), "; - } - pLegendMC->AddEntry(gr, msg.str().c_str(), "l"); - } - { - pCanv->cd(2); - auto* gr = new TGraph(nPoints, trkPointZ.data(), trkPointY.data()); - gr->SetMarkerStyle(1); - gr->SetMarkerColor(currColor); - gr->SetLineColor(currColor); - gr->SetLineStyle(currStyle); - gr->SetLineWidth(kLineWidth); - gr->Draw("LSAME"); - } - } + // ** MC track distributions ** - // Draw MC points - int nPoints = fMCData.GetNofPoints(); - for (int iP = 0; iP < nPoints; ++iP) { - const auto& point = fMCData.GetPoint(iP); - bool bHasHit = point.GetHitIndexes().size() > 0; - Marker_t style = bHasHit ? kMarkerPointWHit : kMarkerPointWOHit; - Color_t color = mMCtrkColors.at(point.GetTrackId()); - { - pCanv->cd(1); - auto* marker = new TMarker(point.GetZ(), point.GetX(), style); - marker->SetMarkerColor(color); - marker->Draw("same"); + // MC momentum + auto* pc_mc_pMC = + MakeQaObject<CbmQaCanvas>("mc_pMC", "MC reconstructable track MC momentum", kCXSIZEPX * 3, kCYSIZEPX * 2); + DrawTrackDistributions(pc_mc_pMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_mc_pMC; }); - auto* pText = new TText(point.GetZ() + 2., point.GetX() + 2., Form("%d", point.GetStationId())); - pText->SetTextColor(color); - pText->SetTextSize(kFontSize); - pText->Draw("same"); - } - { - pCanv->cd(2); - auto* marker = new TMarker(point.GetZ(), point.GetY(), style); - marker->SetMarkerColor(color); - marker->Draw("same"); + // MC rapidity + auto* pc_mc_yMC = + MakeQaObject<CbmQaCanvas>("mc_yMC", "MC reconstructable track MC rapidity", kCXSIZEPX * 3, kCYSIZEPX * 2); + DrawTrackDistributions(pc_mc_yMC, [&](ETrackType t) -> TH1F* { return fvpTrackHistograms[t]->fph_mc_yMC; }); - auto* pText = new TText(point.GetZ() + 2., point.GetY() + 2., Form("%d", point.GetStationId())); - pText->SetTextColor(color); - pText->SetTextSize(kFontSize); - pText->Draw("same"); - } - } + // MC rapidity vs. MC momentum + // auto* pc_mc_pMC_yMC = + MakeQaObject<CbmQaCanvas>("mc_pMC_yMC", "MC track MC mom. vs. rapidity ", kCXSIZEPX * 3, kCYSIZEPX * 2); + DrawSetOf<TH2F>(vCmpTypesGeneral, [&](ETrackType t) -> TH2F* { return fvpTrackHistograms[t]->fph_reco_pMC_yMC; }); - // Draw reconstructed tracks - int nRecoTracks = fvRecoTracks.size(); - std::vector<char> vbHitUsed(fvHits.size()); - std::vector<Color_t> vRecoTrkColors(fvHits.size()); // Reco hit ID vs. color - if (nRecoTracks > 0) { - for (int iTr = 0; iTr < nRecoTracks; ++iTr) { - const auto& trk = fvRecoTracks[iTr]; - Color_t currColor = 1; - Style_t currStyle = trk.IsGhost() ? kLineRecoTrackGhost : kLineRecoTrackNotGhost; - if (trk.IsGhost()) { - currColor = kColorGhost[iColorGhost]; - iColorGhost = (iColorGhost + 1) % static_cast<int>(kColorGhost.size()); - } - else { - int iTmc = trk.GetMatchedMCTrackIndex(); - currColor = iTmc > -1 ? mMCtrkColors[iTmc] : 1; - } + // ** Efficiencies ** - int nHits = trk.GetNofHits(); - std::vector<double> trkHitX(nHits); - std::vector<double> trkHitY(nHits); - std::vector<double> trkHitZ(nHits); - for (int iHLoc = 0; iHLoc < nHits; ++iHLoc) { - int iH = trk.GetHitIndexes()[iHLoc]; - const auto& hit = fvHits[iH]; - vbHitUsed[iH] = true; - trkHitX[iHLoc] = hit.GetX(); - trkHitY[iHLoc] = hit.GetY(); - trkHitZ[iHLoc] = hit.GetZ(); - vRecoTrkColors[iH] = currColor; - } + // MC momentum + auto* pc_eff_pMC = MakeQaObject<CbmQaCanvas>("eff_pMC", "Tracking Eff. vs. MC momentum", kCXSIZEPX * 3, kCYSIZEPX); + DrawTrackEfficiens(pc_eff_pMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_pMC; }); - { - pCanv->cd(1); - auto* gr = new TGraph(nHits, trkHitZ.data(), trkHitX.data()); - gr->SetMarkerStyle(1); - gr->SetMarkerColor(currColor); - gr->SetLineColor(currColor); - gr->SetLineStyle(currStyle); - gr->SetLineWidth(kLineWidth + 2); - gr->Draw("LSAME"); - } - { - pCanv->cd(2); - auto* gr = new TGraph(nHits, trkHitZ.data(), trkHitY.data()); - gr->SetMarkerStyle(1); - gr->SetMarkerColor(currColor); - gr->SetLineColor(currColor); - gr->SetLineStyle(currStyle); - gr->SetLineWidth(kLineWidth + 2); - gr->Draw("LSAME"); - std::stringstream msg; - msg << "ID=" << trk.index << ", "; - msg << "MC_IDs=("; - for (int iTmc : trk.GetMCTrackIndexes()) { - msg << iTmc << ' '; - } - msg << "), "; - msg << "purity=" << trk.GetMaxPurity(); - pLegendReco->AddEntry(gr, msg.str().c_str(), "l"); - } - } - } - // Draw hits - int nHits = fvHits.size(); - if (nHits > 0) { - for (int iH = 0; iH < nHits; ++iH) { - const auto& hit = fvHits[iH]; - bool bFake = hit.GetBestMcPointId() == -1; - bool bUsed = vbHitUsed[iH]; + auto* pc_eff_yMC = MakeQaObject<CbmQaCanvas>("eff_yMC", "Tracking Eff. vs. MC rapidity", kCXSIZEPX * 3, kCYSIZEPX); + DrawTrackEfficiens(pc_eff_yMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_yMC; }); - Marker_t style = bFake ? kMarkerHitWOPoint : kMarkerHitWPoint; - Color_t color = bUsed ? vRecoTrkColors[iH] : 1; - { - pCanv->cd(1); - auto* marker = new TMarker(hit.GetZ(), hit.GetX(), style); - marker->SetMarkerColor(color); - marker->Draw("same"); - } - { - pCanv->cd(2); - auto* marker = new TMarker(hit.GetZ(), hit.GetY(), style); - marker->SetMarkerColor(color); - marker->Draw("same"); - } + auto* pc_eff_thetaMC = + MakeQaObject<CbmQaCanvas>("eff_thetaMC", "Tracking Eff. vs. MC polar angle", kCXSIZEPX * 3, kCYSIZEPX); + DrawTrackEfficiens(pc_eff_thetaMC, + [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_thetaMC; }); + + auto* pc_eff_phiMC = + MakeQaObject<CbmQaCanvas>("eff_phiMC", "Tracking Eff. vs. MC azimuthal angle", kCXSIZEPX * 3, kCYSIZEPX); + DrawTrackEfficiens(pc_eff_phiMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_phiMC; }); + + auto* pc_eff_etaMC = + MakeQaObject<CbmQaCanvas>("eff_etaMC", "Tracking Eff. vs. MC pseudorapidity", kCXSIZEPX * 3, kCYSIZEPX); + DrawTrackEfficiens(pc_eff_etaMC, [&](ETrackType t) -> TProfile* { return fvpTrackHistograms[t]->fph_eff_etaMC; }); + + + // ** Pulls and residuals ** + // NOTE: stored in a subdirectory for a given track type and point type + for (int iType = 0; iType < ETrackType::kEND; ++iType) { + if (fvbTrackTypeOn[iType] && fvpTrackHistograms[iType]->IsMCUsed()) { + fvpTrackHistograms[iType]->fpFitQaFirstHit->CreateResidualPlot(); + fvpTrackHistograms[iType]->fpFitQaFirstHit->CreatePullPlot(); + fvpTrackHistograms[iType]->fpFitQaLastHit->CreateResidualPlot(); + fvpTrackHistograms[iType]->fpFitQaLastHit->CreatePullPlot(); } } } - pCanv->cd(1); - pLegendReco->Draw(); - pCanv->cd(2); - pLegendMC->Draw(); + + return kSUCCESS; } // --------------------------------------------------------------------------------------------------------------------- // -InitStatus OutputQa::InitTimeSlice() +InitStatus OutputQa::InitDataBranches() { - int nMCTracks = 0; - int nMCPoints = 0; - int nHits = 0; - int nRecoTracks = 0; - // Read reconstructed input - fpTSReader->ReadEvent(this->GetCurrentEvent()); - nHits = fvHits.size(); - nRecoTracks = fvRecoTracks.size(); + LOG_IF(fatal, !fpParameters.get()) + << fName << ": CA parameters object is not defined. Please, provide initializer or read parameters from binary " + << "via OutputQa::ReadParameters(filename) from the qa macro"; - fMonitor.IncrementCounter(EMonitorKey::kTrack, nRecoTracks); - fMonitor.IncrementCounter(EMonitorKey::kHit, nHits); + // Turn off detectors that are not used in the reconstruction setup - // Match tracks and points - // Read MC tracks and points + fbUseMvd = fbUseMvd && (fpParameters->GetNstationsActive(ca::EDetectorID::kMvd) > 0); + fbUseSts = fbUseSts && (fpParameters->GetNstationsActive(ca::EDetectorID::kSts) > 0); + fbUseMuch = fbUseMuch && (fpParameters->GetNstationsActive(ca::EDetectorID::kMuch) > 0); + fbUseTrd = fbUseTrd && (fpParameters->GetNstationsActive(ca::EDetectorID::kTrd) > 0); + fbUseTof = fbUseTof && (fpParameters->GetNstationsActive(ca::EDetectorID::kTof) > 0); + + // Turn off detectors, which hits are not presented in input tree + + auto* fairManager = FairRootManager::Instance(); + fbUseMvd = fbUseMvd && fairManager->GetObject("MvdHit"); + fbUseSts = fbUseSts && fairManager->GetObject("StsHit"); + fbUseMuch = fbUseMuch && fairManager->GetObject("MuchPixelHit"); + fbUseTrd = fbUseTrd && fairManager->GetObject("TrdHit"); + fbUseTof = fbUseTof && fairManager->GetObject("TofHit"); + + + LOG(info) << fName << ": Detector subsystems used:"; + LOG(info) << "\tMVD: " << (fbUseMvd ? "ON" : "OFF"); + LOG(info) << "\tSTS: " << (fbUseSts ? "ON" : "OFF"); + LOG(info) << "\tMuCh: " << (fbUseMuch ? "ON" : "OFF"); + LOG(info) << "\tTRD: " << (fbUseTrd ? "ON" : "OFF"); + LOG(info) << "\tTOF: " << (fbUseTof ? "ON" : "OFF"); + + LOG(info) << fName << ": Initializing data branches"; + + // Initialize IO data manager + if (!fpDataManager.get()) { + fpDataManager = std::make_shared<ca::DataManager>(); + } + + // Initialize time slice reader instance + fpTSReader->SetTrackingMode(fTrackingMode); + fpTSReader->SetDetector(ca::EDetectorID::kMvd, fbUseMvd); + fpTSReader->SetDetector(ca::EDetectorID::kSts, fbUseSts); + fpTSReader->SetDetector(ca::EDetectorID::kMuch, fbUseMuch); + fpTSReader->SetDetector(ca::EDetectorID::kTrd, fbUseTrd); + fpTSReader->SetDetector(ca::EDetectorID::kTof, fbUseTof); + + fpTSReader->RegisterParameters(fpParameters); + fpTSReader->RegisterTracksContainer(fvRecoTracks); + fpTSReader->RegisterQaHitContainer(fvHits); + fpTSReader->RegisterHitIndexContainer(fvHitIds); + fpTSReader->SetSortQaHits(true); + if (!fpTSReader->InitRun()) { + return kFATAL; + } + + // Initialize MC module if (IsMCUsed()) { - fpMCModule->InitEvent(this->GetCurrentEvent()); - nMCPoints = fMCData.GetNofPoints(); - nMCTracks = fMCData.GetNofTracks(); - fpMCModule->MatchHits(); - fpMCModule->MatchTracks(); - fMonitor.IncrementCounter(EMonitorKey::kMcPoint, nMCPoints); - fMonitor.IncrementCounter(EMonitorKey::kMcTrack, nMCTracks); + fpMCModule = std::make_shared<MCModule>(fVerbose, fPerformanceMode); + fpMCModule->SetDetector(ca::EDetectorID::kMvd, fbUseMvd); + fpMCModule->SetDetector(ca::EDetectorID::kSts, fbUseSts); + fpMCModule->SetDetector(ca::EDetectorID::kMuch, fbUseMuch); + fpMCModule->SetDetector(ca::EDetectorID::kTrd, fbUseTrd); + fpMCModule->SetDetector(ca::EDetectorID::kTof, fbUseTof); - if (fbDrawEvents && nMCPoints > std::max(0, fEvtDisplayMinNofPoints)) { - DrawEvent(); + fpMCModule->RegisterMCData(fMCData); + fpMCModule->RegisterRecoTrackContainer(fvRecoTracks); + fpMCModule->RegisterHitIndexContainer(fvHitIds); + fpMCModule->RegisterQaHitContainer(fvHits); + fpMCModule->RegisterParameters(fpParameters); + fpMCModule->RegisterFirstHitIndexes(fpTSReader->GetHitFirstIndexDet()); + if (!fpMCModule->InitRun()) { + return kFATAL; } } - LOG_IF(info, fVerbose > 1) << fName << ": Data sample consists of " << nHits << " hits, " << nRecoTracks - << " reco tracks, " << nMCTracks << " MC tracks, " << nMCPoints << " MC points"; + // Initialize monitor + fMonitor.SetCounterName(EMonitorKey::kEvent, "N events"); + fMonitor.SetCounterName(EMonitorKey::kTrack, "N reco tracks"); + fMonitor.SetCounterName(EMonitorKey::kHit, "N hits"); + fMonitor.SetCounterName(EMonitorKey::kMcTrack, "N MC tracks"); + fMonitor.SetCounterName(EMonitorKey::kMcPoint, "N MC points"); + fMonitor.SetRatioKeys({EMonitorKey::kEvent, EMonitorKey::kTrack}); return kSUCCESS; } -// --------------------------------------------------------------------------------------------------------------------- -// -void OutputQa::ReadParameters(const char* filename) -{ - ca::InitManager manager; - manager.ReadParametersObject(filename); - fpParameters = std::make_shared<ca::Parameters>(manager.TakeParameters()); - - LOG(info) << '\n' << fpParameters->ToString(1); -} - // --------------------------------------------------------------------------------------------------------------------- // void OutputQa::InitDrawingAttributes() @@ -1032,3 +893,142 @@ void OutputQa::InitDrawingAttributes() fvTrackDrawAtts[ETrackType::kSecMUP] = {kMagenta - 6, 26}; fvTrackDrawAtts[ETrackType::kSecMUM] = {kMagenta - 10, 32}; } + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus OutputQa::InitHistograms() +{ + auto RegisterTrackQa = [&](const char* typeName, const char* title, ETrackType type, bool bSuppressMC = false) { + if (!fvbTrackTypeOn[type]) { + return; + } + bool bUseMC = IsMCUsed() && !bSuppressMC; + fvsTrackTypeName[type] = typeName; + fvpTrackHistograms[type] = std::make_unique<TrackTypeQa>(typeName, fsPrefix.Data(), bUseMC, fpvObjList); + fvpTrackHistograms[type]->SetRootFolderName(fsRootFolderName + "/" + typeName); + fvpTrackHistograms[type]->SetTitle(title); + fvpTrackHistograms[type]->RegisterParameters(fpParameters); + fvpTrackHistograms[type]->RegisterRecoHits(fvHits); + fvpTrackHistograms[type]->RegisterRecoTracks(fvRecoTracks); + fvpTrackHistograms[type]->RegisterMCData(fMCData); + fvpTrackHistograms[type]->SetDrawAtt(fvTrackDrawAtts[type].fColor, fvTrackDrawAtts[type].fMarker); + fvpTrackHistograms[type]->Init(); + }; + + //for (int i = 0; i < ETrackType::kEND; ++i) { + //LOG(info) << i << ' ' << fvpTrackHistograms[i].get() << ' ' << fvbTrackTypeOn[i]; + //} + + // TODO: Replace these parameters into the AddTrackType method!!! + RegisterTrackQa("all", "all", ETrackType::kAll); + if (IsMCUsed()) { + RegisterTrackQa("prim_long_fast", "primary long fast", ETrackType::kPrimLongFast); + RegisterTrackQa("prim_long", "primary long", ETrackType::kPrimLong); + RegisterTrackQa("ghost", "ghost", ETrackType::kGhost, true); + RegisterTrackQa("prim", "primary", ETrackType::kPrim); + RegisterTrackQa("prim_fast", "primary fast", ETrackType::kPrimFast); + RegisterTrackQa("sec", "secondary", ETrackType::kSec); + RegisterTrackQa("sec_fast", "secondary fast", ETrackType::kSecFast); + RegisterTrackQa("all_pi", "all #pi^{#pm}", ETrackType::kAllPI); + RegisterTrackQa("prim_pi", "primary #pi^{#pm}", ETrackType::kPrimPI); + RegisterTrackQa("prim_pip", "primary #pi^{#plus}", ETrackType::kPrimPIP); + RegisterTrackQa("prim_pim", "primary #pi^{#minus}", ETrackType::kPrimPIM); + RegisterTrackQa("sec_pi", "secondary #pi^{#pm}", ETrackType::kSecPI); + RegisterTrackQa("sec_pip", "secondary #pi^{#plus}", ETrackType::kSecPIP); + RegisterTrackQa("sec_pim", "secondary #pi^{#minus}", ETrackType::kSecPIM); + RegisterTrackQa("all_e", "all e^{#pm}", ETrackType::kAllE); + RegisterTrackQa("prim_e", "primary e^{#pm}", ETrackType::kPrimE); + RegisterTrackQa("prim_ep", "primary e^{#plus}", ETrackType::kPrimEP); + RegisterTrackQa("prim_em", "primary e^{#minus}", ETrackType::kPrimEM); + RegisterTrackQa("sec_e", "secondary e^{#pm}", ETrackType::kSecE); + RegisterTrackQa("sec_ep", "secondary e^{#plus}", ETrackType::kSecEP); + RegisterTrackQa("sec_em", "secondary e^{#minus}", ETrackType::kSecEM); + RegisterTrackQa("all_mu", "all #mu^{#pm}", ETrackType::kAllMU); + RegisterTrackQa("prim_mu", "primary #mu^{#pm}", ETrackType::kPrimMU); + RegisterTrackQa("prim_mup", "primary #mu^{#plus}", ETrackType::kPrimMUP); + RegisterTrackQa("prim_mum", "primary #mu^{#minus}", ETrackType::kPrimMUM); + RegisterTrackQa("sec_mu", "secondary #mu^{#pm}", ETrackType::kSecMU); + RegisterTrackQa("sec_mup", "secondary #mu^{#plus}", ETrackType::kSecMUP); + RegisterTrackQa("sec_mum", "secondary #mu^{#minus}", ETrackType::kSecMUM); + RegisterTrackQa("all_k", "all K^{#pm}", ETrackType::kAllK); + RegisterTrackQa("prim_k", "primary K^{#pm}", ETrackType::kPrimK); + RegisterTrackQa("prim_kp", "primary K^{#plus}", ETrackType::kPrimKP); + RegisterTrackQa("prim_km", "primary K^{#minus}", ETrackType::kPrimKM); + RegisterTrackQa("sec_k", "secondary K^{#pm}", ETrackType::kSecK); + RegisterTrackQa("sec_kp", "secondary K^{#plus}", ETrackType::kSecKP); + RegisterTrackQa("sec_km", "secondary K^{#minus}", ETrackType::kSecKM); + RegisterTrackQa("all_ppbar", "all p/#bar{p}", ETrackType::kAllPPBAR); + RegisterTrackQa("prim_ppbar", "primary p/#bar{p}", ETrackType::kPrimPPBAR); + RegisterTrackQa("prim_p", "primary p", ETrackType::kPrimP); + RegisterTrackQa("prim_pbar", "primary #bar{p}", ETrackType::kPrimPBAR); + RegisterTrackQa("sec_ppbar", "secondary p/#bar{p}", ETrackType::kSecPPBAR); + RegisterTrackQa("sec_p", "secondary p", ETrackType::kSecP); + RegisterTrackQa("sec_pbar", "secondary #bar{p}", ETrackType::kSecPBAR); + } + + // Init default track types for the summary table + if (!fmSummaryTableEntries.size()) { + // clang-format off + fmSummaryTableEntries = { + ETrackType::kPrimLongFast, + ETrackType::kPrimLong, + ETrackType::kPrimFast, + ETrackType::kAll, + ETrackType::kPrim, + ETrackType::kSec + }; + // clang-format on + } + + return kSUCCESS; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus OutputQa::InitTimeSlice() +{ + int nMCTracks = 0; + int nMCPoints = 0; + int nHits = 0; + int nRecoTracks = 0; + + // Read reconstructed input + fpTSReader->ReadEvent(this->GetCurrentEvent()); + nHits = fvHits.size(); + nRecoTracks = fvRecoTracks.size(); + + fMonitor.IncrementCounter(EMonitorKey::kTrack, nRecoTracks); + fMonitor.IncrementCounter(EMonitorKey::kHit, nHits); + + // Match tracks and points + // Read MC tracks and points + if (IsMCUsed()) { + fpMCModule->InitEvent(this->GetCurrentEvent()); + nMCPoints = fMCData.GetNofPoints(); + nMCTracks = fMCData.GetNofTracks(); + fpMCModule->MatchHits(); + fpMCModule->MatchTracks(); + fMonitor.IncrementCounter(EMonitorKey::kMcPoint, nMCPoints); + fMonitor.IncrementCounter(EMonitorKey::kMcTrack, nMCTracks); + + if (fbDrawEvents && nMCPoints > std::max(0, fEvtDisplayMinNofPoints)) { + DrawEvent(); + } + } + LOG_IF(info, fVerbose > 1) << fName << ": Data sample consists of " << nHits << " hits, " << nRecoTracks + << " reco tracks, " << nMCTracks << " MC tracks, " << nMCPoints << " MC points"; + + + return kSUCCESS; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void OutputQa::ReadParameters(const char* filename) +{ + ca::InitManager manager; + manager.ReadParametersObject(filename); + fpParameters = std::make_shared<ca::Parameters>(manager.TakeParameters()); + + LOG(info) << '\n' << fpParameters->ToString(1); +}