From 3bdd9fbe56191c0e708206310ff1052b989e0e19 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Wed, 9 Nov 2022 20:27:24 +0100
Subject: [PATCH] L1: estimated search windows are introduced to tracking
---
macro/run/run_reco.C | 3 +-
reco/L1/CMakeLists.txt | 2 +
reco/L1/CbmL1.cxx | 25 +-
reco/L1/CbmL1.h | 25 +-
reco/L1/CbmL1MCTrack.cxx | 4 +-
reco/L1/CbmL1MCTrack.h | 34 +--
reco/L1/CbmL1Performance.cxx | 89 ++++--
reco/L1/CbmL1ReadEvent.cxx | 5 +-
reco/L1/L1Algo/L1Constants.h | 4 +
reco/L1/L1Algo/L1InitManager.cxx | 42 +++
reco/L1/L1Algo/L1InitManager.h | 18 +-
reco/L1/L1Algo/L1Parameters.cxx | 3 +
reco/L1/L1Algo/L1Parameters.cxx.orig | 360 ++++++++++++++++++++++++
reco/L1/L1Algo/L1Parameters.h | 30 +-
reco/L1/L1Algo/L1SearchWindow.cxx | 112 ++++++++
reco/L1/L1Algo/L1SearchWindow.h | 167 +++++++++++
reco/L1/L1LinkDef.h | 1 +
reco/L1/catools/CaToolsWindowFinder.cxx | 84 ++++++
reco/L1/catools/CaToolsWindowFinder.h | 106 +++++++
19 files changed, 1047 insertions(+), 67 deletions(-)
create mode 100644 reco/L1/L1Algo/L1Parameters.cxx.orig
create mode 100644 reco/L1/L1Algo/L1SearchWindow.cxx
create mode 100644 reco/L1/L1Algo/L1SearchWindow.h
create mode 100644 reco/L1/catools/CaToolsWindowFinder.cxx
create mode 100644 reco/L1/catools/CaToolsWindowFinder.h
diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C
index 5c3a90adc3..2ed01f00b1 100644
--- a/macro/run/run_reco.C
+++ b/macro/run/run_reco.C
@@ -383,7 +383,8 @@ void run_reco(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice =
auto l1 = (debugWithMC) ? new CbmL1("L1", 2, 3, 0) : new CbmL1("L1", 0);
// L1 configuration file (optional)
//l1->SetInputConfigName(TString(gSystem->Getenv("VMCWORKDIR")) + "/reco/L1/L1Algo/L1ConfigExample.yaml");
- if (debugWithMC) { l1->SetOutputMcTracksNtupleFilename("mcTracksNtuple"); }
+ // 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 b4e445b438..3f45f45072 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -71,6 +71,7 @@ set(SRCS
catools/CaToolsMcDataManager.cxx
catools/CaToolsPerformance.cxx
catools/CaToolsDebugger.cxx
+ catools/CaToolsWindowFinder.cxx
ParticleFinder/CbmL1PFFitter.cxx
ParticleFinder/CbmL1PFMCParticle.cxx
@@ -94,6 +95,7 @@ set(HEADERS
CbmL1Vtx.h
L1Algo/L1Def.h
L1Algo/L1Vector.h
+ catools/CaToolsWindowFinder.h
)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index ca0d359f2e..2cc894d87b 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -825,6 +825,19 @@ InitStatus CbmL1::Init()
//fInitManager.PushBackCAIteration(trackingIterAllSecJump);
}
+ // *******************************
+ // ** Initialize search windows **
+ // *******************************
+
+ if (fsInputSearchWindowsFilename.size()) {
+ fInitManager.DevSetIsParSearchWUsed();
+ fInitManager.ReadSearchWindows(fsInputSearchWindowsFilename);
+ }
+ else {
+ fInitManager.DevSetIsParSearchWUsed(false);
+ }
+
+
// **********************
// ** Set special cuts **
// **********************
@@ -1181,7 +1194,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
EfficienciesPerformance();
HistoPerformance();
TrackFitPerformance();
- if (fsMcTracksOutputFilename.size()) { DumpMCTracksToNtuple(); }
+ if (fsMcTripletsOutputFilename.size()) { DumpMCTripletsToTree(); }
// TimeHist();
/// WriteSIMDKFData();
}
@@ -1230,11 +1243,11 @@ void CbmL1::Finish()
outfile->Delete();
}
- if (fpMcTracksOutFile) {
- fpMcTracksOutFile->cd();
- fpMcTracksTree->Write();
- fpMcTracksOutFile->Close();
- fpMcTracksOutFile->Delete();
+ if (fpMcTripletsOutFile) {
+ fpMcTripletsOutFile->cd();
+ fpMcTripletsTree->Write();
+ fpMcTripletsOutFile->Close();
+ fpMcTripletsOutFile->Delete();
}
gFile = currentFile;
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index f3b5edc187..8219400432 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -263,14 +263,20 @@ public:
/// \param filename Name of the input tracking configuration file
void SetInputConfigName(const char* filename);
+ /// \brief Sets a name for the input search window file
+ /// If the filename is not defined, a default track finder with Kalman filter is used. Otherwise, an experimental
+ /// version of track finder with estimated hit search windows is utilised
+ /// \param filename Name of the input search window file
+ void SetInputSearchWindowFilename(const char* filename) { fsInputSearchWindowsFilename = filename; }
+
/// \brief Sets a name for the output configuration file
/// \param filename Name of the input tracking configuration file
void SetOutputConfigName(const char* filename) { fInitManager.SetOutputConfigName(std::string(filename)); }
- /// \brief Sets output file for MC tracks ntuple
- /// If the filename is empty string, ntuple is not filled
+ /// \brief Sets output file for MC triplets tree
+ /// If the filename is empty string, tree is not filled
/// \param filename Name of the output file name
- void SetOutputMcTracksNtupleFilename(const char* filename) { fsMcTracksOutputFilename = std::string(filename); }
+ void SetOutputMcTripletsTreeFilename(const char* filename) { fsMcTripletsOutputFilename = std::string(filename); }
/// Sets flag: to correct input hits on MC or not
/// \param flag: true - hits will be corrected on MC information
@@ -424,9 +430,9 @@ private:
/// Fills performance histograms
void HistoPerformance();
- /// Writes MC tracks to ntuple
+ /// Writes MC triplets to tree
/// \note Executed only if the filename for MC tracks ntuple output is defined
- void DumpMCTracksToNtuple();
+ void DumpMCTripletsToTree();
// ** STandAlone Package service-functions **
@@ -471,7 +477,8 @@ private:
static void writedir2current(TObject* obj); // help procedure
private:
- std::string fInputDataFilename = ""; ///< File name to read/write input hits
+ std::string fInputDataFilename = ""; ///< File name to read/write input hits
+ std::string fsInputSearchWindowsFilename = ""; ///< File name to read search windows
// ***************************
// ** Member variables list **
@@ -660,9 +667,9 @@ private:
static const int fNGhostHistos = 9;
TH1F* fGhostHisto[fNGhostHistos] {nullptr};
- TFile* fpMcTracksOutFile = nullptr; ///< File to save MC-tracks ntuple
- TNtuple* fpMcTracksTree = nullptr; ///< Ntuple to save MC-tracks
- std::string fsMcTracksOutputFilename = ""; ///< Name of file to save MC-tracks ntuple
+ TFile* fpMcTripletsOutFile = nullptr; ///< File to save MC-triplets tree
+ TTree* fpMcTripletsTree = nullptr; ///< Tree to save MC-triplets
+ std::string fsMcTripletsOutputFilename = ""; ///< Name of file to save MC-triplets tree
int fFindParticlesMode {0}; // 0 - don't run FindParticles
// 1 - run, all MC particle is reco-able
diff --git a/reco/L1/CbmL1MCTrack.cxx b/reco/L1/CbmL1MCTrack.cxx
index 46e57059b3..00c6385bae 100644
--- a/reco/L1/CbmL1MCTrack.cxx
+++ b/reco/L1/CbmL1MCTrack.cxx
@@ -26,7 +26,8 @@
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1Hit.h"
-CbmL1MCTrack::CbmL1MCTrack(double mass_, double q_, TVector3 vr, TLorentzVector vp, int _ID, int _mother_ID, int _pdg)
+CbmL1MCTrack::CbmL1MCTrack(double mass_, double q_, TVector3 vr, TLorentzVector vp, int _ID, int _mother_ID, int _pdg,
+ unsigned int _process_ID)
: mass(mass_)
, q(q_)
, p(vp.P())
@@ -39,6 +40,7 @@ CbmL1MCTrack::CbmL1MCTrack(double mass_, double q_, TVector3 vr, TLorentzVector
, ID(_ID)
, mother_ID(_mother_ID)
, pdg(_pdg)
+ , process_ID(_process_ID)
{
}
diff --git a/reco/L1/CbmL1MCTrack.h b/reco/L1/CbmL1MCTrack.h
index 59c820b832..007b07bbd3 100644
--- a/reco/L1/CbmL1MCTrack.h
+++ b/reco/L1/CbmL1MCTrack.h
@@ -38,7 +38,8 @@ public:
CbmL1MCTrack(int _ID) : ID(_ID) {};
- CbmL1MCTrack(double mass, double q, TVector3 vr, TLorentzVector vp, int ID, int mother_ID, int pdg);
+ CbmL1MCTrack(double mass, double q, TVector3 vr, TLorentzVector vp, int ID, int mother_ID, int pdg,
+ unsigned int procID);
// CbmL1MCTrack(TmpMCPoints &mcPoint, TVector3 vr, TLorentzVector vp, int ID, int mother_ID);
bool IsPrimary() const { return mother_ID < 0; }
@@ -72,21 +73,22 @@ private:
void CalculateIsReconstructable();
public:
- double mass = 0.;
- double q = 0.;
- double p = 0.;
- double x = 0.;
- double y = 0.;
- double z = 0.;
- double px = 0.;
- double py = 0.;
- double pz = 0.;
- double time = 0.;
- int ID = -1;
- int iFile = -1;
- int iEvent = -1;
- int mother_ID = -1;
- int pdg = -1;
+ double mass = 0.;
+ double q = 0.;
+ double p = 0.;
+ double x = 0.;
+ double y = 0.;
+ double z = 0.;
+ double px = 0.;
+ double py = 0.;
+ double pz = 0.;
+ double time = 0.;
+ int ID = -1;
+ int iFile = -1;
+ int iEvent = -1;
+ int mother_ID = -1;
+ int pdg = -1;
+ unsigned int process_ID = (unsigned int) -1;
bool isSignal {0};
L1Vector<int> Points {"CbmL1MCTrack::Points"}; // indices of pints in CbmL1::fvMCPoints
L1Vector<int> Hits {"CbmL1MCTrack::Hits"}; // indices of hits in algo->vHits or L1::vHits
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 3bfd392054..99d15a41f4 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -2251,9 +2251,9 @@ void CbmL1::InputPerformance()
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmL1::DumpMCTracksToNtuple()
+void CbmL1::DumpMCTripletsToTree()
{
- if (!fpMcTracksTree) {
+ if (!fpMcTripletsTree) {
auto* currentDir = gDirectory;
auto* currentFile = gFile;
@@ -2262,16 +2262,16 @@ void CbmL1::DumpMCTracksToNtuple()
boost::filesystem::path p = (FairRunAna::Instance()->GetUserOutputFileName()).Data();
std::string dir = p.parent_path().string();
if (dir.empty()) dir = ".";
- std::string filename = dir + "/" + fsMcTracksOutputFilename + "." + p.filename().string();
+ std::string filename = dir + "/" + fsMcTripletsOutputFilename + "." + p.filename().string();
std::cout << "\033[1;31mSAVING A TREE: " << filename << "\033[0m\n";
- fpMcTracksOutFile = new TFile(filename.c_str(), "RECREATE");
- fpMcTracksOutFile->cd();
- fpMcTracksTree = new TNtuple("t", "McTracks", "motherId:pdg:p:q:zv:s:x0:y0:z0:x1:y1:z1:x2:y2:z2");
- // motherID - id of mother particle (< 0, if primary)
+ fpMcTripletsOutFile = new TFile(filename.c_str(), "RECREATE");
+ fpMcTripletsOutFile->cd();
+ fpMcTripletsTree = new TTree("t", "MC Triplets");
+ // motherId - id of mother particle (< 0, if primary)
// pdg - PDG code of particle
- // p - absolute value of track momentum [GeV/c]
- // q - charge of track [e]
+ // processId - id of the process (ROOT::TMCProcessID)
+ // pq - absolute value of track momentum [GeV/c], divided by charge [e]
// zv - z-component of track vertex [cm]
// s - global index of station
// x0, y0, z0 - position of the left MC point in a triplet [cm]
@@ -2282,6 +2282,43 @@ void CbmL1::DumpMCTracksToNtuple()
gDirectory = currentDir;
}
+ // Variables for tree branches
+ int brMotherId = -1; ///< mother ID of track
+ int brPdg = -1; ///< PDG code of track
+ unsigned int brProcId = (unsigned int) -1; ///< process ID (see ROOT::TMCProcessID)
+ float brP = 0.f; ///< abs. momentum of track [GeV/c]
+ int brQ = 0; ///< charge of track [e]
+ float brVertexZ = 0.f; ///< z-component of the MC track vertex [cm]
+ int brStation = -1; ///< global index of the left MC point station
+
+ float brX0 = 0.f; ///< x-component of the left MC point in a triplet [cm]
+ float brY0 = 0.f; ///< y-component of the left MC point in a triplet [cm]
+ float brZ0 = 0.f; ///< z-component of the left MC point in a triplet [cm]
+ float brX1 = 0.f; ///< x-component of the middle MC point in a triplet [cm]
+ float brY1 = 0.f; ///< y-component of the middle MC point in a triplet [cm]
+ float brZ1 = 0.f; ///< z-component of the middle MC point in a triplet [cm]
+ float brX2 = 0.f; ///< x-component of the right MC point in a triplet [cm]
+ float brY2 = 0.f; ///< y-component of the right MC point in a triplet [cm]
+ float brZ2 = 0.f; ///< z-component of the right MC point in a triplet [cm]
+
+ // Register branches
+ fpMcTripletsTree->Branch("brMotherId", &brMotherId, "motherId/I");
+ fpMcTripletsTree->Branch("brPdg", &brPdg, "pdg/I");
+ fpMcTripletsTree->Branch("brProcId", &brProcId, "processId/i");
+ fpMcTripletsTree->Branch("brP", &brP, "p/F");
+ fpMcTripletsTree->Branch("brQ", &brQ, "q/I");
+ fpMcTripletsTree->Branch("brVertexZ", &brVertexZ, "zv/F");
+ fpMcTripletsTree->Branch("brStation", &brStation, "s/I");
+ fpMcTripletsTree->Branch("brX0", &brX0, "x0/F");
+ fpMcTripletsTree->Branch("brY0", &brY0, "y0/F");
+ fpMcTripletsTree->Branch("brZ0", &brZ0, "z0/F");
+ fpMcTripletsTree->Branch("brX1", &brX1, "x1/F");
+ fpMcTripletsTree->Branch("brY1", &brY1, "y1/F");
+ fpMcTripletsTree->Branch("brZ1", &brZ1, "z1/F");
+ fpMcTripletsTree->Branch("brX2", &brX2, "x2/F");
+ fpMcTripletsTree->Branch("brY2", &brY2, "y2/F");
+ fpMcTripletsTree->Branch("brZ2", &brZ2, "z2/F");
+
struct ReducedMcPoint {
int s; ///< global active index of tracking station
float x; ///< x-component of point position [cm]
@@ -2306,21 +2343,25 @@ void CbmL1::DumpMCTracksToNtuple()
// Condition to collect only triplets without gaps in stations
// TODO: SZh 20.10.2022 Add cases for jump iterations
if (vPoints[i + 1].s == vPoints[i].s + 1 && vPoints[i + 2].s == vPoints[i].s + 2) {
- fpMcTracksTree->Fill(tr.mother_ID, ///< index of mother particle
- tr.pdg, ///< PDG code
- tr.p, ///< absolute value of track momentum [GeV/c]
- tr.q, ///< charge of track [e]
- tr.z, ///< z-position of track vertex [cm]
- vPoints[i].s, ///< global index of active tracking station
- vPoints[i].x, ///< x-component of the left MC point in a triplet [cm]
- vPoints[i].y, ///< y-component of the left MC point in a triplet [cm]
- vPoints[i].z, ///< z-component of the left MC point in a triplet [cm]
- vPoints[i + 1].x, ///< x-component of the middle MC point in a triplet [cm]
- vPoints[i + 1].y, ///< y-component of the middle MC point in a triplet [cm]
- vPoints[i + 1].z, ///< z-component of the middle MC point in a triplet [cm]
- vPoints[i + 2].x, ///< x-component of the right MC point in a triplet [cm]
- vPoints[i + 2].y, ///< y-component of the right MC point in a triplet [cm]
- vPoints[i + 2].z); ///< z-component of the right MC point in a triplet [cm]
+ // Fill MC-triplets tree
+ brMotherId = tr.mother_ID;
+ brPdg = tr.pdg;
+ brProcId = tr.process_ID;
+ brP = tr.p;
+ brQ = tr.q;
+ brVertexZ = tr.z;
+ brStation = vPoints[i].s;
+ brX0 = vPoints[i].x;
+ brY0 = vPoints[i].y;
+ brZ0 = vPoints[i].z;
+ brX1 = vPoints[i + 1].x;
+ brY1 = vPoints[i + 1].y;
+ brZ1 = vPoints[i + 1].z;
+ brX2 = vPoints[i + 2].x;
+ brY2 = vPoints[i + 2].y;
+ brZ2 = vPoints[i + 2].z;
+
+ fpMcTripletsTree->Fill();
}
}
}
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 5ac7aac611..c77594bf4d 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -1368,7 +1368,8 @@ void CbmL1::Fill_vMCTracks()
TLorentzVector vp;
MCTrack->GetStartVertex(vr);
MCTrack->Get4Momentum(vp);
- Int_t pdg = MCTrack->GetPdgCode();
+ Int_t pdg = MCTrack->GetPdgCode();
+ unsigned int processID = MCTrack->GetGeantProcessId();
Double_t q = 0, mass = 0.;
if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0;
@@ -1379,7 +1380,7 @@ void CbmL1::Fill_vMCTracks()
//cout << "mc track " << iMCTrack << " pdg " << pdg << " z " << vr.Z() << endl;
Int_t iTrack = fvMCTracks.size(); //or iMCTrack?
- CbmL1MCTrack T(mass, q, vr, vp, iTrack, mother_ID, pdg);
+ CbmL1MCTrack T(mass, q, vr, vp, iTrack, mother_ID, pdg, processID);
// CbmL1MCTrack T(mass, q, vr, vp, iMCTrack, mother_ID, pdg);
T.time = MCTrack->GetStartT();
T.iFile = iFile;
diff --git a/reco/L1/L1Algo/L1Constants.h b/reco/L1/L1Algo/L1Constants.h
index 77d3a4c03c..6d8e3cbeed 100644
--- a/reco/L1/L1Algo/L1Constants.h
+++ b/reco/L1/L1Algo/L1Constants.h
@@ -37,6 +37,10 @@ namespace L1Constants
constexpr int kMaxNthreads = 1u << kThreadBits; ///< Max number of threads, 2^6 = 64
constexpr int kMaxNtriplets = 1u << kTripletBits; ///< Max number of triplets, 2^20 = 1,048,576
+ /// Max number of track groups
+ /// NOTE: For a "track group" definition see L1Parameters.h, GetSearchWindow function
+ constexpr int kMaxNtrackGroups = 4;
+
// TODO: Clarify the meaning of these coefficients
constexpr int kCoeff = 64 / 4; ///< TODO:
constexpr int kSingletPortionSize = 1024 / kCoeff; ///< portion of left hits
diff --git a/reco/L1/L1Algo/L1InitManager.cxx b/reco/L1/L1Algo/L1InitManager.cxx
index bd84d8effe..9da08993de 100644
--- a/reco/L1/L1Algo/L1InitManager.cxx
+++ b/reco/L1/L1Algo/L1InitManager.cxx
@@ -164,6 +164,8 @@ bool L1InitManager::FormParametersContainer()
}
}
+ if (!fParameters.fDevIsParSearchWUsed) { fInitController.SetFlag(EInitKey::kSearchWindows, true); }
+
// Check initialization
this->CheckInit();
@@ -272,6 +274,46 @@ 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);
+ if (!ifs) { LOG(fatal) << "L1InitManager: search window file \"" << fileName << "\" was not found"; }
+
+ try {
+ boost::archive::binary_iarchive ia(ifs);
+ int nWindows = -1;
+ ia >> nWindows;
+ std::stringstream errMsg;
+ for (int iW = 0; iW < nWindows; ++iW) {
+ L1SearchWindow swBuffer;
+ ia >> swBuffer;
+ int iStationID = swBuffer.GetStationID();
+ int iTrackGrID = swBuffer.GetTrackGroupID();
+ if (iStationID < 0 || iStationID > L1Constants::size::kMaxNstations) {
+ errMsg << "\t- wrong station id for entry " << iW << ": " << iStationID << " (should be between 0 and "
+ << L1Constants::size::kMaxNstations << ")\n";
+ }
+ if (iTrackGrID < 0 || iTrackGrID > L1Constants::size::kMaxNtrackGroups) {
+ errMsg << "\t- wrong track group id for entry " << iW << ": " << iTrackGrID << " (should be between 0 and "
+ << L1Constants::size::kMaxNtrackGroups << ")\n";
+ }
+ fParameters.fSearchWindows[iTrackGrID * L1Constants::size::kMaxNstations + iStationID] = swBuffer;
+ }
+ if (errMsg.str().size()) {
+ LOG(fatal) << "L1InitManager: some errors occurred while reading search windows: " << errMsg.str();
+ }
+ }
+ catch (const std::exception&) {
+ LOG(fatal) << "L1InitManager: search windows file \"" << fileName
+ << "\" has incorrect data format or was corrupted";
+ }
+
+ fInitController.SetFlag(EInitKey::kSearchWindows, true);
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
bool L1InitManager::SendParameters(L1Algo* pAlgo)
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index 35eeb04f17..796877548c 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -83,10 +83,11 @@ private:
kStationsInfo, ///< 5) If all the planned stations were added to the manager
kCAIterationsNumberCrosscheck, ///< 6) If the number of CA track finder is initialized
kCAIterations, ///< 7) If the CA track finder iterations were initialized
- kTrackingLevel, ///< 8)
- kGhostSuppression, ///< 9)
- kMomentumCutOff, ///< 10)
- kEnd ///< 11) [technical] number of entries in the enumeration
+ kSearchWindows, ///< 8) If the hit search windows were initialized
+ kTrackingLevel, ///< 9)
+ kGhostSuppression, ///< 10)
+ kMomentumCutOff, ///< 11)
+ kEnd ///< 12) [technical] number of entries in the enumeration
};
using L1DetectorIDIntMap_t = std::unordered_map<L1DetectorID, int, L1Utils::EnumClassHash>;
@@ -193,6 +194,10 @@ public:
/// \param fileName Name of input file
void ReadParametersObject(const std::string& fileName);
+ /// Reads search windows from file
+ /// \param fileName Name of input file
+ void ReadSearchWindows(const std::string& fileName);
+
/// Sets a set of active tracking detector IDs
void SetActiveDetectorIDs(const L1DetectorIDSet_t& detectorIDs);
@@ -261,6 +266,11 @@ public:
/// Flag to match triplets using Mc information
void DevSetIsMatchTripletsViaMc(bool value = true) { fParameters.fDevIsMatchTripletsViaMc = value; }
+ /// Flag to use estimated hit search windows
+ /// true: estimated search windows will be used in track finder
+ /// false: Kalman filter will be used in track finder
+ void DevSetIsParSearchWUsed(bool value = true) { fParameters.fDevIsParSearchWUsed = value; }
+
private:
// *********************
// ** Private methods **
diff --git a/reco/L1/L1Algo/L1Parameters.cxx b/reco/L1/L1Algo/L1Parameters.cxx
index d2c190ee76..876a9bc5da 100644
--- a/reco/L1/L1Algo/L1Parameters.cxx
+++ b/reco/L1/L1Algo/L1Parameters.cxx
@@ -278,6 +278,9 @@ std::string L1Parameters::ToString(int verbosity, int indentLevel) const
aStream << indent << indentChar << "Non-approx. field is used: " << fDevIsMatchDoubletsViaMc << '\n';
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
+
aStream << indent
<< "--------------------------------------------------------------------------------------------------\n";
return aStream.str();
diff --git a/reco/L1/L1Algo/L1Parameters.cxx.orig b/reco/L1/L1Algo/L1Parameters.cxx.orig
new file mode 100644
index 0000000000..768656ee98
--- /dev/null
+++ b/reco/L1/L1Algo/L1Parameters.cxx.orig
@@ -0,0 +1,360 @@
+/* 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 4fa9318884..a86d09947c 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -18,6 +18,7 @@
#include "L1CAIteration.h"
#include "L1Constants.h"
#include "L1Material.h"
+#include "L1SearchWindow.h"
#include "L1Station.h"
#include "L1Utils.h"
#include "L1Vector.h"
@@ -131,6 +132,17 @@ public:
/// \param iStation Index of station in the active stations container
const L1Station& GetStation(int iStation) const { return fStations[iStation]; }
+ /// Gets a search window for a selected station and track group
+ /// \note For a particular track finder iteration one can select a track group, which is defined by the minimal
+ /// momentum of tracks (fast tracks, all tracks), their vertex (primary or secondary tracks), or by particle
+ /// (electrons, muons, hadrons, etc.)
+ /// \param iStation Global index of active station
+ /// \param iTrackGr Index of a track group
+ const L1SearchWindow& GetSearchWindow(int iStation, int iTrackGr) const
+ {
+ return fSearchWindows[iTrackGr * L1Constants::size::kMaxNstations + iStation];
+ }
+
/// Gets reference to the array of station thickness map
const L1MaterialContainer_t& GetThicknessMaps() const { return fThickMap; }
@@ -238,6 +250,13 @@ private:
/// active index: 0 -1 1 2 -1 3 4 5 6 7 0 0 0 0
alignas(L1Constants::misc::kAlignment) std::array<int, L1Constants::size::kMaxNstations> fActiveStationGlobalIDs {};
+ /// Map of search windows vs. active station global index and tracks group
+ /// The tracks group can be defined by minimum momentum (fast/all tracks), origin (primary/secondary) and particle type
+ /// (electron, muon, all). Other options also can be added
+ alignas(L1Constants::misc::kAlignment)
+ std::array<L1SearchWindow, L1Constants::size::kMaxNstations* L1Constants::size::kMaxNtrackGroups> fSearchWindows =
+ {};
+
int fTrackingLevel = 0; ///< tracking level
int fGhostSuppression = 0; ///< flag: if true, ghost tracks are suppressed
float fMomentumCutOff = 0; ///< minimum momentum of tracks
@@ -246,10 +265,12 @@ private:
// ** Flags for development **
// ***************************
- bool fDevIsIgnoreHitSearchAreas {false}; ///< Process all hits on the station ignoring hit search area
- bool fDevIsUseOfOriginalField {false}; ///< Force use of original field
- bool fDevIsMatchDoubletsViaMc {false}; ///< Flag to match doublets using MC information
- bool fDevIsMatchTripletsViaMc {false}; ///< Flag to match triplets using Mc information
+ bool fDevIsIgnoreHitSearchAreas = false; ///< Process all hits on the station ignoring hit search area
+ bool fDevIsUseOfOriginalField = false; ///< Force use of original field
+ bool fDevIsMatchDoubletsViaMc = false; ///< Flag to match doublets using MC information
+ bool fDevIsMatchTripletsViaMc = false; ///< Flag to match triplets using MC information
+ bool fDevIsParSearchWUsed = false; ///< Flag: when true, the parametrized search windows are used in track
+ ///< finder; when false, the Kalman filter is used instead
/// Serialization function
friend class boost::serialization::access;
@@ -273,6 +294,7 @@ private:
ar& fNstationsGeometry;
ar& fNstationsActive;
ar& fActiveStationGlobalIDs;
+ ar& fSearchWindows;
ar& fTrackingLevel;
ar& fGhostSuppression;
diff --git a/reco/L1/L1Algo/L1SearchWindow.cxx b/reco/L1/L1Algo/L1SearchWindow.cxx
new file mode 100644
index 0000000000..c7418a58ef
--- /dev/null
+++ b/reco/L1/L1Algo/L1SearchWindow.cxx
@@ -0,0 +1,112 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+#include "L1SearchWindow.h"
+
+#include <iomanip>
+#include <sstream>
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+L1SearchWindow::L1SearchWindow(int stationID, int trackGrID) : fStationID(stationID), fTrackGroupID(trackGrID)
+{
+ assert(stationID > -1);
+ assert(trackGrID > -1);
+
+ // Case for constant windows (TEMPORARY: we should add selection of different windows):
+ static_assert(kNpars == 1);
+}
+
+// TODO: SZh 08.11.2022: Probably, we should have the assertions in the InitManager and remove them from here, since
+// this class is supposed to be used inside the algorithm core
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetParamDxMaxVsX0(int id, float val)
+{
+ assert(id > -1 && id < kNpars);
+ fvParams[kDxMaxVsX0 * kNpars + id] = val;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetParamDxMinVsX0(int id, float val)
+{
+ assert(id > -1 && id < kNpars);
+ fvParams[kDxMinVsX0 * kNpars + id] = val;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetParamDxMaxVsY0(int id, float val)
+{
+ assert(id > -1 && id < kNpars);
+ fvParams[kDxMaxVsY0 * kNpars + id] = val;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetParamDxMinVsY0(int id, float val)
+{
+ assert(id > -1 && id < kNpars);
+ fvParams[kDxMinVsY0 * kNpars + id] = val;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetParamDyMaxVsX0(int id, float val)
+{
+ assert(id > -1 && id < kNpars);
+ fvParams[kDyMaxVsX0 * kNpars + id] = val;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetParamDyMinVsX0(int id, float val)
+{
+ assert(id > -1 && id < kNpars);
+ fvParams[kDyMinVsX0 * kNpars + id] = val;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetParamDyMaxVsY0(int id, float val)
+{
+ assert(id > -1 && id < kNpars);
+ fvParams[kDyMaxVsY0 * kNpars + id] = val;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetParamDyMinVsY0(int id, float val)
+{
+ assert(id > -1 && id < kNpars);
+ fvParams[kDyMinVsY0 * kNpars + id] = val;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+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 << "\tparameters:\n";
+ msg << "\t\t" << std::setw(6) << std::setfill(' ') << "No." << ' ';
+ msg << std::setw(12) << std::setfill(' ') << "dx_max(x0)" << ' ';
+ msg << std::setw(12) << std::setfill(' ') << "dx_min(x0)" << ' ';
+ msg << std::setw(12) << std::setfill(' ') << "dx_max(y0)" << ' ';
+ msg << std::setw(12) << std::setfill(' ') << "dx_min(y0)" << ' ';
+ msg << std::setw(12) << std::setfill(' ') << "dy_max(x0)" << ' ';
+ msg << std::setw(12) << std::setfill(' ') << "dy_min(x0)" << ' ';
+ msg << std::setw(12) << std::setfill(' ') << "dy_max(y0)" << ' ';
+ msg << std::setw(12) << std::setfill(' ') << "dy_min(y0)" << '\n';
+ for (int iPar = 0; iPar < kNpars; ++iPar) {
+ msg << "\t\t" << std::setw(6) << std::setfill(' ') << iPar << ' ';
+ for (int iDep = 0; iDep < kNdeps; ++iDep) {
+ msg << std::setw(12) << std::setfill(' ') << fvParams[iDep * kNpars + iPar] << ' ';
+ }
+ }
+ return msg.str();
+}
diff --git a/reco/L1/L1Algo/L1SearchWindow.h b/reco/L1/L1Algo/L1SearchWindow.h
new file mode 100644
index 0000000000..13dc50c0d9
--- /dev/null
+++ b/reco/L1/L1Algo/L1SearchWindow.h
@@ -0,0 +1,167 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file L1SearchWindow.h
+/// \brief Provides parameterisation for hits searching window in the CA tracking (header)
+/// \date 08.11.2022
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#ifndef L1SearchWindow_h
+#define L1SearchWindow_h 1
+
+#include <boost/serialization/array.hpp>
+
+#include <array>
+#include <string>
+
+/// TODO: SZh 8.11.2022: add selection of parameterisation
+
+/// \brief Class L1SearchWindow defines a parameterisation of hits search window for CA tracking algorithm
+/// TODO: SZh 8.11.2022: add description
+class L1SearchWindow {
+public:
+ /// Constructor
+ /// \param stationID Global index of active station
+ /// \param trackGrID Track group ID
+ L1SearchWindow(int stationID, int trackGrID);
+
+ /// Default constructor
+ L1SearchWindow() = default;
+
+ /// Destructor
+ ~L1SearchWindow() = default;
+
+ /// Copy constructor
+ L1SearchWindow(const L1SearchWindow& other) = default;
+
+ /// Copy assignment operator
+ L1SearchWindow& operator=(const L1SearchWindow& other) = default;
+
+ /// Move constructor
+ L1SearchWindow(L1SearchWindow&& other) = default;
+
+ /// Move assignment operator
+ L1SearchWindow& operator=(L1SearchWindow&& other) = default;
+
+ /// Parameterisation function for dx_max(x0)
+ float DxMaxVsX0(float /*x*/) const { return fvParams[kDxMaxVsX0 /*+ 0*/]; }
+
+ /// Parameterisation function for dx_min(x0)
+ float DxMinVsX0(float /*x*/) const { return fvParams[kDxMinVsX0 /*+ 0*/]; }
+
+ /// Parameterisation function for dx_max(y0)
+ float DxMaxVsY0(float /*x*/) const { return fvParams[kDxMaxVsY0 /*+ 0*/]; }
+
+ /// Parameterisation function for dx_min(y0)
+ float DxMinVsY0(float /*x*/) const { return fvParams[kDxMinVsY0 /*+ 0*/]; }
+
+ /// Parameterisation function for dy_max(x0)
+ float DyMaxVsX0(float /*x*/) const { return fvParams[kDyMaxVsX0 /*+ 0*/]; }
+
+ /// Parameterisation function for dy_min(x0)
+ float DyMinVsX0(float /*x*/) const { return fvParams[kDyMinVsX0 /*+ 0*/]; }
+
+ /// Parameterisation function for dy_max(y0)
+ float DyMaxVsY0(float /*y*/) const { return fvParams[kDyMaxVsY0 /*+ 0*/]; }
+
+ /// Parameterisation function for dy_min(y0)
+ float DyMinVsY0(float /*y*/) const { return fvParams[kDyMinVsY0 /*+ 0*/]; }
+
+
+ /// Gets station id
+ int GetStationID() const { return fStationID; }
+
+ /// Gets track group id
+ int GetTrackGroupID() const { return fTrackGroupID; }
+
+
+ /// TODO: SZh 08.11.2022: Implement variadic template function
+ /// TODO: SZh 08.11.2022: Try to reduce number of functions
+ /// Sets parameters for dx_max(x0)
+ /// \param id Parameter index
+ /// \param val Parameter value
+ void SetParamDxMaxVsX0(int id, float val);
+
+ /// Sets parameters for dx_min(x0)
+ /// \param id Parameter index
+ /// \param val Parameter value
+ void SetParamDxMinVsX0(int id, float val);
+
+ /// Sets parameters for dx_max(y0)
+ /// \param id Parameter index
+ /// \param val Parameter value
+ void SetParamDxMaxVsY0(int id, float val);
+
+ /// Sets parameters for dx_min(y0)
+ /// \param id Parameter index
+ /// \param val Parameter value
+ void SetParamDxMinVsY0(int id, float val);
+
+ /// Sets parameters for dy_max(x0)
+ /// \param id Parameter index
+ /// \param val Parameter value
+ void SetParamDyMaxVsX0(int id, float val);
+
+ /// Sets parameters for dy_min(x0)
+ /// \param id Parameter index
+ /// \param val Parameter value
+ void SetParamDyMinVsX0(int id, float val);
+
+ /// Sets parameters for dy_max(y0)
+ /// \param id Parameter index
+ /// \param val Parameter value
+ void SetParamDyMaxVsY0(int id, float val);
+
+ /// Sets parameters for dy_min(y0)
+ /// \param id Parameter index
+ /// \param val Parameter value
+ void SetParamDyMinVsY0(int id, float val);
+
+ /// String representation of the contents
+ std::string ToString() const;
+
+private:
+ static constexpr unsigned char kNpars = 1; ///< Max number of parameters for one dependency
+ static constexpr unsigned char kNdeps = 8; ///< Number of the dependencies
+
+ /// Enumeration for dependencies stored
+ enum EDependency
+ {
+ kDxMaxVsX0,
+ kDxMinVsX0,
+ kDxMaxVsY0,
+ kDxMinVsY0,
+ kDyMaxVsX0,
+ kDyMinVsX0,
+ kDyMaxVsY0,
+ kDyMinVsY0
+ };
+
+ /// Search window parameter array containing parameters from
+ /// - dx_max vs x0 - indexes [0 .. kNpars - 1]
+ /// - dx_min vs x0 - indexes [kNpars .. (2 * kNpars - 1)]
+ /// - dx_max vs y0 - indexes [2 * kNpars .. (3 * kNpars - 1)]
+ /// - dx_min vs y0 - indexes [3 * kNpars .. (4 * kNpars - 1)]
+ /// - dy_max vs y0 - indexes [4 * kNpars .. (5 * kNpars - 1)]
+ /// - dy_min vs y0 - indexes [5 * kNpars .. (6 * kNpars - 1)]
+ /// - dy_max vs y0 - indexes [6 * kNpars .. (7 * kNpars - 1)]
+ /// - dy_min vs y0 - indexes [7 * kNpars .. (8 * kNpars - 1)]
+ std::array<float, kNdeps* kNpars> fvParams = {0};
+
+ 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?)
+
+ /// Serialization function
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int)
+ {
+ ar& fvParams;
+ ar& fStationID;
+ ar& fTrackGroupID;
+ }
+};
+
+#endif // L1SearchWindow_h
diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h
index 152819d5b8..7a1cccc9b4 100644
--- a/reco/L1/L1LinkDef.h
+++ b/reco/L1/L1LinkDef.h
@@ -30,5 +30,6 @@
#pragma link C++ class CbmTrackingInputQaSts + ;
#pragma link C++ class CbmTrackerInputQaTrd + ;
#pragma link C++ class CbmTrackerInputQaTof + ;
+#pragma link C++ class ca::tools::WindowFinder + ;
#endif
diff --git a/reco/L1/catools/CaToolsWindowFinder.cxx b/reco/L1/catools/CaToolsWindowFinder.cxx
new file mode 100644
index 0000000000..d663c7cc5b
--- /dev/null
+++ b/reco/L1/catools/CaToolsWindowFinder.cxx
@@ -0,0 +1,84 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file CaToolsWindowFinder.cxx
+/// \brief Framework for CA tracking hit-finder window estimation from MC (implementation)
+/// \since 14.10.2022
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#include "CaToolsWindowFinder.h"
+
+#include <Logger.h>
+
+#include "TChain.h"
+
+#include "L1Constants.h"
+
+using namespace ca::tools;
+using namespace L1Constants::clrs; // for colored logs
+
+ClassImp(ca::tools::WindowFinder);
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+WindowFinder::WindowFinder() : fpMcTripletsTree(new TChain(kTreeName))
+{
+ LOG(info) << kGNb << "Call" << kCLb << "WindowFinder constructor\n" << kCL;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::AddInputFile(const char* filename)
+{
+ assert(fpMcTripletsTree);
+ auto addingRes = fpMcTripletsTree->AddFile(filename);
+ if (addingRes != 1) {
+ LOG(error) << "WindowFinder::AddInputFile: File \"" << filename << "\" cannot be added to the MC triplets chain, "
+ << "some errors occurred\n";
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::Process(EBuildingMode) { fpMcTripletsTree->Draw("zv"); }
+
+
+// ****************************
+// ** Setters implementation **
+// ****************************
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::SetBinning(int nBinsX, int nBinsY)
+{
+ assert(nBinsX > 0);
+ assert(nBinsY > 0);
+ fNbinsX = nBinsX;
+ fNbinsY = nBinsY;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::SetEpsilon(float eps)
+{
+ assert(eps > 0. && eps <= 1.);
+ fEps = eps;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::SetNslices(int nSlices)
+{
+ assert(nSlices > 0);
+ fNslices = nSlices;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void WindowFinder::SetTarget(double x, double y, double z)
+{
+ fTargetPos[0] = x;
+ fTargetPos[1] = y;
+ fTargetPos[2] = z;
+}
diff --git a/reco/L1/catools/CaToolsWindowFinder.h b/reco/L1/catools/CaToolsWindowFinder.h
new file mode 100644
index 0000000000..1ea7e53314
--- /dev/null
+++ b/reco/L1/catools/CaToolsWindowFinder.h
@@ -0,0 +1,106 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file CaToolsWindowFinder.h
+/// \brief Framework for CA tracking hit-finder window estimation from MC (header)
+/// \since 14.10.2022
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#ifndef CaToolsWindowFinder_h
+#define CaToolsWindowFinder_h 1
+
+#include "TObject.h"
+
+// TODO: Replace tmp asserts with exceptions
+
+class TChain;
+
+namespace ca
+{
+ namespace tools
+ {
+
+ /// Enumeration of parameterisations used to describe a search window
+ enum class EWindowForm
+ {
+ kConstant,
+ kLinear,
+ kParabolic,
+ kElliptic
+ };
+
+ /// Enumeration
+ enum class EBuildingMode
+ {
+ kDouplets, ///< Windows are built from doublets and a target
+ kTriplets ///< Windows are built from triplets without using target
+ };
+
+ /// TODO: ... write an instruction ...
+ class WindowFinder : public TObject {
+ public:
+ // TODO: TEMPORARY CONSTANT EXPRESSIONS (TO BE MOVED TO A SEPARATE HEADER)
+ static constexpr const char* kTreeName = "t"; ///< Name of the input MC triplets tree
+
+ public:
+ /// Default constructor
+ WindowFinder();
+
+ virtual ~WindowFinder() = default;
+
+ /// Copy and move are forbidden
+ WindowFinder(const WindowFinder&) = delete;
+ WindowFinder(WindowFinder&&) = delete;
+ WindowFinder& operator=(const WindowFinder&) = delete;
+ WindowFinder& operator=(WindowFinder&&) = delete;
+
+ /// Adds an input file with a tree object of MC triplets
+ /// \note Multiple file input is possible
+ void AddInputFile(const char* filename);
+
+ /// Process a tree (TEST)
+ void Process(EBuildingMode buildingMode);
+
+ /// Sets binning of the dx (dy) distribution
+ /// \param nBinsX Number of bins for the x-axis
+ /// \param nBinsY Number of bins for the y-axis
+ void SetBinning(int nBinsX, int nBinsY);
+
+ /// 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
+ /// 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 components of the target center position
+ /// \param x Position x-component [cm]
+ /// \param y Position y-component [cm]
+ /// \param z Position z-component [cm]
+ void SetTarget(double x, double y, double z);
+
+
+ // ******************************
+ // ** Member variables **
+ // ******************************
+
+ private:
+ std::array<double, 3> fTargetPos = {0}; ///< Target position {x, y, z} [cm]
+
+ // Window extraction settings
+
+ 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
+ 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
+
+ ClassDef(WindowFinder, 0);
+ };
+ } // namespace tools
+} // namespace ca
+
+#endif // CaToolsWindowsFinder_h
--
GitLab