From e667cbbab2f13daec69b5af1e6a4920a8d8fef84 Mon Sep 17 00:00:00 2001
From: Valentina Akishina <v.akishina@gsi.de>
Date: Mon, 16 May 2022 18:02:18 +0000
Subject: [PATCH] L1: skip stations before the target
---
reco/L1/CbmL1.cxx | 143 ++++++++++++++++++++---
reco/L1/CbmL1.h | 9 ++
reco/L1/CbmL1ReadEvent.cxx | 39 +++++--
reco/L1/L1Algo/L1Algo.cxx | 4 +-
reco/L1/ParticleFinder/CbmL1PFFitter.cxx | 14 ++-
5 files changed, 176 insertions(+), 33 deletions(-)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 440f4f8ce6..f99495a506 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -498,9 +498,6 @@ InitStatus CbmL1::Init()
NStsStations = CbmStsSetup::Instance()->GetNofStations();
NStation = NMvdStations + NStsStations + NMuchStations + NTrdStations + NTOFStation;
- geo.push_back(NStation);
- geo.push_back(NMvdStations);
- geo.push_back(NStsStations);
// field
const int M = 5; // polinom order
@@ -559,7 +556,124 @@ InitStatus CbmL1::Init()
//
// } // target field
+ int NMvdStationsTemp = NMvdStations;
+ int NStsStationsTemp = NStsStations;
+ int NTrdStationsTemp = NTrdStations;
+ int NMuchStationsTemp = NMuchStations;
+ int NTOFStationTemp = NTOFStation;
+ int NStationTemp = NStation;
+
+ vector<int> StationIdxConverter;
+
+ auto& target = CbmKF::Instance()->vTargets[0];
+
+ for (Int_t ist = 0; ist < NStation; ist++) {
+ if (ist < NMvdStations) {
+ CbmMvdDetector* mvdDetector = CbmMvdDetector::Instance();
+ CbmMvdStationPar* mvdStationPar = mvdDetector->GetParameterFile();
+ CbmKFTube& t = CbmKF::Instance()->vMvdMaterial[ist];
+
+ if (target.z < t.z) {
+ StationIdxConverter.push_back(ist);
+ StationIdxReverseConverter.push_back(StationIdxConverter.size() - 1);
+ }
+ else {
+ NMvdStationsTemp--;
+ NStationTemp--;
+ StationIdxReverseConverter.push_back(-1);
+ }
+ }
+
+ if (ist >= NMvdStations && ist < (NMvdStations + NStsStations)) {
+ CbmStsStation* station = CbmStsSetup::Instance()->GetStation(ist - NMvdStations);
+ if (target.z < station->GetZ()) {
+ StationIdxConverter.push_back(ist);
+ StationIdxReverseConverter.push_back(StationIdxConverter.size() - 1);
+ }
+ else {
+ NStsStationsTemp--;
+ NStationTemp--;
+ StationIdxReverseConverter.push_back(-1);
+ }
+ }
+
+ if ((ist < (NMvdStations + NStsStations + NMuchStations)) && (ist >= (NMvdStations + NStsStations))) {
+
+ int iStation = (ist - NMvdStations - NStsStations) / 3;
+ CbmMuchStation* station = (CbmMuchStation*) fGeoScheme->GetStation(iStation);
+ CbmMuchLayer* layer = station->GetLayer((ist - NMvdStations - NStsStations) % 3);
+ if (target.z < layer->GetZ()) {
+ StationIdxConverter.push_back(ist);
+ StationIdxReverseConverter.push_back(StationIdxConverter.size() - 1);
+ }
+ else {
+ NTrdStationsTemp--;
+ NStationTemp--;
+ StationIdxReverseConverter.push_back(-1);
+ }
+ }
+
+ if ((ist < (NMvdStations + NStsStations + NTrdStations + NMuchStations))
+ && (ist >= (NMvdStations + NStsStations + NMuchStations))) {
+
+ int num = ist - NMvdStations - NStsStations - NMuchStations;
+ int skip = num;
+
+ if ((fTrackingMode == L1Algo::TrackingMode::kMcbm) && (fMissingHits))
+ if (num > 0) skip++;
+ int ModuleId = fTrdDigiPar->GetModuleId(skip);
+
+ CbmTrdParModDigi* module = (CbmTrdParModDigi*) fTrdDigiPar->GetModulePar(ModuleId);
+
+ if (target.z < module->GetZ()) {
+ StationIdxConverter.push_back(ist);
+ StationIdxReverseConverter.push_back(StationIdxConverter.size() - 1);
+ }
+ else {
+ NMuchStationsTemp--;
+ NStationTemp--;
+ StationIdxReverseConverter.push_back(-1);
+ }
+ }
+
+ if ((ist < (NMvdStations + NStsStations + NTrdStations + NMuchStations + NTOFStation))
+ && (ist >= (NMvdStations + NStsStations + NMuchStations + NTrdStations))) {
+
+ int num = ist - (NMvdStations + NStsStations + NTrdStations + NMuchStations);
+
+ if (target.z < TofStationZ[num]) {
+ StationIdxConverter.push_back(ist);
+ StationIdxReverseConverter.push_back(StationIdxConverter.size() - 1);
+ }
+ else {
+ NTOFStationTemp--;
+ NStationTemp--;
+ StationIdxReverseConverter.push_back(-1);
+ }
+ }
+ }
+
+ NMvdStationsGeom = NMvdStations;
+ NStsStationsGeom = NStsStations;
+ NTrdStationsGeom = NTrdStations;
+ NMuchStationsGeom = NMuchStations;
+ NTOFStationGeom = NTOFStation;
+ NStationGeom = NStation;
+
+ NMvdStations = NMvdStationsTemp;
+ NStsStations = NStsStationsTemp;
+ NTrdStations = NTrdStationsTemp;
+ NMuchStations = NMuchStationsTemp;
+ NTOFStation = NTOFStationTemp;
+ NStation = NStationTemp;
+
+ geo.push_back(NStation);
+ geo.push_back(NMvdStations);
+ geo.push_back(NStsStations);
+
+
for (Int_t ist = 0; ist < NStation; ist++) {
+
double z = 0;
double Xmax = 0, Ymax = 0;
if (ist < NMvdStations) {
@@ -567,7 +681,7 @@ InitStatus CbmL1::Init()
CbmMvdDetector* mvdDetector = CbmMvdDetector::Instance();
CbmMvdStationPar* mvdStationPar = mvdDetector->GetParameterFile();
- CbmKFTube& t = CbmKF::Instance()->vMvdMaterial[ist];
+ CbmKFTube& t = CbmKF::Instance()->vMvdMaterial[StationIdxConverter[ist]];
geo.push_back(1);
geo.push_back(t.z);
geo.push_back(t.dz);
@@ -590,7 +704,8 @@ InitStatus CbmL1::Init()
if (ist >= NMvdStations && ist < (NMvdStations + NStsStations)) {
- CbmStsStation* station = CbmStsSetup::Instance()->GetStation(ist - NMvdStations);
+
+ CbmStsStation* station = CbmStsSetup::Instance()->GetStation(StationIdxConverter[ist] - NMvdStationsGeom);
geo.push_back(0);
geo.push_back(station->GetZ());
geo.push_back(station->GetSensorD());
@@ -626,12 +741,12 @@ InitStatus CbmL1::Init()
if ((ist < (NMvdStations + NStsStations + NMuchStations)) && (ist >= (NMvdStations + NStsStations))) {
- int iStation = (ist - NMvdStations - NStsStations) / 3;
+ int iStation = (StationIdxConverter[ist] - NMvdStationsGeom - NStsStationsGeom) / 3;
CbmMuchStation* station = (CbmMuchStation*) fGeoScheme->GetStation(iStation);
- CbmMuchLayer* layer = station->GetLayer((ist - NMvdStations - NStsStations) % 3);
+ CbmMuchLayer* layer = station->GetLayer((StationIdxConverter[ist] - NMvdStationsGeom - NStsStationsGeom) % 3);
z = layer->GetZ();
@@ -662,12 +777,7 @@ InitStatus CbmL1::Init()
if ((ist < (NMvdStations + NStsStations + NTrdStations + NMuchStations))
&& (ist >= (NMvdStations + NStsStations + NMuchStations))) {
- int num = ist - NMvdStations - NStsStations - NMuchStations;
-
- // if (num == 0) continue;//true_station = 0;
- // if (!true_station) continue;
-
- // Int_t nrModules = fTrdDigiPar->GetNrOfModules();
+ int num = StationIdxConverter[ist] - NMvdStationsGeom - NStsStationsGeom - NMuchStationsGeom;
int skip = num;
@@ -703,12 +813,13 @@ InitStatus CbmL1::Init()
LOG(info) << "L1: Trd station " << num << " at z " << stationZ << endl;
}
+
if ((ist < (NMvdStations + NStsStations + NTrdStations + NMuchStations + NTOFStation))
&& (ist >= (NMvdStations + NStsStations + NMuchStations + NTrdStations))) {
geo.push_back(4);
- int num = ist - (NMvdStations + NStsStations + NTrdStations + NMuchStations);
+ int num = StationIdxConverter[ist] - (NMvdStationsGeom + NStsStationsGeom + NTrdStationsGeom + NMuchStationsGeom);
geo.push_back(TofStationZ[num]);
@@ -728,6 +839,7 @@ InitStatus CbmL1::Init()
LOG(info) << "L1: Tof station " << num << " at z " << z << endl;
}
+
double dx = 1.; // step for the field approximation
double dy = 1.;
@@ -746,6 +858,7 @@ InitStatus CbmL1::Init()
b0(i) = b1(i) = b2(i) = 0.;
}
+
if (CbmKF::Instance()->GetMagneticField()) {
for (double x = -Xmax; x <= Xmax; x += dx)
for (double y = -Ymax; y <= Ymax; y += dy) {
@@ -784,6 +897,7 @@ InitStatus CbmL1::Init()
C[2][i] = c2(i);
}
}
+
geo.push_back(N);
for (int k = 0; k < 3; k++) {
for (int j = 0; j < N; j++) {
@@ -817,7 +931,6 @@ InitStatus CbmL1::Init()
//algo->SetParameters(fParameters);
-
#ifdef FEATURING_L1ALGO_INIT
/********************************************************************************************************************
* EXPERIMENTAL FEATURE: usage of L1InitManager for L1Algo initialization *
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 23e2d41bcf..5cfef845d7 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -275,6 +275,15 @@ private:
int NMuchStations {0}; ///< number of much stations
int NTrdStations {0}; ///< number of trd stations
int NTOFStation {0}; ///< number of tof stations
+ ///
+ int NMvdStationsGeom;
+ int NStsStationsGeom;
+ int NTrdStationsGeom;
+ int NMuchStationsGeom;
+ int NTOFStationGeom;
+ int NStationGeom;
+
+ vector<int> StationIdxReverseConverter;
Int_t fPerformance {0}; // 0 - w\o perf. 1 - L1-Efficiency definition. 2 - QA-Eff.definition
double fTrackingTime {0.}; // time of track finding
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index c5ab9e67de..5c8dcd469f 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -223,7 +223,6 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
Int_t iFile = set_it->first;
Int_t iEvent = set_it->second;
-
if (fMvdPoints) {
Int_t nMvdPointsInEvent = fMvdPoints->Size(iFile, iEvent);
double maxDeviation = 0;
@@ -463,8 +462,12 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.ExtIndex = -(1 + j);
th.id = tmpHits.size();
th.iStation = mh->GetStationNr();
- th.iStripF = firstDetStrip + j;
- th.iStripB = th.iStripF;
+
+ int stIdx = StationIdxReverseConverter[mh->GetStationNr()];
+ if (stIdx == -1) continue;
+ th.iStation = stIdx; //mh->GetStationNr() - 1;
+ th.iStripF = firstDetStrip + j;
+ th.iStripB = th.iStripF;
if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
TVector3 pos, err;
@@ -553,10 +556,15 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
CbmStsHit* mh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndexSort));
th.ExtIndex = hitIndexSort;
th.Det = 1;
- th.iStation =
- NMvdStations + CbmStsSetup::Instance()->GetStationNumber(mh->GetAddress()); //mh->GetStationNr() - 1;
- th.iStripF = firstDetStrip + mh->GetFrontClusterId();
- th.iStripB = firstDetStrip + mh->GetBackClusterId();
+ int stIdx =
+ StationIdxReverseConverter[CbmStsSetup::Instance()->GetStationNumber(mh->GetAddress()) + NMvdStationsGeom];
+
+ if (stIdx == -1) continue;
+
+
+ th.iStation = stIdx; //mh->GetStationNr() - 1;
+ th.iStripF = firstDetStrip + mh->GetFrontClusterId();
+ th.iStripB = firstDetStrip + mh->GetBackClusterId();
if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
if (NStrips <= th.iStripB) { NStrips = th.iStripB + 1; }
@@ -702,8 +710,10 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
int DetId = stationNumber * 3 + layerNumber;
+ int stIdx = StationIdxReverseConverter[DetId + NMvdStationsGeom + NStsStationsGeom];
+ if (stIdx == -1) continue;
+ th.iStation = stIdx; //mh->GetStationNr() - 1;
- th.iStation = DetId + NMvdStations + NStsStations;
//Get time
th.time = mh->GetTime() - 14.5;
th.dt = mh->GetTimeError();
@@ -798,9 +808,13 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
int sta = mh->GetPlaneId();
+ int stIdx =
+ StationIdxReverseConverter[mh->GetPlaneId() + NMvdStationsGeom + NStsStationsGeom + NMuchStationsGeom];
+ if (stIdx == -1) continue;
+
if ((fTrackingMode == L1Algo::TrackingMode::kMcbm) && (sta > 1) && (fMissingHits)) { sta = sta - 1; }
- th.iStation = NMvdStations + sta + NStsStations + NMuchStations;
+ th.iStation = stIdx;
// if (mh->GetPlaneId()==0) continue;
// if (mh->GetPlaneId()==1) continue;
@@ -979,7 +993,10 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
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;
+ int stIdx =
+ StationIdxReverseConverter[sttof + NMvdStationsGeom + NStsStationsGeom + NMuchStationsGeom + NTrdStationsGeom];
+ if (stIdx == -1) continue;
+ th.iStation = stIdx;
L1Station& st = algo->vStations[th.iStation];
th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
@@ -1427,6 +1444,7 @@ void CbmL1::HitMatch()
Int_t iEvent = link.GetEntry();
Int_t iIndex = link.GetIndex();
+
if (fLegacyEventMode) {
iFile = vFileEvent.begin()->first;
iEvent = vFileEvent.begin()->second;
@@ -1435,6 +1453,7 @@ void CbmL1::HitMatch()
Double_t dpnt = dFEI(iFile, iEvent, nMvdPoints + iIndex);
DFEI2I::iterator pnt_it = dFEI2vMCPoints.find(dpnt);
+
assert(pnt_it != dFEI2vMCPoints.end());
totalWeight += link.GetWeight();
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 81b39ed6ad..9f29749eb0 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -180,8 +180,6 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
//std::cout << "cos_b "<<c_b<< " sin_b " << s_b << std::endl;
//std::cout << "cos_f "<<c_f<< " sin_f " << s_f << std::endl;
-
-
int N = static_cast<int>(geo[ind++]);
for (int iC = 0; iC < N; iC++)
st.fieldSlice.cx[iC] = geo[ind++];
@@ -211,7 +209,7 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
Sy += dz * sy + dz * dz * By / 2.;
sy += dz * By;
//st.Sy = Sy; // commented, because is not used in the code (S.Zharko)
- z0 = st.z;
+ z0 = st.z;
}
}
// for( int iS = 0; iS < NStations; ++iS ) { /// Grid is created for each station with the same step: xStep,yStep
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index fac70237a7..d3f080747c 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -186,7 +186,9 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
posx = hit->GetX();
posy = hit->GetY();
posz = hit->GetZ();
- ista = hit->GetStationNr();
+ // ista = hit->GetStationNr();
+ ista = CbmL1::Instance()->StationIdxReverseConverter[hit->GetStationNr()];
+ if (ista == -1) continue;
}
else {
if (!listStsHits) continue;
@@ -196,8 +198,12 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
posx = hit->GetX();
posy = hit->GetY();
posz = hit->GetZ();
- ista = CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress())
- + NMvdStations; //hit->GetStationNr() - 1 + NMvdStations;
+ // ista = CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress())
+ // + NMvdStations; //hit->GetStationNr() - 1 + NMvdStations;
+ ista =
+ CbmL1::Instance()->StationIdxReverseConverter[CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress())
+ + CbmL1::Instance()->NMvdStationsGeom];
+ if (ista == -1) continue;
}
w[ista][iVec] = 1.f;
@@ -409,7 +415,6 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<L1FieldRe
L1FieldRegion fld _fvecalignment;
L1FieldValue fB[3], fB_temp _fvecalignment;
fvec zField[3];
-
FairRootManager* fManger = FairRootManager::Instance();
TClonesArray* listStsHits = (TClonesArray*) fManger->GetObject("StsHit");
TClonesArray* listMvdHits = 0;
@@ -519,7 +524,6 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<L1FieldRe
fvec chi = sqrt(fabs(0.5 * (dx * dx * c[0] - 2 * dx * dy * c[1] + dy * dy * c[2]) / d));
fvec isNull = fvec(fabs(d) < 1.e-20);
chi = fvec(fvec(!isNull) & chi) + fvec(isNull & fvec(0));
-
for (int iVec = 0; iVec < nTracks_SIMD; iVec++)
chiToVtx.push_back(chi[iVec]);
--
GitLab