diff --git a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx
index 219f32f1978cf30036dd8ac44faa59c473b24a44..025d21f88f2d55d459c550e52c0714c070f0751d 100644
--- a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx
+++ b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx
@@ -1,9 +1,10 @@
-/* Copyright (C) 2020-2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
+/* Copyright (C) 2020-2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Pierre-Alain Loizeau [committer], Andreas Redelbach */
+ Authors: Pierre-Alain Loizeau [committer], Andreas Redelbach, Alexandru Bercuci */
#include "CbmMcbmCheckTimingAlgo.h"
+#include "CbmBmonDigi.h"
#include "CbmDigiManager.h"
#include "CbmFlesHistosTools.h"
#include "CbmMuchBeamTimeDigi.h"
@@ -13,16 +14,16 @@
#include "CbmTofDigi.h"
#include "CbmTrdDigi.h"
-#include "FairRootManager.h"
-#include "FairRunOnline.h"
+#include <FairRootManager.h>
+#include <FairRunOnline.h>
#include <Logger.h>
-#include "TF1.h"
-#include "TH1.h"
-#include "TH2.h"
-#include "THttpServer.h"
#include <TDirectory.h>
+#include <TF1.h>
#include <TFile.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <THttpServer.h>
#include <iomanip>
#include <iostream>
@@ -82,7 +83,7 @@ void CbmMcbmCheckTimingAlgo::CheckDataPresence(CheckTimingDetector detToCheck)
/// Handle special case for Bmon as not yet supported in DigiManager
if (ECbmModuleId::kBmon == detToCheck.detId) {
// Get a pointer to the previous already existing data level
- fpBmonDigiVec = ioman->InitObjectAs<std::vector<CbmTofDigi> const*>("BmonDigi");
+ fpBmonDigiVec = ioman->InitObjectAs<std::vector<CbmBmonDigi> const*>("BmonDigi");
if (!fpBmonDigiVec) {
LOG(fatal) << "No storage with Bmon digis found while it should be used. "
"Stopping there!";
@@ -210,7 +211,7 @@ void CbmMcbmCheckTimingAlgo::ProcessTs()
break;
} // case ECbmModuleId::kPsd:
case ECbmModuleId::kBmon: {
- CheckInterSystemOffset<CbmTofDigi>();
+ CheckInterSystemOffset<CbmBmonDigi>();
break;
} // case ECbmModuleId::kBmon:
default: {
@@ -317,7 +318,7 @@ void CbmMcbmCheckTimingAlgo::CheckInterSystemOffset()
break;
} // case ECbmModuleId::kPsd:
case ECbmModuleId::kBmon: {
- FillTimeOffsetHistos<CbmTofDigi>(dRefTime, dRefCharge, uDetIdx);
+ FillTimeOffsetHistos<CbmBmonDigi>(dRefTime, dRefCharge, uDetIdx);
break;
} // case ECbmModuleId::kBmon:
default: {
@@ -427,7 +428,7 @@ void CbmMcbmCheckTimingAlgo::FillTimeOffsetHistos(const Double_t dRefTime, const
fvDets[uDetIdx].iPrevRefFirstDigi = uFirstDigiInWin;
}
-// ---- GetDigiTime ------------------------------------------------------
+// ---- GetDigiInfo ------------------------------------------------------
template<class Digi>
uint CbmMcbmCheckTimingAlgo::GetDigiInfo(UInt_t uDigi, std::vector<std::tuple<double, double, uint>>* vec, ECbmModuleId)
{
@@ -437,22 +438,21 @@ uint CbmMcbmCheckTimingAlgo::GetDigiInfo(UInt_t uDigi, std::vector<std::tuple<do
vec->push_back(std::make_tuple(digi->GetTime(), digi->GetCharge(), digi->GetAddress()));
return 1;
}
+// --- STS specific digi info ---------------------
template<>
-uint CbmMcbmCheckTimingAlgo::GetDigiInfo<CbmTofDigi>(UInt_t uDigi, std::vector<std::tuple<double, double, uint>>* vec,
- ECbmModuleId detId)
+uint CbmMcbmCheckTimingAlgo::GetDigiInfo<CbmStsDigi>(UInt_t uDigi, std::vector<std::tuple<double, double, uint>>* vec,
+ ECbmModuleId)
{
- /** Template specialization for ToF in order to distinguish Bmon digis.
- */
+ const CbmStsDigi* digi = fDigiMan->Get<CbmStsDigi>(uDigi);
vec->clear();
- const CbmTofDigi* pDigi(nullptr);
- if (detId == ECbmModuleId::kBmon) pDigi = &fpBmonDigiVec->at(uDigi);
- else
- pDigi = fDigiMan->Get<CbmTofDigi>(uDigi);
-
- if (pDigi == nullptr) return 0;
- vec->push_back(std::make_tuple(pDigi->GetTime(), pDigi->GetCharge(), pDigi->GetAddress()));
+ if (digi == nullptr) return 0;
+ int32_t add = digi->GetAddress(), uId(CbmStsAddress::GetElementId(add, EStsElementLevel::kStsUnit)),
+ lId(CbmStsAddress::GetElementId(add, EStsElementLevel::kStsLadder)),
+ mId(CbmStsAddress::GetElementId(add, EStsElementLevel::kStsModule));
+ vec->push_back(std::make_tuple(digi->GetTime(), digi->GetCharge(), uId * 10 + lId * 3 + mId));
return 1;
}
+// --- TRD(1D,2D) specific digi info ---------------------
template<>
uint CbmMcbmCheckTimingAlgo::GetDigiInfo<CbmTrdDigi>(UInt_t uDigi, std::vector<std::tuple<double, double, uint>>* vec,
ECbmModuleId)
@@ -491,29 +491,40 @@ uint CbmMcbmCheckTimingAlgo::GetDigiInfo<CbmTrdDigi>(UInt_t uDigi, std::vector<s
}
return n;
}
+// --- ToF specific digi info ---------------------
+template<>
+uint CbmMcbmCheckTimingAlgo::GetDigiInfo<CbmTofDigi>(UInt_t uDigi, std::vector<std::tuple<double, double, uint>>* vec,
+ ECbmModuleId)
+{
+ const CbmTofDigi* digi = fDigiMan->Get<CbmTofDigi>(uDigi);
+ vec->clear();
+ if (digi == nullptr) return 0;
+ int32_t add = digi->GetAddress(), smID(CbmTofAddress::GetSmId(add)), smTyp(CbmTofAddress::GetSmType(add)),
+ rpcID(CbmTofAddress::GetRpcId(add));
+ // chID(CbmTofAddress::GetChannelId(add))
+ // chSd(CbmTofAddress::GetChannelSide(add));
+ vec->push_back(std::make_tuple(digi->GetTime(), digi->GetCharge(), smTyp * 100 + smID * 10 + rpcID));
+ return 1;
+}
// ---- GetViewId ------------------------------------------------------
template<class Digi>
-int CbmMcbmCheckTimingAlgo::GetViewId(CheckTimingDetector, std::tuple<double, double, uint>)
+int CbmMcbmCheckTimingAlgo::GetViewId(CheckTimingDetector det, std::tuple<double, double, uint> info)
{
- return 0;
-}
-template<>
-int CbmMcbmCheckTimingAlgo::GetViewId<CbmTrdDigi>(CheckTimingDetector det, std::tuple<double, double, uint> info)
-{
- uint moduleId = std::get<2>(info);
+ if (det.vName.size() == 1) return 0;
+
+ uint modId = std::get<2>(info);
uint iview = 0;
for (auto view : det.vName) {
- if (view.compare(std::to_string(moduleId)) == 0) return iview;
+ if (view.compare(std::to_string(modId)) == 0) return iview;
iview++;
}
- std::string sFullModId = det.sName + " module " + std::to_string(moduleId);
- if (0 == fUnimplementedTrdViewWarned.count(sFullModId)) {
- LOG(warning) << "Trd condition not implemented for " << sFullModId << ". Skipping it from now on.";
- fUnimplementedTrdViewWarned.insert(sFullModId);
+ std::string sFullId = det.sName + " mod " + std::to_string(modId);
+ if (0 == fUnimplementedView[det.detId].count(sFullId)) {
+ LOG(warning) << det.sName << " condition not implemented for " << sFullId << ". Skipping it from now on.";
+ fUnimplementedView[det.detId].insert(sFullId);
}
-
return -1;
}
diff --git a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h
index 13e4c0b52a934f20c590f78b221418564439840b..7aebc36b614e9acbf45388a2248c01902101c66b 100644
--- a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h
+++ b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h
@@ -1,6 +1,6 @@
-/* Copyright (C) 2020-2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
+/* Copyright (C) 2020-2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Pierre-Alain Loizeau [committer], Andreas Redelbach */
+ Authors: Pierre-Alain Loizeau [committer], Andreas Redelbach, Alexandru Bercuci */
#ifndef CBMMCBMCHECKTIMINGALGO_H
#define CBMMCBMCHECKTIMINGALGO_H
@@ -24,6 +24,7 @@ class CbmTrdDigi;
class CbmTofDigi;
class CbmRichDigi;
class CbmPsdDigi;
+class CbmBmonDigi;
class CheckTimingDetector {
public:
CheckTimingDetector() { ; }
@@ -108,7 +109,8 @@ private:
void CheckInterSystemOffset();
template<class Digi>
void FillTimeOffsetHistos(const Double_t dRefTime, const Double_t dRefCharge, UInt_t uDetIdx);
- /** @brief Retrieve digi (time,charge,addres) info. Use template specialization for special cases (e.g. Bmon, TRD, PSD)
+ /** @brief Retrieve digi (time,charge,addres) info.
+ * SHOULD BE IMPLEMENTED BY DETECTORS IF MORE DIFFERENTIAL STUDIES ARE NEEDED
* @param iDigi digi index in the digi array
* @param vec on return contains the signal(s), time(s) and address pairs. Should be allocated by the user
* @param detId if needed spec
@@ -117,9 +119,9 @@ private:
template<class Digi>
UInt_t GetDigiInfo(UInt_t iDigi, std::vector<std::tuple<double, double, UInt_t>>* vec,
ECbmModuleId detId = ECbmModuleId::kNotExist);
- /** @brief Retrieve the detector view corresponding to the digi data (@see CheckTimingDetector::vName)
+ /** @brief Retrieve the detector view corresponding to the digi data (@see CheckTimingDetector::vName).
* @param det detector definitions
- * @param info tuple of digi info (time, charge address)
+ * @param info tuple of digi info (time, charge, address or more specific info)
* @return the view index for the curent data or -1 if there is none
*/
template<class Digi>
@@ -129,7 +131,7 @@ private:
CbmDigiManager* fDigiMan = nullptr; //!
/** Bmon is not included in CbmDigiManager, so add it explicitly here **/
- const std::vector<CbmTofDigi>* fpBmonDigiVec = nullptr; //!
+ const std::vector<CbmBmonDigi>* fpBmonDigiVec = nullptr; //!
/** @brief Pointer to the Timeslice start time used to write it to the output tree
@remark since we hand this to the FairRootManager it also wants to delete it and we do not have to take care of deletion
@@ -164,7 +166,7 @@ private:
Double_t fRichPeakWidthNs = 40.;
Double_t fPsdPeakWidthNs = 20.;
- std::unordered_set<std::string> fUnimplementedTrdViewWarned = {};
+ std::map<ECbmModuleId, std::unordered_set<std::string>> fUnimplementedView = {};
ClassDefNV(CbmMcbmCheckTimingAlgo, 1);
};
diff --git a/macro/beamtime/mcbm2024/CMakeLists.txt b/macro/beamtime/mcbm2024/CMakeLists.txt
index 56a2a539ab95c2c89735c29663c9eae51a8f17bd..a24c93a4f6d220a94678994f0e3742002824b02d 100644
--- a/macro/beamtime/mcbm2024/CMakeLists.txt
+++ b/macro/beamtime/mcbm2024/CMakeLists.txt
@@ -8,19 +8,19 @@ Install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
PATTERN "*.sh"
)
# # RICH calibration file, resolve symlink to get the full file
-# Get_filename_component(_resolvedRichFile ${CMAKE_CURRENT_SOURCE_DIR}/icd_offset_it_0.data REALPATH)
+# get_filename_component(_resolvedRichFile ${CMAKE_CURRENT_SOURCE_DIR}/icd_offset_it_0.data REALPATH)
# Install(FILES ${_resolvedRichFile}
# DESTINATION share/cbmroot/macro/beamtime/mcbm2022
# )
-#
-# # SLURM scripts, bash scripts
-# Install(DIRECTORY online
-# DESTINATION share/cbmroot/macro/beamtime/mcbm2022
-# FILES_MATCHING PATTERN "*.sbatch"
-# PATTERN "*.sh"
-# )
-#
-# # Just the empty folder for output
-# Install(DIRECTORY data
-# DESTINATION share/cbmroot/macro/beamtime/mcbm2022
-# PATTERN "*" EXCLUDE)
+
+# SLURM scripts, bash scripts
+Install(DIRECTORY online
+ DESTINATION share/cbmroot/macro/beamtime/mcbm2024
+ FILES_MATCHING PATTERN "*.sbatch"
+ PATTERN "*.sh"
+ )
+
+# Just the empty folder for output
+Install(DIRECTORY data
+ DESTINATION share/cbmroot/macro/beamtime/mcbm2024
+ PATTERN "*" EXCLUDE)
diff --git a/macro/beamtime/mcbm2024/check_timing_any.C b/macro/beamtime/mcbm2024/check_timing_any.C
new file mode 100644
index 0000000000000000000000000000000000000000..f2d458bf9b96979c695e7a37171662cd572c0de6
--- /dev/null
+++ b/macro/beamtime/mcbm2024/check_timing_any.C
@@ -0,0 +1,101 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau [committer], Alexandru Bercuci*/
+
+void check_timing_any(TString fileName, UInt_t uRunId = 0, Int_t nEvents = 0, TString outDir = "data/")
+{
+
+ // ========================================================================
+ // Adjust this part according to your requirements
+
+ // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
+ Int_t iVerbose = 3;
+
+ // MC file
+
+ TString srcDir = gSystem->Getenv("VMCWORKDIR");
+ TString outFileName = fileName(0, fileName.Index(".digi.root"));
+
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // ------------------------------------------------------------------------
+
+ // ----- Analysis run --------------------------------------------------
+ FairRunOnline* fRun = new FairRunOnline();
+ fRun->ActivateHttpServer(100, 8080); // refresh each 100 events
+ fRun->SetSink(new FairRootFileSink(Form("%s/%s.sink.root", outDir.Data(), outFileName.Data())));
+ FairFileSource* inputSource = new FairFileSource(fileName);
+ fRun->SetSource(inputSource);
+ // ------------------------------------------------------------------------
+
+ CbmMcbmCheckTimingTask* timeChecker = new CbmMcbmCheckTimingTask();
+ /// Default is using Bmon as reference
+ timeChecker->SetReferenceDetector(ECbmModuleId::kBmon, "Bmon", -1000., 1000., 320., 182, 190);
+ /// Remove detectors not present in 03-2024
+ timeChecker->RemoveCheckDetector(ECbmModuleId::kMuch);
+ timeChecker->RemoveCheckDetector(ECbmModuleId::kPsd);
+ /// Add detectors with wider range
+ timeChecker->AddCheckDetector(ECbmModuleId::kSts, "Sts", -100., 150., 2500);
+ timeChecker->AddCheckDetector(ECbmModuleId::kRich, "Rich", -40., 40., 800);
+ timeChecker->AddCheckDetector(ECbmModuleId::kTrd, "Trd", -300., 300., 150);
+ timeChecker->AddCheckDetector(ECbmModuleId::kTof, "Tof", -40, 40, 800);
+ /// Add extra differential analysis on specific detectors
+ /// The modifications have to be accompanied by a parallel implementation of
+ /// the template function CbmMcbmCheckTimingAlgo::GetDigiInfo<CbmDetDigi>
+ /// TODO Should be implemented as function pointer defined in this macro
+ timeChecker->SetDetectorDifferential(ECbmModuleId::kSts,
+ {"0", "1", "3", "4", "10", "11", "13", "14", "16", "17",
+ "18"}); // select STS modules based on general expression U*10 + L*3 + M
+ timeChecker->SetDetectorDifferential(ECbmModuleId::kTrd, {"5", "21", "37"}); // select modules
+ timeChecker->SetDetectorDifferential(
+ ECbmModuleId::kTof,
+ {"0", "1", "2", "3", "4", "10", "11", "12", "13", "14", "20", "21", "22",
+ "23", "24", "30", "31", "32", "33", "34", "40", "41", "42", "43", "44", "200",
+ "201", "202", "203", "204", "600", "601", "900", "901", "910", "911"}); // select ToF RPCs based on the general expression Typ*100 + SmId * 10 + RpcId
+ if (0 < uRunId) timeChecker->SetOutFilename(Form("%s/%s.tck.root", outDir.Data(), outFileName.Data()));
+ fRun->AddTask(timeChecker);
+
+
+ // ----- Intialise and run --------------------------------------------
+ fRun->Init();
+ cout << "Starting run" << endl;
+ if (nEvents < 0) {
+ fRun->Run(0, 0); // run until end of input file
+ }
+ else {
+ fRun->Run(0, nEvents); // process N Events
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Finish -------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ cout << endl << endl;
+ cout << "Macro finished successfully." << endl;
+ cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
+ cout << endl;
+ // ------------------------------------------------------------------------
+
+ // 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();
+ cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+ cout << maxMemory;
+ cout << "</DartMeasurement>" << endl;
+
+ Float_t cpuUsage = ctime / rtime;
+ cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+ cout << cpuUsage;
+ cout << "</DartMeasurement>" << endl;
+ /*
+ FairMonitor* tempMon = FairMonitor::GetMonitor();
+ tempMon->Print();
+*/
+ // RemoveGeoManager();
+ cout << " Test passed" << endl;
+ cout << " All ok " << endl;
+}
diff --git a/macro/beamtime/mcbm2024/mcbm_event_reco_L1.C b/macro/beamtime/mcbm2024/mcbm_event_reco_L1.C
new file mode 100644
index 0000000000000000000000000000000000000000..69714abdcb84706b2d06011247433f635871dbfa
--- /dev/null
+++ b/macro/beamtime/mcbm2024/mcbm_event_reco_L1.C
@@ -0,0 +1,746 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau, Adrian Weber [committer], Alexandru Bercuci*/
+
+// --------------------------------------------------------------------------
+//
+// Macro for reconstruction of mcbm data (2024)
+// Combined reconstruction (event builder + cluster + hit finder + CA) for different subsystems.
+//
+// --------------------------------------------------------------------------
+
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+/// FIXME: Disable clang formatting to keep easy parameters overview
+/* clang-format off */
+Bool_t mcbm_event_reco_L1(UInt_t uRunId = 3105,
+ Int_t nTimeslices = 10,
+ TString sInpDir = "./data/",
+ TString sOutDir = "./data/",
+ Int_t iUnpFileIndex = -1,
+ Bool_t bMVD = kFALSE,
+ Bool_t bSTS = kTRUE,
+ Bool_t bTRD = kTRUE,
+ Bool_t bTRD2d = kTRUE,
+ Bool_t bRICH = kTRUE,
+ Bool_t bMUCH = kTRUE,
+ Bool_t bTOF = kTRUE,
+ Bool_t bTOFtr = kTRUE,
+ Bool_t bPSD = kFALSE,
+ Bool_t bAli = kTRUE,
+ Bool_t bEvB = kTRUE,
+ Bool_t bL1 = kFALSE,
+ Bool_t bQA = kFALSE,
+ Bool_t bFSD = kFALSE,
+ TString sInpFile = ""
+ )
+{
+ /// FIXME: Re-enable clang formatting after parameters initial values setting
+ /* clang-format on */
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "INFO"; //"INFO";
+ TString logVerbosity = "VERYLOW"; //"VERYLOW";
+ // ------------------------------------------------------------------------
+
+ // ----- Environment --------------------------------------------------
+ TString myName = "mcbm_event_reco"; // this macro's name for screen output
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+
+ // ----- In- and output file names ------------------------------------
+ /// Standardized RUN ID
+ TString sRunId = TString::Format("%04u", uRunId);
+
+ gSystem->Exec("mkdir " + sOutDir);
+ gSystem->Exec("cp $VMCWORKDIR/macro/run/.rootrc .");
+
+ TString outFile = sOutDir + "/" + sRunId;
+ /// Initial pattern
+ TString inFile = sInpDir + sRunId + ".digi";
+
+ /// Add index of splitting at unpacking level if needed
+ if (0 <= iUnpFileIndex) {
+ inFile += TString::Format("_%02u", iUnpFileIndex);
+ // the input par file is not split during unpacking!
+ outFile += TString::Format("_%02u", iUnpFileIndex);
+ } // if ( 0 <= uUnpFileIndex )
+ /// Add ROOT file suffix
+ inFile += ".root";
+ TString parFileOut = outFile + ".par.root";
+ TString geoFileOut = outFile + ".geo.root";
+ outFile += ".rec.root";
+
+ if ("" != sInpFile) {
+ inFile = sInpFile;
+ outFile = inFile;
+ outFile.ReplaceAll(".digi.root", "");
+ parFileOut = outFile + ".par.root";
+ geoFileOut = outFile + ".geo.root";
+ outFile += ".rec.root";
+ }
+ // ---------------------------------------------
+
+
+ // ----- TOF defaults ------------------------
+ // Your folder with the Tof Calibration files;
+ TString cFileId = sRunId + ".5.000";
+ TString TofFileFolder = "";
+
+ // ===> PAL 2022/11/04: overwriten by block around l.510!
+ Int_t calMode = 93;
+ Int_t calSel = 1;
+ Int_t calSm = 2;
+ Int_t RefSel = 11;
+ Double_t dDeadtime = 50.;
+ Int_t iSel2 = -1;
+
+ // Tracking
+ Int_t iSel = 22002;
+ Int_t iTrackingSetup = 1; // 2 for checking without beam counter;
+ Int_t iGenCor = 1;
+ Double_t dScalFac = 2.5;
+ Double_t dChi2Lim2 = 4.;
+ Bool_t bUseSigCalib = kFALSE;
+ Int_t iCalOpt = 110; // 0=do not call CbmTofCalibrator
+ Int_t iTrkPar = 0; // 4 for check without beam counter
+ Double_t dTOffScal = 1.;
+ // ------------------------------------------------------------------------
+
+ // ----- TOF Calibration Settings ---------------------------------------
+ TString cCalId = "490.100.5.0";
+ if (uRunId >= 759) cCalId = "759.100.4.0";
+ if (uRunId >= 812) cCalId = "831.100.4.0";
+ if (uRunId >= 1588) cCalId = "1588.50.6.0";
+ if (uRunId >= 2160) cCalId = "2160.50.4.0";
+ if (uRunId >= 2352) cCalId = "2391.5.000";
+ if (uRunId >= 2700) cCalId = "2912.1";
+
+ Int_t iCalSet = 30040500; // calibration settings
+ if (uRunId >= 759) iCalSet = 10020500;
+ if (uRunId >= 812) iCalSet = 10020500;
+ if (uRunId >= 1588) iCalSet = 12002002;
+ if (uRunId >= 2160) iCalSet = 700900500;
+ if (uRunId >= 2352) iCalSet = 22002500;
+ if (uRunId >= 2700) iCalSet = 12032500;
+
+ Double_t Tint = 100.; // coincidence time interval
+ Int_t iTrackMode = 2; // 2 for TofTracker
+ const Int_t iTofCluMode = 1;
+ // ------------------------------------------------------------------------
+
+ // --- Load the geometry setup ----
+ // This is currently only required by the TRD (parameters)
+ cbm::mcbm::ToForceLibLoad dummy; /// Needed to trigger loading of the library as no fct dict in ROOT6 and CLING
+ TString geoSetupTag = "";
+ try {
+ geoSetupTag = cbm::mcbm::GetSetupFromRunId(uRunId);
+ }
+ catch (const std::invalid_argument& e) {
+ std::cout << "Error in mapping from runID to setup name: " << e.what() << std::endl;
+ return kFALSE;
+ }
+
+ TString geoFile = sInpDir + geoSetupTag + ".geo.root";
+ CbmSetup* geoSetup = CbmSetup::Instance();
+ geoSetup->LoadSetup(geoSetupTag);
+
+ // You can modify the pre-defined setup by using
+ geoSetup->SetActive(ECbmModuleId::kMvd, bMVD);
+ geoSetup->SetActive(ECbmModuleId::kSts, bSTS);
+ geoSetup->SetActive(ECbmModuleId::kMuch, bMUCH);
+ geoSetup->SetActive(ECbmModuleId::kRich, bRICH);
+ geoSetup->SetActive(ECbmModuleId::kTrd, bTRD);
+ geoSetup->SetActive(ECbmModuleId::kTrd2d, bTRD2d);
+ geoSetup->SetActive(ECbmModuleId::kTof, bTOF);
+ geoSetup->SetActive(ECbmModuleId::kPsd, bPSD);
+ geoSetup->SetActive(ECbmModuleId::kFsd, bFSD);
+
+
+ //----- Load Parameters --------------------------------------------------
+ TList* parFileList = new TList();
+ TofFileFolder = Form("%s/%s", TofFileFolder.Data(), cCalId.Data());
+
+ // ----- TRD digitisation parameters -------------------------------------
+ TString geoTagTrd;
+ if (geoSetup->IsActive(ECbmModuleId::kTrd)) {
+ if (geoSetup->GetGeoTag(ECbmModuleId::kTrd, geoTagTrd)) {
+ TString paramFilesTrd(Form("%s/parameters/trd/trd_%s", srcDir.Data(), geoTagTrd.Data()));
+ std::vector<TString> paramFilesVecTrd = {"asic", "digi", "gas", "gain"};
+ for (auto parIt : paramFilesVecTrd) {
+ parFileList->Add(new TObjString(Form("%s.%s.par", paramFilesTrd.Data(), parIt.Data())));
+ }
+ }
+ for (auto parFileVecIt : *parFileList) {
+ LOG(debug) << Form("TrdParams - %s - added to parameter file list\n", parFileVecIt->GetName());
+ }
+ }
+ // ----- TOF digitisation parameters -------------------------------------
+ TString geoTagTof;
+ if (geoSetup->IsActive(ECbmModuleId::kTof)) {
+ geoSetup->GetGeoTag(ECbmModuleId::kTof, geoTagTof);
+ TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTagTof + ".digibdf.par");
+ parFileList->Add(tofBdfFile);
+ std::cout << "-I- " << myName << ": Using parameter file " << tofBdfFile->GetString() << std::endl;
+
+ // TFile* fgeo = new TFile(geoFile);
+ // TGeoManager* geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
+ // if (NULL == geoMan) {
+ // cout << "<E> FAIRGeom not found in geoFile " << geoFile.Data() << endl;
+ // return 1;
+ // }
+ }
+ // ------------------------------------------------------------------------
+
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // ------------------------------------------------------------------------
+
+
+ // ----- FairRunAna ---------------------------------------------------
+ FairRunAna* run = new FairRunAna();
+ FairFileSource* inputSource = new FairFileSource(inFile);
+ run->SetSource(inputSource);
+
+ FairRootFileSink* outputSink = new FairRootFileSink(outFile);
+ run->SetSink(outputSink);
+ run->SetGeomFile(geoFile);
+
+ // Define output file for FairMonitor histograms
+ TString monitorFile{outFile};
+ monitorFile.ReplaceAll(".rec.", ".rec.monitor.");
+ FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Logger settings ----------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ //FairLogger::GetLogger()->SetLogScreenLevel("DEBUG");
+ // ------------------------------------------------------------------------
+
+
+ // =========================================================================
+ // === Alignment Correction ===
+ // =========================================================================
+ // (Fairsoft Apr21p2 or newer is needed)
+ if (bAli) {
+
+ TString alignmentMatrixFileName = srcDir + "/parameters/mcbm/AlignmentMatrices_" + geoSetupTag + ".root";
+
+ if (alignmentMatrixFileName.Length() != 0) {
+
+ std::cout << "-I- " << myName << ": Applying alignment for file " << alignmentMatrixFileName << std::endl;
+
+ // Define the basic structure which needs to be filled with information
+ // This structure is stored in the output file and later passed to the
+ // FairRoot framework to do the (miss)alignment
+ std::map<std::string, TGeoHMatrix>* matrices{nullptr};
+
+ // read matrices from disk
+ LOG(info) << "Filename: " << alignmentMatrixFileName;
+ TFile* misalignmentMatrixRootfile = new TFile(alignmentMatrixFileName, "READ");
+ if (misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->GetObject("MisalignMatrices", matrices);
+ misalignmentMatrixRootfile->Close();
+ }
+ else {
+ LOG(error) << "Could not open file " << alignmentMatrixFileName << "\n Exiting";
+ exit(1);
+ }
+
+ if (matrices) {
+ run->AddAlignmentMatrices(*matrices);
+ }
+ else {
+ LOG(error) << "Alignment required but no matrices found."
+ << "\n Exiting";
+ exit(1);
+ }
+ }
+ }
+ // ------------------------------------------------------------------------
+
+
+ // --------------------event builder---------------------------------------
+ // ----- EventBuilder Settings----------------
+ // mCbm track trigger Tof, Bmon & STS (case 4 of mcbm_unp_event.C)
+ const Int_t eb_TriggerMinNumberBmon{1};
+ const Int_t eb_TriggerMaxNumberBmon{2};
+ const Int_t eb_TriggerMinNumberSts{2};
+ const Int_t eb_TriggerMinNumberStsLayers{1};
+ const Int_t eb_TriggerMinNumberMuch{1};
+ const Int_t eb_TriggerMinNumberTof{8};
+ const Int_t eb_TriggerMinNumberTofLayers{4};
+ const Int_t eb_TriggerMinNumberRich{0};
+
+ const Int_t eb_TrigWinMinBmon{-50};
+ const Int_t eb_TrigWinMaxBmon{50};
+ const Int_t eb_TrigWinMinSts{-60};
+ const Int_t eb_TrigWinMaxSts{60};
+ const Int_t eb_TrigWinMinMuch{-100};
+ const Int_t eb_TrigWinMaxMuch{100};
+ const Int_t eb_TrigWinMinTrd1d{-300};
+ const Int_t eb_TrigWinMaxTrd1d{300};
+ const Int_t eb_TrigWinMinTrd2d{-200};
+ const Int_t eb_TrigWinMaxTrd2d{200};
+ const Int_t eb_TrigWinMinTof{-20};
+ const Int_t eb_TrigWinMaxTof{60};
+ const Int_t eb_TrigWinMinRich{-60};
+ const Int_t eb_TrigWinMaxRich{60};
+
+ if (bEvB) {
+ CbmTaskBuildRawEvents* evBuildRaw = new CbmTaskBuildRawEvents();
+ evBuildRaw->SetVerbose(3);
+ //Choose between NoOverlap, MergeOverlap, AllowOverlap
+ evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::AllowOverlap);
+ //For time being it is needed. We will remove CbmMuchBeamTimeDigi.
+ evBuildRaw->ChangeMuchBeamtimeDigiFlag();
+ // Remove detectors where digis not found
+ if (!bRICH || !geoSetup->IsActive(ECbmModuleId::kRich)) evBuildRaw->RemoveDetector(kRawEventBuilderDetRich);
+ if (!bMUCH || !geoSetup->IsActive(ECbmModuleId::kMuch)) evBuildRaw->RemoveDetector(kRawEventBuilderDetMuch);
+ if (!bPSD || !geoSetup->IsActive(ECbmModuleId::kPsd)) evBuildRaw->RemoveDetector(kRawEventBuilderDetPsd);
+ if (!bTRD || !geoSetup->IsActive(ECbmModuleId::kTrd)) evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd);
+ //TODO temporary exclude 2D from event building (asked by M. Beyer)
+ if (!bTRD2d) evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd2D);
+ if (!bSTS || !geoSetup->IsActive(ECbmModuleId::kSts)) evBuildRaw->RemoveDetector(kRawEventBuilderDetSts);
+ if (!bTOF || !geoSetup->IsActive(ECbmModuleId::kTof)) evBuildRaw->RemoveDetector(kRawEventBuilderDetTof);
+ if (!bFSD || !geoSetup->IsActive(ECbmModuleId::kFsd)) evBuildRaw->RemoveDetector(kRawEventBuilderDetFsd);
+
+ // Set Bmon as reference detector.
+ // Select only Bmon1 [newly installed for mcbm2024 data] (AB)
+ evBuildRaw->SetReferenceDetector(kRawEventBuilderDetBmon, {0, 1});
+ // For making MuCh as seed detector
+ // evBuildRaw->SetReferenceDetector(kRawEventBuilderDetMuch);
+
+ // Use sliding window seed builder with STS
+ //evBuildRaw->SetReferenceDetector(kRawEventBuilderDetUndef);
+ //evBuildRaw->AddSeedTimeFillerToList(kRawEventBuilderDetTof);
+ //evBuildRaw->SetSlidingWindowSeedFinder(10, 10, 50);
+ //evBuildRaw->SetSeedFinderQa(true); // optional QA information for seed finder
+ //evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, 1000);
+ //evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1);
+ //evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, -500, 500);
+
+ // void SetTsParameters(double TsStartTime, double TsLength, double TsOverLength):
+ // => TsStartTime=0, TsLength=128 + 1.28 ms in 2022, TsOverLength=1.28 ms (1MS) in mCBM2022
+ evBuildRaw->SetTsParameters(0.0, 1.28e8, 1.28e6);
+
+ if (geoSetup->IsActive(ECbmModuleId::kTof)) {
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kTof, eb_TriggerMinNumberTof);
+ evBuildRaw->SetTriggerMinLayersNumber(ECbmModuleId::kTof, eb_TriggerMinNumberTofLayers);
+ evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kTof, -1);
+ }
+
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kBmon, eb_TriggerMinNumberBmon);
+ evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kBmon, eb_TriggerMaxNumberBmon);
+
+ if (geoSetup->IsActive(ECbmModuleId::kSts)) {
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, eb_TriggerMinNumberSts);
+ evBuildRaw->SetTriggerMinLayersNumber(ECbmModuleId::kSts, eb_TriggerMinNumberStsLayers);
+ evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1);
+ }
+ if (bMUCH && geoSetup->IsActive(ECbmModuleId::kMuch)) {
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kMuch, eb_TriggerMinNumberMuch);
+ evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kMuch, -1);
+ }
+ if (geoSetup->IsActive(ECbmModuleId::kRich)) {
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kRich, eb_TriggerMinNumberRich);
+ evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kRich, -1);
+ }
+
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kBmon, eb_TrigWinMinBmon, eb_TrigWinMaxBmon);
+ if (geoSetup->IsActive(ECbmModuleId::kTof))
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kTof, eb_TrigWinMinTof, eb_TrigWinMaxTof);
+ if (geoSetup->IsActive(ECbmModuleId::kSts))
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, eb_TrigWinMinSts, eb_TrigWinMaxSts);
+ if (bMUCH && geoSetup->IsActive(ECbmModuleId::kMuch))
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kMuch, eb_TrigWinMinMuch, eb_TrigWinMaxMuch);
+ if (geoSetup->IsActive(ECbmModuleId::kTrd))
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kTrd, eb_TrigWinMinTrd1d, eb_TrigWinMaxTrd1d);
+ if (geoSetup->IsActive(ECbmModuleId::kTrd))
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kTrd2d, eb_TrigWinMinTrd2d, eb_TrigWinMaxTrd2d);
+ if (geoSetup->IsActive(ECbmModuleId::kRich))
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kRich, eb_TrigWinMinRich, eb_TrigWinMaxRich);
+
+ run->AddTask(evBuildRaw);
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Reconstruction tasks -----------------------------------------
+
+
+ // =========================================================================
+ // === local STS Reconstruction ===
+ // =========================================================================
+
+ if (bSTS && geoSetup->IsActive(ECbmModuleId::kSts)) {
+ CbmRecoSts* recoSts = new CbmRecoSts();
+ if (bEvB) {
+ recoSts->SetMode(ECbmRecoMode::EventByEvent);
+ }
+ else {
+ recoSts->SetMode(ECbmRecoMode::Timeslice);
+ }
+ recoSts->SetVerbose(3);
+ recoSts->SetTimeCutDigisAbs(20.0); // cluster finder: time cut in ns
+ recoSts->SetTimeCutClustersAbs(20.0); // hit finder: time cut in ns
+
+ // Sensor params
+ CbmStsParSensor sensor6cm(CbmStsSensorClass::kDssdStereo);
+ sensor6cm.SetPar(0, 6.2092); // Extension in x
+ sensor6cm.SetPar(1, 6.2); // Extension in y
+ sensor6cm.SetPar(2, 0.03); // Extension in z
+ sensor6cm.SetPar(3, 5.9692); // Active size in y
+ sensor6cm.SetPar(4, 1024.); // Number of strips front side
+ sensor6cm.SetPar(5, 1024.); // Number of strips back side
+ sensor6cm.SetPar(6, 0.0058); // Strip pitch front side
+ sensor6cm.SetPar(7, 0.0058); // Strip pitch back side
+ sensor6cm.SetPar(8, 0.0); // Stereo angle front side
+ sensor6cm.SetPar(9, 7.5); // Stereo angle back side
+
+ CbmStsParSensor sensor12cm(sensor6cm); // copy all parameters, change then only the y size
+ sensor12cm.SetPar(1, 12.4); // Extension in y
+ sensor12cm.SetPar(3, 12.1692); // Active size in y
+
+ // --- Addresses for sensors
+ // --- They are defined in each station as sensor 1, module 1, halfladderD (2), ladder 1
+ // Int_t GetAddress(UInt_t unit = 0, UInt_t ladder = 0, UInt_t halfladder = 0, UInt_t module = 0, UInt_t sensor = 0,
+ // UInt_t side = 0, UInt_t version = kCurrentVersion);
+
+ // setup available starting with sts_v24b. Available if reconstructing from rra files (AB 14.05.2024)
+ Int_t stsAddress00 = CbmStsAddress::GetAddress(0, 0, 0, 0, 0, 0); // U0 L0 M0 6 cm
+ Int_t stsAddress01 = CbmStsAddress::GetAddress(1, 0, 1, 0, 0, 0); // U1 L0 M0 6 cm
+ Int_t stsAddress02 = CbmStsAddress::GetAddress(1, 0, 1, 1, 0, 0); // U1 L0 M1 6 cm
+ Int_t stsAddress03 = CbmStsAddress::GetAddress(1, 1, 1, 0, 0, 0); // U1 L1 M0 6 cm
+ Int_t stsAddress04 = CbmStsAddress::GetAddress(1, 1, 1, 1, 0, 0); // U1 L1 M1 6 cm
+ Int_t stsAddress05 = CbmStsAddress::GetAddress(2, 0, 1, 0, 0, 0); // U2 L0 M0 6 cm
+ Int_t stsAddress06 = CbmStsAddress::GetAddress(2, 0, 1, 1, 0, 0); // U2 L0 M1 12 cm
+ Int_t stsAddress07 = CbmStsAddress::GetAddress(2, 1, 1, 0, 0, 0); // U2 L1 M0 6 cm
+ Int_t stsAddress08 = CbmStsAddress::GetAddress(2, 1, 1, 1, 0, 0); // U2 L1 M1 12 cm
+ Int_t stsAddress09 = CbmStsAddress::GetAddress(2, 2, 1, 0, 0, 0); // U2 L2 M0 6 cm
+ Int_t stsAddress10 = CbmStsAddress::GetAddress(2, 2, 1, 1, 0, 0); // U2 L2 M1 6 cm
+ Int_t stsAddress11 = CbmStsAddress::GetAddress(2, 2, 1, 2, 0, 0); // U2 L2 M2 6 cm
+
+ // // setup available up to sts_v22f. Available if reconstructing from tsa files (AB 14.05.2024)
+ // Int_t stsAddress00 = 0;
+ // Int_t stsAddress01 = CbmStsAddress::GetAddress(0, 0, 1, 0, 0, 0); // U0 L0 M0 6 cm
+ // Int_t stsAddress02 = CbmStsAddress::GetAddress(0, 0, 1, 1, 0, 0); // U0 L0 M1 6 cm
+ // Int_t stsAddress03 = CbmStsAddress::GetAddress(0, 1, 1, 0, 0, 0); // U0 L1 M0 6 cm
+ // Int_t stsAddress04 = CbmStsAddress::GetAddress(0, 1, 1, 1, 0, 0); // U0 L1 M1 6 cm
+ // Int_t stsAddress05 = CbmStsAddress::GetAddress(1, 0, 1, 0, 0, 0); // U1 L0 M0 6 cm
+ // Int_t stsAddress06 = CbmStsAddress::GetAddress(1, 0, 1, 1, 0, 0); // U1 L0 M1 12 cm
+ // Int_t stsAddress07 = CbmStsAddress::GetAddress(1, 1, 1, 0, 0, 0); // U1 L1 M0 6 cm
+ // Int_t stsAddress08 = CbmStsAddress::GetAddress(1, 1, 1, 1, 0, 0); // U1 L1 M1 12 cm
+ // Int_t stsAddress09 = CbmStsAddress::GetAddress(1, 2, 1, 0, 0, 0); // U1 L2 M0 6 cm
+ // Int_t stsAddress10 = CbmStsAddress::GetAddress(1, 2, 1, 1, 0, 0); // U1 L2 M1 6 cm
+ // Int_t stsAddress11 = CbmStsAddress::GetAddress(1, 2, 1, 2, 0, 0); // U1 L2 M2 6 cm
+
+ std::cout << "STS address00 " << std::dec << stsAddress00 << " " << std::hex << stsAddress00 << std::endl;
+ std::cout << "STS address01 " << std::dec << stsAddress01 << " " << std::hex << stsAddress01 << std::endl;
+ std::cout << "STS address02 " << std::dec << stsAddress02 << " " << std::hex << stsAddress02 << std::endl;
+ std::cout << "STS address03 " << std::dec << stsAddress03 << " " << std::hex << stsAddress03 << std::endl;
+ std::cout << "STS address04 " << std::dec << stsAddress04 << " " << std::hex << stsAddress04 << std::endl;
+ std::cout << "STS address05 " << std::dec << stsAddress05 << " " << std::hex << stsAddress05 << std::endl;
+ std::cout << "STS address06 " << std::dec << stsAddress06 << " " << std::hex << stsAddress06 << std::endl;
+ std::cout << "STS address07 " << std::dec << stsAddress07 << " " << std::hex << stsAddress07 << std::endl;
+ std::cout << "STS address08 " << std::dec << stsAddress08 << " " << std::hex << stsAddress08 << std::endl;
+ std::cout << "STS address09 " << std::dec << stsAddress09 << " " << std::hex << stsAddress09 << std::endl;
+ std::cout << "STS address10 " << std::dec << stsAddress10 << " " << std::hex << stsAddress10 << std::endl;
+ std::cout << "STS address11 " << std::dec << stsAddress11 << " " << std::hex << stsAddress11 << std::endl;
+
+ // --- Now we can define the sensor parameter set and tell recoSts to use it
+ auto sensorParSet = new CbmStsParSetSensor("CbmStsParSetSensor", "STS sensor parameters"
+ "mcbm2024");
+ if (stsAddress00) sensorParSet->SetParSensor(stsAddress00, sensor6cm);
+ sensorParSet->SetParSensor(stsAddress01, sensor6cm);
+ sensorParSet->SetParSensor(stsAddress02, sensor6cm);
+ sensorParSet->SetParSensor(stsAddress03, sensor6cm);
+ sensorParSet->SetParSensor(stsAddress04, sensor6cm);
+ sensorParSet->SetParSensor(stsAddress05, sensor6cm);
+ sensorParSet->SetParSensor(stsAddress06, sensor12cm);
+ sensorParSet->SetParSensor(stsAddress07, sensor6cm);
+ sensorParSet->SetParSensor(stsAddress08, sensor12cm);
+ sensorParSet->SetParSensor(stsAddress09, sensor6cm);
+ sensorParSet->SetParSensor(stsAddress10, sensor6cm);
+ sensorParSet->SetParSensor(stsAddress11, sensor6cm);
+ recoSts->UseSensorParSet(sensorParSet);
+
+ // ASIC params: #ADC channels, dyn. range, threshold, time resol., dead time,
+ // noise RMS, zero-threshold crossing rate
+ auto parAsic = new CbmStsParAsic(128, 31, 75000., 3000., 5., 800., 1000., 3.9789e-3);
+
+ // Module params: number of channels, number of channels per ASIC
+ auto parMod = new CbmStsParModule(2048, 128);
+ parMod->SetAllAsics(*parAsic);
+ recoSts->UseModulePar(parMod);
+
+ // Sensor conditions: full depletion voltage, bias voltage, temperature,
+ // coupling capacitance, inter-strip capacitance
+ auto sensorCond = new CbmStsParSensorCond(70., 140., 268., 17.5, 1.);
+ recoSts->UseSensorCond(sensorCond);
+
+ run->AddTask(recoSts);
+ std::cout << "-I- : Added task " << recoSts->GetName() << std::endl;
+ // ------------------------------------------------------------------------
+ }
+
+ // =========================================================================
+ // === local MUCH Reconstruction ===
+ // =========================================================================
+ if (bMUCH) {
+ TString muchGeoTag;
+ if (geoSetup->GetGeoTag(ECbmModuleId::kMuch, muchGeoTag)) {
+ // --- Parameter file name
+ TString geoTag;
+ geoSetup->GetGeoTag(ECbmModuleId::kMuch, geoTag);
+ Int_t muchFlag = 0;
+ if (geoTag.Contains("mcbm")) muchFlag = 1;
+ TString sectorFile = gSystem->Getenv("VMCWORKDIR");
+ sectorFile += "/parameters/much/much_" + geoTag(0, 4) + "_mcbm_digi_sector.root";
+ //sectorFile += "/parameters/much/much_v22j_mcbm_digi_sector.root";
+ std::cout << "Using parameter file " << sectorFile << std::endl;
+
+ // --- Initialization of the digi scheme
+ auto muchGeoScheme = CbmMuchGeoScheme::Instance();
+ if (!muchGeoScheme->IsInitialized()) {
+ muchGeoScheme->Init(sectorFile.Data(), muchFlag);
+ }
+ // --- Hit finder for GEMs
+ FairTask* muchReco = new CbmMuchFindHitsGem(sectorFile.Data(), muchFlag);
+ run->AddTask(muchReco);
+ std::cout << "-I- " << myName << ": Added task " << muchReco->GetName() << " with parameter file " << sectorFile
+ << std::endl;
+ }
+ }
+ // ------------------------------------------------------------------------
+
+ // =========================================================================
+ // === local TRD Reconstruction ===
+ // =========================================================================
+
+ if (bTRD && geoSetup->IsActive(ECbmModuleId::kTrd)) {
+ CbmTrdClusterFinder* trdCluster;
+ Double_t triggerThreshold = 0.5e-6; // SIS100
+
+ trdCluster = new CbmTrdClusterFinder();
+ trdCluster->SetNeighbourEnable(true, false);
+ trdCluster->SetMinimumChargeTH(triggerThreshold);
+ trdCluster->SetRowMerger(true);
+ CbmTrdClusterFinder::UseRecoHelpers();
+ run->AddTask(trdCluster);
+ std::cout << "-I- : Added task " << trdCluster->GetName() << std::endl;
+
+ CbmTrdHitProducer* trdHit = new CbmTrdHitProducer();
+ //trdHit->SetHitTimeOffset(363); // hit time synchronization for TRD2D determined on run 2391
+ trdHit->SetHitTimeOffset(80); // hit time synchronization for TRD2D determined on run 2914
+ run->AddTask(trdHit);
+ std::cout << "-I- : Added task " << trdHit->GetName() << std::endl;
+ }
+
+ // =========================================================================
+ // === RICH Reconstruction ===
+ // =========================================================================
+
+ if (bRICH && geoSetup->IsActive(ECbmModuleId::kRich)) {
+ // ----- Local reconstruction of RICH Hits ------------------------------
+ CbmRichMCbmHitProducer* hitProd = new CbmRichMCbmHitProducer();
+ hitProd->SetMappingFile(std::string(srcDir.Data())
+ + "/macro/rich/mcbm/beamtime/mRICH_Mapping_vert_20190318_elView.geo");
+ hitProd->SetIcdFilenameBase(std::string(srcDir.Data()) + "/macro/beamtime/mcbm2022/icd_offset_it");
+ hitProd->setToTLimits(23.7, 30.0);
+ hitProd->applyToTCut();
+ hitProd->applyICDCorrection();
+ run->AddTask(hitProd);
+ // ------------------------------------------------------------------------
+
+ // ----- Local reconstruction in RICh -> Finding of Rings ---------------
+ CbmRichReconstruction* richReco = new CbmRichReconstruction();
+ richReco->UseMCbmSetup();
+ run->AddTask(richReco);
+ // ------------------------------------------------------------------------
+ }
+
+ // =========================================================================
+ // === TOF Hitfinding ===
+ // =========================================================================
+
+ TString parPath = srcDir + "/parameters/mcbm/";
+ if (bTOF && geoSetup->IsActive(ECbmModuleId::kTof)) {
+ TString cFname;
+ switch (iTofCluMode) {
+ case 1: {
+ // ----- TOF defaults ------------------------
+ Int_t calMode = 93;
+ Int_t calSel = 1;
+ Int_t calSm = 0;
+ Int_t RefSel = 0;
+ Double_t dDeadtime = 50.;
+ Int_t iSel2 = 500;
+ Bool_t bOut = kTRUE;
+
+ // ------------------------------------------------------------------------
+ gROOT->LoadMacro(srcDir + "/macro/beamtime/mcbm2022/ini_tof_clusterizer.C");
+ Char_t* cCmd =
+ Form("ini_tof_clusterizer(%d,%d,%d,%d,\"%s\",%d,%d,%d,%f,\"%s\",\"%s\")", calMode, calSel, calSm, RefSel,
+ cFileId.Data(), iCalSet, (Int_t) bOut, iSel2, dDeadtime, cCalId.Data(), parPath.Data());
+
+ cout << "<I> " << cCmd << endl;
+ gInterpreter->ProcessLine(cCmd);
+ // disable histogramming
+ CbmTofEventClusterizer* tofClust = CbmTofEventClusterizer::Instance();
+ tofClust->SetDutId(-1); // to disable histogramming
+ } break;
+
+ default: {
+ ;
+ }
+ }
+ // -------------------------------------------------------------------------
+
+ // =========================================================================
+ // === Tof Tracking ===
+ // =========================================================================
+
+ if (bTOFtr) {
+ cout << "<I> Initialize Tof tracker by ini_trks" << endl;
+
+ gROOT->LoadMacro(srcDir + "/macro/beamtime/mcbm2022/ini_tof_trks.C");
+ Char_t* cCmd =
+ Form("ini_tof_trks(%d,%d,%d,%6.2f,%8.1f,\"%s\",%d,%d,%d,%f,\"%s\")", iSel, iTrackingSetup, iGenCor, dScalFac,
+ dChi2Lim2, cCalId.Data(), (Int_t) bUseSigCalib, iCalOpt, iTrkPar, dTOffScal, parPath.Data());
+ cout << "<I> " << cCmd << endl;
+ gInterpreter->ProcessLine(cCmd);
+
+ CbmTofFindTracks* tofFindTracks = CbmTofFindTracks::Instance();
+ Int_t iNStations = tofFindTracks->GetNStations();
+ }
+ }
+
+
+ //Constant Field
+ // CbmConstField *fMagField = new CbmConstField();
+ // fMagField->SetFieldXYZ(0, 0 ,0 ); // values are in kG
+ // fMagField->SetFieldRegions(-10, -10 ,-10 , 10, 10 , 10 );
+ // run->SetField(fMagField);
+
+ // =========================================================================
+ // === L1 ===
+ // =========================================================================
+ if (bL1) {
+ run->AddTask(new CbmTrackingDetectorInterfaceInit());
+
+ CbmL1* l1 = new CbmL1("L1");
+ l1->SetMcbmMode();
+ l1->SetVerbose(3);
+ run->AddTask(l1);
+
+ CbmL1GlobalTrackFinder* globalTrackFinder = new CbmL1GlobalTrackFinder();
+ FairTask* globalFindTracks = new CbmL1GlobalFindTracksEvents(globalTrackFinder);
+ run->AddTask(globalFindTracks);
+ }
+ // =========================================================================
+ // === QA ===
+ // =========================================================================
+ if (bQA) {
+ // e.g for RICH:
+ CbmRichMCbmQaReal* qaTask = new CbmRichMCbmQaReal();
+ Int_t taskId = 1;
+ if (taskId < 0) {
+ qaTask->SetOutputDir(Form("result_run%d", uRunId));
+ }
+ else {
+ qaTask->SetOutputDir(Form("result_run%d_%05d", uRunId, taskId));
+ }
+ //qaTask->XOffsetHistos(+25.0);
+ qaTask->XOffsetHistos(-4.1);
+ if (uRunId > 2351) qaTask->XOffsetHistos(0.0);
+ qaTask->SetMaxNofDrawnEvents(100);
+ qaTask->SetTotRich(23.7, 30.0);
+ qaTask->SetTriggerRichHits(eb_TriggerMinNumberRich);
+ qaTask->SetTriggerTofHits(0); // eb_TriggerMinNumberTof);
+ qaTask->SetSEDisplayRingOnly();
+ run->AddTask(qaTask);
+ }
+ // ------------------------------------------------------------------------
+
+ // ----- Parameter database --------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Set runtime DB" << std::endl;
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ FairParRootFileIo* parIo1 = new FairParRootFileIo();
+ FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+ FairParRootFileIo* parIo3 = new FairParRootFileIo();
+ //parIo1->open(parFileIn.Data(), "READ");
+ //rtdb->setFirstInput(parIo1);
+ parIo2->open(parFileList, "in");
+ rtdb->setSecondInput(parIo2);
+ parIo3->open(parFileOut.Data(), "RECREATE");
+ // ------------------------------------------------------------------------
+ rtdb->setOutput(parIo3);
+
+ // ----- 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, nTimeslices);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Finish -------------------------------------------------------
+
+ rtdb->print();
+ rtdb->saveOutput();
+ run->CreateGeometryFile(geoFileOut);
+
+ 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 " << outFile << std::endl;
+ std::cout << "Parameter file is " << parFileOut << std::endl;
+ std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << std::endl;
+ std::cout << 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();
+ // ------------------------------------------------------------------------
+
+ /// --- Screen output for automatic tests
+ std::cout << " Test passed" << std::endl;
+ std::cout << " All ok " << std::endl;
+
+ return kTRUE;
+}
diff --git a/macro/qa/run_recoQa.C b/macro/qa/run_recoQa.C
index 5243739d3705bbe9d4b9de6cdc754133997ef31d..b93f19c8b6e2106c34af7d810c34fc09f35d3069 100644
--- a/macro/qa/run_recoQa.C
+++ b/macro/qa/run_recoQa.C
@@ -36,7 +36,7 @@
#include <TStopwatch.h>
#endif
-void run_recoQa(Int_t nEvents = -1, TString dataset = "2391_node8_0_0000",
+void run_recoQa(Int_t nEvents = -1, TString recFile = "2391_node8_0_0000.rec.root",
TString setupName = "mcbm_beam_2022_05_23_nickel", Bool_t bUseMC = kFALSE)
{
@@ -44,7 +44,7 @@ void run_recoQa(Int_t nEvents = -1, TString dataset = "2391_node8_0_0000",
// Adjust this part according to your requirements
// ----- Logger settings ----------------------------------------------
- FairLogger::GetLogger()->SetLogScreenLevel("DEBUG");
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
// ------------------------------------------------------------------------
@@ -55,10 +55,11 @@ void run_recoQa(Int_t nEvents = -1, TString dataset = "2391_node8_0_0000",
// ------------------------------------------------------------------------
// ----- In- and output file names ------------------------------------
+ TString dataset = recFile;
+ dataset.ReplaceAll(".rec.root", "");
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";
@@ -69,6 +70,8 @@ void run_recoQa(Int_t nEvents = -1, TString dataset = "2391_node8_0_0000",
CbmSetup* setup = CbmSetup::Instance();
setup->LoadSetup(setupName);
+ //setup->SetActive(ECbmModuleId::kSts, false);
+ //setup->SetActive(ECbmModuleId::kTof, false);
// In general, the following parts need not be touched
// ========================================================================
@@ -82,6 +85,15 @@ void run_recoQa(Int_t nEvents = -1, TString dataset = "2391_node8_0_0000",
gDebug = 0;
// ------------------------------------------------------------------------
+ std::map<std::string, TGeoHMatrix>* matrices{nullptr};
+ TFile* misalignmentMatrixRootfile = TFile::Open(Form("AlignmentMatrices_%s.root", setupName.Data()), "READ");
+ if (misalignmentMatrixRootfile && misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->GetObject("MisalignMatrices", matrices);
+ misalignmentMatrixRootfile->Close();
+ }
+ else
+ LOG(warning) << "Could not find usable alignment file. Skip.";
+
// ----- FairRunAna ---------------------------------------------------
FairFileSource* inputSource = new FairFileSource(recFile);
if (bUseMC) {
@@ -102,6 +114,8 @@ void run_recoQa(Int_t nEvents = -1, TString dataset = "2391_node8_0_0000",
FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
// ------------------------------------------------------------------------
+ if (matrices) run->AddAlignmentMatrices(*matrices);
+
// ----- MCDataManager (legacy mode) -----------------------------------
if (bUseMC) {
std::cout << "-I- " << myName << ": Adding MC manager and MC to reco matching tasks\n";
@@ -141,10 +155,21 @@ void run_recoQa(Int_t nEvents = -1, TString dataset = "2391_node8_0_0000",
run->AddTask(l1);
// ----- Reco QA --------------------------------------------
+ CbmRecoQaTask::Detector* det(nullptr);
+ CbmRecoQaTask::Detector::View* view(nullptr);
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});
+ // recoQa->SetSetupClass(CbmRecoQaTask::kMcbm22);
+ recoQa->SetSetupClass(CbmRecoQaTask::kMcbm24);
+ recoQa->AddEventFilter(CbmRecoQaTask::EventFilter::eEventCut::kMultTrk)->SetFilter({1, 2});
+ recoQa->AddTrackFilter(CbmRecoQaTask::TrackFilter::eTrackCut::kSts)->SetFilter({3});
+ recoQa->AddTrackFilter(CbmRecoQaTask::TrackFilter::eTrackCut::kTrd)->SetFilter({2});
+ recoQa->AddTrackFilter(CbmRecoQaTask::TrackFilter::eTrackCut::kTof)->SetFilter({1});
+ // if ((det = recoQa->GetDetector(ECbmModuleId:kSts))) {
+ // view = det->GetView("U0L0M0");
+ // view->SetProjection();
+ // }
+ // 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);
diff --git a/macro/run/run_unpack_tsa.C b/macro/run/run_unpack_tsa.C
index afe3c100314468967e4a676c92a44adc090e4224..4072f0dd45bf8c5d37cf26e6fb5d8868dc542cc8 100644
--- a/macro/run/run_unpack_tsa.C
+++ b/macro/run/run_unpack_tsa.C
@@ -315,6 +315,24 @@ void run_unpack_tsa(std::vector<std::string> infile = {"test.tsa"}, UInt_t runid
if (2350 <= runid) {
trd1Dconfig->SetSystemTimeOffset(1300); // [ns] value to be updated
}
+ if (runid >= 2700) //mCBM2024
+ {
+ for (int elinkId = 0; elinkId < 36; ++elinkId) {
+ trd1Dconfig->SetElinkTimeOffset(20736, elinkId, -36);
+ trd1Dconfig->SetElinkTimeOffset(20737, elinkId, -37);
+ trd1Dconfig->SetElinkTimeOffset(20738, elinkId, -65);
+ trd1Dconfig->SetElinkTimeOffset(20739, elinkId, -52);
+ trd1Dconfig->SetElinkTimeOffset(20740, elinkId, -50);
+ trd1Dconfig->SetElinkTimeOffset(20741, elinkId, -49);
+
+ //trd1Dconfig->SetElinkTimeOffset(20992, elinkId, 0); //no correlation
+ trd1Dconfig->SetElinkTimeOffset(20993, elinkId, -50);
+ trd1Dconfig->SetElinkTimeOffset(20994, elinkId, -57);
+ trd1Dconfig->SetElinkTimeOffset(20995, elinkId, -52);
+ trd1Dconfig->SetElinkTimeOffset(20996, elinkId, -52);
+ trd1Dconfig->SetElinkTimeOffset(20997, elinkId, -53);
+ }
+ }
}
// -------------
diff --git a/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.cxx b/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.cxx
index 51cffc37a3a6203c1f30b5211a6c6aca530feb80..b687c94fc406232355892efb3ec29059df218783 100644
--- a/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.cxx
+++ b/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.cxx
@@ -1,6 +1,6 @@
-/* Copyright (C) 2020-2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
+/* Copyright (C) 2020-2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Pierre-Alain Loizeau [committer], Dominik Smith */
+ Authors: Pierre-Alain Loizeau [committer], Dominik Smith, Alexandru Bercuci */
#include "CbmAlgoBuildRawEvents.h"
@@ -14,8 +14,8 @@
#include "CbmRichDigi.h"
#include "CbmStsDigi.h"
#include "CbmTofDigi.h"
+#include "CbmTrdAddress.h"
#include "CbmTrdDigi.h"
-
#include "TimesliceMetaData.h"
/// FAIRROOT headers
@@ -26,7 +26,7 @@
/// FAIRSOFT headers (geant, boost, ...)
#include <TCanvas.h>
#include <TClonesArray.h>
-#include <TFolder.h>
+#include <TDirectoryFile.h>
#include <TH1.h>
#include <TH2.h>
#include <THttpServer.h>
@@ -36,6 +36,8 @@
/// C/C++ headers
#include <algorithm>
+#define VERBOSE 0
+
template<>
void CbmAlgoBuildRawEvents::LoopOnSeeds<Double_t>();
@@ -284,7 +286,10 @@ void CbmAlgoBuildRawEvents::LoopOnSeeds<Double_t>()
Double_t dTime = fSeedTimes->at(uSeed);
/// Check if seed in acceptance window (is this needed here?)
- if (dTime < fdSeedWindowBeg) { continue; }
+ if (dTime < fdSeedWindowBeg) {
+ LOG(debug4) << Form("CbmAlgoBuildRawEvents :: reject for time %f < %f\n", dTime, fdSeedWindowBeg);
+ continue;
+ }
else if (fdSeedWindowEnd < dTime) {
break;
}
@@ -305,10 +310,17 @@ void CbmAlgoBuildRawEvents::LoopOnSeeds()
for (UInt_t uDigi = 0; uDigi < uNbRefDigis; ++uDigi) {
LOG(debug) << Form("Checking seed %6u / %6u", uDigi, uNbRefDigis);
const DigiSeed* pDigi = GetDigi<DigiSeed>(uDigi);
+ // hack for mCBM2024 data and Bmon station selection
+ if (fRefDet.detId == ECbmModuleId::kBmon && !filterBmon(pDigi->GetAddress())) continue;
+
const Double_t dTime = pDigi->GetTime();
+ //printf("time = %f %d %d\n", dTime, CbmTofAddress::GetChannelSide(add), CbmTofAddress::GetChannelId(add));
/// Check if seed in acceptance window
- if (dTime < fdSeedWindowBeg) { continue; }
+ if (dTime < fdSeedWindowBeg) {
+ LOG(debug4) << Form("CbmAlgoBuildRawEvents :: reject for time %f < %f\n", dTime, fdSeedWindowBeg);
+ continue;
+ }
else if (fdSeedWindowEnd < dTime) {
break;
}
@@ -329,6 +341,8 @@ void CbmAlgoBuildRawEvents::CheckSeed(Double_t dSeedTime, UInt_t uSeedDigiIdx)
{
/// If previous event valid and event overlap not allowed, check if we are in overlap
/// and react accordingly
+
+ LOG(debug4) << Form("CbmAlgoBuildRawEvents :: CheckSeed(%f, %d)\n", dSeedTime, uSeedDigiIdx);
if (nullptr != fCurrentEvent
&& (EOverlapModeRaw::AllowOverlap != fOverMode || dSeedTime - fdPrevEvtTime < GetSeedTimeWinRange())
&& dSeedTime - fdPrevEvtTime < fdWidestTimeWinRange) {
@@ -552,7 +566,7 @@ void CbmAlgoBuildRawEvents::SearchMatches(Double_t dSeedTime, RawEventBuilderDet
/// This algo relies on time sorted vectors for the selected detectors
UInt_t uLocalIndexStart = detMatch.fuStartIndex;
UInt_t uLocalIndexEnd = detMatch.fuStartIndex;
-
+ LOG(debug4) << Form("CbmAlgoBuildRawEvents :: SearchMatches(%f, %s)\n", dSeedTime, detMatch.sName.data());
/// Check the Digis until out of window
const UInt_t uNbSelDigis = GetNofDigis(detMatch.detId);
/// Loop on size of vector
@@ -561,6 +575,31 @@ void CbmAlgoBuildRawEvents::SearchMatches(Double_t dSeedTime, RawEventBuilderDet
const Double_t dTime = pDigi->GetTime();
const Double_t dTimeDiff = dTime - dSeedTime;
LOG(debug4) << detMatch.sName << Form(" => Checking match %6u / %6u, dt %f", uDigi, uNbSelDigis, dTimeDiff);
+ // int32_t add = pDigi->GetAddress(), dId[3] = {-1, -1, -1};
+ // switch (detMatch.detId) {
+ // case ECbmModuleId::kBmon:
+ // dId[0] = CbmTofAddress::GetChannelSide(add);
+ // dId[1] = CbmTofAddress::GetChannelId(add);
+ // break;
+ // case ECbmModuleId::kSts:
+ // dId[0] = CbmStsAddress::GetElementId(add, EStsElementLevel::kStsUnit);
+ // dId[1] = CbmStsAddress::GetElementId(add, EStsElementLevel::kStsLadder);
+ // dId[2] = CbmStsAddress::GetElementId(add, EStsElementLevel::kStsModule);
+ // break;
+ // case ECbmModuleId::kTrd:
+ // case ECbmModuleId::kTrd2d:
+ // dId[0] = CbmTrdAddress::GetLayerId(add);
+ // dId[1] = CbmTrdAddress::GetModuleId(add);
+ // dId[2] = CbmTrdAddress::GetModuleAddress(add);
+ // break;
+ // case ECbmModuleId::kTof:
+ // dId[0] = CbmTofAddress::GetSmId(add);
+ // dId[1] = CbmTofAddress::GetSmType(add);
+ // dId[2] = CbmTofAddress::GetRpcId(add);
+ // break;
+ // default: break;
+ // }
+ // LOG(debug4) << Form("CbmAlgoBuildRawEvents :: Checking match %6u / %6u, dt=%f %s[%d %d %d]\n", uDigi, uNbSelDigis, dTimeDiff, detMatch.sName.data(), dId[0], dId[1], dId[2]);
/// Check if within time window, update start/stop indices if needed
if (dTimeDiff < detMatch.fdTimeWinBeg) {
@@ -597,15 +636,18 @@ void CbmAlgoBuildRawEvents::SearchMatches(Double_t dSeedTime, RawEventBuilderDet
/// Update the StartIndex and EndIndex for the next event seed
detMatch.fuStartIndex = uLocalIndexStart;
detMatch.fuEndIndex = uLocalIndexEnd;
+ LOG(debug4) << Form("CbmAlgoBuildRawEvents :: SearchMatches(%d, %d)\n", uLocalIndexStart, uLocalIndexEnd);
}
void CbmAlgoBuildRawEvents::AddDigiToEvent(const RawEventBuilderDetector& det, Int_t _entry)
{
+ LOG(debug4) << Form("CbmAlgoBuildRawEvents :: :: AddDigiToEvent(%s)\n", det.sName.data());
fCurrentEvent->AddData(det.dataType, _entry);
}
void CbmAlgoBuildRawEvents::CheckTriggerCondition(Double_t dSeedTime)
{
+ LOG(debug4) << Form("CbmAlgoBuildRawEvents :: :: CheckTriggerCondition(%f)\n", dSeedTime);
/// Check if event is filling trigger conditions and clear it if not
if (HasTrigger(fCurrentEvent)) {
fdPrevEvtTime = dSeedTime;
@@ -642,51 +684,65 @@ Bool_t CbmAlgoBuildRawEvents::HasTrigger(CbmEvent* event)
return kTRUE;
}
-void CbmAlgoBuildRawEvents::SetBmonEventTime(CbmEvent* event)
+bool CbmAlgoBuildRawEvents::SetBmonEventTime(CbmEvent* event)
{
const int32_t iNbDigis = event->GetNofData(ECbmDataType::kBmonDigi);
+ double eventTime(0.);
+ bool timeSet(false);
- if (0 < iNbDigis) {
- uint idx = event->GetIndex(ECbmDataType::kBmonDigi, 0);
+ LOG(debug4) << Form("CbmAlgoBuildRawEvents :: SetBmonEventTime(%p, %d)\n", (void*) event, iNbDigis);
+ for (int idigi = 0; idigi < iNbDigis; ++idigi) {
+ uint idx = event->GetIndex(ECbmDataType::kBmonDigi, idigi);
const CbmBmonDigi* pDigi = GetDigi<CbmBmonDigi>(idx);
- double eventTime = pDigi->GetTime();
+ if (nullptr == pDigi) continue;
+ if (!filterBmon(pDigi->GetAddress())) continue;
- for (int idigi = 1; idigi < iNbDigis; ++idigi) {
- idx = event->GetIndex(ECbmDataType::kBmonDigi, idigi);
- pDigi = GetDigi<CbmBmonDigi>(idx);
- if (nullptr == pDigi) continue;
- eventTime = std::min(pDigi->GetTime(), eventTime);
+ if (!timeSet) {
+ eventTime = pDigi->GetTime();
+ timeSet = true;
}
- event->SetStartTime(eventTime);
+ else
+ eventTime = std::min(pDigi->GetTime(), eventTime);
}
+
+ if (timeSet)
+ event->SetStartTime(eventTime);
+ else
+ LOG(warning) << "CbmAlgoBuildRawEvents::SetBmonEventTime : failed";
+
+ return timeSet;
}
Bool_t CbmAlgoBuildRawEvents::CheckTriggerConditions(CbmEvent* event, const RawEventBuilderDetector& det)
{
+ LOG(debug4) << Form("CbmAlgoBuildRawEvents :: :: CheckTriggerConditions(%p, %s)\n", (void*) event, det.sName.data());
/// Check if both Trigger conditions disabled for this detector
if (0 == det.fuTriggerMinDigis && det.fiTriggerMaxDigis < 0) { return kTRUE; }
/// Check if detector present
if (!CheckDataAvailable(det.detId)) {
- LOG(warning) << "Event does not have digis storage for " << det.sName
- << " while the following trigger min/max are defined: " << det.fuTriggerMinDigis << " "
- << det.fiTriggerMaxDigis;
+ LOG(debug2) << "Event does not have digis storage for " << det.sName
+ << " while the following trigger min/max are defined: " << det.fuTriggerMinDigis << " "
+ << det.fiTriggerMaxDigis;
return kFALSE;
}
/// Check trigger rejection by minimal/maximal number or absence, if enabled/requested
- int32_t iNbDigis = event->GetNofData(det.dataType);
+ int32_t iNbDigis = event->GetNofData(det.dataType);
+ int32_t iNbFilteredDigis = (det.detId == ECbmModuleId::kBmon ? getNofFilteredBmonDigis(event) : iNbDigis);
/// a.Check trigger rejection by minimal number (if enabled)
- if (0 < det.fuTriggerMinDigis && ((-1 == iNbDigis) || (static_cast<UInt_t>(iNbDigis) < det.fuTriggerMinDigis))) {
- LOG(debug2) << "Event does not have enough digis: " << iNbDigis << " vs " << det.fuTriggerMinDigis << " for "
- << det.sName;
+ if (0 < det.fuTriggerMinDigis
+ && ((-1 == iNbFilteredDigis) || (static_cast<UInt_t>(iNbFilteredDigis) < det.fuTriggerMinDigis))) {
+ LOG(debug2) << "Event does not have enough digis: " << iNbFilteredDigis << " vs " << det.fuTriggerMinDigis
+ << " for " << det.sName;
return kFALSE;
}
/// b.Check trigger rejection by maximal number (if enabled)
- if (0 <= det.fiTriggerMaxDigis && det.fiTriggerMaxDigis < iNbDigis) {
- LOG(debug2) << "Event Has too many digis: " << iNbDigis << " vs " << det.fiTriggerMaxDigis << " for " << det.sName;
+ if (0 <= det.fiTriggerMaxDigis && det.fiTriggerMaxDigis < iNbFilteredDigis) {
+ LOG(debug2) << "Event Has too many digis: " << iNbFilteredDigis << " vs " << det.fiTriggerMaxDigis << " for "
+ << det.sName;
return kFALSE;
}
@@ -959,7 +1015,7 @@ uint64_t CbmAlgoBuildRawEvents::GetSizeFromDigisNb(ECbmModuleId detId, uint64_t
//----------------------------------------------------------------------
void CbmAlgoBuildRawEvents::CreateHistograms()
{
- outFolder = new TFolder("AlgoBuildRawEventsHist", " AlgoBuildRawEvents Histos");
+ outFolder = new TDirectoryFile("AlgoBuildRawEventsHist", " AlgoBuildRawEvents Histos");
outFolder->Clear();
fhEventTime = new TH1F("hEventTime", "seed time of the events; Seed time within TS [s]; Events", 100000, 0, 0.2);
@@ -1587,7 +1643,7 @@ void CbmAlgoBuildRawEvents::AddDetector(ECbmModuleId selDet, ECbmDataType dataTy
fdTimeWinBegIn, fdTimeWinEndIn));
}
-void CbmAlgoBuildRawEvents::SetReferenceDetector(RawEventBuilderDetector refDetIn)
+void CbmAlgoBuildRawEvents::SetReferenceDetector(RawEventBuilderDetector refDetIn, std::vector<bool> select)
{
/// Loop on selection detectors
for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
@@ -1627,6 +1683,21 @@ void CbmAlgoBuildRawEvents::SetReferenceDetector(RawEventBuilderDetector refDetI
UpdateWidestTimeWinRange();
/// Detect usage of BMon to set Event StartTime
CheckBmonInUse();
+
+ if (fbBmonInUse) {
+ if (select.size()) {
+ LOG(info) << "CbmAlgoBuildRawEvents::SetReferenceDetector => Detected Bmon station selection.";
+ fUseBmonMap.assign(select.size(), false);
+ int it0(0);
+ for (auto t0 : select)
+ SwitchBmonStation(it0++, (t0 > 0));
+ }
+ }
+ else {
+ if (select.size())
+ LOG(warning) << "CbmAlgoBuildRawEvents::SetReferenceDetector => Detected use of selector\nfor a reference "
+ "detector which does not support this option. Skip selection for the moment.";
+ }
}
void CbmAlgoBuildRawEvents::AddDetector(RawEventBuilderDetector selDet)
@@ -1841,4 +1912,49 @@ void CbmAlgoBuildRawEvents::UpdateWidestTimeWinRange()
}
}
+void CbmAlgoBuildRawEvents::SwitchBmonStation(int id, bool on)
+{
+ if (id < 0 || id >= (int) fUseBmonMap.size()) {
+ LOG(warning) << "CbmAlgoBuildRawEvents::SwitchBmonStation: Bmon station id outside range. Skip.";
+ return;
+ }
+ LOG(info) << "CbmAlgoBuildRawEvents::SwitchBmonStation: Bmon" << id << " station switched " << (on ? "ON" : "OFF")
+ << " for trigger.";
+ fUseBmonMap[id] = on;
+}
+
+bool CbmAlgoBuildRawEvents::filterBmon(int32_t add)
+{
+ // skip any filtering if not explicitly requested by user.
+ if (!fUseBmonMap.size()) return true;
+
+ // consistency check on the Bmon detector
+ if (CbmTofAddress::GetSmType(add) != 5) {
+ LOG(warning) << "CbmAlgoBuildRawEvents::filterBmon: Bmon digi with wrong address [GetSmType() != 5]. Skip.";
+ return false;
+ }
+ size_t mod = (size_t) CbmTofAddress::GetChannelSide(add);
+ // select Bmon detector
+ if (!fUseBmonMap[mod]) {
+ LOG(debug2) << "CbmAlgoBuildRawEvents::filterBmon : reject seed from Bmon" << mod;
+ return false;
+ }
+ return true;
+}
+
+int32_t CbmAlgoBuildRawEvents::getNofFilteredBmonDigis(CbmEvent* event)
+{
+ // skip any filtering if not explicitly requested by user.
+ if (!fUseBmonMap.size()) return event->GetNofData(ECbmDataType::kBmonDigi);
+
+ int32_t ndigi(0);
+ for (size_t idigi = 0; idigi < event->GetNofData(ECbmDataType::kBmonDigi); ++idigi) {
+ uint idx = event->GetIndex(ECbmDataType::kBmonDigi, idigi);
+ const CbmBmonDigi* pDigi = GetDigi<CbmBmonDigi>(idx);
+ if (!filterBmon(pDigi->GetAddress())) continue;
+ ndigi++;
+ }
+ return ndigi;
+}
+
ClassImp(CbmAlgoBuildRawEvents)
diff --git a/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.h b/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.h
index 8f49b27b93ab768e6259ffffbfac9766ab883f6b..bf474c867a1cb3932dc73362255932f1a0a0482b 100644
--- a/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.h
+++ b/reco/eventbuilder/digis/CbmAlgoBuildRawEvents.h
@@ -1,6 +1,6 @@
-/* Copyright (C) 2020-2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
+/* Copyright (C) 2020-2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Pierre-Alain Loizeau [committer], Dominik Smith */
+ Authors: Pierre-Alain Loizeau [committer], Dominik Smith, Alexandru Bercuci*/
#ifndef CBMALGOBUILDRAWEVENTS_H
#define CBMALGOBUILDRAWEVENTS_H
@@ -40,7 +40,7 @@ class TProfile;
class TNamed;
class TStopwatch;
class TCanvas;
-class TFolder;
+class TDirectoryFile;
enum class EOverlapModeRaw
{
@@ -153,7 +153,7 @@ public:
void AddDetector(ECbmModuleId selDet, ECbmDataType dataTypeIn, std::string sNameIn, UInt_t uTriggerMinDigisIn = 0,
Int_t iTriggerMaxDigisIn = -1, Double_t fdTimeWinBegIn = -100, Double_t fdTimeWinEndIn = 100);
- void SetReferenceDetector(RawEventBuilderDetector refDetIn);
+ void SetReferenceDetector(RawEventBuilderDetector refDetIn, std::vector<bool> select = {});
void AddDetector(RawEventBuilderDetector selDet);
void RemoveDetector(RawEventBuilderDetector selDet);
@@ -162,7 +162,6 @@ public:
void SetTriggerMinLayersNumber(ECbmModuleId selDet, UInt_t uVal);
void SetTriggerWindow(ECbmModuleId selDet, Double_t dWinBeg, Double_t dWinEnd);
void SetHistogramMaxDigiNb(ECbmModuleId selDet, Double_t dDigiNbMax);
-
void SetTsParameters(Double_t dTsStartTime, Double_t dTsLength, Double_t dTsOverLength)
{
fdTsStartTime = dTsStartTime;
@@ -225,7 +224,7 @@ public:
}
// Output folder for histograms
- TFolder* GetOutFolder() { return outFolder; }
+ TDirectoryFile* GetOutFolder() { return outFolder; }
/// Data output access
std::vector<CbmEvent*>& GetEventVector() { return fEventVector; }
@@ -257,11 +256,19 @@ private:
void UpdateTimeWinBoundariesExtrema();
void UpdateWidestTimeWinRange();
- void SetBmonEventTime(CbmEvent* event);
+ void SwitchBmonStation(int id, bool on = true);
+ bool SetBmonEventTime(CbmEvent* event);
void CheckBmonInUse();
- TFolder* outFolder; // oputput folder to store histograms
+ /** \brief Filter Bmon stations. Hack added for the mCBM2024 data (AB)
+ * \param[in] add address of the Bmon digi
+ * \return true if selected
+ * \sa fuUseBmonMap, SelectBmonStations(), getNofFilteredDigis() */
+ bool filterBmon(int32_t add);
+ int32_t getNofFilteredBmonDigis(CbmEvent* ev);
+
+ TDirectoryFile* outFolder; // oputput folder to store histograms
/// Constants
static constexpr Double_t kdDefaultTimeWinBeg = -100.0;
@@ -274,6 +281,8 @@ private:
Bool_t fbUseMuchBeamtimeDigi = kTRUE; //! Switch between MUCH digi classes
Bool_t fbGetTimings = kFALSE; //! Measure CPU time using stopwatch
Bool_t fbUseTsMetaData = kTRUE; //! Read Ts Parameters from input tree
+ std::vector<bool> fUseBmonMap = {}; //! bit map for Bmon trigger. Defined by user
+
/// Event building mode and detectors selection
EOverlapModeRaw fOverMode {EOverlapModeRaw::AllowOverlap};
@@ -372,7 +381,7 @@ private:
Double_t fdPrevEvtTime = 0.; //! Save previous time information
Double_t fdPrevEvtEndTime = 0.; //! Save previous event last digi time information
- ClassDefNV(CbmAlgoBuildRawEvents, 1);
+ ClassDefNV(CbmAlgoBuildRawEvents, 2);
};
#endif // CBMALGOBUILDRAWEVENTS_H
diff --git a/reco/eventbuilder/digis/CbmTaskBuildRawEvents.h b/reco/eventbuilder/digis/CbmTaskBuildRawEvents.h
index 00d6d79a2cef6bc419efef8ed68e5ed80cb10bdd..7b5ab4214fac5d9cf256a47ede680d3464a00816 100644
--- a/reco/eventbuilder/digis/CbmTaskBuildRawEvents.h
+++ b/reco/eventbuilder/digis/CbmTaskBuildRawEvents.h
@@ -75,9 +75,9 @@ public:
fbFillHistos = bFlag;
if (nullptr != fpAlgo) fpAlgo->SetFillHistos(fbFillHistos);
}
- void SetReferenceDetector(RawEventBuilderDetector refDet)
+ void SetReferenceDetector(RawEventBuilderDetector refDet, std::vector<bool> select = {})
{
- if (nullptr != fpAlgo) fpAlgo->SetReferenceDetector(refDet);
+ if (nullptr != fpAlgo) fpAlgo->SetReferenceDetector(refDet, select);
}
void AddDetector(RawEventBuilderDetector selDet)
{
diff --git a/reco/qa/CMakeLists.txt b/reco/qa/CMakeLists.txt
index f31759495be590b7ca18051c896642ecba9ff413..8dc5d62e6ef29ce756bcc3db179010735c90c315 100644
--- a/reco/qa/CMakeLists.txt
+++ b/reco/qa/CMakeLists.txt
@@ -13,7 +13,7 @@ set(SRCS
set(LIBRARY_NAME CbmRecoQa)
set(LINKDEF RecoQaLinkDef.h)
set(PUBLIC_DEPENDENCIES
- #Algo
+ CbmSimSteer
CbmBase
CbmData
FairRoot::Base
diff --git a/reco/qa/CbmRecoQaTask.cxx b/reco/qa/CbmRecoQaTask.cxx
index 259210dc0bc8d04b1caff856e418143974fb52bb..17c82881d75aaf27ca2bdbe7b30645b50703b555 100644
--- a/reco/qa/CbmRecoQaTask.cxx
+++ b/reco/qa/CbmRecoQaTask.cxx
@@ -5,21 +5,92 @@
// CBM headers
#include "CbmRecoQaTask.h"
+#include "../../core/detectors/rich/utils/CbmRichUtil.h"
+#include "CbmDefs.h"
+#include "CbmEvent.h"
+#include "CbmFsdHit.h"
#include "CbmGlobalTrack.h"
+#include "CbmMvdHit.h"
+#include "CbmPsdHit.h"
+#include "CbmRichHit.h"
+#include "CbmSetup.h"
+#include "CbmStsAddress.h"
+#include "CbmStsHit.h"
#include "CbmTimeSlice.h"
#include "CbmTofAddress.h"
+#include "CbmTofHit.h"
#include "CbmTrdAddress.h"
+#include "CbmTrdHit.h"
// FAIR headers
-#include "FairRootManager.h"
-#include "FairSink.h"
+#include <FairRootManager.h>
+#include <FairSink.h>
// ROOT headers
-#include "TClonesArray.h"
-#include "TH2D.h"
+#include <TClonesArray.h>
+#include <TGeoManager.h>
+#include <TGeoNode.h>
+#include <TH2D.h>
using namespace std;
+std::bitset<kRecoQaNConfigs> CbmRecoQaTask::fuRecoConfig = {};
+
//_____________________________________________________________________
CbmRecoQaTask::CbmRecoQaTask() : FairTask("RecoQA", 0) {}
+//_____________________________________________________________________
+CbmRecoQaTask::Detector* CbmRecoQaTask::AddDetector(ECbmModuleId id)
+{
+ /* Interface to the user. Check that the det defined in the outside world fulfills the quality criteria.
+ */
+ switch (id) {
+ case ECbmModuleId::kMvd:
+ case ECbmModuleId::kSts:
+ case ECbmModuleId::kMuch:
+ case ECbmModuleId::kRich:
+ case ECbmModuleId::kTrd:
+ case ECbmModuleId::kTrd2d:
+ case ECbmModuleId::kTof:
+ case ECbmModuleId::kFsd:
+ case ECbmModuleId::kPsd:
+ if (GetDetector(id)) {
+ LOG(warn) << GetName() << "::AddDetector : det " << ToString(id)
+ << " already registered. Using \"CbmRecoQaTask::GetDetector(ECbmModuleId).\"";
+ return GetDetector(id);
+ }
+ LOG(debug) << GetName() << "::AddDetector(" << ToString(id) << ").";
+ fDetQa.emplace(id, id);
+ break;
+ default: LOG(warn) << GetName() << "::AddDetector : unsupported det " << ToString(id); return nullptr;
+ }
+
+ return &fDetQa[id];
+}
+
+//_____________________________________________________________________
+CbmRecoQaTask::Detector* CbmRecoQaTask::GetDetector(ECbmModuleId id)
+{
+ if (fDetQa.find(id) == fDetQa.end()) return nullptr;
+ return &fDetQa[id];
+}
+
+//_____________________________________________________________________
+const CbmTrack* CbmRecoQaTask::GetTrack(ECbmModuleId did, int id) const
+{
+ if (fTracks.find(did) == fTracks.end()) {
+ LOG(warning) << GetName() << " missing tracks for " << ToString(did);
+ return nullptr;
+ }
+ if (id < 0) {
+ LOG(warning) << GetName() << " negative trk idx for " << ToString(did);
+ return nullptr;
+ }
+ const CbmTrack* trk = (const CbmTrack*) (fTracks.at(did))->At(id);
+ if (!trk) {
+ LOG(debug1) << GetName() << " missing trk_" << id << " for " << ToString(did);
+ return nullptr;
+ }
+ return trk;
+}
+
//_____________________________________________________________________
InitStatus CbmRecoQaTask::Init()
{
@@ -34,301 +105,388 @@ InitStatus CbmRecoQaTask::Init()
return kERROR;
}
- fTracks = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrack"));
- if (!fTracks) {
- LOG(fatal) << GetName() << "::Init: Global track array not found!";
- return kERROR;
- }
+ fGTracks = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrack"));
+ if (!fGTracks)
+ LOG(warn) << GetName() << "::Init: Global track array not found!";
+ else
+ fuRecoConfig.set(kRecoTracks);
- 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.";
- }
+ fEvents = static_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
+ if (!fEvents)
+ LOG(warn) << GetName() << "::Init: No event found. Some results will be missing.";
+ else
+ fuRecoConfig.set(kRecoEvents);
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();
+ LOG(info) << GetName() << "::Init: Registered projections " << fProjs.size();
- int itrkLy(0), iprjLy(0);
- TFolder *modDir(nullptr), *smDir(nullptr);
- fOutFolder.Clear();
switch (fSetupClass) {
case kMcbm22:
LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2022\".";
- // STS U0 - station 0
- modDir = fOutFolder.AddFolder("STS", "STS reziduals to CA");
- 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++;
-
- // ===============================================================================
- // ToF Sm0 & Sm3 (Typ0) - station 5
- modDir = fOutFolder.AddFolder("ToF", "ToF reziduals to CA");
- 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);
- }
- itrkLy++;
-
- // ===============================================================================
- // TRG - upstream projections
- modDir = fOutFolder.AddFolder("Target", "Target tomography with CA");
- iprjLy = itrkLy;
- for (auto zprj : fProjView) {
- LOG(debug) << GetName() << "::Init: Project trks @ " << zprj.first;
- fModView[iprjLy].push_back(ModuleView());
- fModView[iprjLy][0].hdXx.first = 1;
- fModView[iprjLy][0].hdXx.second =
- new TH2D(Form("hxy%d", iprjLy - itrkLy), Form("Project @ z = %+5.2f; X_{trk} (cm); Y_{trk} (cm)", zprj.first),
- 700, -25, 10, 100, -5, 5);
- modDir->Add(fModView[iprjLy][0].hdXx.second);
- fModView[iprjLy][0].hdYy.first = 1;
- fModView[iprjLy][0].hdYy.second =
- new TH2D(Form("hxt%d", iprjLy), Form("Project @ z = %+5.2f; X_{trk} (cm); T_{trk} - T_{ev} (ns)", zprj.first),
- 700, -25, 10, 350, -15.5, 5.5);
- modDir->Add(fModView[iprjLy][0].hdYy.second);
- iprjLy++;
- }
+ InitMcbm22();
+ break;
+ case kMcbm24:
+ LOG(info) << GetName() << "::Init: Setup config for \"mCBM 2024\". ";
+ InitMcbm24();
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;
}
+ fOutFolder.Clear();
+ for (auto& detp : fDetQa) { // SYSTEM LOOP
+ auto& det = detp.second;
+
+ fHits[det.id] = static_cast<TClonesArray*>(ioman->GetObject(det.hit.name.data()));
+ if (!fHits[det.id]) {
+ LOG(warn) << GetName() << "::Init: Hit array for " << ToString(det.id) << " not found!";
+ continue;
+ }
+
+ if (fGTracks && det.trk.id != ECbmDataType::kUnknown) {
+ fTracks[det.id] = static_cast<TClonesArray*>(ioman->GetObject(det.trk.name.data()));
+ if (!fTracks[det.id]) LOG(warn) << GetName() << "::Init: Track array for " << ToString(det.id) << " not found!";
+ }
+ if (!det.Init(&fOutFolder)) continue;
+ }
+
+ // TDirectoryFile* trkDir = (TDirectoryFile*)fOutFolder.mkdir("track", "CA QA");
+ // for (auto& trkly : fTrkView) {
+ // for (auto& view : trkly.second) view.Register(trkDir);
+ // }
+
+ TDirectoryFile* trgDir = (TDirectoryFile*) fOutFolder.mkdir("target", "Target tomography with CA");
+ for (auto& prj : fProjs)
+ prj.second.Register(trgDir);
+
return kSUCCESS;
}
+/// Hit classification on system and view
+template<class Hit>
+bool CbmRecoQaTask::Detector::View::HasAddress(const CbmHit* h, double&, double&) const
+{
+ LOG(error) << "Unprocessed hit in view " << name;
+ cout << h->ToString();
+
+ return false;
+}
+// STS specialization
+template<>
+bool CbmRecoQaTask::Detector::View::HasAddress<CbmStsHit>(const CbmHit* h, double& x, double& y) const
+{
+ int32_t a = h->GetAddress();
+ uint32_t sel[] = {CbmStsAddress::GetElementId(a, EStsElementLevel::kStsUnit),
+ CbmStsAddress::GetElementId(a, EStsElementLevel::kStsLadder),
+ CbmStsAddress::GetElementId(a, EStsElementLevel::kStsModule)};
+ uint idx(0);
+ for (auto ii : fSelector) {
+ if (uint(ii) != sel[idx]) break;
+ idx++;
+ }
+ if (idx != fSelector.size())
+ return false;
+ else
+ LOG(debug4) << "Accept Sts hit for " << sel[0] << " " << sel[1] << " " << sel[2];
+
+ const CbmStsHit* h0 = dynamic_cast<const CbmStsHit*>(h);
+ if (!h0) {
+ LOG(error) << "Failed loading STS hit in view " << name;
+ return false;
+ }
+ x = h0->GetX();
+ y = h0->GetY();
+ return true;
+}
+// Rich specialization
+template<>
+bool CbmRecoQaTask::Detector::View::HasAddress<CbmRichHit>(const CbmHit* h, double& x, double& y) const
+{
+ int16_t uId = CbmRichUtil::GetDirichId(h->GetAddress());
+ // TODO implement at RichUtil level
+ int8_t modId = (uId >> 8) & 0xF;
+ if (modId >= fSelector[0]) return false;
+
+ const CbmRichHit* h0 = dynamic_cast<const CbmRichHit*>(h);
+ if (!h0) {
+ LOG(error) << "Failed loading RICH hit in view " << name;
+ return false;
+ }
+ x = h0->GetX();
+ y = h0->GetY();
+ return true;
+}
+// TRD specialization
+template<>
+bool CbmRecoQaTask::Detector::View::HasAddress<CbmTrdHit>(const CbmHit* h, double& x, double& y) const
+{
+ int16_t uId = CbmTrdAddress::GetModuleAddress(h->GetAddress());
+ if (uId != fSelector[0])
+ return false;
+ else
+ LOG(debug4) << "Accept Trd hit for " << uId;
+
+
+ const CbmTrdHit* h0 = dynamic_cast<const CbmTrdHit*>(h);
+ if (!h0) {
+ LOG(error) << "Failed loading TRD hit in view " << name;
+ return false;
+ }
+ x = h0->GetX();
+ y = h0->GetY();
+ return true;
+}
+// ToF specialization
+template<>
+bool CbmRecoQaTask::Detector::View::HasAddress<CbmTofHit>(const CbmHit* h, double& x, double& y) const
+{
+ int32_t a = h->GetAddress();
+ int32_t sel[] = {CbmTofAddress::GetSmId(a), CbmTofAddress::GetSmType(a), CbmTofAddress::GetRpcId(a)};
+ uint idx(0);
+ for (auto ii : fSelector) {
+ if (ii != sel[idx]) break;
+ idx++;
+ }
+ if (idx != fSelector.size())
+ return false;
+ else
+ LOG(debug4) << "Accept Tof hit for " << sel[0] << " " << sel[1] << " " << sel[2];
+
+ const CbmTofHit* h0 = dynamic_cast<const CbmTofHit*>(h);
+ if (!h0) {
+ LOG(error) << "Failed loading ToF hit in view " << name;
+ return false;
+ }
+ x = h0->GetX();
+ y = h0->GetY();
+ return true;
+}
+template<class Hit>
+bool CbmRecoQaTask::Detector::View::Load(const CbmHit* h, const CbmEvent* ev)
+{
+ double x(0), y(0);
+ if (!HasAddress<Hit>(h, x, y)) return false;
+
+ //printf("View[%s] z(%d)[%f]\n", name.data(), (int)h->GetType(), h->GetZ());
+ int scale(0);
+ TH2* hh(nullptr);
+ if (ev) {
+ if (fProjection.find(eViewProjection::kChdT) != fProjection.end()) {
+ scale = get<0>(fProjection[eViewProjection::kChdT]);
+ hh = get<2>(fProjection[eViewProjection::kChdT]);
+ hh->Fill(x, scale * (h->GetTime() - ev->GetStartTime()));
+ }
+ if (fProjection.find(eViewProjection::kXYh) != fProjection.end()) {
+ hh = get<2>(fProjection[eViewProjection::kXYh]);
+ hh->Fill(x, y);
+ }
+ }
+ if (fProjection.find(eViewProjection::kXYs) != fProjection.end()) {
+ hh = get<2>(fProjection[eViewProjection::kXYs]);
+ hh->Fill(x, y);
+ }
+ return true;
+}
+
+bool CbmRecoQaTask::Detector::View::Load(const CbmKfTrackFitter::FitNode* n)
+{
+ 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]);
+ // pullt = dt / sqrt(n->fMt.Dt2()[0] + n->fTrack.GetCovariance(5, 5)[0]);
+
+ for (auto& projection : fProjection) {
+ int scale = get<0>(projection.second);
+ TH2* hh = get<2>(projection.second);
+ if (!hh) continue;
+
+ switch (projection.first) {
+ case eViewProjection::kXYa: hh->Fill(n->fMxy.X()[0], n->fMxy.Y()[0]); break;
+ case eViewProjection::kXYp: hh->Fill(n->fTrack.X()[0], n->fTrack.Y()[0]); break;
+ case eViewProjection::kXdX: hh->Fill(n->fMxy.X()[0], scale * dx); break;
+ case eViewProjection::kYdY: hh->Fill(n->fMxy.Y()[0], scale * dy); break;
+ case eViewProjection::kWdT: hh->Fill(n->fMxy.X()[0], scale * dt); break;
+ case eViewProjection::kXpX: hh->Fill(n->fMxy.X()[0], pullx); break;
+ case eViewProjection::kYpY: hh->Fill(n->fMxy.Y()[0], pully); break;
+ default: break;
+ }
+ }
+ return true;
+}
+
//_____________________________________________________________________
void CbmRecoQaTask::Exec(Option_t*)
{
- LOG(debug) << GetName() << "::Exec : Evs[" << (fEvents ? fEvents->GetEntriesFast() : 0) << "] Trks["
- << fTracks->GetEntriesFast() << "]";
- int itrack(0), nnodes(0);
- for (auto trackObj : *fTracks) {
- //printf("GlobalTrack %d\n", itrack);
+ LOG(info) << GetName() << "::Exec : Evs[" << (fEvents ? fEvents->GetEntriesFast() : 0) << "] Trks["
+ << (fGTracks ? fGTracks->GetEntriesFast() : 0) << "]";
+
+ int iev(0), itrack(0), nnodes(0);
+ if (fEvents) {
+ for (auto evObj : *fEvents) { // EVENT LOOP
+ auto ev = dynamic_cast<CbmEvent*>(evObj);
+ // track QA filtering
+ if (!FilterEvent(ev)) continue;
+
+ for (auto detp : fDetQa) { // SYSTEM LOOP
+ auto det = detp.second;
+ // check SYSTEM @ hit level
+ int nh = max(int(0), int(ev->GetNofData(det.hit.id)));
+
+ if (!fHits[det.id]) {
+ LOG(error) << GetName() << "::Exec() : Hits for " << ToString(det.id) << " not available. Skip.";
+ continue;
+ }
+ for (int ih(0); ih < nh; ih++) { // HIT LOOP
+ const int jh = ev->GetIndex(det.hit.id, ih);
+ const CbmHit* hit = dynamic_cast<CbmHit*>(fHits[det.id]->At(jh));
+ // printf("%s hit_%2d dt=%f time_event=%f time_hit=%f\n", ToString(det.id).data(), ih, ev->GetStartTime() - hit->GetTime(), ev->GetStartTime(), hit->GetTime());
+
+ bool ret(false);
+ for (auto view : det.fViews) { // VIEW LOOP
+ switch (det.id) {
+ case ECbmModuleId::kMvd: ret = view.Load<CbmMvdHit>(hit, ev); break;
+ case ECbmModuleId::kSts: ret = view.Load<CbmStsHit>(hit, ev); break;
+ case ECbmModuleId::kMuch: ret = view.Load<CbmMuchPixelHit>(hit, ev); break;
+ case ECbmModuleId::kRich: ret = view.Load<CbmRichHit>(hit, ev); break;
+ case ECbmModuleId::kTrd: ret = view.Load<CbmTrdHit>(hit, ev); break;
+ case ECbmModuleId::kTrd2d: ret = view.Load<CbmTrdHit>(hit, ev); break;
+ case ECbmModuleId::kTof: ret = view.Load<CbmTofHit>(hit, ev); break;
+ case ECbmModuleId::kFsd: ret = view.Load<CbmFsdHit>(hit, ev); break;
+ case ECbmModuleId::kPsd: ret = view.Load<CbmPsdHit>(hit, ev); break;
+ default: LOG(fatal) << GetName() << "::Exec : unsupported det " << ToString(det.id); break;
+ }
+ if (ret) break;
+ } // END VIEW LOOP
+ } // END HIT LOOP
+ } // END SYSTEM LOOP
+
+ if (!fGTracks) continue;
+
+ int ntrk = max(int(0), int(ev->GetNofData(ECbmDataType::kGlobalTrack)));
+ for (int itrk(0); itrk < ntrk; itrk++) { // TRACK LOOP
+ auto trackObj = (CbmGlobalTrack*) (*fGTracks)[ev->GetIndex(ECbmDataType::kGlobalTrack, itrk)];
+
+ // track QA filtering
+ if (!FilterTrack(trackObj)) continue;
+
+ 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;
+ }
+ trkKf.fMsQp0 = 1. / 0.5;
+ trkKf.fIsMsQp0Set = true;
+
+ // minimalistic QA tool for tracks used for target projection
+ int nhit[(size_t) ECbmModuleId::kLastModule] = {0};
+ for (auto& n : trkKf.fNodes) {
+ // check if hit is well defined
+ if (n.fHitSystemId == ECbmModuleId::kNotExist) continue;
+ // check if detector is registered
+ if (fDetQa.find(n.fHitSystemId) == fDetQa.end()) continue;
+
+ auto det = fDetQa[n.fHitSystemId];
+ // det.Print();
+ auto view = det.FindView(n.fMxy.X()[0], n.fMxy.Y()[0], n.fZ);
+ if (!view) {
+ LOG(debug) << GetName() << "::Exec: view for tracking layer " << ToString(det.id) << " not defined.";
+ continue;
+ }
+
+ n.fIsTimeSet = false;
+ n.fIsXySet = false;
+ trkKf.MakeConsistent();
+ fFitter.FitTrack(trkKf);
+ view->Load(&n);
+ nhit[(int) n.fHitSystemId]++;
+ n.fIsTimeSet = true;
+ n.fIsXySet = true;
+ }
+
+ // Fit track for intermediate projections
+ int iprj(0);
+ if (!nnodes) nnodes = (int) trkKf.fNodes.size();
+ for (auto zprj : fProjs) {
+ 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;
+ if (fProjs.find(n.fReference1 + nnodes) == fProjs.end()) {
+ LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fReference1 + nnodes << " not defined.";
+ continue;
+ }
+
+ // View* mod = &fProjs[n.fReference1 + nnodes][0];
+ // mod->Load(&n);
+ //
+ // if (mod->hdXx.second)
+ // mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X()[0], mod->hdXx.first * n.fTrack.Y()[0]);
+ // if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first * n.fTrack.X()[0], n.fTrack.Time()[0]);
+ }
+ itrack++;
+ } // END TRACK LOOP
+ iev++;
+ } // END EVENT LOOP
+
+ LOG(info) << GetName() << "::Exec : Evs(%)[" << 1.e2 * iev / fEvents->GetEntriesFast() << "] Trks(%)["
+ << (fGTracks ? 1.e3 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
+ return; // END EbyE QA
+ }
+
+ // FULL TS LOOP
+ for (auto detp : fDetQa) { // SYSTEM LOOP
+ auto det = detp.second;
+ // check SYSTEM @ hit level
+ if (!fHits[det.id]) {
+ LOG(error) << GetName() << "::Exec() : Hits for " << ToString(det.id) << " not available. Skip.";
+ continue;
+ }
+ int nh = fHits[det.id]->GetEntriesFast();
+
+ for (int ih(0); ih < nh; ih++) { // HIT LOOP
+ const CbmHit* hit = dynamic_cast<CbmHit*>(fHits[det.id]->At(ih));
+ // printf("%s hit_%2d dt=%f time_event=%f time_hit=%f\n", ToString(det.id).data(), ih, ev->GetStartTime() - hit->GetTime(), ev->GetStartTime(), hit->GetTime());
+
+ bool ret(false);
+ for (auto view : det.fViews) { // VIEW LOOP
+ switch (det.id) {
+ case ECbmModuleId::kMvd: ret = view.Load<CbmMvdHit>(hit); break;
+ case ECbmModuleId::kSts: ret = view.Load<CbmStsHit>(hit); break;
+ case ECbmModuleId::kMuch: ret = view.Load<CbmMuchPixelHit>(hit); break;
+ case ECbmModuleId::kRich: ret = view.Load<CbmRichHit>(hit); break;
+ case ECbmModuleId::kTrd: ret = view.Load<CbmTrdHit>(hit); break;
+ case ECbmModuleId::kTrd2d: ret = view.Load<CbmTrdHit>(hit); break;
+ case ECbmModuleId::kTof: ret = view.Load<CbmTofHit>(hit); break;
+ case ECbmModuleId::kFsd: ret = view.Load<CbmFsdHit>(hit); break;
+ case ECbmModuleId::kPsd: ret = view.Load<CbmPsdHit>(hit); break;
+ default: LOG(fatal) << GetName() << "::Exec : unsupported det " << ToString(det.id); break;
+ }
+ if (ret) break;
+ } // END VIEW LOOP
+ } // END HIT LOOP
+ } // END SYSTEM LOOP
+ if (!fGTracks) {
+ LOG(info) << GetName() << "::Exec : TS local reco only.";
+ return; // END TS QA (no tracking)
+ }
+
+ for (auto trackObj : *fGTracks) { // TRACK LOOP
auto track = dynamic_cast<const CbmGlobalTrack*>(trackObj);
+
+ // track QA filtering
+ if (!FilterTrack(track)) continue;
+
CbmKfTrackFitter::Track trkKf;
if (!fFitter.CreateGlobalTrack(trkKf, *track)) {
LOG(fatal) << GetName() << "::Exec: can not create the track for the fit! ";
@@ -341,66 +499,32 @@ void CbmRecoQaTask::Exec(Option_t*)
// minimalistic QA tool for tracks used for target projection
int nhit[(size_t) ECbmModuleId::kLastModule] = {0};
for (auto& n : trkKf.fNodes) {
+ // check if hit is well defined
if (n.fHitSystemId == ECbmModuleId::kNotExist) continue;
+ // check if detector is registered
+ if (fDetQa.find(n.fHitSystemId) == fDetQa.end()) continue;
+
+ auto det = fDetQa[n.fHitSystemId];
+ // det.Print();
+ auto view = det.FindView(n.fMxy.X()[0], n.fMxy.Y()[0], n.fZ);
+ if (!view) {
+ LOG(debug) << GetName() << "::Exec: view for tracking layer " << ToString(det.id) << " not defined.";
+ continue;
+ }
+
n.fIsTimeSet = false;
n.fIsXySet = false;
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]);
- // pullt = dt / sqrt(n.fMt.Dt2()[0] + n.fTrack.GetCovariance(5, 5)[0]);
-
- // add granularity for specific systems
- int modId = 0;
- switch (n.fHitSystemId) {
- case ECbmModuleId::kTof:
- modId = CbmTofAddress::GetRpcId(n.fHitAddress);
- nhit[(int) n.fHitSystemId]++;
- break;
- default: nhit[(int) n.fHitSystemId]++; break;
- }
-
- // check that we have correctly created all views
- if (fModView.find(n.fMaterialLayer) == fModView.end()) {
- LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fMaterialLayer << " not defined.";
- continue;
- }
- if (fModView[n.fMaterialLayer].size() <= (size_t) modId) {
- LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fMaterialLayer << " @ " << modId
- << " not defined.";
- continue;
- }
- ModuleView* mod = &fModView[n.fMaterialLayer][modId];
- if (!mod) {
- LOG(warn) << GetName() << "::Exec: can not access view for tracking layer " << n.fMaterialLayer << " @ "
- << 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);
+ view->Load(&n);
+ nhit[(int) n.fHitSystemId]++;
n.fIsTimeSet = true;
n.fIsXySet = true;
}
- // track QA filtering
- if (!FilterTrack(nhit)) continue;
-
// Fit track with all hits ON for vx projection
int iprj(0);
- for (auto zprj : fProjView) {
+ for (auto zprj : fProjs) {
CbmKfTrackFitter::FitNode n = CbmKfTrackFitter::FitNode();
n.fZ = zprj.first;
n.fIsXySet = true;
@@ -411,49 +535,1041 @@ void CbmRecoQaTask::Exec(Option_t*)
fFitter.FitTrack(trkKf);
for (auto& n : trkKf.fNodes) {
if (n.fReference1 < 0) continue;
- if (fModView.find(n.fReference1 + nnodes) == fModView.end()) {
+ if (fProjs.find(n.fReference1 + nnodes) == fProjs.end()) {
LOG(debug) << GetName() << "::Exec: view for tracking layer " << n.fReference1 + nnodes << " not defined.";
continue;
}
- ModuleView* mod = &fModView[n.fReference1 + nnodes][0];
- if (mod->hdXx.second)
- mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X()[0], mod->hdXx.first * n.fTrack.Y()[0]);
- if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first * n.fTrack.X()[0], n.fTrack.Time()[0]);
+ // View* mod = &fProjs[n.fReference1 + nnodes][0];
+ // mod->Load(&n);
+ // if (mod->hdXx.second)
+ // mod->hdXx.second->Fill(mod->hdXx.first * n.fTrack.X()[0], mod->hdXx.first * n.fTrack.Y()[0]);
+ // if (mod->hdYy.second) mod->hdYy.second->Fill(mod->hdYy.first * n.fTrack.X()[0], n.fTrack.Time()[0]);
}
itrack++;
- }
+ } // END TRACK LOOP
+ LOG(info) << GetName() << "::Exec : Trks(%)[" << (fGTracks ? 1.e3 * itrack / fGTracks->GetEntriesFast() : 0) << "]";
+
+ return; // END TS QA (tracking)
}
//_____________________________________________________________________
void CbmRecoQaTask::SetProjections(std::vector<double> vzpj)
{
- fProjView.clear();
+ fProjs.clear();
for (auto zprj : vzpj) {
- fProjView[zprj] = new TH2D(Form("hxx"), Form("; X (cm); Y (cm)"), 450, -25, 20, 450, -25, 20);
+ // fProjView[zprj] = new TH2D(Form("hxx"), Form("; X (cm); Y (cm)"), 450, -25, 20, 450, -25, 20);
+ fProjs.emplace(zprj, Detector::View(Form("z%.2f", zprj), "", {}));
}
}
//_____________________________________________________________________
-bool CbmRecoQaTask::FilterTrack(int* ptr)
+void CbmRecoQaTask::Finish()
{
- bool ret = false;
- switch (fSetupClass) {
- case kMcbm22: {
- int* nhit = ptr;
- if (nhit[int(ECbmModuleId::kSts)] == 2 && nhit[int(ECbmModuleId::kTrd2d)] && nhit[int(ECbmModuleId::kTof)])
- ret = true;
+ FairSink* sink = FairRootManager::Instance()->GetSink();
+ sink->WriteObject(&fOutFolder, nullptr);
+}
+
+//_____________________________________________________________________
+bool CbmRecoQaTask::FilterEvent(const CbmEvent* ptr)
+{
+ for (auto evCut : fFilterEv) {
+ if (!evCut.Accept(ptr, this)) return false;
+ }
+ return true;
+}
+
+//_____________________________________________________________________
+bool CbmRecoQaTask::FilterTrack(const CbmGlobalTrack* ptr)
+{
+ for (auto trkCut : fFilterTrk) {
+ if (!trkCut.Accept(ptr, this)) return false;
+ }
+ return true;
+}
+
+//_____________________________________________________________________
+void CbmRecoQaTask::InitMcbm22()
+{
+ CbmSetup* setup = CbmSetup::Instance();
+ if (!setup) {
+ LOG(fatal) << GetName() << "::InitMcbm22() : Missing setup definition.";
+ return;
+ }
+
+ TString dtag;
+ Detector::View* v(nullptr);
+ if (setup->IsActive(ECbmModuleId::kSts)) {
+ dtag = "v22f_mcbm";
+ Detector* sts = AddDetector(ECbmModuleId::kSts);
+ // U0 L0 M0
+ v = sts->AddView(
+ "U0L0M0",
+ Form("/cave_1/sts_%s_0/Station01_1/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", dtag.Data()),
+ {0, 0, 0});
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eViewProjection::kYdY, 4, "mm");
+ v = sts->AddView(
+ "U0L0M1",
+ Form("/cave_1/sts_%s_0/Station01_1/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", dtag.Data()),
+ {0, 0, 1});
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eViewProjection::kYdY, 4, "mm");
+ // U0 L1
+ v = sts->AddView(
+ "U0L1M0",
+ Form("/cave_1/sts_%s_0/Station01_1/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", dtag.Data()),
+ {0, 1, 0});
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eViewProjection::kYdY, 4, "mm");
+ v = sts->AddView(
+ "U0L1M1",
+ Form("/cave_1/sts_%s_0/Station01_1/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", dtag.Data()),
+ {0, 1, 1});
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v->SetProjection(eViewProjection::kYdY, 4, "mm");
+ // U1 L0
+ v = sts->AddView(
+ "U1L0M0",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module03_1/Sensor03_1", dtag.Data()),
+ {1, 0, 0});
+ v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
+ v->SetProjection(eViewProjection::kYdY, 3, "mm");
+ v = sts->AddView(
+ "U1L0M1",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module04_2/Sensor04_1", dtag.Data()),
+ {1, 0, 1});
+ v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
+ v->SetProjection(eViewProjection::kYdY, 3, "mm");
+ // U1 L1
+ v = sts->AddView(
+ "U1L1M0",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module03_1/Sensor03_1", dtag.Data()),
+ {1, 1, 0});
+ v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
+ v->SetProjection(eViewProjection::kYdY, 3, "mm");
+ v = sts->AddView(
+ "U1L1M1",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module04_2/Sensor04_1", dtag.Data()),
+ {1, 1, 1});
+ v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
+ v->SetProjection(eViewProjection::kYdY, 3, "mm");
+ // U1 L2
+ v = sts->AddView(
+ "U1L2M0",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_1/Sensor03_1", dtag.Data()),
+ {1, 2, 0});
+ v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
+ v->SetProjection(eViewProjection::kYdY, 3, "mm");
+ v = sts->AddView(
+ "U1L2M1",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_2/Sensor03_1", dtag.Data()),
+ {1, 2, 1});
+ v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
+ v->SetProjection(eViewProjection::kYdY, 3, "mm");
+ v = sts->AddView(
+ "U1L2M2",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_3/Sensor03_1", dtag.Data()),
+ {1, 2, 2});
+ v->SetProjection(eViewProjection::kXdX, 0.75, "mm");
+ v->SetProjection(eViewProjection::kYdY, 3, "mm");
+ }
+ if (setup->IsActive(ECbmModuleId::kTrd)) {
+ dtag = "v22h_mcbm";
+ Detector* trd = AddDetector(ECbmModuleId::kTrd2d);
+ // Trd2D
+ v = trd->AddView("2D", Form("/cave_1/trd_%s_0/layer01_20101/module9_101001001/gas_0", dtag.Data()), {5});
+ v->SetProjection(eViewProjection::kChdT, 400, "ns");
+ v->SetProjection(eViewProjection::kXdX, 10, "mm");
+ v->SetProjection(eViewProjection::kYdY, 20, "mm");
+ // Trd1Dx
+ trd = AddDetector(ECbmModuleId::kTrd);
+ v = trd->AddView("1Dx", Form("/cave_1/trd_%s_0/layer02_10202/module8_101002001", dtag.Data()), {21});
+ v->SetProjection(eViewProjection::kChdT, 400, "ns");
+ v->SetProjection(eViewProjection::kXdX, 1.5, "cm");
+ // Trd1Dy
+ v = trd->AddView("1Dy", Form("/cave_1/trd_%s_0/layer03_11303/module8_101303001", dtag.Data()), {37});
+ v->SetProjection(eViewProjection::kChdT, 400, "ns");
+ }
+ if (setup->IsActive(ECbmModuleId::kTof)) {
+ dtag = "v21k_mcbm";
+ Detector* tof = AddDetector(ECbmModuleId::kTof);
+ tof->hit.name = "TofCalHit";
+ vector<int> tofSelect(3);
+ // add type 0
+ tofSelect[1] = 0; // ToF type
+ for (int ism(0); ism < 4; ism++) {
+ tofSelect[0] = ism;
+ for (int irpc(0); irpc < 5; irpc++) {
+ tofSelect[2] = irpc;
+ v = tof->AddView(Form("Sm%dRpc%d", tofSelect[0], tofSelect[2]),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(),
+ dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
+ tofSelect);
+ v->SetProjection(eViewProjection::kChdT, 15, "ns");
+ }
+ }
+ // add type 2
+ tofSelect[1] = 2; // ToF type
+ tofSelect[0] = 0; // Tof SM id
+ for (int irpc(0); irpc < 5; irpc++) {
+ tofSelect[2] = irpc;
+ tof->AddView(Form("Sm_2%dRpc%d", tofSelect[0], tofSelect[2]),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(), dtag.Data(),
+ tofSelect[1], tofSelect[0], tofSelect[2]),
+ tofSelect);
+ }
+ }
+
+ // ===============================================================================
+ // TRG - upstream projections
+ // int iprjLy = 0;
+ // for (auto prj : fProjs) {
+ // LOG(debug) << GetName() << "::Init: Project trks @ " << prj.first;
+ // prj.second.name = Form("Z=%+4.1f", prj.first);
+ // prj.second.hdXx.first = 1;
+ // prj.second.hdXx.second =
+ // new TH2D(Form("Hxy%d", iprjLy), Form("Project @ z = %+5.2f; X_{trk} (cm); Y_{trk} (cm)", prj.first),
+ // 700, -25, 10, 100, -5, 5);
+ // prj.second.hdT.first = 1;
+ // prj.second.hdT.second =
+ // new TH2D(Form("Hxt%d", iprjLy), Form("Project @ z = %+5.2f; X_{trk} (cm); T_{trk} - T_{ev} (ns)", prj.first),
+ // 700, -25, 10, 350, -15.5, 5.5);
+ // iprjLy++;
+ // }
+}
+
+//_____________________________________________________________________
+void CbmRecoQaTask::InitMcbm24()
+{
+ CbmSetup* setup = CbmSetup::Instance();
+ if (!setup) {
+ LOG(fatal) << GetName() << "::InitMcbm24() : Missing setup definition.";
+ return;
+ }
+ TString dtag;
+ Detector::View* v(nullptr);
+ if (setup->IsActive(ECbmModuleId::kSts)) {
+ // setup->GetGeoTag(ECbmModuleId::kSts, dtag);
+ // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
+ dtag = "v24c_mcbm";
+ Detector* sts = AddDetector(ECbmModuleId::kSts);
+ // U0 L0 M0
+ v = sts->AddView(
+ "U0L0M0",
+ Form("/cave_1/sts_%s_0/Station01_1/Ladder13_1/HalfLadder13u_1/HalfLadder13u_Module03_1/Sensor03_1", dtag.Data()),
+ {0, 0, 0});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ // U1 L0 M0
+ v = sts->AddView(
+ "U1L0M0",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", dtag.Data()),
+ {1, 0, 0});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v = sts->AddView(
+ "U1L0M1",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder09_1/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", dtag.Data()),
+ {1, 0, 1});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ // U1 L1
+ v = sts->AddView(
+ "U1L1M0",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_1/Sensor03_1", dtag.Data()),
+ {1, 1, 0});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v = sts->AddView(
+ "U1L1M1",
+ Form("/cave_1/sts_%s_0/Station02_2/Ladder09_2/HalfLadder09d_2/HalfLadder09d_Module03_2/Sensor03_1", dtag.Data()),
+ {1, 1, 1});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ // U2 L0
+ v = sts->AddView(
+ "U2L0M0",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module03_1/Sensor03_1", dtag.Data()),
+ {2, 0, 0});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v = sts->AddView(
+ "U2L0M1",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder10_1/HalfLadder10d_2/HalfLadder10d_Module04_2/Sensor04_1", dtag.Data()),
+ {2, 0, 1});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ // U2 L1
+ v = sts->AddView(
+ "U2L1M0",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module03_1/Sensor03_1", dtag.Data()),
+ {2, 1, 0});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v = sts->AddView(
+ "U2L1M1",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder12_2/HalfLadder12d_2/HalfLadder12d_Module04_2/Sensor04_1", dtag.Data()),
+ {2, 1, 1});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ // U2 L2
+ v = sts->AddView(
+ "U2L2M0",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_1/Sensor03_1", dtag.Data()),
+ {2, 2, 0});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ v = sts->AddView(
+ "U2L2M1",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_2/Sensor03_1", dtag.Data()),
+ {2, 2, 1});
+ v = sts->AddView(
+ "U2L2M2",
+ Form("/cave_1/sts_%s_0/Station03_3/Ladder11_3/HalfLadder11d_2/HalfLadder11d_Module03_3/Sensor03_1", dtag.Data()),
+ {2, 2, 2});
+ v->SetProjection(eViewProjection::kYdY, 5, "mm");
+ v->SetProjection(eViewProjection::kXdX, 1, "mm");
+ }
+ if (setup->IsActive(ECbmModuleId::kTrd)) {
+ // setup->GetGeoTag(ECbmModuleId::kTrd, dtag);
+ // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
+ dtag = "v24e_mcbm";
+ Detector* trd = AddDetector(ECbmModuleId::kTrd2d);
+ // Trd2D
+ v = trd->AddView("2D", Form("/cave_1/trd_%s_0/layer01_20101/module9_101001001/gas_0", dtag.Data()), {5});
+ v->SetProjection(eViewProjection::kChdT, 300, "ns");
+ // Trd1Dx
+ trd = AddDetector(ECbmModuleId::kTrd);
+ v = trd->AddView("1Dx", Form("/cave_1/trd_%s_0/layer02_10202/module5_101002001", dtag.Data()), {21});
+ v->SetProjection(eViewProjection::kChdT, 300, "ns");
+ v->SetProjection(eViewProjection::kXdX, 7, "cm");
+ // Trd1Dy
+ v = trd->AddView("1Dy", Form("/cave_1/trd_%s_0/layer03_11303/module5_101103001", dtag.Data()), {37});
+ v->SetProjection(eViewProjection::kChdT, 300, "ns");
+ v->SetProjection(eViewProjection::kXdX, 7, "cm");
+ }
+ if (setup->IsActive(ECbmModuleId::kTof)) {
+ // setup->GetGeoTag(ECbmModuleId::kTof, dtag);
+ // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
+ dtag = "v24d_mcbm";
+ Detector* tof = AddDetector(ECbmModuleId::kTof);
+ vector<int> tofSelect(3);
+ // add type 0
+ for (int ism(0); ism < 6; ism++) {
+ tofSelect[0] = ism;
+ for (int irpc(0); irpc < 5; irpc++) {
+ tofSelect[2] = irpc;
+ v = tof->AddView(Form("Sm%dRpc%d", ism, irpc),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(),
+ dtag.Data(), tofSelect[1], tofSelect[0], tofSelect[2]),
+ tofSelect);
+ v->SetProjection(eViewProjection::kXdX, 7, "cm");
+ v->SetProjection(eViewProjection::kYdY, 7, "cm");
+ }
+ }
+ // add type 6 (Buch)
+ tofSelect[0] = 0;
+ tofSelect[1] = 6;
+ for (int irpc(0); irpc < 2; irpc++) {
+ tofSelect[2] = irpc;
+ tof->AddView(Form("BuchRpc%d", irpc),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(), dtag.Data(),
+ tofSelect[1], tofSelect[0], tofSelect[2]),
+ tofSelect);
+ }
+ // add type 9
+ tofSelect[1] = 9;
+ for (int ism(0); ism < 2; ism++) {
+ tofSelect[0] = ism;
+ for (int irpc(0); irpc < 2; irpc++) {
+ tofSelect[2] = irpc;
+ tof->AddView(Form("Test%dRpc%d", ism, irpc),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(), dtag.Data(),
+ tofSelect[1], tofSelect[0], tofSelect[2]),
+ tofSelect);
+ }
+ }
+ // add type 2
+ tofSelect[1] = 2;
+ for (int ism(0); ism < 2; ism++) {
+ tofSelect[0] = ism;
+ for (int irpc(0); irpc < 5; irpc++) {
+ tofSelect[2] = irpc;
+ tof->AddView(Form("Sm2%dRpc%d", ism, irpc),
+ Form("/cave_1/tof_%s_0/tof_%sStand_1/module_%d_%d/gas_box_0/counter_%d", dtag.Data(), dtag.Data(),
+ tofSelect[1], tofSelect[0], tofSelect[2]),
+ tofSelect);
+ }
+ }
+ }
+ if (setup->IsActive(ECbmModuleId::kRich)) {
+ // setup->GetGeoTag(ECbmModuleId::kRich, dtag);
+ // LOG(debug) << GetName() << "::InitMcbm24() : found " << dtag;
+ dtag = "v24a_mcbm";
+ Detector* rich = AddDetector(ECbmModuleId::kRich);
+ // Aerogel 1
+ rich->AddView("Aerogel", Form("/cave_1/rich_%s_0/box_1/Gas_1", dtag.Data()), {4});
+ }
+}
+
+//========= DETECTOR ======================
+CbmRecoQaTask::Detector::Detector(ECbmModuleId did)
+{
+ id = did;
+ switch (id) {
+ case ECbmModuleId::kMvd: new (&hit) Data(ECbmDataType::kMvdHit, "MvdHit"); break;
+ case ECbmModuleId::kSts:
+ new (&hit) Data(ECbmDataType::kStsHit, "StsHit");
+ new (&trk) Data(ECbmDataType::kStsTrack, "StsTrack");
+ break;
+ case ECbmModuleId::kMuch: new (&hit) Data(ECbmDataType::kMuchPixelHit, "MuchHit"); break;
+ case ECbmModuleId::kRich: new (&hit) Data(ECbmDataType::kRichHit, "RichHit"); break;
+ case ECbmModuleId::kTrd:
+ case ECbmModuleId::kTrd2d:
+ new (&hit) Data(ECbmDataType::kTrdHit, "TrdHit");
+ new (&trk) Data(ECbmDataType::kTrdTrack, "TrdTrack");
+ break;
+ case ECbmModuleId::kTof:
+ new (&hit) Data(ECbmDataType::kTofHit, "TofHit");
+ new (&trk) Data(ECbmDataType::kTofTrack, "TofTrack");
break;
+ case ECbmModuleId::kFsd: new (&hit) Data(ECbmDataType::kFsdHit, "FsdHit"); break;
+ case ECbmModuleId::kPsd: new (&hit) Data(ECbmDataType::kPsdHit, "PsdHit"); break;
+ default:
+ LOG(warn) << "QA unsupported for Detector=" << ToString(did);
+ id = ECbmModuleId::kNotExist;
+ break;
+ }
+}
+
+//+++++++++++++++++++++ Detector ++++++++++++++++++++++++++++++++
+//_________________________________________________________________
+CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::AddView(const char* n, const char* p, std::vector<int> set)
+{
+ View* v(nullptr);
+ fViews.emplace_back(n, p, set);
+ v = &fViews.back();
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kXYs);
+
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kChdT);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kXYh);
+
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kXYa);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kXYp);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kXdX);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kYdY);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kWdT);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kXpX);
+ v->AddProjection(CbmRecoQaTask::eViewProjection::kYpY);
+
+ return v;
+}
+
+CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::GetView(const char* n)
+{
+ for (auto& view : fViews)
+ if (view.name.compare(n) == 0) return &view;
+
+ return nullptr;
+}
+
+CbmRecoQaTask::Detector::View* CbmRecoQaTask::Detector::FindView(double x, double y, double z)
+{
+ for (auto& v : fViews) {
+ // printf("Looking for p[%f %f %f] in V[%s-%s] @ %f %f %f (%f %f %f)\n", x, y, z, ToString(id).data(), v.name.data(), v.pos[0], v.pos[1], v.pos[2], v.size[0], v.size[1], v.size[2]);
+ if (abs(v.pos[0] - x) > 0.5 * v.size[0]) continue;
+ if (abs(v.pos[1] - y) > 0.5 * v.size[1]) continue;
+ if (abs(v.pos[2] - z) > 0.5 * v.size[2]) continue;
+ return &v;
+ }
+
+ return nullptr;
+}
+
+bool CbmRecoQaTask::Detector::Init(TDirectoryFile* fOut)
+{
+ /* Query the geometry and trigger detector views init
+ */
+ if (!gGeoManager) {
+ LOG(fatal) << "CbmRecoQaTask::Detector::Init() " << ToString(id) << " missing geometry.";
+ return false;
+ }
+ LOG(debug) << "CbmRecoQaTask::Detector::Init() : " << ToString(id);
+
+ TDirectoryFile* modDir =
+ (TDirectoryFile*) fOut->mkdir(ToString(id).data(), Form("Reco QA for %s", ToString(id).data()));
+ bool ret(true);
+ for (auto& view : fViews) {
+ bool vret = view.Init(ToString(id).data());
+ view.Register(modDir);
+ ret = ret && vret;
+ }
+ return ret;
+}
+
+void CbmRecoQaTask::Detector::Print() const
+{
+ stringstream s;
+ s << "D[" << ToString(id) << "] views[" << fViews.size() << "]\n";
+ for (auto v : fViews)
+ s << v.ToString();
+ cout << s.str();
+}
+
+//+++++++++++++++++++++ Detector::View ++++++++++++++++++++++++++++++++
+//_______________________________________________________________________
+bool CbmRecoQaTask::Detector::View::AddProjection(eViewProjection prj, float range, const char* unit)
+{
+ if (fProjection.find(prj) != fProjection.end()) {
+ LOG(warn) << "Projection " << ToString(prj) << " already registered";
+ return false;
+ }
+ int scale(1);
+ if (strcmp(unit, "mm") == 0)
+ scale = 10;
+ else if (strcmp(unit, "um") == 0)
+ scale = 10000;
+ else if (strcmp(unit, "ps") == 0)
+ scale = 1000;
+ else if (strcmp(unit, "eV") == 0)
+ scale = 1000;
+ else if (strcmp(unit, "cm") == 0)
+ scale = 1;
+ else if (strcmp(unit, "ns") == 0)
+ scale = 1;
+ else if (strcmp(unit, "keV") == 0)
+ scale = 1;
+ else if (strcmp(unit, "a.u.") == 0)
+ scale = 1;
+ else
+ LOG(warn) << "Projection units " << unit << " not registered. Natural units will be used.";
+
+ fProjection[prj] = make_tuple(scale, range, nullptr);
+ return true;
+}
+
+bool CbmRecoQaTask::Detector::View::SetProjection(eViewProjection prj, float range, const char* unit)
+{
+ if (fProjection.find(prj) == fProjection.end()) {
+ LOG(warn) << "Projection " << ToString(prj)
+ << " not initialized. Calling \"CbmRecoQaTask::Detector::View::AddProjection()\"";
+ return AddProjection(prj, range, unit);
+ }
+ int scale(1);
+ if (strcmp(unit, "mm") == 0)
+ scale = 10;
+ else if (strcmp(unit, "um") == 0)
+ scale = 10000;
+ else if (strcmp(unit, "ps") == 0)
+ scale = 1000;
+ else if (strcmp(unit, "eV") == 0)
+ scale = 1000;
+ else if (strcmp(unit, "cm") == 0)
+ scale = 1;
+ else if (strcmp(unit, "ns") == 0)
+ scale = 1;
+ else if (strcmp(unit, "keV") == 0)
+ scale = 1;
+ else if (strcmp(unit, "a.u.") == 0)
+ scale = 1;
+ else
+ LOG(warn) << "Projection units " << unit << " not registered. Natural units will be used.";
+
+ get<0>(fProjection[prj]) = scale;
+ get<1>(fProjection[prj]) = range;
+ return true;
+}
+
+bool CbmRecoQaTask::Detector::View::Init(const char* dname)
+{
+ if (!gGeoManager->cd(path.data())) return false;
+ TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
+ const double* tr(m->GetTranslation());
+ TGeoVolume* v = gGeoManager->GetCurrentVolume();
+ TGeoShape* bb = v->GetShape();
+ double w_lo, w_hi;
+ for (int i(0); i < 3; i++) { // loop over the spatial dimensions
+ size[i] = bb->GetAxisRange(i + 1, w_lo, w_hi);
+ pos[i] = 0.5 * (w_lo + w_hi) + tr[i];
+ }
+
+ LOG(info) << "CbmRecoQaTask::Detector(" << dname << ")::View(" << name << ")::Init() : size [" << size[0] << " x "
+ << size[1] << " x " << size[2] << "].";
+
+ // TODO short-hand for XY display resolution
+ int dscale = 10;
+ if (strcmp(dname, "Tof") == 0 || strcmp(dname, "Rich") == 0) dscale = 1;
+
+ int nbinsx, nbinsy;
+ string unit;
+ for (auto& projection : fProjection) {
+ int scale = get<0>(projection.second);
+ float yrange = get<1>(projection.second);
+ char xy_id = 0;
+ double xlo = pos[0] - 0.5 * size[0], xhi = pos[0] + 0.5 * size[0], ylo = pos[1] - 0.5 * size[1],
+ yhi = pos[1] + 0.5 * size[1];
+ switch (projection.first) {
+ case eViewProjection::kXdX:
+ nbinsx = 10 * ceil(size[0]); // mm binning
+ unit = makeYrange(scale, yrange);
+ nbinsy = 200; //* ceil(yrange);
+ get<2>(projection.second) = new TH2D(Form("hxx_%s%s", dname, name.data()),
+ Form("X resolution %s [%s]; X_{%s-%s} (cm); #Delta X (%s)", name.data(),
+ dname, dname, name.data(), unit.data()),
+ nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
+ break;
+
+ case eViewProjection::kYdY:
+ nbinsx = 10 * ceil(size[1]); // mm binning
+ unit = makeYrange(scale, yrange);
+ nbinsy = 200; //* ceil(yrange);
+ get<2>(projection.second) = new TH2D(Form("hyy_%s%s", dname, name.data()),
+ Form("Y resolution %s [%s]; Y_{%s-%s} (cm); #Delta Y (%s)", name.data(),
+ dname, dname, name.data(), unit.data()),
+ nbinsx, ylo, yhi, nbinsy, -yrange, yrange);
+ break;
+
+ case eViewProjection::kXpX:
+ nbinsx = 10 * ceil(size[0]); // mm binning
+ if (yrange < 0) {
+ LOG(debug) << "ProjectionP[" << name << "] using default range.";
+ yrange = 5; // +- 5 sigma range
+ }
+ nbinsy = 200; //* ceil(yrange);
+ get<2>(projection.second) =
+ new TH2D(Form("hpx_%s%s", dname, name.data()),
+ Form("X pulls %s [%s]; X_{%s-%s} (cm); pull(X)", name.data(), dname, dname, name.data()), nbinsx,
+ xlo, xhi, nbinsy, -yrange, yrange);
+ break;
+
+ case eViewProjection::kYpY:
+ nbinsx = 10 * ceil(size[1]); // mm binning
+ if (yrange < 0) {
+ LOG(debug) << "ProjectionP[" << name << "] using default range.";
+ yrange = 5; // +- 5 sigma range
+ }
+ nbinsy = 200; //* ceil(yrange);
+ get<2>(projection.second) =
+ new TH2D(Form("hpy_%s%s", dname, name.data()),
+ Form("Y pulls %s [%s]; Y_{%s-%s} (cm); pull(Y)", name.data(), dname, dname, name.data()), nbinsx,
+ ylo, yhi, nbinsy, -yrange, yrange);
+ break;
+
+ case eViewProjection::kWdT:
+ nbinsx = 10 * ceil(size[0]); // mm binning
+ unit = makeTrange(scale, yrange);
+ nbinsy = 2 * ceil(yrange);
+ get<2>(projection.second) = new TH2D(Form("hxt_%s%s", dname, name.data()),
+ Form("Hit_Trk %s [%s]; X_{%s-%s} (cm); #Delta time_{TRK} (%s)",
+ name.data(), dname, dname, name.data(), unit.data()),
+ nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
+ break;
+
+ case eViewProjection::kChdT:
+ nbinsx = dscale * ceil(size[0]); // mm binning
+ unit = makeTrange(scale, yrange);
+ nbinsy = 10 * ceil(yrange);
+ get<2>(projection.second) = new TH2D(Form("hct_%s%s", dname, name.data()),
+ Form("Hit_Event %s [%s]; X_{%s-%s} (cm); #Delta time_{EV} (%s)",
+ name.data(), dname, dname, name.data(), unit.data()),
+ nbinsx, xlo, xhi, nbinsy, -yrange, yrange);
+ break;
+
+ case eViewProjection::kXYa:
+ xy_id = 'a';
+ // fall through
+ case eViewProjection::kXYs:
+ if (!xy_id) xy_id = 's';
+ // fall through
+ case eViewProjection::kXYh:
+ if (!xy_id) xy_id = 'h';
+ nbinsx = dscale * ceil(size[0]); // mm binning
+ nbinsy = dscale * ceil(size[1]);
+ get<2>(projection.second) =
+ new TH2D(Form("hxy%c_%s%s", xy_id, dname, name.data()),
+ Form("Hit_{%s} %s [%s]; X_{%s-%s} (cm); Y_{%s-%s} (cm)", (xy_id == 'h' ? "reco" : "attach"),
+ name.data(), dname, dname, name.data(), dname, name.data()),
+ nbinsx, xlo, xhi, nbinsy, ylo, yhi);
+ break;
+ case eViewProjection::kXYp:
+ nbinsx = dscale * ceil(size[0] * 2.); // mm binning
+ nbinsy = dscale * ceil(size[1] * 2.);
+ xlo = pos[0] - size[0];
+ xhi = pos[0] + size[0];
+ ylo = pos[1] - size[1];
+ yhi = pos[1] + size[1];
+ get<2>(projection.second) = new TH2D(Form("hxyp_%s%s", dname, name.data()),
+ Form("Trk_{proj.} %s [%s]; X_{%s-%s} (cm); Y_{%s-%s} (cm)", name.data(),
+ dname, dname, name.data(), dname, name.data()),
+ nbinsx, xlo, xhi, nbinsy, ylo, yhi);
+ break;
+ }
+ }
+ return true;
+}
+
+string CbmRecoQaTask::Detector::View::ToString() const
+{
+ stringstream s;
+ s << "V[" << name << "] path[" << path << "] @ [" << pos[0] << ", " << pos[1] << ", " << pos[2] << "] size["
+ << size[0] << ", " << size[1] << ", " << size[2] << "] projections[" << fProjection.size() << "] sel[";
+ for (auto ii : fSelector)
+ s << ii << " ";
+ s << "]\n";
+ return s.str();
+}
+
+uint CbmRecoQaTask::Detector::View::Register(TDirectoryFile* fOut)
+{
+ uint n(0);
+ TDirectoryFile* lDir = (TDirectoryFile*) fOut->mkdir(name.data(), Form("QA for vol[%s]", path.data()));
+ TDirectoryFile* sDir =
+ (!CbmRecoQaTask::fuRecoConfig[kRecoEvents] ? (TDirectoryFile*) lDir->mkdir("TS", "TimeSlice QA") : nullptr);
+ TDirectoryFile* eDir =
+ (CbmRecoQaTask::fuRecoConfig[kRecoEvents] ? (TDirectoryFile*) lDir->mkdir("event", "Local Reco QA") : nullptr);
+ TDirectoryFile* tDir =
+ (CbmRecoQaTask::fuRecoConfig[kRecoTracks] ? (TDirectoryFile*) lDir->mkdir("track", "CA QA") : nullptr);
+ LOG(debug) << "CbmRecoQaTask::Detector::View[" << name << "]::Register() : "
+ << "\n\t level \"TS \" [" << (sDir ? "ON" : "OFF") << "]"
+ << "\n\t level \"EVT\" [" << (eDir ? "ON" : "OFF") << "]"
+ << "\n\t level \"TRK\" [" << (tDir ? "ON" : "OFF") << "]";
+
+ for (auto& projection : fProjection) {
+ TH2* hh = get<2>(projection.second);
+ if (!hh) continue;
+ switch (projection.first) {
+ case eViewProjection::kXYa:
+ case eViewProjection::kXYp:
+ case eViewProjection::kXdX:
+ case eViewProjection::kYdY:
+ case eViewProjection::kWdT:
+ case eViewProjection::kXpX:
+ case eViewProjection::kYpY:
+ if (tDir) tDir->Add(hh);
+ break;
+ case eViewProjection::kChdT:
+ case eViewProjection::kXYh:
+ if (eDir) eDir->Add(hh);
+ break;
+ case eViewProjection::kXYs:
+ if (sDir) sDir->Add(hh);
+ break;
}
- case kMcbm24: break;
- case kDefault: break;
+ n++;
}
+ LOG(debug) << "CbmRecoQaTask::Detector::View[" << name << "]::Register() : " << n << " projections for "
+ << fOut->GetName();
+ return n;
+}
+
+string CbmRecoQaTask::Detector::View::ToString(eViewProjection prj)
+{
+ string s;
+ switch (prj) {
+ case eViewProjection::kXYa: s = "x-y (attach)"; break;
+ case eViewProjection::kXYp: s = "x-y (track)"; break;
+ case eViewProjection::kXdX: s = "dx-x"; break;
+ case eViewProjection::kYdY: s = "dy-y"; break;
+ case eViewProjection::kXpX: s = "pull(x)-x"; break;
+ case eViewProjection::kYpY: s = "pull(y)-y"; break;
+ case eViewProjection::kWdT: s = "dt(trk)-w"; break;
+ case eViewProjection::kChdT: s = "dt(ev)-w"; break;
+ case eViewProjection::kXYh: s = "x-y (hit)"; break;
+ case eViewProjection::kXYs: s = "x-y (TS)"; break;
+ default: LOG(error) << "Detector::View::ToString() : Unknown projection " << int(prj); break;
+ }
+ return s;
+}
+
+string CbmRecoQaTask::Detector::View::makeYrange(const int scale, float& yrange)
+{
+ bool kDefaultRange = false;
+ string unit = "cm";
+ if (yrange < 0) {
+ LOG(debug) << "ProjectionY[" << name << "] using default range.";
+ kDefaultRange = true;
+ yrange = 3; // +- 6 cm default range
+ }
+ if (scale == 10) { // mm resolution
+ unit = "mm";
+ if (kDefaultRange) yrange = 10; // +- 20 mm default range
+ }
+ else if (scale == 10000) { // micron resolution
+ unit = "#mu m";
+ if (kDefaultRange) yrange = 150; // +- 300 um default range
+ }
+ return unit;
+}
+
+string CbmRecoQaTask::Detector::View::makeTrange(const int scale, float& yrange)
+{
+ bool kDefaultRange = false;
+ string unit = "ns";
+ if (yrange < 0) {
+ LOG(debug) << "ProjectionT[" << name << "] using default range.";
+ kDefaultRange = true;
+ yrange = 40; // +- 80 ns default range
+ }
+ if (scale == 1000) { // ps resolution
+ unit = "ps";
+ if (kDefaultRange) yrange = 300; // +- 600 ps default range
+ }
+ return unit;
+}
+
+CbmRecoQaTask::EventFilter* CbmRecoQaTask::AddEventFilter(CbmRecoQaTask::EventFilter::eEventCut type)
+{
+ for (auto cut : fFilterEv) {
+ if (cut.fType == type) {
+ LOG(warning) << GetName() << "::AddEventFilter event filter : " << cut.ToString() << " already on the list.";
+ return nullptr;
+ }
+ }
+ fFilterEv.emplace_back(type);
+ return &fFilterEv.back();
+}
+
+CbmRecoQaTask::TrackFilter* CbmRecoQaTask::AddTrackFilter(CbmRecoQaTask::TrackFilter::eTrackCut type)
+{
+ for (auto cut : fFilterTrk) {
+ if (cut.fType == type) {
+ LOG(warning) << GetName() << "::AddTrackFilter track filter : " << cut.ToString() << " already on the list.";
+ return nullptr;
+ }
+ }
+ fFilterTrk.emplace_back(type);
+ return &fFilterTrk.back();
+}
+
+//+++++++++++++++++++++ EventFilter ++++++++++++++++++++++++++++++++
+//____________________________________________________________________
+bool CbmRecoQaTask::EventFilter::Accept(const CbmEvent* ptr, const CbmRecoQaTask* /* lnk*/)
+{
+ stringstream ss;
+ bool ret(true);
+ size_t val(0);
+ switch (fType) {
+ case eEventCut::kMultTrk:
+ val = ptr->GetNofData(ECbmDataType::kGlobalTrack);
+ if (fMinTrack > 0 && val < size_t(fMinTrack)) {
+ ss << "NofTrack[" << val << "] < min[" << fMinTrack << "].";
+ ret = false;
+ break;
+ }
+ if (ret && fMaxTrack > 0 && val > size_t(fMaxTrack)) {
+ ss << "NofTrack[" << val << "] > max[" << fMaxTrack << "].";
+ ret = false;
+ }
+ break;
+ case eEventCut::kMultHit:
+ val = ptr->GetNofData(ECbmDataType::kStsHit);
+ if (fMultHit[0] > 0 && val > size_t(fMultHit[0])) {
+ ss << "Sts hits [" << val << "] > max[" << fMultHit[0] << "].";
+ ret = false;
+ break;
+ }
+ val = ptr->GetNofData(ECbmDataType::kTrdHit);
+ if (fMultHit[1] > 0 && val > size_t(fMultHit[1])) {
+ ss << "Trd hits [" << val << "] > max[" << fMultHit[1] << "].";
+ ret = false;
+ break;
+ }
+ val = ptr->GetNofData(ECbmDataType::kTofHit);
+ if (fMultHit[2] > 0 && val > size_t(fMultHit[2])) {
+ ss << "Tof hits [" << val << "] > max[" << fMultHit[2] << "].";
+ ret = false;
+ }
+ break;
+ default: break;
+ }
+ if (!ret) LOG(debug2) << "Event reject for : " << ss.str();
return ret;
}
+bool CbmRecoQaTask::EventFilter::SetFilter(std::vector<float> cuts)
+{
+ switch (fType) {
+ case eEventCut::kMultTrk:
+ if (cuts.size() < 2 || cuts[1] < cuts[0]) {
+ LOG(warning) << "Improper definition for event filter :\n\t" << ToString() << endl;
+ HelpMess();
+ return false;
+ }
+ fMinTrack = int(cuts[0]);
+ fMaxTrack = int(cuts[1]);
+ break;
+ case eEventCut::kMultHit:
+ if (cuts.size() < 3) {
+ LOG(warning) << "Improper definition for event filter :\n\t" << ToString() << endl;
+ HelpMess();
+ return false;
+ }
+ for (int i(0); i < int(CbmRecoQaTask::EventFilter::eEventDef::kNofDetHit); i++)
+ fMultHit[i] = int(cuts[i]);
+ break;
+ case eEventCut::kTrigger: break;
+ case eEventCut::kVertex: break;
+ default: break;
+ }
+ return true;
+}
+std::string CbmRecoQaTask::EventFilter::ToString() const
+{
+ stringstream ss;
+ switch (fType) {
+ case eEventCut::kMultTrk: ss << "kMultTrk : \"cut on track multiplicity"; break;
+ case eEventCut::kMultHit: ss << "kMultHit : \"cut on hit multiplicity"; break;
+ case eEventCut::kTrigger: ss << "kTrigger : \"cut on trigger conditions"; break;
+ case eEventCut::kVertex: ss << "kVertex : \"cut on vertex definition"; break;
+ default: break;
+ }
+ return ss.str();
+}
+void CbmRecoQaTask::EventFilter::HelpMess() const
+{
+ LOG(info) << "CbmRecoQaTask::EventFilter : Usage";
+ switch (fType) {
+ case eEventCut::kMultTrk:
+ LOG(info) << ToString();
+ LOG(info) << "\tDepends one two values : min and max of the no of global tracks / event.";
+ LOG(info) << "\tValue should follow the relation max >= min.";
+ break;
+ case eEventCut::kMultHit:
+ LOG(info) << ToString();
+ LOG(info) << "\tDepends one three values : total hits in the following systems (in this order) : STS TRD TOF.";
+ LOG(info) << "\tNegative values are interpreted as no-cut.";
+ break;
+ default:
+ LOG(info) << ToString();
+ LOG(info) << "\tNo user info.";
+ break;
+ }
+}
+//+++++++++++++++++++++ TrackFilter ++++++++++++++++++++++++++++++++
+//____________________________________________________________________
+bool CbmRecoQaTask::TrackFilter::Accept(const CbmGlobalTrack* ptr, const CbmRecoQaTask* lnk)
+{
+ bool ret(true);
+ stringstream ss;
+ int idx(-1);
+ const CbmTrack* trk(nullptr);
+ switch (fType) {
+ case eTrackCut::kSts:
+ if (fNSts <= 0 && !fStsHits.size()) break;
+ if ((idx = ptr->GetStsTrackIndex()) < 0) {
+ ss << "Sts trk index missing";
+ ret = false;
+ break;
+ }
+ if (!(trk = lnk->GetTrack(ECbmModuleId::kSts, idx))) {
+ ss << "Sts trk nullptr";
+ ret = false;
+ break;
+ }
+ // cut on the total no of STS hits/trk
+ if (trk->GetNofHits() < fNSts) {
+ ss << "Sts hits/trk [" << trk->GetNofHits() << "] < min[" << fNSts << "].";
+ ret = false;
+ break;
+ }
+ // cut on the distribution of STS hits/trk/station
+ for (auto stsStation : fStsHits) {
+ if (!stsStation) continue;
+ // check hit in this particular station
+ }
+ break;
+ case eTrackCut::kTrd:
+ if (fNTrd <= 0) break;
+ if ((idx = ptr->GetTrdTrackIndex()) < 0) {
+ ss << "Trd trk index missing";
+ ret = false;
+ break;
+ }
+ if (!(trk = lnk->GetTrack(ECbmModuleId::kTrd, idx))) {
+ ss << "Trd trk nullptr";
+ ret = false;
+ break;
+ }
+ // cut on the total no of TRD hits/trk
+ if (trk->GetNofHits() < fNTrd) {
+ ss << "Trd` hits/trk [" << trk->GetNofHits() << "] < min[" << fNTrd << "].";
+ ret = false;
+ break;
+ }
+ break;
+ case eTrackCut::kTof:
+ if (fNTof <= 0) break;
+ if ((idx = ptr->GetTofTrackIndex()) < 0) {
+ ss << "Tof trk index missing";
+ ret = false;
+ break;
+ }
+ if (!(trk = lnk->GetTrack(ECbmModuleId::kTof, idx))) {
+ ss << "Tof trk nullptr";
+ ret = false;
+ break;
+ }
+ // cut on the total no of ToF hits/trk
+ if (trk->GetNofHits() < fNTof) {
+ ss << "Tof` hits/trk [" << trk->GetNofHits() << "] < min[" << fNTof << "].";
+ ret = false;
+ break;
+ }
+ break;
+ default: break;
+ }
-//_____________________________________________________________________
-void CbmRecoQaTask::Finish()
+ if (!ret) LOG(debug2) << "Track reject for : " << ss.str();
+ return ret;
+}
+bool CbmRecoQaTask::TrackFilter::SetFilter(std::vector<float> cuts)
{
- FairSink* sink = FairRootManager::Instance()->GetSink();
- sink->WriteObject(&fOutFolder, nullptr);
+ switch (fType) {
+ case eTrackCut::kSts:
+ if (cuts.size() < 1) {
+ LOG(warning) << "Improper definition for track filter :\n\t" << ToString() << endl;
+ //HelpMess();
+ return false;
+ }
+ fNSts = int(cuts[0]);
+ break;
+ case eTrackCut::kTrd:
+ if (cuts.size() < 1) {
+ LOG(warning) << "Improper definition for track filter :\n\t" << ToString() << endl;
+ //HelpMess();
+ return false;
+ }
+ fNTrd = int(cuts[0]);
+ break;
+ case eTrackCut::kTof:
+ if (cuts.size() < 1) {
+ LOG(warning) << "Improper definition for track filter :\n\t" << ToString() << endl;
+ //HelpMess();
+ return false;
+ }
+ fNTof = int(cuts[0]);
+ break;
+ default: break;
+ }
+ return true;
+}
+std::string CbmRecoQaTask::TrackFilter::ToString() const
+{
+ stringstream ss;
+ switch (fType) {
+ case eTrackCut::kSts: ss << "kSts : \"cut on no of STS hits / track\""; break;
+ case eTrackCut::kMuch: ss << "kMuch : \"cut on no of Much hits / track\""; break;
+ case eTrackCut::kRich: ss << "kRich : \"cut on no of Rich hits / track\""; break;
+ case eTrackCut::kTrd: ss << "kTrd : \"cut on no of TRD hits / track\""; break;
+ case eTrackCut::kTof: ss << "kTof : \"cut on no of Tof hits / track\""; break;
+ default: break;
+ }
+ return ss.str();
}
+/* clang-format off */
+ClassImp(CbmRecoQaTask)
+ClassImp(CbmRecoQaTask::Detector)
+ClassImp(CbmRecoQaTask::Detector::View)
+ClassImp(CbmRecoQaTask::TrackFilter)
+ClassImp(CbmRecoQaTask::EventFilter)
+ /* clang-format on */
diff --git a/reco/qa/CbmRecoQaTask.h b/reco/qa/CbmRecoQaTask.h
index 44dca4e5d08b3c6959cb9e3de45e1d81aa2747cc..2f219191cea9137862edb379716d0a8924bdfdf6 100644
--- a/reco/qa/CbmRecoQaTask.h
+++ b/reco/qa/CbmRecoQaTask.h
@@ -6,66 +6,244 @@
#define CBMRECOQATASK_H 1
#include "CbmKfTrackFitter.h"
-// #include "CbmQaHist.h"
+#include "CbmTrack.h"
-#include "FairTask.h"
-#include "TFolder.h"
+#include <FairTask.h>
+#include <TDirectoryFile.h>
+
+#include <bitset>
#include <map>
#include <vector>
+#define kRecoQaNConfigs 8
+
+class CbmEvent;
+class CbmHit;
class CbmTimeSlice;
class TClonesArray;
+class TH2;
class TH2D;
class CbmRecoQaTask : public FairTask {
public:
- enum Setup
+ enum eRecoConfig
+ {
+ kRecoEvents = 0, /// has events reconstructed (CbmEvent branch)
+ kRecoTracks, /// has tracks reconstructed (GlobalTrack branch)
+ kStsHits, /// has STS hits (StsHit branch)
+ kTrdHits, /// has TRD` hits (TrdHit branch)
+ kTofHits, /// has ToF hits (TofHit branch)
+ kRichHits, /// has Rich hits (RichHit branch)
+ kMuchHits, /// has Much hits (MuchHit branch)
+ };
+ enum eSetup
{
kMcbm22 = 0,
kMcbm24,
kDefault
};
+ enum class eViewProjection : int
+ {
+ kXYa = 0, /// X-Y hit coorelation in track filtered data
+ kXYp, /// X-Y track projections
+ kXdX, /// X to TRK residuals as function of local X in view
+ kYdY, /// Y to TRK residuals as function of local Y in view
+ kWdT, /// Time to TRK residuals as function of local high resolution coordinate in view
+ kXpX, /// X to TRK pulls as function of local X in view
+ kYpY, /// Y to TRK pulls as function of local Y in view
+ kChdT, /// Time to EV residuals as function of coordinate in view
+ kXYh, /// X-Y hit coorelation in local view
+ kXYs /// X-Y hit coorelation in time slice
+ };
+ struct Detector {
+ struct Data {
+ Data() = default;
+ Data(ECbmDataType i, const char* n) : id(i), name(n) { ; }
+ ECbmDataType id = ECbmDataType::kUnknown;
+ std::string name = "nn";
+ };
+ struct View {
+ friend struct Detector;
+ friend class CbmRecoQaTask;
+
+ View() = default;
+ View(const char* n, const char* p, std::vector<int> set) : name(n), path(p), fSelector(set) { ; }
+ virtual ~View() = default;
+ bool SetProjection(eViewProjection prj, float range, const char* unit);
+ template<class Hit>
+ bool HasAddress(const CbmHit* h, double& x, double& y) const;
+ template<class Hit>
+ bool Load(const CbmHit* h, const CbmEvent* ev = nullptr);
+ bool Load(const CbmKfTrackFitter::FitNode* n);
+ std::string ToString() const;
+ static std::string ToString(eViewProjection prj);
+
+ std::string name = ""; /// name describing the module
+ std::string path = ""; /// path to the geo volume describing the module
+ double size[3] = {0.}; /// detection element geometrical dx dy dz dimmensions
+ double pos[3] = {0.}; /// detection element center x0 y0 z0
+ std::vector<int> fSelector = {}; /// defining subset of the address set for this view
+ std::map<eViewProjection, std::tuple<int, float, TH2*>> fProjection =
+ {}; /// map of projections indexed by their type. Each projection contains also the relative scale [int] wrt to default units (ns, cm, keV) and the range [float].
+ protected:
+ bool AddProjection(eViewProjection prj, float range = -1, const char* unit = "cm");
+ /** \brief convert geo address into pointer to geometry.
+ * By the time of the call the geometry has to be available in memory.
+ * Failing to identify the named physical node will rezult in error */
+ bool Init(const char* dname);
+ /** \brief build directory structure for all projections of current view.*/
+ uint Register(TDirectoryFile* f);
+ /** helper functions to estimate the representation (y) axis
+ * \param[in] scale read-only unit defining parameter
+ * \param[in] range read-write full range on the ordinate
+ * \return units for the ordinate and the range value
+ */
+ std::string makeTrange(const int scale, float& range);
+ std::string makeYrange(const int scale, float& range);
+
+ ClassDef(CbmRecoQaTask::Detector::View, 1); // A QA significant correlation
+ };
+ Detector(ECbmModuleId did = ECbmModuleId::kNotExist);
+ virtual ~Detector() = default;
+ View* AddView(const char* n, const char* p, std::vector<int> set);
+ View* GetView(const char* n);
+ View* FindView(double x, double y, double z);
+ /** \brief Check geometry and trigger Init() for all registered views. Build main directory outut for the current detector.
+ * Failing to identify the geometry will rezult in fatal error.
+ * \return true If ALL the subsequent calls to Init result in a true */
+ bool Init(TDirectoryFile* f);
+ void Print() const;
+
+ ECbmModuleId id;
+ Data hit;
+ Data trk;
+ std::vector<View> fViews = {};
+
+ ClassDef(CbmRecoQaTask::Detector, 1); // QA representation of a detector unit
+ };
+
+ struct TrackFilter {
+ enum class eTrackCut : int
+ {
+ kSts = 0, /// cut on no of STS hits / track
+ kMuch, /// cut on no of Much hits / track
+ kRich, /// cut on no of Rich hits / track
+ kTrd, /// cut on no of TRD hits / track
+ kTof, /// cut on no of Tof hits / track
+ kNone /// no cut
+ };
+ TrackFilter(eTrackCut typ = eTrackCut::kNone) : fType(typ) { ; }
+ virtual ~TrackFilter() = default;
+ bool Accept(const CbmGlobalTrack* ptr, const CbmRecoQaTask* lnk);
+ bool SetFilter(std::vector<float> cuts);
+ std::string ToString() const;
+
+ eTrackCut fType = eTrackCut::kNone;
+
+ private:
+ // STS definition of track
+ int fNSts = -1;
+ std::vector<bool> fStsHits = {};
+ // TRD definition of track
+ int fNTrd = -1;
+ // ToF definition of track
+ int fNTof = -1;
+ ClassDef(CbmRecoQaTask::TrackFilter, 1); // Track cut implementation
+ };
+
+ struct EventFilter {
+ enum class eEventCut : int
+ {
+ kMultTrk = 0, /// cut on track multiplicity
+ kMultHit, /// cut on hit multiplicity
+ kTrigger, /// cut on trigger conditions
+ kVertex, /// cut on vertex definition
+ kNone /// no cut
+ };
+ enum class eEventDef : int
+ {
+ kNofDetHit = 3 /// no of detectors considered for the hit cut [STS TRD ToF]
+ };
+ EventFilter(eEventCut typ = eEventCut::kNone) : fType(typ) { ; }
+ virtual ~EventFilter() = default;
+ bool Accept(const CbmEvent* ptr, const CbmRecoQaTask* lnk);
+ bool SetFilter(std::vector<float> cuts);
+ std::string ToString() const;
+
+ eEventCut fType = eEventCut::kNone;
+
+ private:
+ // track multiplicity cut definition
+ int fMinTrack = 0;
+ int fMaxTrack = 0;
+ // hit multiplicity cut definition
+ int fMultHit[(int) eEventDef::kNofDetHit] = {0}; /// max no of hits/ev for the systems [STS TRD ToF]
+
+ /** \brief Helper function : display usage message for each filter case*/
+ void HelpMess() const;
+
+ ClassDef(CbmRecoQaTask::EventFilter, 1); // Event cut implementation
+ };
public:
CbmRecoQaTask();
virtual ~CbmRecoQaTask() { ; }
-
+ /** Copy the qa test defined for detector det from the steering macro to the current class */
+ virtual Detector* AddDetector(ECbmModuleId did);
+ virtual EventFilter* AddEventFilter(EventFilter::eEventCut cut);
+ virtual TrackFilter* AddTrackFilter(TrackFilter::eTrackCut cut);
+ virtual Detector* GetDetector(ECbmModuleId did);
+ /** \brief Retrieve detector specific track by index.
+ * \param[in] did system Identifier
+ * \param[in] id track id in the TClonesArray
+ * \return pointer to track or nullptr if track/system not available
+ */
+ virtual const CbmTrack* GetTrack(ECbmModuleId did, int id) const;
+ /** \brief Perform initialization of data sources and projections */
virtual InitStatus Init();
/** \brief Executed task **/
virtual void Exec(Option_t* option);
- /** \brief Filter tracks for further use (e.g. track projections)
- * \param[in] ptr information on which to filter (e.g. nhits from detectors)
- * \return true if track accepted
- */
- virtual bool FilterTrack(int* ptr);
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; }
+ void SetSetupClass(CbmRecoQaTask::eSetup setup) { fSetupClass = setup; }
+
+ protected:
+ static std::bitset<kRecoQaNConfigs> fuRecoConfig;
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;
- };
+ /** \brief Filter events for QA use (e.g. event multiplicity)
+ * \param[in] ptr cbm event
+ * \return true if event accepted
+ */
+ virtual bool FilterEvent(const CbmEvent* ptr);
+ /** \brief Filter tracks for further use (e.g. track projections)
+ * \param[in] ptr global track
+ * \return true if track accepted
+ */
+ virtual bool FilterTrack(const CbmGlobalTrack* ptr);
+ void InitMcbm22();
+ void InitMcbm24();
+
CbmKfTrackFitter fFitter;
- TClonesArray* fTracks = nullptr; //! reconstructed global tracks / event
+ TClonesArray* fGTracks = nullptr; //! reconstructed global tracks / event
+ std::map<ECbmModuleId, TClonesArray*> fTracks = {}; //! 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::map<ECbmModuleId, TClonesArray*> fHits = {}; //! reconstructed hits
+ std::map<ECbmModuleId, Detector> fDetQa = {};
+ std::vector<EventFilter> fFilterEv = {};
+ std::vector<TrackFilter> fFilterTrk = {};
+
+ TDirectoryFile fOutFolder = {"RecoQA", "CA track driven reco QA"}; //!
+ eSetup fSetupClass = eSetup::kMcbm24;
+ std::map<double, Detector::View> fProjs = {}; //!
ClassDef(CbmRecoQaTask, 1); // Reconstruction QA analyzed from CA perspective
};
diff --git a/reco/qa/RecoQaLinkDef.h b/reco/qa/RecoQaLinkDef.h
index ecbf6054cf0f27e7c3f030e6bf13845e335dd5b6..f4e30615242b36ac0c72d130db7179e8f68d37ae 100644
--- a/reco/qa/RecoQaLinkDef.h
+++ b/reco/qa/RecoQaLinkDef.h
@@ -13,5 +13,9 @@
#pragma link C++ class CbmRecoQa + ;
#pragma link C++ class CbmTrackingTrdQa + ;
#pragma link C++ class CbmRecoQaTask + ;
+#pragma link C++ class CbmRecoQaTask::Detector + ;
+#pragma link C++ class CbmRecoQaTask::Detector::View + ;
+#pragma link C++ class CbmRecoQaTask::EventFilter + ;
+#pragma link C++ class CbmRecoQaTask::TrackFilter + ;
#endif