From e01908fe6f8b6b4bba469ec292f14d9a305d1b0f Mon Sep 17 00:00:00 2001
From: Valentina <v.akishina@gsi.de>
Date: Tue, 23 Nov 2021 12:16:05 +0100
Subject: [PATCH] L1: add broken triplets for mCBM
---
macro/mcbm/mcbm_reco_event.C | 118 ++++++++++++++---------------
reco/L1/CbmL1.cxx | 7 +-
reco/L1/CbmL1.h | 2 +
reco/L1/CbmL1ReadEvent.cxx | 17 ++---
reco/L1/L1Algo/L1Algo.cxx | 3 +-
reco/L1/L1Algo/L1Algo.h | 9 ++-
reco/L1/L1Algo/L1CATrackFinder.cxx | 39 +++++-----
7 files changed, 102 insertions(+), 93 deletions(-)
diff --git a/macro/mcbm/mcbm_reco_event.C b/macro/mcbm/mcbm_reco_event.C
index 440da89178..7fe3ef71de 100644
--- a/macro/mcbm/mcbm_reco_event.C
+++ b/macro/mcbm/mcbm_reco_event.C
@@ -359,73 +359,71 @@ void mcbm_reco_event(Int_t nEvents = 10, TString dataset = "data/test", const ch
// ------------------------------------------------------------------------
// -------- L1 CA Track Finder ---------------------------------------
- if (strcmp(setupName, "mcbm_beam_2020_03") == 0) {
-
- CbmKF* kalman = new CbmKF();
- run->AddTask(kalman);
- CbmL1* l1 = new CbmL1();
- l1->SetLegacyEventMode(1);
- l1->SetMcbmMode();
- l1->SetUseHitErrors(1);
-
- // --- Material budget file names
- TString mvdGeoTag;
- if (setup->GetGeoTag(ECbmModuleId::kMvd, mvdGeoTag)) {
- TString parFile = gSystem->Getenv("VMCWORKDIR");
- parFile = parFile + "/parameters/mvd/mvd_matbudget_" + mvdGeoTag + ".root";
- std::cout << "Using material budget file " << parFile << std::endl;
- l1->SetMvdMaterialBudgetFileName(parFile.Data());
- }
- TString stsGeoTag;
- if (setup->GetGeoTag(ECbmModuleId::kSts, stsGeoTag)) {
- TString parFile = gSystem->Getenv("VMCWORKDIR");
- parFile = parFile + "/parameters/sts/sts_matbudget_v19a.root";
- std::cout << "Using material budget file " << parFile << std::endl;
- l1->SetStsMaterialBudgetFileName(parFile.Data());
- }
-
- TString muchGeoTag;
- if (setup->GetGeoTag(ECbmModuleId::kMuch, muchGeoTag)) {
+ CbmKF* kalman = new CbmKF();
+ run->AddTask(kalman);
+ CbmL1* l1 = new CbmL1();
+ l1->SetLegacyEventMode(1);
+ l1->SetMcbmMode();
+ l1->SetUseHitErrors(1);
+ if (strcmp(setupName, "mcbm_beam_2021_07_surveyed") == 0) l1->SetMissingHits(1);
+
+ // --- Material budget file names
+ TString mvdGeoTag;
+ if (setup->GetGeoTag(ECbmModuleId::kMvd, mvdGeoTag)) {
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile = parFile + "/parameters/mvd/mvd_matbudget_" + mvdGeoTag + ".root";
+ std::cout << "Using material budget file " << parFile << std::endl;
+ l1->SetMvdMaterialBudgetFileName(parFile.Data());
+ }
+ TString stsGeoTag;
+ if (setup->GetGeoTag(ECbmModuleId::kSts, stsGeoTag)) {
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile = parFile + "/parameters/sts/sts_matbudget_v19a.root";
+ std::cout << "Using material budget file " << parFile << std::endl;
+ l1->SetStsMaterialBudgetFileName(parFile.Data());
+ }
- // --- Parameter file name
- TString geoTag;
- setup->GetGeoTag(ECbmModuleId::kMuch, geoTag);
- Int_t muchFlag = 0;
- if (geoTag.Contains("mcbm")) muchFlag = 1;
+ TString muchGeoTag;
+ if (setup->GetGeoTag(ECbmModuleId::kMuch, muchGeoTag)) {
- TString parFile = gSystem->Getenv("VMCWORKDIR");
- parFile = parFile + "/parameters/much/much_" + geoTag(0, 4) + "_digi_sector.root";
- std::cout << "L1: Using parameter file " << parFile << std::endl;
- l1->SetMuchPar(parFile);
+ // --- Parameter file name
+ TString geoTag;
+ setup->GetGeoTag(ECbmModuleId::kMuch, geoTag);
+ Int_t muchFlag = 0;
+ if (geoTag.Contains("mcbm")) muchFlag = 1;
- TString parFile2 = gSystem->Getenv("VMCWORKDIR");
- parFile2 = parFile2 + "/parameters/much/much_matbudget_" + geoTag + ".root ";
- std::cout << "Using material budget file " << parFile2 << std::endl;
- l1->SetMuchMaterialBudgetFileName(parFile2.Data());
- }
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile = parFile + "/parameters/much/much_" + geoTag(0, 4) + "_digi_sector.root";
+ std::cout << "L1: Using parameter file " << parFile << std::endl;
+ l1->SetMuchPar(parFile);
+
+ TString parFile2 = gSystem->Getenv("VMCWORKDIR");
+ parFile2 = parFile2 + "/parameters/much/much_matbudget_" + geoTag + ".root ";
+ std::cout << "Using material budget file " << parFile2 << std::endl;
+ l1->SetMuchMaterialBudgetFileName(parFile2.Data());
+ }
- TString trdGeoTag;
- if (setup->GetGeoTag(ECbmModuleId::kTrd, trdGeoTag)) {
- TString parFile = gSystem->Getenv("VMCWORKDIR");
- parFile = parFile + "/parameters/trd/trd_matbudget_" + trdGeoTag + ".root ";
- std::cout << "Using material budget file " << parFile << std::endl;
- l1->SetTrdMaterialBudgetFileName(parFile.Data());
- }
+ TString trdGeoTag;
+ if (setup->GetGeoTag(ECbmModuleId::kTrd, trdGeoTag)) {
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile = parFile + "/parameters/trd/trd_matbudget_" + trdGeoTag + ".root ";
+ std::cout << "Using material budget file " << parFile << std::endl;
+ l1->SetTrdMaterialBudgetFileName(parFile.Data());
+ }
- TString tofGeoTag;
- if (setup->GetGeoTag(ECbmModuleId::kTof, tofGeoTag)) {
- TString parFile = gSystem->Getenv("VMCWORKDIR");
- parFile = parFile + "/parameters/tof/tof_matbudget_" + tofGeoTag + ".root ";
- std::cout << "Using material budget file " << parFile << std::endl;
- l1->SetTofMaterialBudgetFileName(parFile.Data());
- }
+ TString tofGeoTag;
+ if (setup->GetGeoTag(ECbmModuleId::kTof, tofGeoTag)) {
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile = parFile + "/parameters/tof/tof_matbudget_" + tofGeoTag + ".root ";
+ std::cout << "Using material budget file " << parFile << std::endl;
+ l1->SetTofMaterialBudgetFileName(parFile.Data());
+ }
- run->AddTask(l1);
+ run->AddTask(l1);
- CbmL1GlobalTrackFinder* globalTrackFinder = new CbmL1GlobalTrackFinder();
- FairTask* globalFindTracks = new CbmL1GlobalFindTracksEvents(globalTrackFinder);
- run->AddTask(globalFindTracks);
- }
+ CbmL1GlobalTrackFinder* globalTrackFinder = new CbmL1GlobalTrackFinder();
+ FairTask* globalFindTracks = new CbmL1GlobalFindTracksEvents(globalTrackFinder);
+ run->AddTask(globalFindTracks);
// ----- Parameter database --------------------------------------------
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 6491893064..3ac2a458b4 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -744,7 +744,7 @@ InitStatus CbmL1::Init()
LOG(error) << "-E- CbmL1: Read geometry from file " << fSTAPDataDir + "geo_algo.txt was NOT successful.";
}
- algo->Init(geo, fUseHitErrors, fTrackingMode);
+ algo->Init(geo, fUseHitErrors, fTrackingMode, fMissingHits);
geo.clear();
algo->fRadThick.reset(algo->NStations);
@@ -919,7 +919,10 @@ InitStatus CbmL1::Init()
for (int j = 1, iSta = (NStsStations + NMvdStations + NMuchStations);
iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations); iSta++, j++) {
TString stationNameSts = stationName;
- stationNameSts += j;
+ int skipStation = j;
+ if (skipStation == 2) skipStation = 3;
+ if (skipStation == 3) skipStation = 4;
+ stationNameSts += skipStation;
TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts);
if (!hStaRadLen) {
cout << "L1: incorrect " << fTrdMatBudgetFileName << " file. No " << stationNameSts << "\n";
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index c7b40e570d..053f37eb32 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -155,6 +155,7 @@ public:
void SetLegacyEventMode(bool b) { fLegacyEventMode = b; }
void SetMuchPar(TString fileName) { fMuchDigiFile = fileName; }
void SetUseHitErrors(bool value) { fUseHitErrors = value; }
+ void SetMissingHits(bool value) { fMissingHits = value; }
void SetStsOnlyMode() { fTrackingMode = L1Algo::TrackingMode::kSts; }
void SetMcbmMode() { fTrackingMode = L1Algo::TrackingMode::kMcbm; }
void SetGlobalMode() { fTrackingMode = L1Algo::TrackingMode::kGlobal; }
@@ -232,6 +233,7 @@ public:
TString fMuchDigiFile {}; // Much digitization file name
bool fUseHitErrors {false};
+ bool fMissingHits {false};
L1Algo::TrackingMode fTrackingMode {L1Algo::TrackingMode::kSts};
L1Vector<CbmL1Track> vRTracks {"CbmL1::vRTracks"}; // reconstructed tracks
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index e5ad92f54d..62dddcfa32 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -685,10 +685,10 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.iStation = NMvdStations + sta + NStsStations + NMuchStations;
- // if (mh->GetPlaneId()==0) continue;
- // if (mh->GetPlaneId()==1) continue;
- // if (mh->GetPlaneId()==2) continue;
- // if (mh->GetPlaneId()==3) continue;
+ // if (mh->GetPlaneId()==0) continue;
+ // if (mh->GetPlaneId()==1) continue;
+ // if (mh->GetPlaneId()==2) continue;
+ // if (mh->GetPlaneId()==3) continue;
th.time = mh->GetTime();
@@ -797,13 +797,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
// if (fTofDigiBdfPar->GetTrackingStation(mh) == 4) sttof = 2;
// if (fTofDigiBdfPar->GetTrackingStation(mh) == 5) sttof = 2;
- if ((th.x > 20) && (th.z > 270) && (fTofDigiBdfPar->GetTrackingStation(mh) == 1)) sttof = 2;
-
- th.iStation = sttof + NMvdStations + NStsStations + NMuchStations + NTrdStations;
th.time = mh->GetTime();
-
- th.dt = mh->GetTimeError();
+ th.dt = mh->GetTimeError();
th.dx = mh->GetDx();
th.dy = mh->GetDy();
@@ -825,8 +821,11 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.y = pos.Y();
th.z = pos.Z();
+ if ((th.x > 20) && (th.z > 270) && (fTofDigiBdfPar->GetTrackingStation(mh) == 1)) sttof = 2;
if (th.z > 400) continue;
+ th.iStation = sttof + NMvdStations + NStsStations + NMuchStations + NTrdStations;
+
L1Station& st = algo->vStations[th.iStation];
th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 6328897c8c..8486c58e5d 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -59,7 +59,7 @@ void L1Algo::SetNThreads(unsigned int n)
}
-void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const TrackingMode mode)
+void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const TrackingMode mode, const bool MissingHits)
{
for (int iProc = 0; iProc < 4; iProc++) {
@@ -75,6 +75,7 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
fUseHitErrors = UseHitErrors;
fTrackingMode = mode;
+ fMissingHits = MissingHits;
//lxir039
// for (int i=0; i<8; i++){
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 92506cbb76..7cd96beee2 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -184,7 +184,7 @@ public:
kMcbm
};
- void Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const TrackingMode mode);
+ void Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const TrackingMode mode, const bool MissingHits);
void SetData(L1Vector<L1Hit>& StsHits_, int nStsStrips_, L1Vector<unsigned char>& SFlag_,
const THitI* StsHitsStartIndex_, const THitI* StsHitsStopIndex_);
@@ -207,9 +207,9 @@ public:
void SetNThreads(unsigned int n);
- int NStations {0}; // number of all detector stations
- int NMvdStations {0}; // number of mvd stations
- int NStsStations {0}; // number of sts stations
+ int NStations {0}; // number of all detector stations
+ int NMvdStations {0}; // number of mvd stations
+ int NStsStations {0}; // number of sts stations
int fNfieldStations {0}; // number of stations in the field region
L1Station vStations[fkMaxNstations] _fvecalignment; // station info
@@ -257,6 +257,7 @@ public:
int fNThreads {0};
bool fUseHitErrors {0};
+ bool fMissingHits {0};
TrackingMode fTrackingMode {kSts};
fvec EventTime[fkMaxNthreads][fkMaxNthreads] {{0}};
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index df8c15fdef..f785fccb0a 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -1768,7 +1768,7 @@ void L1Algo::CATrackFinder()
for (isec = 0; isec < fNFindIterations; ++isec) // all finder
{
if (fTrackingMode == kMcbm) {
- if (isec > 1) { continue; }
+ if (isec > 3) { continue; }
}
// n_g1.assign(n_g1.size(), Portion);
@@ -1848,10 +1848,13 @@ void L1Algo::CATrackFinder()
MaxInvMom = 1.0 / 0.5; // max considered q/p
- if (fTrackingMode == kMcbm) MaxInvMom = 1.5 / 0.1; // max considered q/p
+ if (fTrackingMode == kMcbm) MaxInvMom = 1 / 0.3; // max considered q/p
if ((isec == kAllPrimJumpIter) || (isec == kAllSecIter) || (isec == kAllSecJumpIter)) MaxInvMom = 1.0 / 0.1;
if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllSecEIter)) MaxInvMom = 1. / 0.05;
+ if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllSecEIter))
+ if (fTrackingMode == kMcbm) MaxInvMom = 1 / 0.1; // max considered q/p
+
MaxSlopePV = 1.1;
if ( // (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) ||
(isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
@@ -2030,32 +2033,34 @@ void L1Algo::CATrackFinder()
// output
);
-
- if ((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecJumpIter)) {
+ if ((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecJumpIter) || (fMissingHits)) {
Tindex nG_2;
hitsmG_2.clear();
i1G_2.clear();
- DupletsStaPort( // input
- istal, istal + 2, ip, fDupletPortionSize, fDupletPortionStopIndex,
- // output
- TG_1, fldG_1, hitslG_1,
+ if ((fMissingHits && ((istal == 0) || (istal == 1))) || !fMissingHits)
+ DupletsStaPort( // input
+ istal, istal + 2, ip, fDupletPortionSize, fDupletPortionStopIndex,
+ // output
+ TG_1, fldG_1, hitslG_1,
- lmDupletsG[istal],
+ lmDupletsG[istal],
- nG_2, i1G_2, hitsmG_2);
+ nG_2, i1G_2, hitsmG_2);
- TripletsStaPort( // input
- istal, istal + 1, istal + 3, nstaltriplets, T_1, fld_1, hitsl_1,
+ if ((fMissingHits && (istal == 0)) || !fMissingHits)
+ TripletsStaPort( // input
+ istal, istal + 1, istal + 3, nstaltriplets, T_1, fld_1, hitsl_1,
- n_2, i1_2, hitsm_2, lmDupletsG[istal + 1]);
+ n_2, i1_2, hitsm_2, lmDupletsG[istal + 1]);
- TripletsStaPort( // input
- istal, istal + 2, istal + 3, nstaltriplets, TG_1, fldG_1, hitslG_1,
+ if ((fMissingHits && (istal == 1)) || !fMissingHits)
+ TripletsStaPort( // input
+ istal, istal + 2, istal + 3, nstaltriplets, TG_1, fldG_1, hitslG_1,
- nG_2, i1G_2, hitsmG_2, lmDuplets[istal + 2]
+ nG_2, i1G_2, hitsmG_2, lmDuplets[istal + 2]
- );
+ );
}
} //
}
--
GitLab