diff --git a/macro/qa/run_recoQa.C b/macro/qa/run_recoQa.C
new file mode 100644
index 0000000000000000000000000000000000000000..709753cc8116fe89800d1c260032f13de2b36f09
--- /dev/null
+++ b/macro/qa/run_recoQa.C
@@ -0,0 +1,209 @@
+/* Copyright (C) 2024 Hulubei National Institute of Physics and Nuclear Engineering - Horia Hulubei, Bucharest
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Alexandru Bercuci [committer] */
+
+// --------------------------------------------------------------------------
+//
+// Macro for simulation & reconstruction QA
+//
+// The following naming conventions are assumed:
+// Raw data file: [dataset].event.raw.root
+// Transport file: [dataset].tra.root
+// Parameter file: [dataset].par.root
+// Reconstruction file: [dataset].rec.root
+//
+// S. Gorbunov 28/09/2020
+//
+// --------------------------------------------------------------------------
+
+// Includes needed for IDE
+#if !defined(__CLING__)
+#include "CbmDefs.h"
+#include "CbmMCDataManager.h"
+#include "CbmMuchTransportQa.h"
+#include "CbmQaManager.h"
+#include "CbmSetup.h"
+
+#include <FairFileSource.h>
+#include <FairMonitor.h>
+#include <FairParAsciiFileIo.h>
+#include <FairParRootFileIo.h>
+#include <FairRootFileSink.h>
+#include <FairRunAna.h>
+#include <FairRuntimeDb.h>
+#include <FairSystemInfo.h>
+
+#include <TStopwatch.h>
+#endif
+
+void run_recoQa(Int_t nEvents = -1, TString dataset = "2391_node8_0_0000",
+ TString setupName = "mcbm_beam_2022_05_23_nickel", Bool_t bUseMC = kFALSE)
+{
+
+ // ========================================================================
+ // Adjust this part according to your requirements
+
+ // ----- Logger settings ----------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel("DEBUG");
+ FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+ // ------------------------------------------------------------------------
+
+ // ----- Environment --------------------------------------------------
+ int verbose = 6; // verbose level
+ TString myName = "recoQA"; // this macro's name for screen output
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+ // ----- In- and output file names ------------------------------------
+ TString traFile = dataset + ".tra.root";
+ //TString rawFile = dataset + ".event.raw.root";
+ TString parFile = dataset + ".par.root";
+ TString recFile = dataset + ".rec.root";
+ TString sinkFile = dataset + ".rqa.root";
+ TString geoFile = setupName + ".geo.root";
+
+ // ------------------------------------------------------------------------
+
+ // ----- Load the geometry setup -------------------------------------
+ std::cout << '\n';
+
+ CbmSetup* setup = CbmSetup::Instance();
+ setup->LoadSetup(setupName);
+
+ // In general, the following parts need not be touched
+ // ========================================================================
+
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // ------------------------------------------------------------------------
+
+ // ---- Debug option -------------------------------------------------
+ gDebug = 0;
+ // ------------------------------------------------------------------------
+
+ // ----- FairRunAna ---------------------------------------------------
+ FairFileSource* inputSource = new FairFileSource(recFile);
+ if (bUseMC) {
+ inputSource->AddFriend(traFile);
+ }
+ //inputSource->AddFriend(rawFile);
+
+ FairRunAna* run = new FairRunAna();
+ run->SetSource(inputSource);
+ run->SetGenerateRunInfo(kFALSE);
+ run->SetGeomFile(geoFile);
+
+ FairRootFileSink* sink = new FairRootFileSink(sinkFile);
+ run->SetSink(sink);
+
+ TString monitorFile{sinkFile};
+ monitorFile.ReplaceAll("rqa", "rqa.moni");
+ FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
+ // ------------------------------------------------------------------------
+
+ // ----- MCDataManager (legacy mode) -----------------------------------
+ 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();
+ // NOTE: Matching is suppressed, if there are hit and cluster matches already in the tree. If there
+ // are no hit matches, they are produced on this stage.
+ matchRecoToMC->SuppressHitReMatching();
+ run->AddTask(matchRecoToMC);
+ }
+ // ------------------------------------------------------------------------
+
+ // ----- Tracking detector interface --------------------------------------
+ run->AddTask(new CbmTrackingDetectorInterfaceInit());
+ // ------------------------------------------------------------------------
+
+
+ // ------------------------------------------------------------------------
+
+ // TODO: read tracking parameters from the file like CaOutputQa does
+ //TODO: instead of initializing the tracker
+ //TString caParFile = recFile;
+ //caParFile.ReplaceAll(".root", ".L1Parameters.dat");
+ //auto* pCaOutputQa = new cbm::ca::OutputQa(verbose, bUseMC);
+ //pCaOutputQa->ReadParameters(caParFile.Data());
+
+ // L1 CA track finder setup
+ auto l1 = new CbmL1("CA");
+ l1->SetMcbmMode();
+
+ // User configuration example for CA:
+ //l1->SetConfigUser(srcDir + "/macro/L1/configs/ca_params_user_example.yaml");
+ run->AddTask(l1);
+
+ // ----- Reco QA --------------------------------------------
+ CbmRecoQaTask* recoQa = new CbmRecoQaTask();
+ recoQa->SetSetupClass(CbmRecoQaTask::kMcbm22);
+ double cos25 = TMath::Cos(TMath::DegToRad() * 25.);
+ recoQa->SetProjections({15.1 * cos25, 0., -20 * cos25, -38 * cos25, -50.5 * cos25});
+ run->AddTask(recoQa);
+ run->SetGenerateRunInfo(kFALSE);
+
+ // ------------------------------------------------------------------------
+ //----- Load Parameters --------------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Set runtime DB" << std::endl;
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ FairParRootFileIo* parIo1 = new FairParRootFileIo();
+ parIo1->open(parFile.Data(), "in");
+ rtdb->setFirstInput(parIo1);
+
+ // ------------------------------------------------------------------------
+ // ----- Run initialisation -------------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+
+ rtdb->print();
+
+ // ----- Start run ----------------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(0, nEvents);
+ // ------------------------------------------------------------------------
+
+ // ----- Finish -------------------------------------------------------
+ timer.Stop();
+
+
+ FairMonitor::GetMonitor()->Print();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished successfully." << std::endl;
+ std::cout << "Output file is " << sinkFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << std::endl;
+ std::cout << std::endl;
+ std::cout << " Test passed" << std::endl;
+ std::cout << " All ok " << std::endl;
+ // ------------------------------------------------------------------------
+
+ // ----- Resource monitoring ------------------------------------------
+ // Extract the maximal used memory an add is as Dart measurement
+ // This line is filtered by CTest and the value send to CDash
+ FairSystemInfo sysInfo;
+ Float_t maxMemory = sysInfo.GetMaxMemory();
+ std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+ std::cout << maxMemory;
+ std::cout << "</DartMeasurement>" << std::endl;
+
+ Float_t cpuUsage = ctime / rtime;
+ std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+ std::cout << cpuUsage;
+ std::cout << "</DartMeasurement>" << std::endl;
+ // ------------------------------------------------------------------------
+
+ // ----- Function needed for CTest runtime dependency -----------------
+ RemoveGeoManager();
+ // ------------------------------------------------------------------------
+}
diff --git a/reco/qa/CMakeLists.txt b/reco/qa/CMakeLists.txt
index 136abb007695c24efdbc9a2d1faf215712835a11..85ce9991cef57402c931747c7bffa8fbd8d89fa6 100644
--- a/reco/qa/CMakeLists.txt
+++ b/reco/qa/CMakeLists.txt
@@ -4,9 +4,9 @@ set(INCLUDE_DIRECTORIES
set(SRCS
- CbmRecoQaTask.cxx
CbmRecoQa.cxx
CbmTrackingTrdQa.cxx
+ CbmRecoQaTask.cxx
)
set(LIBRARY_NAME CbmRecoQa)
diff --git a/reco/qa/CbmRecoQaTask.cxx b/reco/qa/CbmRecoQaTask.cxx
index 26bb61d2ee985c5c83d2181587a2970bc78af247..3388c58dc07bb1e93f49cfe55b5c8dda5b6f63b4 100644
--- a/reco/qa/CbmRecoQaTask.cxx
+++ b/reco/qa/CbmRecoQaTask.cxx
@@ -3,23 +3,22 @@
Authors: Alexandru Bercuci [committer] */
// CBM headers
-#include "CbmGlobalTrack.h"
#include "CbmRecoQaTask.h"
+
+#include "CbmGlobalTrack.h"
+#include "CbmTimeSlice.h"
+#include "CbmTofAddress.h"
+#include "CbmTrdAddress.h"
// FAIR headers
#include "FairRootManager.h"
+#include "FairSink.h"
// ROOT headers
#include "TClonesArray.h"
-//#include "TObject.h"
+#include "TH2D.h"
+using namespace std;
//_____________________________________________________________________
-CbmRecoQaTask::CbmRecoQaTask() : FairTask("RecoQA", 0)
-{
-}
-
-
-//_____________________________________________________________________
-CbmRecoQaTask::~CbmRecoQaTask()
-{}
+CbmRecoQaTask::CbmRecoQaTask() : FairTask("RecoQA", 0) {}
//_____________________________________________________________________
InitStatus CbmRecoQaTask::Init()
@@ -31,19 +30,273 @@ InitStatus CbmRecoQaTask::Init()
FairRootManager* ioman = FairRootManager::Instance();
if (!ioman) {
- LOG(error) << "CbmRecoQaTask::Init :: RootManager not instantiated!";
+ LOG(fatal) << GetName() << "::Init :: RootManager not instantiated!";
return kERROR;
}
- fGlobalTracks = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrack"));
- if (!fGlobalTracks) {
- LOG(error) << "CbmRecoQaTask::Init: Global track array not found!";
+ fTracks = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrack"));
+ if (!fTracks) {
+ LOG(fatal) << GetName() << "::Init: Global track array not found!";
return kERROR;
}
- fGlobalTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrackMatch"));
- if (!fGlobalTrackMatches) {
- LOG(warn) << "CbmRecoQaTask::Init: MC info for Global track not available !";
+ fTimeSlice = static_cast<CbmTimeSlice*>(ioman->GetObject("TimeSlice."));
+ if (!fTimeSlice) {
+ LOG(warn) << GetName() << "::Init: No time slice object. Checking for events";
+ fEvents = static_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
+ if (!fEvents) LOG(warn) << GetName() << "::Init: No event found. Some results will be missing.";
+ }
+
+ fTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrackMatch"));
+ if (!fTrackMatches) {
+ LOG(warn) << GetName() << "::Init: MC info for Global track not available !";
+ }
+
+ LOG(info) << GetName() << "::Init: Registered projections " << fProjView.size();
+ for (auto zprj : fProjView)
+ LOG(debug) << GetName() << "::Init: Project trks @ " << zprj.first;
+
+ int itrkLy(0);
+ TFolder *modDir(nullptr), *smDir(nullptr);
+ fOutFolder.Clear();
+ switch (fSetupClass) {
+ case kMcbm22:
+ LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2022\".";
+ modDir = fOutFolder.AddFolder("STS", "STS reziduals to CA");
+ // STS U0 - station 0
+ fModView[itrkLy].push_back(ModuleView());
+ fModView[itrkLy][0].hdXx.first = 10;
+ fModView[itrkLy][0].hdXx.second =
+ new TH2D(Form("hxx%d", itrkLy), Form("STS U0 x-dx; X_{STS U0} (cm); #Delta X (mm)"), 160, -8, 8, 400, -2, 2);
+ modDir->Add(fModView[itrkLy][0].hdXx.second);
+ fModView[itrkLy][0].hdYy.first = 10;
+ fModView[itrkLy][0].hdYy.second = new TH2D(
+ Form("hyy%d", itrkLy), Form("STS U0 y-dy; Y_{STS U0} (cm); #Delta Y (mm)"), 160, -8, 8, 350, -3.5, 3.5);
+ modDir->Add(fModView[itrkLy][0].hdYy.second);
+ fModView[itrkLy][0].hdT.first = 1;
+ fModView[itrkLy][0].hdT.second =
+ new TH2D(Form("ht%d", itrkLy), Form("STS U0 x-dt; X_{STS U0} (cm); #Delta T (ns)"), 160, -8, 8, 15, -60, 60);
+ modDir->Add(fModView[itrkLy][0].hdT.second);
+ fModView[itrkLy][0].hpullXphi =
+ new TH2D(Form("hpullx%d", itrkLy), Form("STS U0 pullx-phi; dx/dz_{trk}; pull X"), 160, -0.8, 0.8, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullXphi);
+ fModView[itrkLy][0].hpullYtht =
+ new TH2D(Form("hpully%d", itrkLy), Form("STS U0 pully-tht; dy/dz_{trk}; pull Y"), 160, -0.8, 0.8, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullYtht);
+ itrkLy++;
+ // STS U1 - station 1
+ fModView[itrkLy].push_back(ModuleView());
+ fModView[itrkLy][0].hdXx.first = 10;
+ fModView[itrkLy][0].hdXx.second =
+ new TH2D(Form("hxx%d", itrkLy), Form("STS U1 x-dx; X_{STS U1} (cm); #Delta X (mm)"), 200, -10, 10, 400, -2, 2);
+ modDir->Add(fModView[itrkLy][0].hdXx.second);
+ fModView[itrkLy][0].hdYy.first = 10;
+ fModView[itrkLy][0].hdYy.second = new TH2D(
+ Form("hyy%d", itrkLy), Form("STS U1 y-dy; Y_{STS U1} (cm); #Delta Y (mm)"), 200, -10, 10, 350, -3.5, 3.5);
+ modDir->Add(fModView[itrkLy][0].hdYy.second);
+ fModView[itrkLy][0].hdT.first = 1;
+ fModView[itrkLy][0].hdT.second =
+ new TH2D(Form("ht%d", itrkLy), Form("STS U1 x-dt; X_{STS U1} (cm); #Delta T (ns)"), 200, -10, 10, 15, -60, 60);
+ modDir->Add(fModView[itrkLy][0].hdT.second);
+ fModView[itrkLy][0].hpullXphi =
+ new TH2D(Form("hpullx%d", itrkLy), Form("STS U1 pullx-phi; dx/dz_{trk}; pull X"), 160, -0.8, 0.8, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullXphi);
+ fModView[itrkLy][0].hpullYtht =
+ new TH2D(Form("hpully%d", itrkLy), Form("STS U1 pully-tht; dy/dz_{trk}; pull Y"), 160, -0.8, 0.8, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullYtht);
+ itrkLy++;
+ // ===============================================================================
+ modDir = fOutFolder.AddFolder("TRD", "TRD reziduals to CA");
+ // TRD 2D - station 2
+ fModView[itrkLy].push_back(ModuleView());
+ fModView[itrkLy][0].hdXx.first = 10;
+ fModView[itrkLy][0].hdXx.second =
+ new TH2D(Form("hxx%d", itrkLy), Form("TRD2D x-dx; X_{TRD2D} (cm); #Delta X (mm)"), 500, -25, 25, 240, -6, 6);
+ modDir->Add(fModView[itrkLy][0].hdXx.second);
+ fModView[itrkLy][0].hdYy.first = 1;
+ fModView[itrkLy][0].hdYy.second =
+ new TH2D(Form("hyy%d", itrkLy), Form("TRD2D y-dy; Y_{TRD2D} (cm); #Delta Y (cm)"), 300, -30, 30, 400, -2, 2);
+ modDir->Add(fModView[itrkLy][0].hdYy.second);
+ fModView[itrkLy][0].hdT.first = 1;
+ fModView[itrkLy][0].hdT.second =
+ new TH2D(Form("ht%d", itrkLy), Form("TRD2D x-dt; X_{TRD2D} (cm); #Delta T (ns)"), 500, -25, 25, 25, -100, 100);
+ modDir->Add(fModView[itrkLy][0].hdT.second);
+ fModView[itrkLy][0].hpullXphi =
+ new TH2D(Form("hpullx%d", itrkLy), Form("TRD2D pullx-phi; dx/dz_{trk}; pull X"), 160, -0.8, 0.8, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullXphi);
+ fModView[itrkLy][0].hpullYtht =
+ new TH2D(Form("hpully%d", itrkLy), Form("TRD2D pully-tht; dy/dz_{trk}; pull Y"), 160, -0.8, 0.8, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullYtht);
+ itrkLy++;
+ // TRD 1Dx - station 3
+ fModView[itrkLy].push_back(ModuleView());
+ fModView[itrkLy][0].hdXx.first = 10;
+ fModView[itrkLy][0].hdXx.second = new TH2D(
+ Form("hxx%d", itrkLy), Form("TRD1Dx x-dx; X_{TRD1Dx} (cm); #Delta X (mm)"), 450, -45, 45, 200, -10, 10);
+ modDir->Add(fModView[itrkLy][0].hdXx.second);
+ fModView[itrkLy][0].hdYy.first = 1.;
+ fModView[itrkLy][0].hdYy.second =
+ new TH2D(Form("hyy%d", itrkLy), Form("TRD1Dx y-dy; Y_{TRD1Dx} (cm); #Delta Y (cm)"), 12, -45, 45, 400, -20, 20);
+ modDir->Add(fModView[itrkLy][0].hdYy.second);
+ fModView[itrkLy][0].hdT.first = 1;
+ fModView[itrkLy][0].hdT.second = new TH2D(
+ Form("ht%d", itrkLy), Form("TRD1Dx x-dt; X_{TRD1Dx} (cm); #Delta T (ns)"), 450, -45, 45, 25, -100, 100);
+ modDir->Add(fModView[itrkLy][0].hdT.second);
+ fModView[itrkLy][0].hpullXphi =
+ new TH2D(Form("hpullx%d", itrkLy), Form("TRD1Dx pullx-phi; dx/dz_{trk}; pull X"), 100, -0.5, 0.5, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullXphi);
+ fModView[itrkLy][0].hpullYtht =
+ new TH2D(Form("hpully%d", itrkLy), Form("TRD1Dx pully-tht; dy/dz_{trk}; pull Y"), 100, -0.5, 0.5, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullYtht);
+ itrkLy++;
+ // TRD 1Dy - station 4
+ fModView[itrkLy].push_back(ModuleView());
+ fModView[itrkLy][0].hdXx.first = 1;
+ fModView[itrkLy][0].hdXx.second =
+ new TH2D(Form("hxx%d", itrkLy), Form("TRD1Dy x-dx; X_{TRD1Dy} (cm); #Delta X (cm)"), 12, -45, 45, 400, -20, 20);
+ modDir->Add(fModView[itrkLy][0].hdXx.second);
+ fModView[itrkLy][0].hdYy.first = 10.;
+ fModView[itrkLy][0].hdYy.second = new TH2D(
+ Form("hyy%d", itrkLy), Form("TRD1Dy y-dy; Y_{TRD1Dy} (cm); #Delta Y (mm)"), 500, -50, 50, 400, -20, 20);
+ modDir->Add(fModView[itrkLy][0].hdYy.second);
+ fModView[itrkLy][0].hdT.first = 1;
+ fModView[itrkLy][0].hdT.second = new TH2D(
+ Form("ht%d", itrkLy), Form("TRD1Dy y-dt; Y_{TRD1Dy} (cm); #Delta T (ns)"), 500, -50, 50, 25, -100, 100);
+ modDir->Add(fModView[itrkLy][0].hdT.second);
+ fModView[itrkLy][0].hpullXphi =
+ new TH2D(Form("hpullx%d", itrkLy), Form("TRD1Dy pullx-phi; dx/dz_{trk}; pull X"), 100, -0.5, 0.5, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullXphi);
+ fModView[itrkLy][0].hpullYtht =
+ new TH2D(Form("hpully%d", itrkLy), Form("TRD1Dy pully-tht; dy/dz_{trk}; pull Y"), 100, -0.5, 0.5, 100, -5, 5);
+ modDir->Add(fModView[itrkLy][0].hpullYtht);
+ itrkLy++;
+ // ===============================================================================
+ modDir = fOutFolder.AddFolder("ToF", "ToF reziduals to CA");
+ // ToF Sm0 & Sm3 (Typ0) - station 5
+ smDir = modDir->AddFolder("SM0_0", "ToF reziduals to CA for Sm0 and Sm3 type 0");
+ for (int irpc(0); irpc < 5; irpc++) {
+ fModView[itrkLy].push_back(ModuleView());
+ ModuleView* cdet = &fModView[itrkLy][irpc];
+ cdet->hdXx.first = 1;
+ cdet->hdXx.second =
+ new TH2D(Form("hxx%d%d", itrkLy, irpc), Form("ToF Sm0/0 Rpc%d x-dx; X_{ToF} (cm); #Delta X (cm)", irpc), 900,
+ -25, 65, 60, -3, 3);
+ smDir->Add(cdet->hdXx.second);
+ cdet->hdYy.first = 1;
+ cdet->hdYy.second =
+ new TH2D(Form("hyy%d%d", itrkLy, irpc), Form("ToF Sm0/0 Rpc%d y-dy; Y_{ToF} (cm); #Delta Y (cm)", irpc), 160,
+ -80, 80, 60, -3, 3);
+ smDir->Add(cdet->hdYy.second);
+ cdet->hdT.first = 1.;
+ cdet->hdT.second =
+ new TH2D(Form("ht%d%d", itrkLy, irpc), Form("ToF Sm0/0 Rpc%d x-dt; X_{ToF} (cm); #Delta T (ns)", irpc), 900,
+ -25, 65, 70, -20, 50);
+ smDir->Add(cdet->hdT.second);
+ cdet->hpullXphi =
+ new TH2D(Form("hpullx%d%d", itrkLy, irpc), Form("ToF Sm0/0 Rpc%d pullx-phi; dx/dz_{ToF}; pull X", irpc), 160,
+ -0.3, 0.5, 100, -5, 5);
+ smDir->Add(cdet->hpullXphi);
+ cdet->hpullYtht =
+ new TH2D(Form("hpully%d%d", itrkLy, irpc), Form("ToF Sm0/0 Rpc%d pully-tht; dy/dz_{ToF}; pull Y", irpc), 180,
+ -0.45, 0.45, 100, -5, 5);
+ smDir->Add(cdet->hpullYtht);
+ }
+ itrkLy++;
+
+ // ToF Sm1 (Typ0) - station 6
+ smDir = modDir->AddFolder("SM1_0", "ToF reziduals to CA for Sm1 type 0");
+ for (int irpc(0); irpc < 5; irpc++) {
+ fModView[itrkLy].push_back(ModuleView());
+ ModuleView* cdet = &fModView[itrkLy][irpc];
+ cdet->hdXx.first = 1;
+ cdet->hdXx.second =
+ new TH2D(Form("hxx%d%d", itrkLy, irpc), Form("ToF Sm1/0 Rpc%d x-dx; X_{ToF} (cm); #Delta X (cm)", irpc), 500,
+ -25, 25, 60, -3, 3);
+ smDir->Add(cdet->hdXx.second);
+ cdet->hdYy.first = 1;
+ cdet->hdYy.second =
+ new TH2D(Form("hyy%d%d", itrkLy, irpc), Form("ToF Sm1/0 Rpc%d y-dy; Y_{ToF} (cm); #Delta Y (cm)", irpc), 160,
+ -80, 80, 60, -3, 3);
+ smDir->Add(cdet->hdYy.second);
+ cdet->hdT.first = 1.;
+ cdet->hdT.second =
+ new TH2D(Form("ht%d%d", itrkLy, irpc), Form("ToF Sm1/0 Rpc%d x-dt; X_{ToF} (cm); #Delta T (ns)", irpc), 500,
+ -25, 25, 70, -20, 50);
+ smDir->Add(cdet->hdT.second);
+ cdet->hpullXphi =
+ new TH2D(Form("hpullx%d%d", itrkLy, irpc), Form("ToF Sm1/0 Rpc%d pullx-phi; dx/dz_{ToF}; pull X", irpc), 120,
+ -0.3, 0.3, 100, -5, 5);
+ smDir->Add(cdet->hpullXphi);
+ cdet->hpullYtht =
+ new TH2D(Form("hpully%d%d", itrkLy, irpc), Form("ToF Sm1/0 Rpc%d pully-tht; dy/dz_{ToF}; pull Y", irpc), 180,
+ -0.45, 0.45, 100, -5, 5);
+ smDir->Add(cdet->hpullYtht);
+ }
+ itrkLy++;
+
+ // ToF Sm0 (Typ2) - station 7
+ smDir = modDir->AddFolder("SM0_2", "ToF reziduals to CA for Sm0 type 2");
+ for (int irpc(0); irpc < 5; irpc++) {
+ fModView[itrkLy].push_back(ModuleView());
+ ModuleView* cdet = &fModView[itrkLy][irpc];
+ cdet->hdXx.first = 1;
+ cdet->hdXx.second =
+ new TH2D(Form("hxx%d%d", itrkLy, irpc), Form("ToF Sm0/2 Rpc%d x-dx; X_{ToF} (cm); #Delta X (cm)", irpc), 500,
+ -50, 50, 60, -3, 3);
+ smDir->Add(cdet->hdXx.second);
+ cdet->hdYy.first = 1;
+ cdet->hdYy.second =
+ new TH2D(Form("hyy%d%d", itrkLy, irpc), Form("ToF Sm0/2 Rpc%d y-dy; Y_{ToF} (cm); #Delta Y (cm)", irpc), 160,
+ -80, 80, 60, -3, 3);
+ smDir->Add(cdet->hdYy.second);
+ cdet->hdT.first = 1.;
+ cdet->hdT.second =
+ new TH2D(Form("ht%d%d", itrkLy, irpc), Form("ToF Sm0/2 Rpc%d x-dt; X_{ToF} (cm); #Delta T (ps)", irpc), 500,
+ -50, 50, 100, -500, 500);
+ smDir->Add(cdet->hdT.second);
+ cdet->hpullXphi =
+ new TH2D(Form("hpullx%d%d", itrkLy, irpc), Form("ToF Sm0/2 Rpc%d pullx-phi; dx/dz_{ToF}; pull X", irpc), 120,
+ -0.3, 0.3, 100, -5, 5);
+ smDir->Add(cdet->hpullXphi);
+ cdet->hpullYtht =
+ new TH2D(Form("hpully%d%d", itrkLy, irpc), Form("ToF Sm0/2 Rpc%d pully-tht; dy/dz_{ToF}; pull Y", irpc), 180,
+ -0.45, 0.45, 100, -5, 5);
+ smDir->Add(cdet->hpullYtht);
+ }
+ itrkLy++;
+
+ // ToF Sm2 (Typ0) - station 8
+ smDir = modDir->AddFolder("SM2_0", "ToF reziduals to CA for Sm2 type 0");
+ for (int irpc(0); irpc < 5; irpc++) {
+ fModView[itrkLy].push_back(ModuleView());
+ ModuleView* cdet = &fModView[itrkLy][irpc];
+ cdet->hdXx.first = 1;
+ cdet->hdXx.second =
+ new TH2D(Form("hxx%d%d", itrkLy, irpc), Form("ToF Sm2/0 Rpc%d x-dx; X_{ToF} (cm); #Delta X (cm)", irpc), 600,
+ -30, 30, 60, -3, 3);
+ smDir->Add(cdet->hdXx.second);
+ cdet->hdYy.first = 1;
+ cdet->hdYy.second =
+ new TH2D(Form("hyy%d%d", itrkLy, irpc), Form("ToF Sm2/0 Rpc%d y-dy; Y_{ToF} (cm); #Delta Y (cm)", irpc), 160,
+ -80, 80, 60, -3, 3);
+ smDir->Add(cdet->hdYy.second);
+ cdet->hdT.first = 1.;
+ cdet->hdT.second =
+ new TH2D(Form("ht%d%d", itrkLy, irpc), Form("ToF Sm2/0 Rpc%d x-dt; X_{ToF} (cm); #Delta T (ps)", irpc), 600,
+ -30, 30, 100, -500, 500);
+ smDir->Add(cdet->hdT.second);
+ cdet->hpullXphi =
+ new TH2D(Form("hpullx%d%d", itrkLy, irpc), Form("ToF Sm2/0 Rpc%d pullx-phi; dx/dz_{ToF}; pull X", irpc), 120,
+ -0.3, 0.3, 100, -5, 5);
+ smDir->Add(cdet->hpullXphi);
+ cdet->hpullYtht =
+ new TH2D(Form("hpully%d%d", itrkLy, irpc), Form("ToF Sm2/0 Rpc%d pully-tht; dy/dz_{ToF}; pull Y", irpc), 180,
+ -0.45, 0.45, 100, -5, 5);
+ smDir->Add(cdet->hpullYtht);
+ }
+ break;
+ case kMcbm24: LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2024\". Not implemented."; break;
+ case kDefault: LOG(info) << GetName() << "::Init: Setup config for \"DEFAULT\". Not implemented."; break;
}
return kSUCCESS;
@@ -52,60 +305,99 @@ InitStatus CbmRecoQaTask::Init()
//_____________________________________________________________________
void CbmRecoQaTask::Exec(Option_t*)
{
- for (int iTr = 0; iTr < fGlobalTracks->GetEntriesFast(); iTr++) {
-// if (static_cast<int>(fTracks.size()) >= fMaxNtracks) {
-// break;
-// }
- const CbmGlobalTrack* globalTrack = dynamic_cast<const CbmGlobalTrack*>(fGlobalTracks->At(iTr));
- if (!globalTrack) {
- LOG(fatal) << "RecoQA: null pointer to the global track!";
+ LOG(debug) << GetName() << "::Exec : Evs[" << (fEvents ? fEvents->GetEntriesFast() : 0) << "] Trks["
+ << fTracks->GetEntriesFast() << "]";
+ int itrack(0);
+ for (auto trackObj : *fTracks) {
+ //printf("GlobalTrack %d\n", itrack);
+ auto track = dynamic_cast<const CbmGlobalTrack*>(trackObj);
+ CbmKfTrackFitter::Track trkKf;
+ if (!fFitter.CreateGlobalTrack(trkKf, *track)) {
+ LOG(fatal) << GetName() << "::Exec: can not create the track for the fit! ";
break;
}
- CbmKfTrackFitter::Track trk;
- if (!fFitter.CreateGlobalTrack(trk, *globalTrack)) {
- LOG(fatal) << "RecoQA: can not create the global track for the fit! ";
- break;
- }
- trk.fNodes.insert(trk.fNodes.begin(), CbmKfTrackFitter::FitNode());
- trk.fMsQp0 = 1. / 0.5;
- trk.fIsMsQp0Set = true;
-
- // try first fit with all hits ON
- trk.MakeConsistent();
- fFitter.FitTrack(trk);
- CbmKfTrackFitter::FitNode &node = trk.fNodes.front();
- printf("val @ vx %f fit %f\n", node.fMxy.X()[0], node.fTrack.X()[0]);
-
-
- for (auto& n : trk.fNodes) {
+ trkKf.fMsQp0 = 1. / 0.5;
+ trkKf.fIsMsQp0Set = true;
+
+ for (auto& n : trkKf.fNodes) {
+ if (n.fHitSystemId == ECbmModuleId::kNotExist) continue;
+ //printf(" %d z[%6.2f]\n", n.fMaterialLayer, n.fZ);
+ //printf(" D[%d] hit[%d %d] Flags[%c %c %c %c]\n",
+ // (int)n.fHitSystemId, n.fHitAddress, n.fHitIndex,
+ // (n.fIsXySet ? 'y' : 'n'), (n.fIsTimeSet ? 'y' : 'n'), (n.fIsFitted ? 'y' : 'n'), (n.fIsRadThickSet ? 'y' : 'n'));
+ //printf(" hit x[%+6.2f] y[%+6.2f] t[%f]\n", (double)n.fMxy.X()[0], (double)n.fMxy.Y()[0], (double)n.fMt.T()[0]);
+
n.fIsTimeSet = false;
n.fIsXySet = false;
- trk.MakeConsistent();
- fFitter.FitTrack(trk);
- double dx = n.fMxy.X()[0] - n.fTrack.X()[0],
- dy = n.fMxy.Y()[0] - n.fTrack.Y()[0],
+ trkKf.MakeConsistent();
+ fFitter.FitTrack(trkKf);
+ double dx = n.fMxy.X()[0] - n.fTrack.X()[0], dy = n.fMxy.Y()[0] - n.fTrack.Y()[0],
+ dt = n.fMt.T()[0] - n.fTrack.Time()[0],
pullx = dx / sqrt(n.fMxy.Dx2()[0] + n.fTrack.GetCovariance(0, 0)[0]),
- pully = dy / sqrt(n.fMxy.Dy2()[0] + n.fTrack.GetCovariance(1, 1)[0]);
+ pully = dy / sqrt(n.fMxy.Dy2()[0] + n.fTrack.GetCovariance(1, 1)[0]),
+ pullt = dt / sqrt(n.fMt.Dt2()[0] + n.fTrack.GetCovariance(5, 5)[0]);
+ //printf(" res dx[%+6.4f] dy[%+6.4f] dt[%+6.4f]\n pull x[%+6.4f] y[%+6.4f] t[%+6.4f]\n", dx, dy, dt, pullx, pully, pullt);
+
+ // add granularity for specific systems
+ int modId = 0;
+ if (n.fHitSystemId == ECbmModuleId::kTof) modId = CbmTofAddress::GetRpcId(n.fHitAddress);
+
+ // check that we have correctly created all views
+ if (fModView.find(n.fMaterialLayer) == fModView.end()) continue;
+ if (fModView[n.fMaterialLayer].size() <= (size_t)modId) continue;
+ ModuleView* mod = &fModView[n.fMaterialLayer][modId];
+ if (!mod) {
+ LOG(warn) << GetName() << "::Exec: can not access mod view for " << modId;
+ continue;
+ }
+ if (mod->hdXx.second) mod->hdXx.second->Fill(n.fMxy.X()[0], mod->hdXx.first * dx);
+ if (mod->hdYy.second) mod->hdYy.second->Fill(n.fMxy.Y()[0], mod->hdYy.first * dy);
+ if (mod->hdT.second) {
+ if (n.fHitSystemId == ECbmModuleId::kTrd
+ && CbmTrdAddress::GetModuleAddress(n.fHitAddress)
+ == 37) { // TODO mark generically the low resolution coordinate of TRD
+ mod->hdT.second->Fill(n.fMxy.Y()[0], mod->hdT.first * dt);
+ }
+ else
+ mod->hdT.second->Fill(n.fMxy.X()[0], mod->hdT.first * dt);
+ }
+ if (mod->hpullXphi) mod->hpullXphi->Fill(n.fTrack.Tx()[0], pullx);
+ if (mod->hpullYtht) mod->hpullYtht->Fill(n.fTrack.Ty()[0], pully);
+ n.fIsTimeSet = true;
+ n.fIsXySet = true;
}
-// //if (t.fNstsHits < 1) continue;
-// //if (t.fNtrdHits < 2) continue;
-// if (t.fNstsHits + t.fNmuchHits + t.fNtrd1dHits + t.fNtrd2dHits + t.fNtofHits < fNtrackingStations - 1) continue;
-// //if (t.fNtrdHits < 3) continue;
-// if (t.fNstsHits < 2) continue;
-// if (t.fNtofHits < 2) continue;
-//
- fFitter.FitTrack(trk);
- fTracks.push_back(trk);
+ // Fit track with all hits ON for vx projection
+ int iprj(0);
+ for (auto zprj : fProjView) {
+ CbmKfTrackFitter::FitNode n = CbmKfTrackFitter::FitNode();
+ n.fZ = zprj.first;
+ n.fIsXySet = true;
+ n.fReference1 = iprj++;
+ trkKf.fNodes.push_back(n);
+ }
+ trkKf.MakeConsistent();
+ fFitter.FitTrack(trkKf);
+ for (auto& n : trkKf.fNodes) {
+ if (n.fReference1 < 0) continue;
+ //printf(" prj[%d] z[%+6.4f] x[%+6.4f] y[%+6.4f] t[%+6.4f]\n", n.fReference1, n.fZ, (double)n.fTrack.X()[0], (double)n.fTrack.Y()[0], (double)n.fTrack.Time()[0]);
+ }
+ itrack++;
}
}
//_____________________________________________________________________
-void CbmRecoQaTask::FinishEvent()
+void CbmRecoQaTask::SetProjections(std::vector<double> vzpj)
{
+ fProjView.clear();
+ for (auto zprj : vzpj) {
+ fProjView[zprj] = new TH2D(Form("hxx"), Form("; X (cm); Y (cm)"), 450, -25, 20, 450, -25, 20);
+ }
}
//_____________________________________________________________________
-void CbmRecoQaTask::FinishTask()
+void CbmRecoQaTask::Finish()
{
+ FairSink* sink = FairRootManager::Instance()->GetSink();
+ sink->WriteObject(&fOutFolder, nullptr);
}
diff --git a/reco/qa/CbmRecoQaTask.h b/reco/qa/CbmRecoQaTask.h
index 80fa7dbc6247a567ba4eab04e4ef9f507bbe8b7e..156088158fc38da5371fe2c4017b688caf6cf121 100644
--- a/reco/qa/CbmRecoQaTask.h
+++ b/reco/qa/CbmRecoQaTask.h
@@ -5,34 +5,64 @@
#ifndef CBMRECOQATASK_H
#define CBMRECOQATASK_H 1
+#include "CbmKfTrackFitter.h"
+// #include "CbmQaHist.h"
+
#include "FairTask.h"
+#include "TFolder.h"
-#include "CbmKfTrackFitter.h"
+#include <map>
#include <vector>
+class CbmTimeSlice;
class TClonesArray;
+class TH2D;
class CbmRecoQaTask : public FairTask {
-public:
+ public:
+ enum Setup
+ {
+ kMcbm22 = 0,
+ kMcbm24,
+ kDefault
+ };
+
+ public:
CbmRecoQaTask();
- ~CbmRecoQaTask();
+ virtual ~CbmRecoQaTask() { ; }
virtual InitStatus Init();
/** \brief Executed task **/
virtual void Exec(Option_t* option);
- virtual void FinishEvent();
- virtual void FinishTask();
-
-private:
+ virtual void Finish();
+
+ /** \brief Define the set of extra z positions where the track should be projected in the x-y plane */
+ void SetProjections(std::vector<double> vzpj);
+ void SetSetupClass(CbmRecoQaTask::Setup setup) { fSetupClass = setup; }
+
+ private:
CbmRecoQaTask(const CbmRecoQaTask&);
CbmRecoQaTask& operator=(const CbmRecoQaTask&);
+ struct ModuleView {
+ std::pair<float, TH2D*> hdXx = std::make_pair(1, nullptr);
+ std::pair<float, TH2D*> hdYy = std::make_pair(1, nullptr);
+ std::pair<float, TH2D*> hdT = std::make_pair(1, nullptr);
+ TH2D* hpullXphi = nullptr;
+ TH2D* hpullYtht = nullptr;
+ };
CbmKfTrackFitter fFitter;
- TClonesArray* fGlobalTracks = nullptr; //! reconstructed global tracks / event
- TClonesArray* fGlobalTrackMatches = nullptr; //! MC info for the global tracks
+ TClonesArray* fTracks = nullptr; //! reconstructed global tracks / event
+ TClonesArray* fTrackMatches = nullptr; //! MC info for the global tracks
+ TClonesArray* fEvents = nullptr; //! reconstructed events
+ CbmTimeSlice* fTimeSlice = nullptr; //! Time slice info
+ //
+ TFolder fOutFolder = {"RecoQA", "CA track driven reco QA"}; //!
+ Setup fSetupClass = Setup::kMcbm22;
+ std::map<int, std::vector<ModuleView>> fModView = {}; //!
+ std::map<double, TH2D*> fProjView = {}; //!
- std::vector<CbmKfTrackFitter::Track> fTracks;
- ClassDef(CbmRecoQaTask, 1);
+ ClassDef(CbmRecoQaTask, 1); // Reconstruction QA analyzed from CA perspective
};
#endif // CBMRECOQATASK_H
diff --git a/reco/qa/RecoQaLinkDef.h b/reco/qa/RecoQaLinkDef.h
index 12944fe046571192684836df6e8eaeb88ef1511b..ecbf6054cf0f27e7c3f030e6bf13845e335dd5b6 100644
--- a/reco/qa/RecoQaLinkDef.h
+++ b/reco/qa/RecoQaLinkDef.h
@@ -10,8 +10,8 @@
#pragma link off all functions;
-#pragma link C++ class CbmRecoQaTask + ;
#pragma link C++ class CbmRecoQa + ;
#pragma link C++ class CbmTrackingTrdQa + ;
+#pragma link C++ class CbmRecoQaTask + ;
#endif