diff --git a/macro/L1/search_window.C b/macro/L1/search_window.C
new file mode 100644
index 0000000000000000000000000000000000000000..767f01340ab9865e97484e57daaa164720cca138
--- /dev/null
+++ b/macro/L1/search_window.C
@@ -0,0 +1,83 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file search_window.C
+/// \brief A macro to run a hit-search window estimator for the CA tracking
+/// \since 8.11.2022
+/// \author S.Zharko <s.zharko@gsi.de>
+
+/// Splits a string to a vector of substrings using a char token
+/// \param origString Original string
+/// \param token Char token
+std::vector<TString> SplitString(const char* aString, char token = ':');
+
+// Instruction
+//
+// 1. Creation of MC-triplets tree input
+// In a reconstruction macro, please, run CbmL1::SetOutputMcTripletsTreeFilename("<prefix>"), where <prefix> is a pre-
+// fix for an output file. NOTE: this feature will not work with the "debugWithMC" flag switched off
+//
+// 2. Initialization of search windows file to CA tracking
+// The file contains an array of serialized L1SearchWindow objects. The file can be generated
+// for a given set of station IDs, MC-triplet trees and a sequence of track finder iteration using
+// this macro. The L1SearchWindow objects are optional. They are initialized as far as a corresponding input file
+// is provided. To do so, please, add this function into your reconstruction macro
+//l1->SetInputSearchWindowFilename(TString(gSystem->Getenv("VMCWORKDIR")) + "/macro/L1/SearchWindows.dat");
+
+
+/// Estimates search windows
+void search_window(
+ const char* yamlConfigFile = "", ///< Name of yaml configuration file containing track finder iterations
+ const char* inputTreeFile = "", ///< Name of ROOT input file with MC triplets
+ const char* swOutputFile = "", ///< Name of L1SearchWindow objects output
+ const char* logOutputName = "" ///< Name of pdf log
+)
+{
+ // ----- Create a window finder object
+ ca::tools::WindowFinder wf;
+
+ // ----- Initialize window finder parameters
+ wf.SetOutputName(swOutputFile);
+ wf.ReadTrackingIterationsFromYAML(yamlConfigFile);
+
+ // Set global indexes of active tracking stations by which the search windows will be estimated
+ std::vector<int> vStationIDs = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
+ wf.SetStationIndexes(vStationIDs);
+
+ // Set target (NOTE: for now x and y of the target should be 0.)
+ wf.SetTarget(0., 0., -44.);
+
+ // Additional cut
+ //TCut extraCut = "abs(pdg) == 211"; // Select pions only
+ //wf.SetExtraCut(extraCut);
+
+ // ----- Define input trees with MC triplets
+ auto vsInputTreeFilenames = SplitString(inputTreeFile);
+ for (const auto& inputTreeFilename : vsInputTreeFilenames) {
+ wf.AddInputFile(inputTreeFilename);
+ }
+
+ // ----- Run the windows estimator
+ wf.Process();
+
+ // ----- Save a log to pdf files
+ wf.DumpCanvasesToPdf(logOutputName);
+}
+
+
+// *********************
+// ** Implementations **
+// *********************
+
+std::vector<TString> SplitString(const char* aString, char token)
+{
+ TString sOrigin = TString(aString);
+ TObjArray* strings = sOrigin.Tokenize(TString(token));
+
+ std::vector<TString> res;
+ for (int i = 0; i < strings->GetEntries(); ++i) {
+ res.push_back(((TObjString*) strings->At(i))->String());
+ }
+ return res;
+}
diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C
index 2ed01f00b1f084253e5ade68e3449b5219a9d821..997d408d42e4ebfab33fcae82e49dbc132bffc4a 100644
--- a/macro/run/run_reco.C
+++ b/macro/run/run_reco.C
@@ -381,10 +381,11 @@ void run_reco(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice =
// L1 tracking
auto l1 = (debugWithMC) ? new CbmL1("L1", 2, 3, 0) : new CbmL1("L1", 0);
+
// L1 configuration file (optional)
+ // At the moment, the YAML configuration file defines different parameters for a sequence of track finder
+ // iterations. The same file should be used in ca::tools::WindowFinder class for hit search window estimation
//l1->SetInputConfigName(TString(gSystem->Getenv("VMCWORKDIR")) + "/reco/L1/L1Algo/L1ConfigExample.yaml");
- // Define filename to save MC triplets
- //if (debugWithMC) { l1->SetOutputMcTripletsTreeFilename("mcTripletsTree"); }
// --- Material budget file names
TString mvdGeoTag;
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 3f45f45072de8e994cad5493970df75023897523..a9e5be41c11c0552a33cb44231c58da9c2ab4dd8 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -64,6 +64,7 @@ set(SRCS
L1Algo/L1IODataManager.cxx
L1Algo/L1CloneMerger.cxx
L1Algo/L1ConfigRW.cxx
+ L1Algo/L1SearchWindow.cxx
L1Algo/utils/L1AlgoDraw.cxx
L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
L1Algo/utils/L1AlgoPulls.cxx
@@ -72,6 +73,7 @@ set(SRCS
catools/CaToolsPerformance.cxx
catools/CaToolsDebugger.cxx
catools/CaToolsWindowFinder.cxx
+ catools/CaToolsWFExpression.cxx
ParticleFinder/CbmL1PFFitter.cxx
ParticleFinder/CbmL1PFMCParticle.cxx
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 2cc894d87bac89bc709f5e784835c049721da45b..a12cd94ff93f0850cd7a244b583977a9ebc4df6e 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -632,14 +632,14 @@ InitStatus CbmL1::Init()
auto trackingIterFastPrim = L1CAIteration(trIterDefault, "FastPrimIter");
trackingIterFastPrim.SetMinNhitsStation0(3);
trackingIterFastPrim.SetPickGather(3.0f);
- trackingIterFastPrim.SetMaxInvMom(1.0 / 0.5);
+ trackingIterFastPrim.SetMaxQp(1.0 / 0.5);
trackingIterFastPrim.SetMaxDZ(0);
trackingIterFastPrim.SetTargetPosSigmaXY(1, 1);
trackingIterFastPrim.SetPrimaryFlag(true);
auto trackingIterAllPrim = L1CAIteration(trIterDefault, "AllPrimIter");
trackingIterAllPrim.SetMinNhitsStation0(3);
- trackingIterAllPrim.SetMaxInvMom(1.0 / 0.05);
+ trackingIterAllPrim.SetMaxQp(1.0 / 0.05);
trackingIterAllPrim.SetTargetPosSigmaXY(1, 1);
trackingIterAllPrim.SetPrimaryFlag(true);
@@ -647,7 +647,7 @@ InitStatus CbmL1::Init()
trackingIterFastPrim2.SetTripletChi2Cut(21.1075f);
trackingIterFastPrim2.SetTripletFinalChi2Cut(21.1075f);
trackingIterFastPrim2.SetPickGather(3.0f);
- trackingIterFastPrim2.SetMaxInvMom(1.0 / 0.5);
+ trackingIterFastPrim2.SetMaxQp(1.0 / 0.5);
trackingIterFastPrim2.SetMaxDZ(0);
trackingIterFastPrim2.SetTargetPosSigmaXY(5, 5);
trackingIterFastPrim2.SetPrimaryFlag(true);
@@ -655,7 +655,7 @@ InitStatus CbmL1::Init()
auto trackingIterAllSec = L1CAIteration(trIterDefault, "AllSecIter");
trackingIterAllSec.SetTripletChi2Cut(18.7560f); // = 6.252 * 3; // prob = 0.1
trackingIterAllSec.SetTripletFinalChi2Cut(18.7560f);
- trackingIterAllSec.SetMaxInvMom(1.0 / 0.1);
+ trackingIterAllSec.SetMaxQp(1.0 / 0.1);
trackingIterAllSec.SetTargetPosSigmaXY(10, 10);
trackingIterAllSec.SetMaxSlopePV(1.5);
trackingIterAllSec.SetPrimaryFlag(false);
@@ -664,7 +664,7 @@ InitStatus CbmL1::Init()
trackingIterFastPrimJump.SetTripletChi2Cut(21.1075f); // prob = 0.01
trackingIterFastPrimJump.SetTripletFinalChi2Cut(21.1075f);
trackingIterFastPrimJump.SetPickGather(3.0f);
- trackingIterFastPrimJump.SetMaxInvMom(1.0 / 0.5);
+ trackingIterFastPrimJump.SetMaxQp(1.0 / 0.5);
trackingIterFastPrimJump.SetMaxDZ(0);
trackingIterFastPrimJump.SetTargetPosSigmaXY(5, 5);
trackingIterFastPrimJump.SetPrimaryFlag(true);
@@ -673,7 +673,7 @@ InitStatus CbmL1::Init()
auto trackingIterAllPrimJump = L1CAIteration(trIterDefault, "AllPrimJumpIter");
trackingIterAllPrimJump.SetTripletChi2Cut(18.7560f);
trackingIterAllPrimJump.SetTripletFinalChi2Cut(18.7560f);
- trackingIterAllPrimJump.SetMaxInvMom(1.0 / 0.1);
+ trackingIterAllPrimJump.SetMaxQp(1.0 / 0.1);
trackingIterAllPrimJump.SetTargetPosSigmaXY(5, 5);
trackingIterAllPrimJump.SetPrimaryFlag(true);
trackingIterAllPrimJump.SetJumpedFlag(true);
@@ -681,7 +681,7 @@ InitStatus CbmL1::Init()
auto trackingIterAllSecJump = L1CAIteration(trIterDefault, "AllSecJumpIter");
trackingIterAllSecJump.SetTripletChi2Cut(21.1075f);
trackingIterAllSecJump.SetTripletFinalChi2Cut(21.1075f);
- trackingIterAllSecJump.SetMaxInvMom(1.0 / 0.1);
+ trackingIterAllSecJump.SetMaxQp(1.0 / 0.1);
trackingIterAllSecJump.SetMaxSlopePV(1.5f);
trackingIterAllSecJump.SetTargetPosSigmaXY(10, 10);
trackingIterAllSecJump.SetPrimaryFlag(false);
@@ -689,7 +689,7 @@ InitStatus CbmL1::Init()
auto trackingIterAllPrimE = L1CAIteration(trIterDefault, "AllPrimEIter");
trackingIterAllPrimE.SetMinNhitsStation0(3);
- trackingIterAllPrimE.SetMaxInvMom(1.0 / 0.05);
+ trackingIterAllPrimE.SetMaxQp(1.0 / 0.05);
trackingIterAllPrimE.SetTargetPosSigmaXY(1, 1);
trackingIterAllPrimE.SetPrimaryFlag(true);
trackingIterAllPrimE.SetElectronFlag(true);
@@ -697,7 +697,7 @@ InitStatus CbmL1::Init()
auto trackingIterAllSecE = L1CAIteration(trIterDefault, "AllSecEIter");
trackingIterAllSecE.SetTripletChi2Cut(18.7560f);
trackingIterAllSecE.SetTripletFinalChi2Cut(18.7560f);
- trackingIterAllSecE.SetMaxInvMom(1.0 / 0.05);
+ trackingIterAllSecE.SetMaxQp(1.0 / 0.05);
trackingIterAllSecE.SetMaxSlopePV(1.5f);
trackingIterAllSecE.SetTargetPosSigmaXY(10, 10);
trackingIterAllSecE.SetPrimaryFlag(false);
@@ -711,17 +711,17 @@ InitStatus CbmL1::Init()
trackingIterAllSecE.SetMinNhits(4);
trackingIterAllSecE.SetMinNhitsStation0(4);
- trackingIterAllPrim.SetMaxInvMom(1. / 0.1);
- trackingIterAllPrimE.SetMaxInvMom(1. / 0.1);
- trackingIterAllSecE.SetMaxInvMom(1. / 0.1);
+ trackingIterAllPrim.SetMaxQp(1. / 0.1);
+ trackingIterAllPrimE.SetMaxQp(1. / 0.1);
+ trackingIterAllSecE.SetMaxQp(1. / 0.1);
trackingIterAllSecJump.SetMinNhits(4);
trackingIterAllSecJump.SetMinNhitsStation0(4);
- trackingIterFastPrim.SetMaxInvMom(1.0 / 0.3);
- trackingIterAllSec.SetMaxInvMom(1.0 / 0.3);
- trackingIterFastPrimJump.SetMaxInvMom(1.0 / 0.3);
- trackingIterAllPrimJump.SetMaxInvMom(1.0 / 0.3);
- trackingIterAllSecJump.SetMaxInvMom(1.0 / 0.3);
+ trackingIterFastPrim.SetMaxQp(1.0 / 0.3);
+ trackingIterAllSec.SetMaxQp(1.0 / 0.3);
+ trackingIterFastPrimJump.SetMaxQp(1.0 / 0.3);
+ trackingIterAllPrimJump.SetMaxQp(1.0 / 0.3);
+ trackingIterAllSecJump.SetMaxQp(1.0 / 0.3);
trackingIterAllSec.SetMinNhits(4);
trackingIterAllSec.SetMinNhitsStation0(4);
@@ -751,7 +751,7 @@ InitStatus CbmL1::Init()
it.SetDoubletChi2Cut(7.56327f); // = 1.3449 * 2.f / 3.f; // prob = 0.1
it.SetPickGather(3.0f);
it.SetTripletLinkChi2(16.0);
- it.SetMaxInvMom(1.0 / 0.1); //(1.0 / 0.5);
+ it.SetMaxQp(1.0 / 0.1); //(1.0 / 0.5);
it.SetMaxSlopePV(.1f);
it.SetMaxSlope(.5f);
it.SetMaxDZ(0.05);
@@ -771,7 +771,7 @@ InitStatus CbmL1::Init()
trd2dIter2.SetDoubletChi2Cut(4. * 7.56327f); // = 1.3449 * 2.f / 3.f; // prob = 0.1
trd2dIter2.SetPickGather(3.0f);
trd2dIter2.SetTripletLinkChi2(16.);
- trd2dIter2.SetMaxInvMom(1.0 / 0.05); //(1.0 / 0.5);
+ trd2dIter2.SetMaxQp(1.0 / 0.05); //(1.0 / 0.5);
trd2dIter2.SetMaxSlopePV(.5f);
trd2dIter2.SetMaxSlope(.5f);
trd2dIter2.SetMaxDZ(0.05);
@@ -788,7 +788,7 @@ InitStatus CbmL1::Init()
trd2dIter3.SetDoubletChi2Cut(15.); // = 1.3449 * 2.f / 3.f; // prob = 0.1
trd2dIter3.SetPickGather(3.0f);
trd2dIter3.SetTripletLinkChi2(200.);
- trd2dIter3.SetMaxInvMom(1.0 / 0.1); //(1.0 / 0.5);
+ trd2dIter3.SetMaxQp(1.0 / 0.1); //(1.0 / 0.5);
trd2dIter3.SetMaxSlopePV(.01);
trd2dIter3.SetMaxSlope(.4); //.5f);
trd2dIter3.SetMaxDZ(0.05);
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 99d15a41f43b6dda25b57d6addbf4a1c9f384d30..5b38e7d06c95d4c7a07220dce7a428977871f1b1 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -2254,8 +2254,8 @@ void CbmL1::InputPerformance()
void CbmL1::DumpMCTripletsToTree()
{
if (!fpMcTripletsTree) {
- auto* currentDir = gDirectory;
- auto* currentFile = gFile;
+ TDirectory* currentDir = gDirectory;
+ TFile* currentFile = gFile;
// Get prefix and directory
@@ -2324,7 +2324,6 @@ void CbmL1::DumpMCTripletsToTree()
float x; ///< x-component of point position [cm]
float y; ///< y-component of point position [cm]
float z; ///< z-component of point position [cm]
- bool operator<(const ReducedMcPoint& other) { return s < other.s; }
};
for (const auto& tr : fvMCTracks) {
@@ -2338,7 +2337,9 @@ void CbmL1::DumpMCTripletsToTree()
vPoints.emplace_back(ReducedMcPoint {point.iStation, float(point.x), float(point.y), float(point.z)});
}
- std::sort(vPoints.begin(), vPoints.end());
+ std::sort(vPoints.begin(), vPoints.end(),
+ [](const ReducedMcPoint& lhs, const ReducedMcPoint& rhs) { return lhs.s < rhs.s; });
+
for (unsigned int i = 0; i + 2 < vPoints.size(); ++i) {
// Condition to collect only triplets without gaps in stations
// TODO: SZh 20.10.2022 Add cases for jump iterations
diff --git a/reco/L1/L1Algo/L1CAIteration.cxx b/reco/L1/L1Algo/L1CAIteration.cxx
index 5fe919b75dc50e1accb3feb78bbfce55e7d95d4b..e70460b1206ee14106f70a978584dababe7716a8 100644
--- a/reco/L1/L1Algo/L1CAIteration.cxx
+++ b/reco/L1/L1Algo/L1CAIteration.cxx
@@ -36,7 +36,7 @@ bool L1CAIteration::Check() const
res = this->CheckValueLimits<float>("doublet_chi2_cut", fDoubletChi2Cut, 0.f, kMaxFloat) && res;
res = this->CheckValueLimits<float>("pick_gather", fPickGather, 0.f, kMaxFloat) && res;
res = this->CheckValueLimits<float>("triplet_link_chi2", fTripletLinkChi2, 0.f, kMaxFloat) && res;
- res = this->CheckValueLimits<float>("min_momentum", fMaxInvMom, 1.f / kMaxFloat, 1.f / 0.001f) && res;
+ res = this->CheckValueLimits<float>("max_qp", fMaxQp, 0.001f, kMaxFloat) && res;
res = this->CheckValueLimits<float>("max_slope_pv", fMaxSlopePV, 0.f, kMaxFloat) && res;
res = this->CheckValueLimits<float>("max_slope", fMaxSlope, 0.f, kMaxFloat) && res;
res = this->CheckValueLimits<float>("max_dz", fMaxDZ, 0.f, kMaxFloat) && res;
@@ -85,7 +85,7 @@ std::string L1CAIteration::ToString(int indentLevel) const
aStream << indent << indCh << "Doublet chi2 cut: " << fDoubletChi2Cut << '\n';
aStream << indent << indCh << "Pick gather: " << fPickGather << '\n';
aStream << indent << indCh << "Triplet link chi2: " << fTripletLinkChi2 << '\n';
- aStream << indent << indCh << "Max inverse momentum: " << fMaxInvMom << '\n';
+ aStream << indent << indCh << "Max q / p: " << fMaxQp << '\n';
aStream << indent << indCh << "Max slope at primary vertex: " << fMaxSlopePV << '\n';
aStream << indent << indCh << "Max slope: " << fMaxSlope << '\n';
aStream << indent << indCh << "Max DZ: " << fMaxDZ << '\n';
diff --git a/reco/L1/L1Algo/L1CAIteration.h b/reco/L1/L1Algo/L1CAIteration.h
index 9892efe44581a14a18d478ffb892cf188736af55..1d61d40a388d168ca6a3bf25874c3d4027926a91 100644
--- a/reco/L1/L1Algo/L1CAIteration.h
+++ b/reco/L1/L1Algo/L1CAIteration.h
@@ -75,7 +75,7 @@ public:
float GetMaxDZ() const { return fMaxDZ; }
/// Gets max considered q/p for tracks
- float GetMaxInvMom() const { return fMaxInvMom; }
+ float GetMaxQp() const { return fMaxQp; }
/// Gets max slope (tx\ty) in 3D hit position of a triplet
float GetMaxSlope() const { return fMaxSlope; }
@@ -151,8 +151,7 @@ public:
void SetMaxDZ(float input) { fMaxDZ = input; }
/// Sets max considered q/p for tracks
- /// TODO: Replace with minimum momentum setter (S.Zharko)
- void SetMaxInvMom(float input) { fMaxInvMom = input; }
+ void SetMaxQp(float input) { fMaxQp = input; }
/// Sets max slope (tx\ty) in 3D hit position of a triplet
void SetMaxSlope(float input) { fMaxSlope = input; }
@@ -224,7 +223,7 @@ private:
float fDoubletChi2Cut = 11.3449 * 2.f / 3.f; ///< Doublet chi2 upper cut
float fPickGather = 3.0; ///< Size of region to attach new hits to the created track
float fTripletLinkChi2 = 25.0; ///< Min value of dp^2/dp_error^2, for which two tiplets are neighbours
- float fMaxInvMom = 1.0 / 0.5; ///< Max considered q/p for tracks
+ float fMaxQp = 1.0 / 0.5; ///< Max considered q/p for tracks
float fMaxSlopePV = 1.1; ///< Max slope (tx\ty) in primary vertex
float fMaxSlope = 2.748; ///< Max slope (tx\ty) in 3D hit position of a triplet
float fMaxDZ = 0.f; ///< Correction for accounting overlaping and iff z [cm]
@@ -264,7 +263,7 @@ private:
ar& fDoubletChi2Cut;
ar& fPickGather;
ar& fTripletLinkChi2;
- ar& fMaxInvMom;
+ ar& fMaxQp;
ar& fMaxSlopePV;
ar& fMaxSlope;
ar& fMaxDZ;
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 08afa4422bf19d4df45e58b042e385dda3d36cae..887b191264fa057e03bef51c10cc7479ebd6bd1c 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -1734,7 +1734,7 @@ void L1Algo::CATrackFinder()
fTripletFinalChi2Cut = caIteration.GetTripletFinalChi2Cut();
fPickGather = caIteration.GetPickGather(); //3.0;
fTripletLinkChi2 = caIteration.GetTripletLinkChi2(); //5.0;
- fMaxInvMom = caIteration.GetMaxInvMom(); //1.0 / 0.5; // max considered q/p
+ fMaxInvMom = caIteration.GetMaxQp(); //1.0 / 0.5; // max considered q/p
fMaxSlopePV = caIteration.GetMaxSlopePV(); //1.1;
fMaxSlope = caIteration.GetMaxSlope(); //2.748; // corresponds to 70 grad
diff --git a/reco/L1/L1Algo/L1ConfigExample.yaml b/reco/L1/L1Algo/L1ConfigExample.yaml
index 4ded80ef8f168d90b8352a320807ca80c8f439d6..dccfe35a95653a7e5e2895d1d5130323e0f18bb3 100644
--- a/reco/L1/L1Algo/L1ConfigExample.yaml
+++ b/reco/L1/L1Algo/L1ConfigExample.yaml
@@ -12,83 +12,91 @@ l1:
# Iterations will be used in the same sequence as they defined below
iterations:
- name: "TESTFastPrimIter"
- track_chi2_cut: 10.
- triplet_chi2_cut: 23.4450
- doublet_chi2_cut: 7.56327
- pick_gather: 3.0
- pick_neighbour: 5.0
- min_momentum: 0.5 # GeV/c
- max_slope_pv: 1.1
- max_slope: 2.748
- max_dz: 0. # cm
- min_start_triplet_lvl: 0
- first_station_index: 0
- target_pos_sigma_x: 1. # cm
- target_pos_sigma_y: 1. # cm
- is_primary: true
- is_electron: false
- is_jumped: false
- if_extend_tracks: true
- if_suppress_ghost: false
- is_track_from_triplets: false
+ min_n_hits: 4
+ min_n_hits_sta_0: 3
+ track_chi2_cut: 10.
+ triplet_chi2_cut: 23.4450
+ doublet_chi2_cut: 7.56327
+ pick_gather: 3.0
+ pick_neighbour: 5.0
+ max_qp: 2. # (GeV/c)^-1, 1 / 0.5
+ max_slope_pv: 1.1
+ max_slope: 2.748
+ max_dz: 0. # cm
+ min_start_triplet_lvl: 0
+ first_station_index: 0
+ target_pos_sigma_x: 1. # cm
+ target_pos_sigma_y: 1. # cm
+ is_primary: true
+ is_electron: false
+ is_jumped: false
+ if_extend_tracks: true
+ if_suppress_ghost: false
+ is_track_from_triplets: false
- name: "AllPrimIter"
- track_chi2_cut: 10.
- triplet_chi2_cut: 23.4450
- doublet_chi2_cut: 7.56327
- pick_gather: 4.0
- pick_neighbour: 5.0
- min_momentum: 0.05 # GeV/c
- max_slope_pv: 1.1
- max_slope: 2.748
- max_dz: 0.1 # cm
- min_start_triplet_lvl: 0
- first_station_index: 0
- target_pos_sigma_x: 1. # cm
- target_pos_sigma_y: 1. # cm
- is_primary: true
- is_electron: false
- is_jumped: false
- if_extend_tracks: true
- if_suppress_ghost: false
- is_track_from_triplets: false
+ min_n_hits: 4
+ min_n_hits_sta_0: 3
+ track_chi2_cut: 10.
+ triplet_chi2_cut: 23.4450
+ doublet_chi2_cut: 7.56327
+ pick_gather: 4.0
+ pick_neighbour: 5.0
+ max_qp: 20. # (GeV/c)^-1, 1 / 0.05
+ max_slope_pv: 1.1
+ max_slope: 2.748
+ max_dz: 0.1 # cm
+ min_start_triplet_lvl: 0
+ first_station_index: 0
+ target_pos_sigma_x: 1. # cm
+ target_pos_sigma_y: 1. # cm
+ is_primary: true
+ is_electron: false
+ is_jumped: false
+ if_extend_tracks: true
+ if_suppress_ghost: false
+ is_track_from_triplets: false
- name: "AllPrimJumpIter"
- track_chi2_cut: 10.
- triplet_chi2_cut: 18.7560
- doublet_chi2_cut: 7.56327
- pick_gather: 4.0
- pick_neighbour: 5.0
- min_momentum: 0.1 # GeV/c
- max_slope_pv: 1.1
- max_slope: 2.748
- max_dz: 0.1 # cm
- min_start_triplet_lvl: 0
- first_station_index: 0
- target_pos_sigma_x: 5. # cm
- target_pos_sigma_y: 5. # cm
- is_primary: true
- is_electron: false
- is_jumped: true
- if_extend_tracks: true
- if_suppress_ghost: true
- is_track_from_triplets: false
+ min_n_hits: 4
+ min_n_hits_sta_0: 4
+ track_chi2_cut: 10.
+ triplet_chi2_cut: 18.7560
+ doublet_chi2_cut: 7.56327
+ pick_gather: 4.0
+ pick_neighbour: 5.0
+ max_qp: 10. # (GeV/c)^-1, 1 / 0.1
+ max_slope_pv: 1.1
+ max_slope: 2.748
+ max_dz: 0.1 # cm
+ min_start_triplet_lvl: 0
+ first_station_index: 0
+ target_pos_sigma_x: 5. # cm
+ target_pos_sigma_y: 5. # cm
+ is_primary: true
+ is_electron: false
+ is_jumped: true
+ if_extend_tracks: true
+ if_suppress_ghost: true
+ is_track_from_triplets: false
- name: "AllSecIter"
- track_chi2_cut: 10.
- triplet_chi2_cut: 18.7560
- doublet_chi2_cut: 7.56327
- pick_gather: 4.0
- pick_neighbour: 5.0
- min_momentum: 0.1 # GeV/c
- max_slope_pv: 1.5
- max_slope: 2.748
- max_dz: 0.1 # cm
- min_start_triplet_lvl: 1
- first_station_index: 0
- target_pos_sigma_x: 10. # cm
- target_pos_sigma_y: 10. # cm
- is_primary: false
- is_electron: false
- is_jumped: false
- if_extend_tracks: true
- if_suppress_ghost: true
- is_track_from_triplets: false
+ min_n_hits: 4
+ min_n_hits_sta_0: 4
+ track_chi2_cut: 10.
+ triplet_chi2_cut: 18.7560
+ doublet_chi2_cut: 7.56327
+ pick_gather: 4.0
+ pick_neighbour: 5.0
+ max_qo: 10. # (GeV/c)^-1, 1 / 0.1
+ max_slope_pv: 1.5
+ max_slope: 2.748
+ max_dz: 0.1 # cm
+ min_start_triplet_lvl: 1
+ first_station_index: 0
+ target_pos_sigma_x: 10. # cm
+ target_pos_sigma_y: 10. # cm
+ is_primary: false
+ is_electron: false
+ is_jumped: false
+ if_extend_tracks: true
+ if_suppress_ghost: true
+ is_track_from_triplets: false
...
diff --git a/reco/L1/L1Algo/L1ConfigRW.cxx b/reco/L1/L1Algo/L1ConfigRW.cxx
index 562df5ef82208f2d557121d5da85c2fae064ac1c..1f4fb02bf48813211b8a72d633610e87558ccbdd 100644
--- a/reco/L1/L1Algo/L1ConfigRW.cxx
+++ b/reco/L1/L1Algo/L1ConfigRW.cxx
@@ -15,7 +15,6 @@
#include <yaml-cpp/yaml.h>
-#include "L1CAIteration.h"
#include "L1InitManager.h"
using namespace std::string_literals;
@@ -46,17 +45,10 @@ std::vector<std::string> L1ConfigRW::GetNodeKeys(const YAML::Node& node) const
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1ConfigRW::ReadCAIterations(const YAML::Node& node)
+std::vector<L1CAIteration> L1ConfigRW::ReadCAIterations(const YAML::Node& node) const
{
+ std::vector<L1CAIteration> vIters;
if (node) {
- if (fVerbose > -1) {
- LOG(info) << "L1 config: Reading CA tracking iterations sequence. Default iterations will be ignored";
- }
- fpInitManager->ClearCAIterations();
- fpInitManager->SetCAIterationsNumberCrosscheck(node.size());
- if (fVerbose > 2) { LOG(info) << "L1 config: " << fVerbose << " CA iterations were recorded"; }
- if (fVerbose > 3) { LOG(info) << "L1 config: Recorded iterations:"; }
-
// Loop over input sub-nodes. Each sub-node corresponds to one L1CAIteration object
for (const auto& input : node) {
try {
@@ -67,7 +59,7 @@ void L1ConfigRW::ReadCAIterations(const YAML::Node& node)
caIter.SetDoubletChi2Cut(input["doublet_chi2_cut"].as<float>(caIter.GetDoubletChi2Cut()));
caIter.SetPickGather(input["pick_gather"].as<float>(caIter.GetPickGather()));
caIter.SetTripletLinkChi2(input["triplet_link_chi2"].as<float>(caIter.GetTripletLinkChi2()));
- caIter.SetMaxInvMom(1. / input["min_momentum"].as<float>(caIter.GetMaxInvMom()));
+ caIter.SetMaxQp(input["max_qp"].as<float>(caIter.GetMaxQp()));
caIter.SetMaxSlopePV(input["max_slope_pv"].as<float>(caIter.GetMaxSlopePV()));
caIter.SetMaxSlope(input["max_slope"].as<float>(caIter.GetMaxSlope()));
caIter.SetMaxDZ(input["max_dz"].as<float>(caIter.GetMaxDZ()));
@@ -79,10 +71,10 @@ void L1ConfigRW::ReadCAIterations(const YAML::Node& node)
caIter.SetTrackFromTripletsFlag(input["is_track_from_triplets"].as<bool>(caIter.GetTrackFromTripletsFlag()));
caIter.SetExtendTracksFlag(input["if_extend_tracks"].as<bool>(caIter.GetExtendTracksFlag()));
caIter.SetJumpedFlag(input["is_jumped"].as<bool>(caIter.GetJumpedFlag()));
- caIter.SetMinNhits(input["min_n_hits"].as<bool>(caIter.GetMinNhits()));
- caIter.SetMinNhitsStation0(input["min_n_hits_sta_0"].as<bool>(caIter.GetMinNhitsStation0()));
+ caIter.SetMinNhits(input["min_n_hits"].as<int>(caIter.GetMinNhits()));
+ caIter.SetMinNhitsStation0(input["min_n_hits_sta_0"].as<int>(caIter.GetMinNhitsStation0()));
if (fVerbose > 3) { LOG(info) << "L1 config:\n" << caIter.ToString(1); }
- fpInitManager->PushBackCAIteration(caIter);
+ vIters.push_back(caIter);
}
catch (const YAML::InvalidNode& exc) {
const auto nodeKeys = this->GetNodeKeys(input);
@@ -94,10 +86,7 @@ void L1ConfigRW::ReadCAIterations(const YAML::Node& node)
}
} // Loop over input sub-nodes: end
}
- else {
- LOG(warn) << "L1 config: CA tracking iterations sequence was not found in the parameters file, default will be "
- << "used. To define iterations please use the following path in the YAML file: l1/algorithm/iterations";
- }
+ return vIters;
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -117,7 +106,31 @@ void L1ConfigRW::ReadYaml(const std::string& filename)
}
// Tracking iterations
- this->ReadCAIterations(config["l1"]["algorithm"]["iterations"]);
+ this->InitCAIterations(config["l1"]["algorithm"]["iterations"]);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void L1ConfigRW::InitCAIterations(const YAML::Node& node)
+{
+ assert(fpInitManager); // Class itself can be used without the init manager, but this function cannot
+ auto vIters = this->ReadCAIterations(node);
+ if (vIters.size()) {
+ if (fVerbose > -1) {
+ LOG(info) << "L1 config: Reading CA tracking iterations sequence. Default iterations will be ignored";
+ }
+ fpInitManager->ClearCAIterations();
+ fpInitManager->SetCAIterationsNumberCrosscheck(vIters.size());
+ if (fVerbose > 2) { LOG(info) << "L1 config: " << fVerbose << " CA iterations were recorded"; }
+ if (fVerbose > 3) { LOG(info) << "L1 config: Recorded iterations:"; }
+ for (const auto& caIter : vIters) {
+ fpInitManager->PushBackCAIteration(caIter);
+ }
+ }
+ else {
+ LOG(warn) << "L1 config: CA tracking iterations sequence was not found in the parameters file, default will be "
+ << "used. To define iterations please use the following path in the YAML file: l1/algorithm/iterations";
+ }
}
// ---------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/L1Algo/L1ConfigRW.h b/reco/L1/L1Algo/L1ConfigRW.h
index accde1d92cdda9e78d943cc742983a4173c0691a..0a1a16f40a45191d096d9d5d2c09c3cc2af17881 100644
--- a/reco/L1/L1Algo/L1ConfigRW.h
+++ b/reco/L1/L1Algo/L1ConfigRW.h
@@ -13,6 +13,8 @@
#include <string>
#include <vector>
+#include "L1CAIteration.h"
+
namespace YAML
{
class Node;
@@ -37,10 +39,15 @@ public:
/// \param filename Name of the configuration file
void WriteYaml(const std::string& filename);
-private:
/// Reads CA track finder iterations from YAML node
- /// \param node YAML node
- void ReadCAIterations(const YAML::Node& node);
+ /// \param node YAML node containing the iterations
+ /// \return A vector of iterations
+ std::vector<L1CAIteration> ReadCAIterations(const YAML::Node& node) const;
+
+private:
+ /// Calls CA track finder iterations reader and initialize the iterations in the init manager
+ /// \param node YAML node containing the iterations
+ void InitCAIterations(const YAML::Node& node);
/// Gets parameters content of the node
/// \param node YAML node
diff --git a/reco/L1/L1Algo/L1InitManager.cxx b/reco/L1/L1Algo/L1InitManager.cxx
index 9da08993deeb53ea96d0a3f0048b61b90a9a0acd..dfd1e5eda6f56113924e872588b38f07d7235ab8 100644
--- a/reco/L1/L1Algo/L1InitManager.cxx
+++ b/reco/L1/L1Algo/L1InitManager.cxx
@@ -10,6 +10,7 @@
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
#include <algorithm>
#include <fstream>
@@ -261,12 +262,12 @@ void L1InitManager::PushBackCAIteration(const L1CAIteration& iteration)
void L1InitManager::ReadParametersObject(const std::string& fileName)
{
// Open input binary file
- std::ifstream ifs(fileName, std::ios::binary);
+ std::ifstream ifs(fileName);
if (!ifs) { LOG(fatal) << "L1InitManager: parameters data file \"" << fileName << "\" was not found"; }
// Get L1InputData object
try {
- boost::archive::binary_iarchive ia(ifs);
+ boost::archive::text_iarchive ia(ifs);
ia >> fParameters;
}
catch (const std::exception&) {
@@ -279,12 +280,15 @@ void L1InitManager::ReadParametersObject(const std::string& fileName)
void L1InitManager::ReadSearchWindows(const std::string& fileName)
{
// Open input binary file
- std::ifstream ifs(fileName, std::ios::binary);
+ std::ifstream ifs(fileName);
if (!ifs) { LOG(fatal) << "L1InitManager: search window file \"" << fileName << "\" was not found"; }
try {
- boost::archive::binary_iarchive ia(ifs);
+ boost::archive::text_iarchive ia(ifs);
+ int nPars = -1;
int nWindows = -1;
+ ia >> nPars;
+ assert(nPars == 1); // Currently only the constant windows are available
ia >> nWindows;
std::stringstream errMsg;
for (int iW = 0; iW < nWindows; ++iW) {
diff --git a/reco/L1/L1Algo/L1Parameters.cxx b/reco/L1/L1Algo/L1Parameters.cxx
index 876a9bc5da0425d9f302c8d6e608ab232478665c..2b98dcc874a827a358ff522e7cbfdd227681de7d 100644
--- a/reco/L1/L1Algo/L1Parameters.cxx
+++ b/reco/L1/L1Algo/L1Parameters.cxx
@@ -279,7 +279,16 @@ std::string L1Parameters::ToString(int verbosity, int indentLevel) const
aStream << indent << indentChar << "Doublets vs. MC matching: " << fDevIsMatchDoubletsViaMc << '\n';
aStream << indent << indentChar << "Triplets vs. MC matching: " << fDevIsMatchTripletsViaMc << '\n';
aStream << indent << indentChar << "Use hit search windows: " << fDevIsParSearchWUsed << '\n';
- // TODO: SZh 09.11.2022: Print search windows here
+
+ if (fDevIsParSearchWUsed) {
+ aStream << indent << "SEARCH WINDOWS:\n";
+ for (int iSt = 1; iSt < fNstationsActiveTotal; ++iSt) {
+ for (int iIter = 0; iIter < (int) fCAIterations.size(); ++iIter) {
+ aStream << indent << "- station: " << iSt << ", iteration: " << iIter << '\n';
+ aStream << GetSearchWindow(iSt, iIter).ToString() << '\n';
+ }
+ }
+ }
aStream << indent
<< "--------------------------------------------------------------------------------------------------\n";
diff --git a/reco/L1/L1Algo/L1Parameters.cxx.orig b/reco/L1/L1Algo/L1Parameters.cxx.orig
deleted file mode 100644
index 768656ee986c272ddc898a13267220a0653fec61..0000000000000000000000000000000000000000
--- a/reco/L1/L1Algo/L1Parameters.cxx.orig
+++ /dev/null
@@ -1,360 +0,0 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov, Sergei Zharko [committer] */
-
-/// @file L1Parameters.cxx
-/// @brief Parameter container for the L1Algo library
-/// @since 10.02.2022
-/// @author S.Zharko <s.zharko@gsi.de>
-
-#include "L1Parameters.h"
-
-#include <FairLogger.h>
-
-#include <iomanip>
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-L1Parameters::L1Parameters()
-{
- fActiveStationGlobalIDs.fill(-1); // by default, all stations are inactive, thus all the IDs must be -1
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-<<<<<<< HEAD
-L1Parameters::~L1Parameters() noexcept {}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-L1Parameters::L1Parameters(const L1Parameters& other) noexcept
- : fMaxDoubletsPerSinglet(other.fMaxDoubletsPerSinglet)
- , fMaxTripletPerDoublets(other.fMaxTripletPerDoublets)
- , fCAIterations(other.fCAIterations)
- , fTargetPos(other.fTargetPos)
- , fVertexFieldValue(other.fVertexFieldValue)
- , fVertexFieldRegion(other.fVertexFieldRegion)
- , fStations(other.fStations)
- , fThickMap(other.fThickMap)
- , fNstationsGeometryTotal(other.fNstationsGeometryTotal)
- , fNstationsActiveTotal(other.fNstationsActiveTotal)
- , fNstationsGeometry(other.fNstationsGeometry)
- , fNstationsActive(other.fNstationsActive)
- , fActiveStationGlobalIDs(other.fActiveStationGlobalIDs)
- , fTrackingLevel(other.fTrackingLevel)
- , fGhostSuppression(other.fGhostSuppression)
- , fMomentumCutOff(other.fMomentumCutOff)
- , fDevIsIgnoreHitSearchAreas(other.fDevIsIgnoreHitSearchAreas)
- , fDevIsUseOfOriginalField(other.fDevIsUseOfOriginalField)
- , fDevIsMatchDoubletsViaMc(other.fDevIsMatchDoubletsViaMc)
- , fDevIsMatchTripletsViaMc(other.fDevIsMatchTripletsViaMc)
-{
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-L1Parameters& L1Parameters::operator=(const L1Parameters& other) noexcept
-{
- if (this != &other) { L1Parameters(other).Swap(*this); }
- return *this;
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-L1Parameters::L1Parameters(L1Parameters&& other) noexcept { this->Swap(other); }
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-L1Parameters& L1Parameters::operator=(L1Parameters&& other) noexcept
-{
- if (this != &other) {
- L1Parameters tmp(std::move(other));
- this->Swap(tmp);
- }
- return *this;
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-void L1Parameters::Swap(L1Parameters& other) noexcept
-{
- std::swap(fMaxDoubletsPerSinglet, other.fMaxDoubletsPerSinglet);
- std::swap(fMaxTripletPerDoublets, other.fMaxTripletPerDoublets);
- std::swap(fCAIterations, other.fCAIterations);
- std::swap(fTargetPos, other.fTargetPos);
- std::swap(fVertexFieldValue, other.fVertexFieldValue);
- std::swap(fVertexFieldRegion, other.fVertexFieldRegion);
- std::swap(fStations, other.fStations);
- std::swap(fThickMap, other.fThickMap);
- std::swap(fNstationsGeometryTotal, other.fNstationsGeometryTotal);
- std::swap(fNstationsActiveTotal, other.fNstationsActiveTotal);
- std::swap(fNstationsGeometry, other.fNstationsGeometry);
- std::swap(fNstationsActive, other.fNstationsActive);
- std::swap(fActiveStationGlobalIDs, other.fActiveStationGlobalIDs);
- std::swap(fTrackingLevel, other.fTrackingLevel);
- std::swap(fGhostSuppression, other.fGhostSuppression);
- std::swap(fMomentumCutOff, other.fMomentumCutOff);
- std::swap(fDevIsIgnoreHitSearchAreas, other.fDevIsIgnoreHitSearchAreas);
- std::swap(fDevIsUseOfOriginalField, other.fDevIsUseOfOriginalField);
- std::swap(fDevIsMatchDoubletsViaMc, other.fDevIsMatchDoubletsViaMc);
- std::swap(fDevIsMatchTripletsViaMc, other.fDevIsMatchTripletsViaMc);
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-=======
->>>>>>> L1: explicit implementation of copy and move constructors and assignment operators were replaced with default for L1Parameters and L1CAIteration classes
-void L1Parameters::CheckConsistency() const
-{
- LOG(info) << "Consistency test for L1 parameters object... ";
- /*
- * Arrays containing numbers of stations
- *
- * In the fNstationsActive and fNstationsGeometry array first L1Constants::size::kMaxNdetectors elements are the numbers of stations
- * in particular detector subsystem. The last element in both arrays corresponds to the total number of stations (geometry or
- * active). This fact applies restriction on the arrays: the last element must be equal to the sum of all previous elements.
- *
- */
-
- if (std::accumulate(fNstationsGeometry.cbegin(), fNstationsGeometry.cend(), 0) != fNstationsGeometryTotal) {
- throw std::logic_error("L1Parameters: invalid object condition: total number of stations provided by geometry "
- "differs from the sum of the station numbers for individual detector subsystems");
- };
-
- if (std::accumulate(fNstationsActive.cbegin(), fNstationsActive.cend(), 0) != fNstationsActiveTotal) {
- throw std::logic_error("L1Parameters: invalid object condition: total number of stations active in tracking "
- "differs from the sum of the station numbers for individual detector subsystems");
- };
-
- /*
- * Array of active station IDs
- *
- * In the array of active station IDs,
- * (i) sum of elements, which are not equal -1, must be equal the number of stations,
- * (ii) subsequence of all the elements, which are not equal -1, must be a gapless subset of integer numbers starting with 0
- */
-
- auto filterInactiveStationIDs = [](int x) { return x != -1; };
- int nStationsCheck =
- std::count_if(fActiveStationGlobalIDs.cbegin(), fActiveStationGlobalIDs.cend(), filterInactiveStationIDs);
- if (nStationsCheck != fNstationsActiveTotal) {
- std::stringstream msg;
- msg << "L1Parameters: invalid object condition: array of active station IDs is not consistent "
- << "with the total number of stations (" << nStationsCheck << " vs. " << fNstationsActiveTotal << ' '
- << "expected)";
- throw std::logic_error(msg.str());
- }
-
- // Check, if the subsequence of all the elements, which are not equal -1, must be a gapless subset of integer numbers
- // starting from 0. If it is, the testValue in the end must be equal the number of non -1 elements
-
- std::vector<int> idsCheck(
- nStationsCheck); // temporary vector containing a sequence of active station IDs without -1 elements
- std::copy_if(fActiveStationGlobalIDs.cbegin(), fActiveStationGlobalIDs.cend(), idsCheck.begin(),
- filterInactiveStationIDs);
- bool isStationIDsOk = true;
- for (int id = 0; id < nStationsCheck; ++id) {
- isStationIDsOk = isStationIDsOk && idsCheck[id] == id;
- }
- if (!isStationIDsOk) {
- std::stringstream msg;
- msg << "L1Parameters: invalid object condition: array of active station IDs is not a gapless subset "
- << "of integer numbers starting from 0:\n\t";
- for (auto id : fActiveStationGlobalIDs) {
- msg << std::setw(3) << std::setfill(' ') << id << ' ';
- }
- throw std::logic_error(msg.str());
- }
-
- /*
- * Check magnetic field flags of the stations
- *
- * In a current version of tracking there are three configurations possible to be proceeded:
- * A. All the stations are inside magnetic field
- * B. There is no magnetic field in a setup
- * C. All the first stations are inside magnetic field, all the last stations are outside the field
- * In all the cases the fieldStatus flags should be sorted containing all non-zero elements in the beginning
- * (representing stations placed into magnetic field) and all zero elements in the end of z-axis.
- */
- bool ifFieldStatusFlagsOk = std::is_sorted(
- fStations.cbegin(), fStations.cbegin() + fNstationsActiveTotal,
- [&](const L1Station& lhs, const L1Station& rhs) { return bool(lhs.fieldStatus) > bool(rhs.fieldStatus); });
-
- if (!ifFieldStatusFlagsOk) {
- std::stringstream msg;
- msg << "L1Parameters: invalid object condition: L1 tracking is impossible for a given field configuration:\n";
- for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
- msg << "- station ID: " << iSt << ", field status: " << fStations[iSt].fieldStatus << '\n';
- }
- throw std::logic_error(msg.str());
- }
-
- /*
- * Check target position SIMD vector
- */
-
- L1Utils::CheckSimdVectorEquality(fTargetPos[0], "L1Parameters: target position x");
- L1Utils::CheckSimdVectorEquality(fTargetPos[1], "L1Parameters: target position y");
- L1Utils::CheckSimdVectorEquality(fTargetPos[2], "L1Parameters: target position z");
-
- /*
- * Check vertex field region and value objects at primary vertex
- */
-
- fVertexFieldValue.CheckConsistency();
- fVertexFieldRegion.CheckConsistency();
-
-
- /*
- * Check if each station object is consistent itself, and if all of them are placed after the target
- * NOTE: If a station was not set up, it is accounted inconsistent (uninitialized). In the array of stations there are uninitialized
- * stations possible (with id > NstationsActiveTotal), thus one should NOT run the loop above over all the stations in array
- * but only until *(fNstationsActive.cend() - 1) (== NstationsActiveTotal).
- * TODO: Probably, we should introduce methods, which check the consistency of fully initialized objects such as L1Station,
- * L1MaterialInfo, etc. (S.Zharko)
- */
-
- for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
- fStations[iSt].CheckConsistency();
- if (fStations[iSt].z[0] < fTargetPos[2][0]) {
- std::stringstream msg;
- msg << "L1Parameters: station with global ID = " << iSt << " is placed before target "
- << "(z_st = " << fStations[iSt].z[0] << " [cm] < z_targ = " << fTargetPos[2][0] << " [cm])";
- throw std::logic_error(msg.str());
- }
- }
-
- /*
- * Check thickness maps
- * NOTE: If a L1Material map was not set up, it is accounted inconsistent (uninitialized). In the array of thickness maps for each
- * there are uninitialized elements possible (with id > NstationsActiveTotal), thus one should NOT run the loop above over
- * all the stations in array but only until *(fNstationsActive.cend() - 1) (== NstationsActiveTotal).
- */
-
- for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
- fThickMap[iSt].CheckConsistency();
- }
-
-
- /*
- * Check iterations sequence
- * 1. Number of iterations should be larger then zero
- * 2. Each iteration should contain values within predefined limits
- * 3. Number of iterations with TrackFromTriplets flag turned on no more then 1
- * 4. If the TrackFromTriplets iteration exists, it should be the last one in the sequence
- */
- {
- int nIterations = fCAIterations.size();
- if (!nIterations) {
- std::stringstream msg;
- msg << "L1Parameters: 0 track finder iterations were found. Please, define at least one iteration";
- throw std::logic_error(msg.str());
- }
-
- std::string names = "";
- for (const auto iter : fCAIterations) {
- if (!iter.Check()) { names += iter.GetName() + " "; }
- }
- if (names.size()) {
- std::stringstream msg;
- msg << "L1Parameters: some parameters are out of range for the following iterations: " << names;
- throw std::logic_error(msg.str());
- }
-
- nIterations = std::count_if(fCAIterations.begin(), fCAIterations.end(),
- [=](const L1CAIteration& it) { return it.GetTrackFromTripletsFlag(); });
- if (nIterations > 1) {
- std::stringstream msg;
- msg << "L1Parameters: found " << nIterations << " iterations with GetTrackFromTripletsFlag() == true:\n";
- for (const auto& iter : fCAIterations) {
- if (iter.GetTrackFromTripletsFlag()) { msg << '\t' << "- " << iter.GetName() << '\n'; }
- }
- msg << "Only the one iteration can have GetTrackFromTripletsFlag() == true, and this iteration should be ";
- msg << "the last one";
- throw std::logic_error(msg.str());
- }
-
- if (nIterations == 1 && !(fCAIterations.end() - 1)->GetTrackFromTripletsFlag()) {
- std::stringstream msg;
- msg << "L1Parameters: iteration with GetTrackFromTripletsFlag() == true is not the last in a sequence. ";
- msg << "The GetTrackFromTripletsFlag() value in the iterations sequence: \n";
- for (const auto& iter : fCAIterations) {
- msg << '\t' << "- " << std::setw(15) << std::setfill(' ') << iter.GetName() << ' ';
- msg << std::setw(6) << std::setfill(' ') << iter.GetTrackFromTripletsFlag() << '\n';
- }
- throw std::logic_error(msg.str());
- }
- }
-
- LOG(info) << "Consistency test for L1 parameters object... \033[1;32mpassed\033[0m";
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-void L1Parameters::Print(int /*verbosityLevel*/) const
-{
- LOG(info) << ToString();
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-std::string L1Parameters::ToString(int verbosity, int indentLevel) const
-{
- // FIXME: SZh: Fill it with other parameters
- std::stringstream aStream {};
- constexpr char indentChar = '\t';
- std::string indent(indentLevel, indentChar);
- aStream << '\n';
- aStream << indent << "--------------------------------------------------------------------------------\n";
- aStream << indent << " L1 parameters list\n";
- aStream << indent << "--------------------------------------------------------------------------------\n";
- aStream << indent << "COMPILE TIME CONSTANTS:\n";
- aStream << indent << indentChar << "Bits to code one station: " << L1Constants::size::kStationBits << '\n';
- aStream << indent << indentChar << "Bits to code one thread: " << L1Constants::size::kThreadBits << '\n';
- aStream << indent << indentChar << "Bits to code one triplet: " << L1Constants::size::kTripletBits << '\n';
- aStream << indent << indentChar << "Max number of stations: " << L1Constants::size::kMaxNstations << '\n';
- aStream << indent << indentChar << "Max number of threads: " << L1Constants::size::kMaxNthreads << '\n';
- aStream << indent << indentChar << "Max number of triplets: " << L1Constants::size::kMaxNtriplets << '\n';
- aStream << indent << "RUNTIME CONSTANTS:\n";
- aStream << indent << indentChar << "Max number of doublets per singlet: " << fMaxDoubletsPerSinglet << '\n';
- aStream << indent << indentChar << "Max number of triplets per doublet: " << fMaxTripletPerDoublets << '\n';
- aStream << indent << "CA TRACK FINDER ITERATIONS:\n";
- for (int idx = 0; idx < static_cast<int>(fCAIterations.size()); ++idx) {
- aStream << idx << ") " << fCAIterations[idx].ToString(indentLevel + 1) << '\n';
- }
- aStream << indent << "GEOMETRY:\n";
-
- aStream << indent << indentChar << "TARGET:\n";
- aStream << indent << indentChar << indentChar << "Position:\n";
- for (int dim = 0; dim < 3 /*nDimensions*/; ++dim) {
- aStream << indent << indentChar << indentChar << indentChar << char(120 + dim) << " = " << fTargetPos[dim][0]
- << " cm\n";
- }
- aStream << indent << indentChar << "NUMBER OF STATIONS:\n";
- aStream << indent << indentChar << indentChar << "Number of stations (Geometry): ";
- for (int idx = 0; idx < int(fNstationsGeometry.size()); ++idx) {
- aStream << std::setw(2) << std::setfill(' ') << fNstationsGeometry[idx] << ' ';
- }
- aStream << " | total = " << std::setw(2) << std::setfill(' ') << fNstationsGeometryTotal;
- aStream << '\n';
- aStream << indent << indentChar << indentChar << "Number of stations (Active): ";
- for (int idx = 0; idx < int(fNstationsActive.size()); ++idx) {
- aStream << std::setw(2) << std::setfill(' ') << fNstationsActive[idx] << ' ';
- }
- aStream << " | total = " << std::setw(2) << std::setfill(' ') << fNstationsActiveTotal;
- aStream << '\n';
- aStream << indent << indentChar << indentChar << "Active station indexes: ";
- for (int idx = 0; idx < *(fNstationsActive.end() - 1); ++idx) {
- aStream << std::setw(3) << std::setfill(' ') << fActiveStationGlobalIDs[idx] << ' ';
- }
- aStream << '\n';
-
- aStream << indent << indentChar << "STATIONS:\n ";
- for (int idx = 0; idx < *(fNstationsActive.end() - 1); ++idx) {
- aStream << indent << indentChar << indentChar << fStations[idx].ToString(verbosity) << '\n';
- }
-
- aStream << indent << "FLAGS:\n";
- aStream << indent << "--------------------------------------------------------------------------------\n";
- return aStream.str();
-}
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
index a86d09947c2a5ec2b8ee1c9005126a43a6c88842..e49462ed83b4f5e9718a5015628b9d838705d95d 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -63,7 +63,7 @@ public:
L1Parameters(L1Parameters&& other) noexcept = default;
/// Move assignment operator
- L1Parameters& operator=(L1Parameters&& other) noexcept = default;
+ L1Parameters& operator=(L1Parameters&& other) = default;
/// Prints configuration
void Print(int verbosityLevel = 0) const;
@@ -140,6 +140,8 @@ public:
/// \param iTrackGr Index of a track group
const L1SearchWindow& GetSearchWindow(int iStation, int iTrackGr) const
{
+ assert(iStation < fNstationsActiveTotal && iStation > 0);
+ assert(iTrackGr < int(fCAIterations.size()));
return fSearchWindows[iTrackGr * L1Constants::size::kMaxNstations + iStation];
}
diff --git a/reco/L1/L1Algo/L1SearchWindow.cxx b/reco/L1/L1Algo/L1SearchWindow.cxx
index c7418a58efa62cf3e8f46f67f9e9cb8b902b3313..ee92ae19ec764e3db6ced0a8fc4b3d9adae1b1c1 100644
--- a/reco/L1/L1Algo/L1SearchWindow.cxx
+++ b/reco/L1/L1Algo/L1SearchWindow.cxx
@@ -4,6 +4,7 @@
#include "L1SearchWindow.h"
+#include <cassert>
#include <iomanip>
#include <sstream>
@@ -22,7 +23,7 @@ L1SearchWindow::L1SearchWindow(int stationID, int trackGrID) : fStationID(statio
// this class is supposed to be used inside the algorithm core
// ---------------------------------------------------------------------------------------------------------------------
//
-void SetParamDxMaxVsX0(int id, float val)
+void L1SearchWindow::SetParamDxMaxVsX0(int id, float val)
{
assert(id > -1 && id < kNpars);
fvParams[kDxMaxVsX0 * kNpars + id] = val;
@@ -30,7 +31,7 @@ void SetParamDxMaxVsX0(int id, float val)
// ---------------------------------------------------------------------------------------------------------------------
//
-void SetParamDxMinVsX0(int id, float val)
+void L1SearchWindow::SetParamDxMinVsX0(int id, float val)
{
assert(id > -1 && id < kNpars);
fvParams[kDxMinVsX0 * kNpars + id] = val;
@@ -38,7 +39,7 @@ void SetParamDxMinVsX0(int id, float val)
// ---------------------------------------------------------------------------------------------------------------------
//
-void SetParamDxMaxVsY0(int id, float val)
+void L1SearchWindow::SetParamDxMaxVsY0(int id, float val)
{
assert(id > -1 && id < kNpars);
fvParams[kDxMaxVsY0 * kNpars + id] = val;
@@ -46,7 +47,7 @@ void SetParamDxMaxVsY0(int id, float val)
// ---------------------------------------------------------------------------------------------------------------------
//
-void SetParamDxMinVsY0(int id, float val)
+void L1SearchWindow::SetParamDxMinVsY0(int id, float val)
{
assert(id > -1 && id < kNpars);
fvParams[kDxMinVsY0 * kNpars + id] = val;
@@ -54,7 +55,7 @@ void SetParamDxMinVsY0(int id, float val)
// ---------------------------------------------------------------------------------------------------------------------
//
-void SetParamDyMaxVsX0(int id, float val)
+void L1SearchWindow::SetParamDyMaxVsX0(int id, float val)
{
assert(id > -1 && id < kNpars);
fvParams[kDyMaxVsX0 * kNpars + id] = val;
@@ -62,7 +63,7 @@ void SetParamDyMaxVsX0(int id, float val)
// ---------------------------------------------------------------------------------------------------------------------
//
-void SetParamDyMinVsX0(int id, float val)
+void L1SearchWindow::SetParamDyMinVsX0(int id, float val)
{
assert(id > -1 && id < kNpars);
fvParams[kDyMinVsX0 * kNpars + id] = val;
@@ -70,7 +71,7 @@ void SetParamDyMinVsX0(int id, float val)
// ---------------------------------------------------------------------------------------------------------------------
//
-void SetParamDyMaxVsY0(int id, float val)
+void L1SearchWindow::SetParamDyMaxVsY0(int id, float val)
{
assert(id > -1 && id < kNpars);
fvParams[kDyMaxVsY0 * kNpars + id] = val;
@@ -78,7 +79,7 @@ void SetParamDyMaxVsY0(int id, float val)
// ---------------------------------------------------------------------------------------------------------------------
//
-void SetParamDyMinVsY0(int id, float val)
+void L1SearchWindow::SetParamDyMinVsY0(int id, float val)
{
assert(id > -1 && id < kNpars);
fvParams[kDyMinVsY0 * kNpars + id] = val;
@@ -91,7 +92,7 @@ std::string L1SearchWindow::ToString() const
std::stringstream msg;
msg << "----- CA hits search window: \n";
msg << "\tstation ID: " << fStationID << '\n';
- msg << "\ttracks group ID: " << fTracksGroupID << '\n';
+ msg << "\ttracks group ID: " << fTrackGroupID << '\n';
msg << "\tparameters:\n";
msg << "\t\t" << std::setw(6) << std::setfill(' ') << "No." << ' ';
msg << std::setw(12) << std::setfill(' ') << "dx_max(x0)" << ' ';
diff --git a/reco/L1/L1Algo/L1SearchWindow.h b/reco/L1/L1Algo/L1SearchWindow.h
index 13dc50c0d910495f3b276b54f9f3ab509e263201..98899d9ce4d85d5201c2a54515adab51ede5b9d0 100644
--- a/reco/L1/L1Algo/L1SearchWindow.h
+++ b/reco/L1/L1Algo/L1SearchWindow.h
@@ -11,6 +11,7 @@
#define L1SearchWindow_h 1
#include <boost/serialization/array.hpp>
+#include <boost/serialization/string.hpp>
#include <array>
#include <string>
@@ -39,7 +40,7 @@ public:
L1SearchWindow& operator=(const L1SearchWindow& other) = default;
/// Move constructor
- L1SearchWindow(L1SearchWindow&& other) = default;
+ L1SearchWindow(L1SearchWindow&& other) noexcept = default;
/// Move assignment operator
L1SearchWindow& operator=(L1SearchWindow&& other) = default;
@@ -75,6 +76,9 @@ public:
/// Gets track group id
int GetTrackGroupID() const { return fTrackGroupID; }
+ /// Sets tag
+ /// A tag can be used for insurance, if this search window corresponds to a desired track finder iteration
+ void SetTag(const char* name) { fsTag = name; }
/// TODO: SZh 08.11.2022: Implement variadic template function
/// TODO: SZh 08.11.2022: Try to reduce number of functions
@@ -151,7 +155,7 @@ private:
int fStationID = -1; ///< Global index of active tracking station
int fTrackGroupID = -1; ///< Index of tracks group
- //std::string fsTag = ""; ///< Tag, containing information on the tracks group (TODO: can be omitted?)
+ std::string fsTag = ""; ///< Tag, containing information on the tracks group (TODO: can be omitted?)
/// Serialization function
friend class boost::serialization::access;
@@ -161,6 +165,7 @@ private:
ar& fvParams;
ar& fStationID;
ar& fTrackGroupID;
+ ar& fsTag;
}
};
diff --git a/reco/L1/catools/CaToolsWFExpression.cxx b/reco/L1/catools/CaToolsWFExpression.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7f9fbd2aa39a5365b9f8578ef0991598def5ff57
--- /dev/null
+++ b/reco/L1/catools/CaToolsWFExpression.cxx
@@ -0,0 +1,211 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+
+#include "CaToolsWFExpression.h"
+
+#include "Logger.h"
+
+#include "TBox.h"
+#include "TF1.h"
+#include "TGraph.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TPad.h"
+#include "TTree.h"
+
+#include <algorithm>
+
+using namespace ca::tools;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+WFExpression::WFExpression(TTree* pTree, const char* exprDist, const char* exprParam)
+ : fpTree(pTree)
+ , fsExprDist(exprDist)
+ , fsExprParam(exprParam)
+{
+ // Assertions
+ // TODO: SZh 11.11.2022: REPLACE ALL THE ASSERTIONS WITH EXCEPTIONS!!!!
+ assert(fsExprDist.Length());
+ assert(fsExprParam.Length());
+ assert(fpTree);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::tuple<std::vector<float>, std::vector<float>> WFExpression::CalculateParameters()
+{
+ // Check, if data were initialized
+ assert(fpPadBase);
+ assert(fpPadSlices);
+
+ // Unique name, used to access a histogram
+ fsName = Form("%s:%s_%s", fsExprDist.Data(), fsExprParam.Data(), fCut.GetTitle());
+ fsName.ReplaceAll(" ", "");
+
+ // Expression to be drawn from the tree
+ TString sExpr = Form("%s:%s>>htemp(%d,,,%d)", fsExprDist.Data(), fsExprParam.Data(), fNbinsParam, fNbinsDist);
+
+ // Get base distribution
+ fpPadBase->cd();
+ fpTree->Draw(sExpr.Data(), fCut, "colz");
+ fpHistBase = (TH2F*) gPad->GetPrimitive("htemp");
+ fpHistBase->Sumw2();
+ fpHistBase->SetStats(kFALSE);
+ fpHistBase->SetName(fsName);
+ fpHistBase->SetTitle(fsTitle);
+
+ // Create slices
+ int sliceSize = fNbinsParam / fNslices;
+ fpPadSlices->cd();
+ fpPadSlices->Divide(fNslices == 1 ? 1 : 2, (fNslices + 1) / 2);
+
+ /// Vectors to store the upper and lower boundaries from each slice
+ for (int iS = 0; iS < fNslices; ++iS) {
+ fpPadSlices->cd(iS + 1);
+ int iBinMin = sliceSize * iS;
+ int iBinMax = sliceSize * (iS + 1) - 1;
+ std::tie(fvLoSBoundaries[iS], fvUpSBoundaries[iS], fvSCenters[iS]) = ProcessSlice(iBinMin, iBinMax);
+ assert(fvLoSBoundaries[iS] <= fvUpSBoundaries[iS]);
+ }
+
+ // Draw edges, obtained for different slices
+ fpPadBase->cd();
+ TGraph* pGrUpper = new TGraph(fNslices, fvSCenters.data(), fvUpSBoundaries.data());
+ pGrUpper->SetMarkerColor(kOrange + 2);
+ pGrUpper->SetMarkerStyle(20);
+ pGrUpper->Draw("PSAME");
+ TGraph* pGrLower = new TGraph(fNslices, fvSCenters.data(), fvLoSBoundaries.data());
+ pGrLower->SetMarkerColor(kOrange + 2);
+ pGrLower->SetMarkerStyle(20);
+ pGrLower->Draw("PSAME");
+
+ // Get window parameters
+ GetConstWindowParams();
+
+ return std::tie(fvLoParams, fvUpParams);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::tuple<float, float, float> WFExpression::ProcessSlice(int iBinMin, int iBinMax)
+{
+ // ----- Create a slice
+ assert(iBinMin <= iBinMax);
+ float lEdge = fpHistBase->GetXaxis()->GetBinLowEdge(iBinMin + 1);
+ float uEdge = fpHistBase->GetXaxis()->GetBinUpEdge(iBinMax + 1);
+ float sliceCenter = (lEdge + uEdge) * 0.5;
+ const char* parTitle = fpHistBase->GetXaxis()->GetTitle();
+ TString sSliceTitle = Form("%.2f < %s < %.2f (bins %d-%d)", lEdge, parTitle, uEdge, iBinMin, iBinMax);
+ TString sSliceName = Form("%s_slice%d-%d", fpHistBase->GetName(), iBinMin, iBinMax);
+ TH1F* pHistSlice = (TH1F*) fpHistBase->ProjectionY(sSliceName, iBinMin, iBinMax);
+ pHistSlice->SetTitle(sSliceTitle);
+
+ // Slice drawing attributes
+ constexpr float maximumScale = 1.1;
+ pHistSlice->SetMaximum(pHistSlice->GetMaximum() * maximumScale);
+ pHistSlice->Draw("");
+
+ // ----- Find selection boundaries
+ // A simple algorithm, which finds boundaries containing eps/2 entries to the left and right from the central region.
+ int contentL = 0;
+ int contentR = 0;
+ float windowMin = 0.f;
+ float windowMax = 0.f;
+ float nEntriesTot = pHistSlice->GetEntries();
+
+ // left side
+ for (int iBin = 1; contentL < nEntriesTot * fEps * 0.5; ++iBin) {
+ windowMin = pHistSlice->GetXaxis()->GetBinLowEdge(iBin);
+ contentL += pHistSlice->GetBinContent(iBin);
+ }
+
+ // right side
+ for (int iBin = pHistSlice->GetXaxis()->GetNbins(); contentR < nEntriesTot * fEps * 0.5; --iBin) {
+ windowMax = pHistSlice->GetXaxis()->GetBinUpEdge(iBin);
+ contentR += pHistSlice->GetBinContent(iBin);
+ }
+
+ // ----- Draw selection box on the histogram
+ TBox* pWindowBox = new TBox(windowMin, pHistSlice->GetMinimum(), windowMax, pHistSlice->GetMaximum());
+ pWindowBox->SetFillColor(kGreen);
+ pWindowBox->Draw("SAME");
+ pHistSlice->Draw("SAME");
+
+ return std::tie(windowMin, windowMax, sliceCenter);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WFExpression::GetConstWindowParams()
+{
+ assert(kNpars == 1); // Cannot be used for a different number of parameters
+ assert(fvUpParams.size() == 1);
+ assert(fvLoParams.size() == 1);
+ // In case of several slices we select the most outlying point as a boundary of the window
+ fvUpParams[0] = *(std::max_element(fvUpSBoundaries.begin(), fvUpSBoundaries.end()));
+ fvLoParams[0] = *(std::min_element(fvLoSBoundaries.begin(), fvLoSBoundaries.end()));
+
+ // ----- Draw the estimated parameterisations
+ fpPadBase->cd();
+ float xMin = fpHistBase->GetXaxis()->GetBinLowEdge(1);
+ float xMax = fpHistBase->GetXaxis()->GetBinUpEdge(fpHistBase->GetXaxis()->GetNbins());
+ TF1* fMin = new TF1("fMin", "[0]", xMin, xMax);
+ fMin->SetParameter(0, fvLoParams[0]);
+ fMin->SetLineColor(kOrange + 2);
+ fMin->Draw("LSAME");
+ TF1* fMax = new TF1("fMax", "[0]", xMin, xMax);
+ fMax->SetParameter(0, fvUpParams[0]);
+ fMax->SetLineColor(kOrange + 2);
+ fMax->Draw("LSAME");
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WFExpression::SetEps(float eps)
+{
+ assert(eps > 0 && eps < 1.0);
+ fEps = eps;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WFExpression::SetNslices(int nSlices)
+{
+ assert(nSlices > 0);
+ fNslices = nSlices;
+ fvUpSBoundaries.clear();
+ fvLoSBoundaries.clear();
+ fvSCenters.clear();
+ fvUpSBoundaries.resize(fNslices);
+ fvLoSBoundaries.resize(fNslices);
+ fvSCenters.resize(fNslices);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WFExpression::SetNbins(int nBinsDist, int nBinsParam)
+{
+ assert(nBinsDist > 0);
+ assert(nBinsParam > 0);
+ fNbinsDist = nBinsDist;
+ fNbinsParam = nBinsParam;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WFExpression::SetPadBase(TPad* pad)
+{
+ assert(pad);
+ fpPadBase = pad;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WFExpression::SetPadSlices(TPad* pad)
+{
+ assert(pad);
+ fpPadSlices = pad;
+}
diff --git a/reco/L1/catools/CaToolsWFExpression.h b/reco/L1/catools/CaToolsWFExpression.h
new file mode 100644
index 0000000000000000000000000000000000000000..99bd7f7e8bb94d89b8d7ffcbb42b4831925763af
--- /dev/null
+++ b/reco/L1/catools/CaToolsWFExpression.h
@@ -0,0 +1,131 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+#ifndef CaToolsWFExpression_h
+#define CaToolsWFExpression_h 1
+
+#include "TCut.h"
+#include "TString.h"
+
+#include <cassert>
+#include <tuple>
+#include <vector>
+
+class TTree;
+class TH2F;
+class TPad;
+
+namespace ca
+{
+ namespace tools
+ {
+ /// A helper class for ca::tools::WindowFinder to handle a particular expression (dx vs. x0 etc.) and all related
+ /// methods to work with it.
+ /// \note DISTANCE (expression, axis) -- dx or dy, signed distance between the extrapolated point and real MC-point
+ /// PARAMETER (expression, axis) -- x0 or y0, a parameter, vs. which the distance is studied
+ class WFExpression {
+ public:
+ static constexpr int kNpars = 1; /// TMP: number of parameters
+
+ /// Default constructor
+ WFExpression() = delete;
+
+ /// Constructor
+ /// \param pChain A pointer to a tree with MC triplets
+ /// \param exprDist Expression for distance
+ /// \param exprParam Expression for parameter
+ WFExpression(TTree* pTree, const char* exprDist, const char* exprParam);
+
+ /// Copy constructor
+ WFExpression(const WFExpression& other) = delete;
+
+ /// Move constructor
+ WFExpression(WFExpression&& other) = delete;
+
+ /// Destructor
+ ~WFExpression() = default;
+
+ /// Copy assignment operator
+ WFExpression& operator=(const WFExpression& other) = delete;
+
+ /// Move assignment operator
+ WFExpression& operator=(WFExpression&& other) = delete;
+
+ /// Calculates parameters
+ /// \param pTree Pointer to a tree with MC triplets
+ /// \return A tuple:
+ /// - vector of parameters for upper bound
+ /// - vector of parameters for lower bound
+ std::tuple<std::vector<float>, std::vector<float>> CalculateParameters();
+
+ /// Sets cut, including information on the station and track finder iteration
+ void SetCut(const TCut& cut) { fCut = cut; }
+
+ /// Sets fraction of the events, which can be left out of boundaries
+ void SetEps(float eps);
+
+ /// Sets number of slices
+ void SetNslices(int nSlices);
+
+ /// Sets number of bins
+ /// \param nBinsDist Number of bins along the distance axis
+ /// \param nBinsParam Number of bins along the parameter axis
+ void SetNbins(int nBinsDist, int nBinsParam);
+
+ /// Sets base pad pointer
+ void SetPadBase(TPad* pad);
+
+ /// Sets slices pad pointer
+ void SetPadSlices(TPad* pad);
+
+ /// Sets title of the histograms
+ void SetTitle(const char* title) { fsTitle = title; }
+
+ private:
+ // *****************************
+ // ** Private class functions **
+ // *****************************
+
+ /// Process slice
+ /// \param iBinMax Max bin of the parameter axis to start a projection (date from the bin are included)
+ /// \return Tuple: lower bound, upper bound, slice center
+ std::tuple<float, float, float> ProcessSlice(int iBinMin, int iBinMax);
+
+ /// Gets window parameterisations assuming there is no dependence from parameter
+ void GetConstWindowParams();
+ // TODO: use other functions for other window shapes: GetParabWindowParams, GetEllipticWindowParams etc.
+
+
+ // *********************
+ // ** Class variables **
+ // *********************
+ TTree* fpTree = nullptr; ///< Tree to be analyzed
+ TH2F* fpHistBase = nullptr; ///< Base histogram (distance vs. parameter (x0 or y0))
+
+ TString fsExprDist = ""; ///< Expression along the distance axis
+ TString fsExprParam = ""; ///< Expression along the parameter axis
+
+ TCut fCut = ""; ///< Cut used to draw and expression
+ int fNslices = 8; ///< Number of slices along the parameter axis
+ float fEps = 0.0005; ///< A fraction of triplets, which can be lost
+
+
+ std::vector<float> fvUpSBoundaries = std::vector<float>(fNslices); ///< Upper boundaries for diff. slices
+ std::vector<float> fvLoSBoundaries = std::vector<float>(fNslices); ///< Lower boundaries for diff. slices
+ std::vector<float> fvSCenters = std::vector<float>(fNslices); ///< Slice centers
+ std::vector<float> fvUpParams = std::vector<float>(kNpars); ///< Parameters for max
+ std::vector<float> fvLoParams = std::vector<float>(kNpars); ///< Parameters for min
+
+ // ----- Plotting options
+ int fNbinsParam = 400; ///< Number of bins along the parameter axis
+ int fNbinsDist = 400; ///< Number of bins along the distance axis
+ TString fsTitle = ""; ///< Title of expression
+ TString fsName = ""; ///< Name of the expression (expr + cut, must be unique!!, TODO: make a check)
+ TPad* fpPadBase = nullptr; ///< Pointer to a pad for base histogram
+ TPad* fpPadSlices = nullptr; ///< Pointer to a pad for slices
+ };
+ } // namespace tools
+} // namespace ca
+
+#endif // CaToolsWFExpression_h
diff --git a/reco/L1/catools/CaToolsWindowFinder.cxx b/reco/L1/catools/CaToolsWindowFinder.cxx
index d663c7cc5bd78b0e2dfd9fc0b21b4c00904cf814..a46a424349644c7f438ef353453de6b98c699b79 100644
--- a/reco/L1/catools/CaToolsWindowFinder.cxx
+++ b/reco/L1/catools/CaToolsWindowFinder.cxx
@@ -9,11 +9,25 @@
#include "CaToolsWindowFinder.h"
-#include <Logger.h>
+#include "Logger.h"
+#include "TCanvas.h"
#include "TChain.h"
+#include "TPad.h"
+#include "TPaveText.h"
+#include "TString.h"
+#include <boost/archive/text_oarchive.hpp>
+
+#include <fstream>
+#include <sstream>
+
+#include <yaml-cpp/yaml.h>
+
+#include "CaToolsWFExpression.h"
+#include "L1ConfigRW.h"
#include "L1Constants.h"
+#include "L1SearchWindow.h"
using namespace ca::tools;
using namespace L1Constants::clrs; // for colored logs
@@ -41,7 +55,237 @@ void WindowFinder::AddInputFile(const char* filename)
// ---------------------------------------------------------------------------------------------------------------------
//
-void WindowFinder::Process(EBuildingMode) { fpMcTripletsTree->Draw("zv"); }
+L1SearchWindow WindowFinder::CreateSW(int iStation, const L1CAIteration& caIter)
+{
+
+ // Get a cut
+ TCut cut = GetTrackSelectionCut(iStation, caIter);
+
+ if (TString(fExtraCut.GetTitle()).Length()) { cut = cut && fExtraCut; }
+
+ // Set a unique to a canvas using the cut expression
+ TString sUniqueName = cut.GetTitle();
+ sUniqueName.ReplaceAll(" ", "");
+
+ // Prepare a canvas and pads to draw the histograms
+ TString canvName = TString("canv_") + sUniqueName;
+ TCanvas* canv = new TCanvas(canvName, canvName, 3000, 2000);
+ fvpCanvases.push_back(canv); // Register canvas
+
+ // Pads for base distributions
+ TPad* padBase[EExpr::kEND]; /// dx vs. x0, dx vs. y0, dy vs. x0, dy vs. y0
+ padBase[EExpr::kDxVsX0] = new TPad("padBase_dx_x0", "padBase_dx_x0", 0.00, 0.50, 0.25, 0.80);
+ padBase[EExpr::kDxVsY0] = new TPad("padBase_dx_y0", "padBase_dx_y0", 0.25, 0.50, 0.50, 0.80);
+ padBase[EExpr::kDyVsX0] = new TPad("padBase_dy_x0", "padBase_dy_x0", 0.50, 0.50, 0.75, 0.80);
+ padBase[EExpr::kDyVsY0] = new TPad("padBase_dy_y0", "padBase_dy_y0", 0.75, 0.50, 1.00, 0.80);
+
+ // Pads for slices
+ TPad* padSlices[EExpr::kEND]; /// projections
+ padSlices[EExpr::kDxVsX0] = new TPad("padSlices_dx_x0", "padSlices_dx_x0", 0.00, 0.00, 0.25, 0.50);
+ padSlices[EExpr::kDxVsY0] = new TPad("padSlices_dx_y0", "padSlices_dx_y0", 0.25, 0.00, 0.50, 0.50);
+ padSlices[EExpr::kDyVsX0] = new TPad("padSlices_dy_x0", "padSlices_dy_x0", 0.50, 0.00, 0.75, 0.50);
+ padSlices[EExpr::kDyVsY0] = new TPad("padSlices_dy_y0", "padSlices_dy_y0", 0.75, 0.00, 1.00, 0.50);
+
+ for (int iExpr = 0; iExpr < EExpr::kEND; ++iExpr) {
+ padSlices[iExpr]->Draw();
+ padBase[iExpr]->Draw();
+ }
+
+ // A pad for text
+ TPad* padDescr = new TPad("padDescr", "padDescr", 0.00, 0.80, 1.00, 1.00); /// Definitions of dx and dy, cut
+ padDescr->Draw();
+
+ PrintCaseInformation(padDescr, iStation, caIter);
+
+
+ // Process expressions
+ WFExpression exprDxVsX0(fpMcTripletsTree, GetDistExpr(EExpr::kDxVsX0), "x0");
+ exprDxVsX0.SetEps(fEps);
+ exprDxVsX0.SetNbins(fNbinsY, fNbinsX);
+ exprDxVsX0.SetNslices(fNslices);
+ exprDxVsX0.SetCut(cut);
+ exprDxVsX0.SetTitle("dx vs. x0;x0 [cm];dx [cm]");
+ exprDxVsX0.SetPadBase(padBase[EExpr::kDxVsX0]);
+ exprDxVsX0.SetPadSlices(padSlices[EExpr::kDxVsX0]);
+
+ WFExpression exprDxVsY0(fpMcTripletsTree, GetDistExpr(EExpr::kDxVsY0), "y0");
+ exprDxVsY0.SetEps(fEps);
+ exprDxVsY0.SetNbins(fNbinsY, fNbinsX);
+ exprDxVsY0.SetNslices(fNslices);
+ exprDxVsY0.SetCut(cut);
+ exprDxVsY0.SetTitle("dx vs. y0;y0 [cm];dx [cm]");
+ exprDxVsY0.SetPadBase(padBase[EExpr::kDxVsY0]);
+ exprDxVsY0.SetPadSlices(padSlices[EExpr::kDxVsY0]);
+
+ WFExpression exprDyVsX0(fpMcTripletsTree, GetDistExpr(EExpr::kDyVsX0), "x0");
+ exprDyVsX0.SetEps(fEps);
+ exprDyVsX0.SetNbins(fNbinsY, fNbinsX);
+ exprDyVsX0.SetNslices(fNslices);
+ exprDyVsX0.SetCut(cut);
+ exprDyVsX0.SetTitle("dy vs. x0;x0 [cm];dy [cm]");
+ exprDyVsX0.SetPadBase(padBase[EExpr::kDyVsX0]);
+ exprDyVsX0.SetPadSlices(padSlices[EExpr::kDyVsX0]);
+
+ WFExpression exprDyVsY0(fpMcTripletsTree, GetDistExpr(EExpr::kDyVsY0), "y0");
+ exprDyVsY0.SetEps(fEps);
+ exprDyVsY0.SetNbins(fNbinsY, fNbinsX);
+ exprDyVsY0.SetNslices(fNslices);
+ exprDyVsY0.SetCut(cut);
+ exprDyVsY0.SetTitle("dy vs. y0;y0 [cm];dy [cm]");
+ exprDyVsY0.SetPadBase(padBase[EExpr::kDyVsY0]);
+ exprDyVsY0.SetPadSlices(padSlices[EExpr::kDyVsY0]);
+
+ auto [parsUpperDxVsX0, parsLowerDxVsX0] = exprDxVsX0.CalculateParameters();
+ auto [parsUpperDxVsY0, parsLowerDxVsY0] = exprDxVsY0.CalculateParameters();
+ auto [parsUpperDyVsX0, parsLowerDyVsX0] = exprDyVsX0.CalculateParameters();
+ auto [parsUpperDyVsY0, parsLowerDyVsY0] = exprDyVsY0.CalculateParameters();
+
+ // Prepare a searching window
+ int iCaIter = &caIter - &fvCaIters[0];
+ L1SearchWindow sw(iStation, iCaIter);
+ sw.SetTag(caIter.GetName().c_str());
+ assert(fNparams == 1); // TMP: At the moment only constant windows available
+ for (int iP = 0; iP < fNparams; ++iP) {
+ sw.SetParamDxMaxVsX0(iP, parsUpperDxVsX0[iP]);
+ sw.SetParamDxMinVsX0(iP, parsLowerDxVsX0[iP]);
+ sw.SetParamDxMaxVsY0(iP, parsUpperDxVsY0[iP]);
+ sw.SetParamDxMinVsY0(iP, parsLowerDxVsY0[iP]);
+ sw.SetParamDyMaxVsX0(iP, parsUpperDyVsX0[iP]);
+ sw.SetParamDyMinVsX0(iP, parsLowerDyVsX0[iP]);
+ sw.SetParamDyMaxVsY0(iP, parsUpperDyVsY0[iP]);
+ sw.SetParamDyMinVsY0(iP, parsLowerDyVsY0[iP]);
+ }
+ return sw;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::DumpCanvasesToPdf(const char* filename) const
+{
+ for (int iC = 0; iC < (int) fvpCanvases.size(); ++iC) {
+ TString sCanvName = Form("%s_%d.pdf", filename, iC);
+ fvpCanvases[iC]->SaveAs(sCanvName);
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+const char* WindowFinder::GetDistExpr(EExpr expr) const
+{
+ assert(fTargetPos[0] == 0);
+ assert(fTargetPos[1] == 0);
+ switch (expr) {
+ case EExpr::kDxVsX0:
+ case EExpr::kDxVsY0: return Form("x1 - x0 * (z1 - (%f)) / (z0 - (%f))", fTargetPos[2], fTargetPos[2]);
+ case EExpr::kDyVsX0:
+ case EExpr::kDyVsY0: return Form("y1 - y0 * (z1 - (%f)) / (z0 - (%f))", fTargetPos[2], fTargetPos[2]);
+ default: return ""; // TODO: throw an exception
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+// TODO: SZh 10.11.2022: Should it be a static function??
+TCut WindowFinder::GetTrackSelectionCut(int iStation, const L1CAIteration& caIter) const
+{
+ TCut cut;
+
+ // Select station
+ // NOTE: iStation stands for global index of an active tracking station for which a searching window is estimated.
+ // In the MC triplets tree the station index ("s") is stored for the left station of the triplet. Thus, for
+ // doublets the station is to be selected as "iStation - 1" (search window is estimated on the middle station),
+ // and as "iStation - 2" for the triplets (search window is estimated on the right station).
+ cut = cut && Form("s==%d", iStation - 1);
+
+ // Select q/p
+ cut = cut && Form("abs(q/p)<%f", caIter.GetMaxQp());
+
+ // Select origin (primary/secondary)
+ cut = cut && (caIter.GetPrimaryFlag() ? "processId==0" : "processId!=0");
+
+ // TODO: SZh 10.11.2022: Try to apply other cuts (electron, slope etc.)
+ return cut;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::PrintCaseInformation(TPad* pPad, int iStation, const L1CAIteration& caIter) const
+{
+ pPad->cd();
+
+ TPaveText* textL = new TPaveText(0.01, 0.01, 0.49, 0.99);
+ textL->SetBorderSize(1);
+ textL->SetTextFont(42);
+ textL->SetFillColor(kWhite);
+ textL->SetTextAlign(11);
+ textL->SetTextSize(0.1);
+ textL->AddText(Form("Global index of active station: %d", iStation));
+ textL->AddText(Form("Track finder iteration: %s", caIter.GetName().c_str()));
+ textL->AddText(Form(" - |q/p| < %.2f e(Gev/c)^{-1}", caIter.GetMaxQp()));
+ textL->AddText(Form(" - primary: %s", caIter.GetPrimaryFlag() ? "yes" : "no"));
+ if (TString(fExtraCut.GetTitle()).Length()) { textL->AddText(Form("Optional cut: %s", fExtraCut.GetTitle())); }
+ textL->Draw();
+
+ TPaveText* textR = new TPaveText(0.51, 0.01, 0.99, 0.99);
+ textR->SetBorderSize(1);
+ textR->SetTextFont(42);
+ textR->SetTextAlign(11);
+ textR->SetTextSize(0.1);
+ textR->SetFillColor(kWhite);
+ textR->AddText("Distance expressions:");
+ textR->AddText(Form(" dx = %s", GetDistExpr(EExpr::kDxVsY0)));
+ textR->AddText(Form(" dy = %s", GetDistExpr(EExpr::kDyVsY0)));
+
+ textR->Draw();
+ //pPad->Draw();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::Process(Option_t* /*opt*/)
+{
+ std::ofstream ofs(fsOutputName);
+ if (!ofs) {
+ LOG(error) << "WindowFinder: failed opening file \"" << fsOutputName << " for writing search windows\"";
+ return;
+ }
+
+ boost::archive::text_oarchive oa(ofs);
+ oa << fNparams;
+ oa << int(fvStationIndexes.size() * fvCaIters.size());
+ for (auto iStation : fvStationIndexes) {
+ for (const auto& iter : fvCaIters) {
+ LOG(info) << "WindowFinder: processing global active station index " << kCLb << iStation << kCL
+ << " and track finder iteration " << kCLb << iter.GetName() << kCL;
+ auto sw = CreateSW(iStation, iter);
+ oa << sw;
+ }
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::ReadTrackingIterationsFromYAML(const char* filename)
+{
+ L1ConfigRW reader(/*pInitManager = */ nullptr, 0);
+ YAML::Node config;
+ try {
+ config = YAML::LoadFile(filename);
+ }
+ catch (const YAML::BadFile& exc) {
+ throw std::runtime_error("file does not exist");
+ }
+ catch (const YAML::ParserException& exc) {
+ throw std::runtime_error("file is formatted improperly");
+ }
+
+ // Get iterations
+ fvCaIters = reader.ReadCAIterations(config["l1"]["algorithm"]["iterations"]);
+
+ // Print iterations:
+ for (const auto& iter : fvCaIters) {
+ LOG(info) << iter.ToString() << '\n';
+ }
+}
// ****************************
@@ -74,6 +318,21 @@ void WindowFinder::SetNslices(int nSlices)
fNslices = nSlices;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::SetStationIndexes(const std::vector<int>& vStationIndexes)
+{
+ fvStationIndexes.clear();
+ for (auto iSt : vStationIndexes) {
+ if (iSt < 1) {
+ LOG(warn) << "WindowFinder::SetStationIndexes: attempt to estimate searching window for station with index "
+ << iSt << " < 1. This index will be omitted";
+ continue;
+ }
+ fvStationIndexes.push_back(iSt);
+ }
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
void WindowFinder::SetTarget(double x, double y, double z)
diff --git a/reco/L1/catools/CaToolsWindowFinder.h b/reco/L1/catools/CaToolsWindowFinder.h
index 1ea7e53314524582ddaeacf09030e04b4b48fb79..9ab2dcc0fb7b4e18f4c2bd051caa38b44a6e0ce4 100644
--- a/reco/L1/catools/CaToolsWindowFinder.h
+++ b/reco/L1/catools/CaToolsWindowFinder.h
@@ -10,31 +10,33 @@
#ifndef CaToolsWindowFinder_h
#define CaToolsWindowFinder_h 1
+#include "TCut.h"
#include "TObject.h"
+#include <array>
+#include <vector>
+
+#include "L1CAIteration.h"
+
// TODO: Replace tmp asserts with exceptions
class TChain;
+class L1SearchWindow;
+class TPad;
+class TCanvas;
namespace ca
{
namespace tools
{
-
- /// Enumeration of parameterisations used to describe a search window
- enum class EWindowForm
- {
- kConstant,
- kLinear,
- kParabolic,
- kElliptic
- };
-
- /// Enumeration
- enum class EBuildingMode
+ /// Enumeration to handle processed expressions
+ enum EExpr
{
- kDouplets, ///< Windows are built from doublets and a target
- kTriplets ///< Windows are built from triplets without using target
+ kDxVsX0,
+ kDxVsY0,
+ kDyVsX0,
+ kDyVsY0,
+ kEND
};
/// TODO: ... write an instruction ...
@@ -47,6 +49,7 @@ namespace ca
/// Default constructor
WindowFinder();
+ /// Destructor
virtual ~WindowFinder() = default;
/// Copy and move are forbidden
@@ -59,8 +62,16 @@ namespace ca
/// \note Multiple file input is possible
void AddInputFile(const char* filename);
+ /// Saves canvases to a set of canvases to pdf
+ void DumpCanvasesToPdf(const char* filename = "WFLog") const;
+
/// Process a tree (TEST)
- void Process(EBuildingMode buildingMode);
+ /// \param opt Define options to process:
+ /// 'T' - triplets are used instead of doublets
+ void Process(Option_t* opt = "");
+
+ /// Reads the iterations from YAML config
+ void ReadTrackingIterationsFromYAML(const char* filename);
/// Sets binning of the dx (dy) distribution
/// \param nBinsX Number of bins for the x-axis
@@ -70,11 +81,26 @@ namespace ca
/// Sets a fraction of triplets (doublets), which can be omitted by the window
void SetEpsilon(float eps);
- /// Sets number of slices along the x-axis
+ /// Sets additional cut (can be useful to reduce distribution contamination by outliers)
+ /// \param cut A cut object
+ void SetExtraCut(const TCut& cut) { fExtraCut = cut; }
+
+ /// Sets number of slices along the x-axiso
/// If the number of slices larger then 1, the window bounds will be fitted with a function (...which TODO).
/// Otherwise, the bounds will be represented with constants (independent from x0 or y0)
void SetNslices(int nSlices);
+ /// Sets name of output file with search windows array
+ /// \note Format of output file:
+ /// \note <number of parameters> <number of search windows stored> <array of serialized L1SearchWindow objects>
+ /// \param filename Name of the file
+ void SetOutputName(const char* filename) { fsOutputName = filename; }
+
+ /// Define indexes of stations for which windows are needed
+ /// \note: A stution
+ /// TODO: Get number of stations from a tree ...
+ void SetStationIndexes(const std::vector<int>& vStationIndexes);
+
/// Sets components of the target center position
/// \param x Position x-component [cm]
/// \param y Position y-component [cm]
@@ -82,22 +108,49 @@ namespace ca
void SetTarget(double x, double y, double z);
- // ******************************
- // ** Member variables **
- // ******************************
-
private:
- std::array<double, 3> fTargetPos = {0}; ///< Target position {x, y, z} [cm]
+ /// Creates a search window for a selected station and iteration
+ L1SearchWindow CreateSW(int iStation, const L1CAIteration& caIter);
- // Window extraction settings
+ /// Returns expression for dx or dy to be drawn in a tree
+ const char* GetDistExpr(EExpr expr) const;
+
+ /// Gets a cut for doublets/triplets defined by station and a tracking iteration
+ /// \param iStation Global index of an active station
+ /// \param caIter CA track finder iteration object
+ TCut GetTrackSelectionCut(int iStation, const L1CAIteration& caIter) const;
+
+ /// Prints information on the dx and dy expression as well as used cuts on the pad
+ /// \param pPad A pad to print the information
+ /// \param iStation Global index of an active station
+ /// \param caIter CA track finder iteration object
+ void PrintCaseInformation(TPad* pPad, int iStation, const L1CAIteration& caIter) const;
+
+
+ // *********************
+ // ** Class variables **
+ // *********************
+ int fNparams = 1; ///< number of parameters of the searching window
- int fNbinsX = 360; ///< Number of bins for the X axis of the distribution
- int fNbinsY = 800; ///< Number of bins for the Y axis of the distribution
+ std::string fsOutputName = "SearchWindows.dat"; ///< Name for output file with estimated search windows
+
+ // ----- Input parameters (iterations and stations)
+ std::vector<L1CAIteration> fvCaIters = {}; ///< Tracking iterations
+ std::vector<int> fvStationIndexes = {}; ///< Global indexes of active stations to find the windows
+ std::array<double, 3> fTargetPos = {0}; ///< Target position {x, y, z} [cm]
+
+ TCut fExtraCut = TCut(""); ///< Optional cut on the triplets/doublets
+
+ // Window extraction settings
+ int fNbinsX = 400; ///< Number of bins for the X axis of the distribution
+ int fNbinsY = 400; ///< Number of bins for the Y axis of the distribution
int fNslices = 1; ///< Number of slices along the X axis
float fEps = 0.001; ///< Fraction of triplets (doublets), which can be omitted
TChain* fpMcTripletsTree = nullptr; ///< Chain of trees containing MC triplets, generated in CbmL1 Performance
+ std::vector<TCanvas*> fvpCanvases; // Vector of pointer to cavnases
+
ClassDef(WindowFinder, 0);
};
} // namespace tools