diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C index d914202eb461efb4b9e078f924589d0c51bafd12..d0bcb8e4f9b77431925ba8d954155a5dfe4ed90b 100644 --- a/macro/mcbm/mcbm_qa.C +++ b/macro/mcbm/mcbm_qa.C @@ -28,6 +28,7 @@ #include "CbmKF.h" #include "CbmMCDataManager.h" #include "CbmMatchRecoToMC.h" +#include "CbmMuchDigitizerQa.h" #include "CbmMuchGeoScheme.h" #include "CbmMuchHitFinderQa.h" #include "CbmMuchTransportQa.h" @@ -61,6 +62,30 @@ void mcbm_qa(Int_t nEvents = 0, TString setupName = "mcbm_beam_2020_03", Bool_t bUseMC = kTRUE, TString config = "", + TString benchmarkInput = ""); +/* clang-format on */ + +/* clang-format off */ +/// \brief Qa macro execution function with independent input names +/// \param nEvents Number of events to proceed +/// \param traColFile Collision transport input +/// \param rawFile Digi input +/// \param recFile Reconstruction input +/// \param sinkFile Output filename +/// \param setupName Name of the setup (TODO: can be a run number) +/// \param bUseMC Flag for MC (simulation) usage +/// \param config QA YAML configuraiton file +/// \param benchmarkInput Path to a benchmark QA output, obtained for a given setup and given number of events +void mcbm_qa(Int_t nEvents = 0, + TString traColFile = "data/mcbm_beam_2022_05_nickel.tra.root", + TString rawFile = "data/mcbm_beam_2022_05_nickel.event.raw.root", + TString recFile = "data/mcbm_beam_2022_05_nickel.rec.root", + TString parFile = "data/mcbm_beam_2022_05_nickel.par.root", + TString geoFile = "data/mcbm_beam_2022_05_nickel.geo.root", + TString sinkFile = "data/mcbm_beam_2022_05_nickel.qa.root", + TString setupName = "mcbm_beam_2022_05_nickel", + TString bUseMC = kTRUE, + TString config = "", TString benchmarkInput = "") /* clang-format on */ { @@ -70,27 +95,21 @@ void mcbm_qa(Int_t nEvents = 0, // ----- Logger settings ---------------------------------------------- FairLogger::GetLogger()->SetLogScreenLevel("INFO"); - FairLogger::GetLogger()->SetLogVerbosityLevel("LOW"); + FairLogger::GetLogger()->SetLogVerbosityLevel("HIGH"); // ------------------------------------------------------------------------ // ----- Environment -------------------------------------------------- - int verbose = 6; // verbose level + int verbose = 2; // verbose level TString myName = "mcbm_qa"; // this macro's name for screen output TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory // ------------------------------------------------------------------------ // ----- In- and output file names ------------------------------------ - TString rawFile = dataset + ".event.raw.root"; - TString traFile = dataset + ".tra.root"; - TString parFile = dataset + ".par.root"; - TString recFile = dataset + ".rec.root"; - TString sinkFile = dataset + ".qa.root"; TString qaConfig = (config.Length() ? config : srcDir + "/macro/qa/configs/qa_tasks_config_" + setupName + ".yaml"); if (gSystem->AccessPathName(qaConfig.Data())) { // file not found, using default one std::cout << "-I- " << myName << ": the QA configuration file " << qaConfig << " not found, using default one\n"; qaConfig = srcDir + "/macro/qa/configs/qa_tasks_config_mcbm.yaml"; } - TString benchmarkOut = sinkFile + ".qa.benchmark.root"; // ------------------------------------------------------------------------ @@ -180,13 +199,15 @@ void mcbm_qa(Int_t nEvents = 0, // ----- FairRunAna --------------------------------------------------- FairFileSource* inputSource = new FairFileSource(rawFile); - inputSource->AddFriend(traFile); + if (bUseMC) { + inputSource->AddFriend(traColFile); + } inputSource->AddFriend(recFile); FairRunAna* run = new FairRunAna(); run->SetSource(inputSource); run->SetGenerateRunInfo(kFALSE); - + run->SetGeomFile(geoFile); FairRootFileSink* sink = new FairRootFileSink(sinkFile); run->SetSink(sink); @@ -213,7 +234,7 @@ void mcbm_qa(Int_t nEvents = 0, std::cout << "-I- " << myName << ": Adding MC manager and MC to reco matching tasks\n"; auto* mcManager = new CbmMCDataManager("MCDataManager", 1); - mcManager->AddFile(traFile); + mcManager->AddFile(traColFile); qaManager->AddTask(mcManager); auto* matchRecoToMC = new CbmMatchRecoToMC(); @@ -373,3 +394,23 @@ void mcbm_qa(Int_t nEvents = 0, RemoveGeoManager(); // ------------------------------------------------------------------------ } + +/* clang-format off */ +/// \brief QA macro execution function +/// \param nEvents Number of events to proceed +/// \param dataset Prefix of the output files files upstream the simulation/reconstruction chain +/// \param setupName Name of the setup +/// \param bUseMC Flag for MC (simulation) usage +/// \param config QA YAML configuraiton file +/// \param benchmarkInput Path to a benchmark QA output, obtained for a given setup and given number of events +void mcbm_qa(Int_t nEvents, TString dataset, TString setupName, Bool_t bUseMC, TString config, TString benchmarkInput) +{ + TString rawFile = dataset + ".event.raw.root"; + TString traFile = dataset + ".tra.root"; + TString parFile = dataset + ".par.root"; + TString geoFile = dataset + ".geo.root"; + TString recFile = dataset + ".rec.root"; + TString sinkFile = dataset + ".qa.root"; + mcbm_qa(nEvents, traFile, rawFile, recFile, parFile, geoFile, sinkFile, setupName, bUseMC, config, benchmarkInput); +} +/* clang-format on */ diff --git a/macro/qa/CMakeLists.txt b/macro/qa/CMakeLists.txt index 3f283fd15726ab8636fa549fe33f843bcad4e853..40d12babebdc070a9266d9fdf1460ba6832638e5 100644 --- a/macro/qa/CMakeLists.txt +++ b/macro/qa/CMakeLists.txt @@ -99,7 +99,7 @@ foreach(setup IN LISTS cbm_setup) endforeach(setup IN LISTS cbm_setup) # ============================================================================ -Install(FILES ../run/.rootrc qa_compare.C ./qa_report_ca_offline.C +Install(FILES ../run/.rootrc qa_compare.C ./qa_report_ca_offline.C ./run_recoQa.C DESTINATION share/cbmroot/macro/qa ) #Install(CODE "FILE(MAKE_DIRECTORY \${CMAKE_INSTALL_PREFIX}/share/cbmroot/macro/run/data)") diff --git a/macro/run/CMakeLists.txt b/macro/run/CMakeLists.txt index bf0a50166a81b5a64148203b659ea3b73a05f2ad..6f12f249203d96b2ad038de8b45c60b72f42096a 100644 --- a/macro/run/CMakeLists.txt +++ b/macro/run/CMakeLists.txt @@ -603,7 +603,7 @@ EndIf() # If(DEFINED ENV{RAW_DATA_PATH} ) ##################### # ============================================================================ -Install(FILES .rootrc run_tra_file.C run_tra_beam.C run_digi.C run_reco.C run_qa.C run_unpack_online.C run_unpack_tsa.C +Install(FILES .rootrc run_tra_file.C run_tra_beam.C run_digi.C run_reco.C run_qa.C run_unpack_online.C run_unpack_tsa.C create_mcbm_geo_setup.C qa_config.cbm.yaml DESTINATION share/cbmroot/macro/run ) diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index 87188873f21b30b6d14cea4dd87d9df385ecd50c..1de14f1fc65fa9589a79f35f449218dfbbf3ae03 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -150,6 +150,7 @@ set(PUBLIC_DEPENDENCIES ROOT::Graf ROOT::Hist ROOT::Physics + fmt::fmt ) set(PRIVATE_DEPENDENCIES diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx index 437f1beb1dffe661f941d42b8600d0ebe01b1185..3c458329d28bca73906b222ff0c5c64d7fac8847 100644 --- a/reco/L1/CbmCaMCModule.cxx +++ b/reco/L1/CbmCaMCModule.cxx @@ -429,7 +429,12 @@ void MCModule::ReadMCPointsForDetector<ca::EDetectorID::kTof>() if (iPext < 0) { continue; } - auto* pExtPoint = static_cast<CbmTofPoint*>(pBrPoints->Get(link)); + auto* pExtPoint = dynamic_cast<CbmTofPoint*>(pBrPoints->Get(link)); + if (!pExtPoint) { + LOG(error) << "ca::MCModule: MC Point with link=" << link.ToString() << " does not exist in branch"; + continue; + } + int trkId = pExtPoint->GetTrackID(); int rpcAddr = pExtPoint->GetDetectorID() << kNofBitsRpcAddress; // FIXME: auto key = std::make_pair(trkId, rpcAddr); @@ -537,12 +542,18 @@ void MCModule::ReadMCTracks() // ----- Total number of tracks int nTracksTot = 0; for (const auto& [iFile, iEvent] : fFileEventIDs) { + if (iFile < 0 || iEvent < 0) { + continue; + } nTracksTot += fpMCTracks->Size(iFile, iEvent); /// iFile, iEvent } fpMCData->ReserveNofTracks(nTracksTot); // ----- Loop over MC events for (const auto& [iFile, iEvent] : fFileEventIDs) { + if (iFile < 0 || iEvent < 0) { + continue; + } auto pEvtHeader = dynamic_cast<FairMCEventHeader*>(fpMCEventHeader->Get(iFile, iEvent)); if (!pEvtHeader) { LOG(fatal) << "cbm::ca::MCModule: event header is not found for file " << iFile << " and event " << iEvent; diff --git a/reco/L1/qa/CbmCaInputQaBase.cxx b/reco/L1/qa/CbmCaInputQaBase.cxx index 613c1c1d91a6c3e86d394ff460dcc9a6f5107afd..823787553c6dc4bbeefbb88a57a87d2c239d8b3b 100644 --- a/reco/L1/qa/CbmCaInputQaBase.cxx +++ b/reco/L1/qa/CbmCaInputQaBase.cxx @@ -321,6 +321,8 @@ void CbmCaInputQaBase<DetID>::FillHistograms() } } + LOG(info) << fName << ": Number of hits: " << nHits; + for (int iH = 0; iH < nHits; ++iH) { fHitQaData.Reset(); fHitQaData.SetHitIndex(iH); diff --git a/reco/L1/qa/CbmCaInputQaBase.h b/reco/L1/qa/CbmCaInputQaBase.h index 9f70ed3fcf9d9729559fc91abfbab79aed541140..610571c0d69abc8413dd3c763dfb7159018e9866 100644 --- a/reco/L1/qa/CbmCaInputQaBase.h +++ b/reco/L1/qa/CbmCaInputQaBase.h @@ -304,6 +304,7 @@ class CbmCaInputQaBase : public CbmQaTask { std::vector<TProfile2D*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC std::vector<TH1F*> fvph_reco_eff; ///< Distribution of hit reconstruction efficiency in bins of fvpe_reco_eff_vs_xy + // FIXME: change to private protected: cbm::ca::HitQaData fHitQaData{}; ///< Current hit QA data object diff --git a/reco/L1/qa/CbmCaInputQaSetup.cxx b/reco/L1/qa/CbmCaInputQaSetup.cxx index 4cec7abf9afafe4a3cc1c2cb66cb3a8f4d32cea2..2513fae945d9bed738217c85382925afd09d68e3 100644 --- a/reco/L1/qa/CbmCaInputQaSetup.cxx +++ b/reco/L1/qa/CbmCaInputQaSetup.cxx @@ -52,7 +52,8 @@ void InputQaSetup::CheckInit() const template<cbm::algo::ca::EDetectorID DetID> void InputQaSetup::FillHistogramsDet() { - using Hit_t = HitTypes_t::at<DetID>; + using Hit_t = HitTypes_t::at<DetID>; + using McPoint_t = PointTypes_t::at<DetID>; int nHits = fvpBrHits[DetID]->GetEntriesFast(); for (int iH = 0; iH < nHits; ++iH) { @@ -97,12 +98,50 @@ void InputQaSetup::FillHistogramsDet() fvZmax[DetID][iStLoc] = zHit; } if (iStGeo >= 0) { + fph_hit_xy[iStGeo]->Fill(xHit, yHit); fph_hit_xz[iStGeo]->Fill(zHit, xHit); fph_hit_yz[iStGeo]->Fill(zHit, yHit); } fph_hit_xz.back()->Fill(zHit, xHit); fph_hit_yz.back()->Fill(zHit, yHit); } // iH + + if (IsMCUsed()) { + int nMCevents = fpMCEventList->GetNofEvents(); + for (int iE = 0; iE < nMCevents; ++iE) { + int iFile = fpMCEventList->GetFileIdByIndex(iE); + int iEvent = fpMCEventList->GetEventIdByIndex(iE); + int nPoints = fvpBrPoints[DetID]->Size(iFile, iEvent); + + for (int iP = 0; iP < nPoints; ++iP) { + const auto* pPoint = dynamic_cast<const McPoint_t*>(fvpBrPoints[DetID]->Get(iFile, iEvent, iP)); + if (!pPoint) { + LOG(error) << fName << ": point with iFile=" << iFile << ", iEvent=" << iEvent << ", iP=" << iP + << " for detector " << kDetName[DetID] << " is not found"; + continue; + } + auto address = pPoint->GetDetectorID(); + + int iStLoc = fvpDetInterface[DetID]->GetTrackingStationIndex(address); + if (iStLoc < 0) { + continue; + } + + int iStGeo = fpParameters->GetStationIndexGeometry(iStLoc, DetID); + auto xPoint = pPoint->FairMCPoint::GetX(); + auto yPoint = pPoint->FairMCPoint::GetY(); + auto zPoint = pPoint->FairMCPoint::GetZ(); + + if (iStGeo >= 0) { + fph_point_xy[iStGeo]->Fill(xPoint, yPoint); + fph_point_xz[iStGeo]->Fill(zPoint, xPoint); + fph_point_yz[iStGeo]->Fill(zPoint, yPoint); + } + fph_point_xz.back()->Fill(zPoint, xPoint); + fph_point_yz.back()->Fill(zPoint, yPoint); + } // iP + } // iE + } } // --------------------------------------------------------------------------------------------------------------------- @@ -223,7 +262,10 @@ try { assert(pFairManager); auto pMcManager = IsMCUsed() ? dynamic_cast<CbmMCDataManager*>(pFairManager->GetObject("MCDataManager")) : nullptr; - fpMCEventHeader = IsMCUsed() ? pMcManager->GetObject("MCEventHeader.") : nullptr; + if (IsMCUsed()) { + fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairManager->GetObject("MCEventList.")); + fpMCEventHeader = pMcManager->GetObject("MCEventHeader."); + } fvpBrPoints.fill(nullptr); fvpBrHits.fill(nullptr); @@ -293,9 +335,19 @@ InitStatus InputQaSetup::InitHistograms() MakeQaDirectory(Form("hit_occupancy/%s", kDetName[iD])); } } + MakeQaDirectory("point_occupancy"); + for (int iD = 0; iD < static_cast<int>(ca::EDetectorID::kEND); ++iD) { + if (fvbUseDet[iD]) { + MakeQaDirectory(Form("point_occupancy/%s", kDetName[iD])); + } + } + fph_hit_xy.resize(nStGeo); fph_hit_xz.resize(nStGeo + 1); fph_hit_yz.resize(nStGeo + 1); + fph_point_xy.resize(nStGeo); + fph_point_xz.resize(nStGeo + 1); + fph_point_yz.resize(nStGeo + 1); /// TEMPORARY constexpr int kNbinsZ = 300; @@ -317,14 +369,38 @@ InitStatus InputQaSetup::InitHistograms() TString sTsuff = (iStGeo == nStGeo) ? "" : Form(" in %s station %d", kDetName[detID], iStLoc); // Hit occupancy - sF = (iStGeo == nStGeo) ? "hit_occupancy" : TString(Form("hit_occupancy/%s", kDetName[detID])); - sN = (TString) "hit_xz" + sNsuff; - sT = (TString) "hit occupancy in xz-plane" + sTsuff + ";z_{hit} [cm];x_{hit} [cm]"; + sF = Form("hit_occupancy%s", ((iStGeo == nStGeo) ? "" : Form("/%s", kDetName[detID]))); + sN = Form("hit_xz%s", sNsuff.Data()); + sT = Form("hit occupancy in xz-plane%s;z_{hit} [cm];x_{hit} [cm]", sTsuff.Data()); fph_hit_xz[iStGeo] = MakeQaObject<TH2F>(sF + "/" + sN, sT, kNbinsZ, kMinZ, kMaxZ, kNbinsX, kMinX, kMaxX); - sN = (TString) "hit_yz" + sNsuff; - sT = (TString) "hit occupancy in yz-plane" + sTsuff + ";z_{hit} [cm];y_{hit} [cm]"; + sN = Form("hit_yz%s", sNsuff.Data()); + sT = Form("hit occupancy in yz-plane%s;z_{hit} [cm];y_{hit} [cm]", sTsuff.Data()); fph_hit_yz[iStGeo] = MakeQaObject<TH2F>(sF + "/" + sN, sT, kNbinsZ, kMinZ, kMaxZ, kNbinsY, kMinY, kMaxY); + + if (iStGeo < nStGeo) { + sN = Form("hit_xy%s", sNsuff.Data()); + sT = Form("hit occupancy in xy-plane%s;x_{hit} [cm];y_{hit} [cm]", sTsuff.Data()); + fph_hit_xy[iStGeo] = MakeQaObject<TH2F>(sF + "/" + sN, sT, kNbinsX, kMinX, kMaxX, kNbinsY, kMinY, kMaxY); + } + + if (IsMCUsed()) { + // Point occupancy + sF = Form("point_occupancy%s", ((iStGeo == nStGeo) ? "" : Form("/%s", kDetName[detID]))); + sN = Form("point_xz%s", sNsuff.Data()); + sT = Form("point occupancy in xz-plane%s;z_{point} [cm];x_{point} [cm]", sTsuff.Data()); + fph_point_xz[iStGeo] = MakeQaObject<TH2F>(sF + "/" + sN, sT, kNbinsZ, kMinZ, kMaxZ, kNbinsX, kMinX, kMaxX); + + sN = Form("point_yz%s", sNsuff.Data()); + sT = Form("point occupancy in yz-plane%s;z_{point} [cm];y_{point} [cm]", sTsuff.Data()); + fph_point_yz[iStGeo] = MakeQaObject<TH2F>(sF + "/" + sN, sT, kNbinsZ, kMinZ, kMaxZ, kNbinsY, kMinY, kMaxY); + + if (iStGeo < nStGeo) { + sN = Form("point_xy%s", sNsuff.Data()); + sT = Form("point occupancy in xy-plane%s;x_{point} [cm];y_{point} [cm]", sTsuff.Data()); + fph_point_xy[iStGeo] = MakeQaObject<TH2F>(sF + "/" + sN, sT, kNbinsX, kMinX, kMaxX, kNbinsY, kMinY, kMaxY); + } + } } return kSUCCESS; diff --git a/reco/L1/qa/CbmCaInputQaSetup.h b/reco/L1/qa/CbmCaInputQaSetup.h index 112cefd6789587182a991717fa11d9be3991cee7..0b7c382405d1f8f7a9b53903925d8d20ee66253c 100644 --- a/reco/L1/qa/CbmCaInputQaSetup.h +++ b/reco/L1/qa/CbmCaInputQaSetup.h @@ -14,6 +14,7 @@ #include "CbmL1DetectorID.h" #include "CbmMCDataArray.h" #include "CbmMCDataObject.h" +#include "CbmMCEventList.h" #include "CbmMuchPixelHit.h" #include "CbmMuchPoint.h" #include "CbmMuchTrackingInterface.h" @@ -90,6 +91,7 @@ namespace cbm::ca DetIdArr_t<CbmMCDataArray*> fvpBrPoints = {{nullptr}}; ///< Input branch for MC points CbmMCDataObject* fpMCEventHeader = nullptr; ///< Pointer to MC event header + CbmMCEventList* fpMCEventList = nullptr; ///< Pointer to MC event list std::unique_ptr<ca::Parameters<ca::fvec>> fpParameters = nullptr; ///< Pointer to CA parameters object std::string fsParametersFilename = ""; ///< Filename for the tracking parameters @@ -104,8 +106,10 @@ namespace cbm::ca // ** Histograms ** // ****************** + std::vector<TH2F*> fph_hit_xy = {}; ///< Hit occupancy in x-y plane vs. station std::vector<TH2F*> fph_hit_xz = {}; ///< Hit occupancy in x-z plane vs. station std::vector<TH2F*> fph_hit_yz = {}; ///< Hit occupancy in y-z plane vs. station + std::vector<TH2F*> fph_point_xy = {}; ///< MC point occupancy in x-z plane vs. station std::vector<TH2F*> fph_point_xz = {}; ///< MC point occupancy in x-z plane vs. station std::vector<TH2F*> fph_point_yz = {}; ///< MC point occupancy in y-z plane vs. station }; diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx index 4083bba43fc31bb3653e081f56f54a5291e6002b..4dd3738fd54b69ddd5272522ceac916d5ff53a1c 100644 --- a/reco/L1/qa/CbmCaInputQaSts.cxx +++ b/reco/L1/qa/CbmCaInputQaSts.cxx @@ -120,7 +120,7 @@ void CbmCaInputQaSts::FillHistogramsPerHit() { const auto* pHit = dynamic_cast<const CbmStsHit*>(fpHits->At(fHitQaData.GetHitIndex())); - { // u-coordinate residual per number of digis + if (IsMCUsed()) { // u-coordinate residual per number of digis const auto* pCluster = dynamic_cast<const CbmStsCluster*>(fpClusters->At(pHit->GetFrontClusterId())); assert(pCluster); int nDigis = pCluster->GetNofDigis(); @@ -130,7 +130,7 @@ void CbmCaInputQaSts::FillHistogramsPerHit() fvph_pull_u_Ndig[nDigis]->Fill(fHitQaData.GetPullU()); } - { // v-coordinate residual per number of digis + if (IsMCUsed()) { // v-coordinate residual per number of digis const auto* pCluster = dynamic_cast<const CbmStsCluster*>(fpClusters->At(pHit->GetBackClusterId())); assert(pCluster); int nDigis = pCluster->GetNofDigis(); @@ -233,7 +233,6 @@ InitStatus CbmCaInputQaSts::InitHistograms() // ----- Initialize histograms, which are use MC-information if (IsMCUsed()) { - // Resize histogram vectors fvph_pull_u_Ndig.resize(fkMaxDigisInClusterForPulls + 1, nullptr); fvph_pull_v_Ndig.resize(fkMaxDigisInClusterForPulls + 1, nullptr); diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx index 815754cd5c887cdc8ecd3039fbba7866c9039193..1e0cba6d711da9da3ead9013913476878eaa338b 100644 --- a/reco/L1/qa/CbmCaOutputQa.cxx +++ b/reco/L1/qa/CbmCaOutputQa.cxx @@ -77,6 +77,8 @@ OutputQa::OutputQa(int verbose, bool isMCUsed, ECbmRecoMode recoMode) // Init track type histograms drawing attributes InitDrawingAttributes(); + + SetProcessFullTs(false); // If CbmEvent branch exists, the routine will run in event-based regime } // ---------------------------------------------------------------------------------------------------------------------