From 1d1a6906a434978ee542e962bdf5a69f9993addb Mon Sep 17 00:00:00 2001 From: "s.zharko@gsi.de" <s.zharko@gsi.de> Date: Thu, 12 May 2022 16:23:11 +0200 Subject: [PATCH] L1Algo: geo vector was removed from init --- reco/L1/CbmL1.cxx | 14 +- reco/L1/CbmL1.h | 7 +- reco/L1/CbmL1Performance.cxx | 14 +- reco/L1/L1Algo/L1Algo.cxx | 315 +++++++++++------------ reco/L1/L1Algo/L1Algo.h | 43 ++-- reco/L1/L1Algo/L1CATrackFinder.cxx | 22 +- reco/L1/L1Algo/L1TrackExtender.cxx | 4 +- reco/L1/L1Algo/L1TrackFitter.cxx | 12 +- reco/L1/ParticleFinder/CbmL1PFFitter.cxx | 8 +- 9 files changed, 217 insertions(+), 222 deletions(-) diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index ec33fcef9d..36dd5b6c96 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -954,7 +954,10 @@ InitStatus CbmL1::Init() // Step 4: initialize IDs of detectors active in tracking // TODO: temporary for tests, must be initialized somewhere in run_reco.C or similar (S.Zh.) - fActiveTrackingDetectorIDs = {L1DetectorID::kMvd, L1DetectorID::kSts}; + if (fUseMVD) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kMvd); } + if (fUseMUCH) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kMuch); } + if (fUseTRD) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kTrd); } + if (fUseTOF) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kTof); } fpInitManager->SetActiveDetectorIDs(fActiveTrackingDetectorIDs); constexpr double PI = 3.14159265358; // TODO: why cmath is not used? (S.Zh.) @@ -1265,7 +1268,7 @@ InitStatus CbmL1::Init() TDirectory* oldDir = gDirectory; TFile* rlFile = new TFile(fMvdMatBudgetFileName); cout << "MVD Material budget file is " << fMvdMatBudgetFileName << ".\n"; - for (int j = 0, iSta = 0; iSta < algo->NMvdStations; iSta++, j++) { + for (int j = 0, iSta = 0; iSta < algo->GetNstationsBeforePipe(); iSta++, j++) { TString stationNameMvd = stationName; stationNameMvd += j; TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameMvd); @@ -1305,7 +1308,7 @@ InitStatus CbmL1::Init() else { LOG(warn) << "No MVD material budget file is found. Homogenious budget " "will be used"; - for (int iSta = 0; iSta < algo->NMvdStations; iSta++) { + for (int iSta = 0; iSta < algo->GetNstationsBeforePipe(); iSta++) { algo->fRadThick[iSta].SetBins(1, 100); // mvd should be in +-100 cm square algo->fRadThick[iSta].table.resize(1); @@ -1320,7 +1323,7 @@ InitStatus CbmL1::Init() TDirectory* oldDir = gDirectory; TFile* rlFile = new TFile(fStsMatBudgetFileName); cout << "STS Material budget file is " << fStsMatBudgetFileName << ".\n"; - for (int j = 1, iSta = algo->NMvdStations; iSta < (algo->NMvdStations + NStsStations); iSta++, j++) { + for (int j = 1, iSta = algo->GetNstationsBeforePipe(); iSta < (algo->GetNstationsBeforePipe() + NStsStations); iSta++, j++) { TString stationNameSts = stationName; stationNameSts += j; TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts); @@ -1355,7 +1358,7 @@ InitStatus CbmL1::Init() else { LOG(warn) << "No STS material budget file is found. Homogenious budget " "will be used"; - for (int iSta = algo->NMvdStations; iSta < (algo->NMvdStations + NStsStations); iSta++) { + for (int iSta = algo->GetNstationsBeforePipe(); iSta < (algo->GetNstationsBeforePipe() + NStsStations); iSta++) { algo->fRadThick[iSta].SetBins(1, 100); algo->fRadThick[iSta].table.resize(1); algo->fRadThick[iSta].table[0].resize(1); @@ -1556,7 +1559,6 @@ void CbmL1::Exec(Option_t* /*option*/) {} void CbmL1::Reconstruct(CbmEvent* event) { - fVerbose = 11; // TODO: Remove it (tmp)! (S.Zharko) static int nevent = 0; vFileEvent.clear(); diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h index 745c94aa85..da09b0569a 100644 --- a/reco/L1/CbmL1.h +++ b/reco/L1/CbmL1.h @@ -94,8 +94,8 @@ enum class L1DetectorID kMvd, kSts, kMuch, - kTof, - kTrd + kTrd, + kTof }; // TODO: insert documentation! (S.Zh.) @@ -294,8 +294,7 @@ private: L1InitManager* fpInitManager {nullptr}; ///< Pointer to L1InitManager object of L1 algorithm core - std::set<L1DetectorID> fActiveTrackingDetectorIDs {L1DetectorID::kMvd, - L1DetectorID::kSts}; ///< Set of detectors active in tracking + std::set<L1DetectorID> fActiveTrackingDetectorIDs {L1DetectorID::kSts}; ///< Set of detectors active in tracking L1AlgoInputData* fData {nullptr}; diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx index 6c602bfa69..bdf485e0ce 100644 --- a/reco/L1/CbmL1Performance.cxx +++ b/reco/L1/CbmL1Performance.cxx @@ -221,7 +221,7 @@ struct TL1PerfEfficiencies : public TL1Efficiencies { mc_length_hits.counters[index] += _mc_length_hits; }; - void PrintEff(const std::string& nameOfTable = "efficiency_table") + void PrintEff(bool ifPrintTableToLog = false, const std::string& nameOfTable = "efficiency_table") { L1_assert(nEvents != 0); @@ -272,7 +272,9 @@ struct TL1PerfEfficiencies : public TL1Efficiencies { aTable->SetCell(iC, 7, mc_length_hits.counters[iC] / double(mc.counters[iC])); aTable->SetCell(iC, 8, mc_length.counters[iC] / double(mc.counters[iC])); } - cout << *aTable; // print a table to log + if (ifPrintTableToLog) { + cout << *aTable; // print a table to log + } cout << "Ghost probability : " << ratio_ghosts << " | " << ghosts << endl; gDirectory = curdir; @@ -464,7 +466,7 @@ void CbmL1::EfficienciesPerformance() } // fVerbose > 1 cout << endl; cout << "L1 ACCUMULATED STAT : " << L1_NEVENTS << " EVENTS " << endl << endl; - L1_NTRA.PrintEff(); + L1_NTRA.PrintEff(/*ifPrintTableToLog = */ true); cout << "MC tracks/event found : " << int(double(L1_NTRA.reco.counters[L1_NTRA.indices["total"]]) / double(L1_NEVENTS)) << endl; cout << endl; @@ -1200,7 +1202,7 @@ void CbmL1::TrackFitPerformance() } if (ih < 3) continue; CbmL1MCPoint& mcP = vMCPoints[mc.Points[0]]; - targB = algo->vtxFieldValue; + targB = algo->GetVtxFieldValue(); fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]); L1Extrapolate(trPar, mcP.zIn, trPar.qp, fld); @@ -1311,7 +1313,7 @@ void CbmL1::TrackFitPerformance() } if (ih < 3) continue; CbmL1MCPoint& mcP = vMCPoints[iMC]; - targB = algo->vtxFieldValue; + targB = algo->GetVtxFieldValue(); fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]); L1Extrapolate(trPar, mcP.zOut, trPar.qp, fld); @@ -1431,7 +1433,7 @@ void CbmL1::TrackFitPerformance() L1FieldValue B[3], targB _fvecalignment; float z[3]; - targB = algo->vtxFieldValue; + targB = algo->GetVtxFieldValue(); int ih = 1; for (unsigned int iHit = 0; iHit < it->StsHits.size(); iHit++) { diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx index f0251aaeee..b93b904e3b 100644 --- a/reco/L1/L1Algo/L1Algo.cxx +++ b/reco/L1/L1Algo/L1Algo.cxx @@ -82,137 +82,142 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra // threadNumberToCpuMap[2*i+1] = 15-(i+8); // } - int ind = 0; - { - L1FieldValue B[3]; - fvec z[3]; - for (int i = 0; i < 3; i++) { - z[i] = geo[ind++]; - B[i].x = geo[ind++]; - B[i].y = geo[ind++]; - B[i].z = geo[ind++]; -#ifndef TBB2 - std::cout << "L1Algo Input Magnetic field:" << z[i][0] << " " << B[i].x[0] << " " << B[i].y[0] << " " << B[i].z[0] - << std::endl; -#endif // TBB2 - } - vtxFieldRegion.Set(B[0], z[0], B[1], z[1], B[2], z[2]); - vtxFieldValue = B[0]; - } + //int ind = 0; + //{ + // L1FieldValue B[3]; + // fvec z[3]; + // for (int i = 0; i < 3; i++) { + // z[i] = geo[ind++]; + // B[i].x = geo[ind++]; + // B[i].y = geo[ind++]; + // B[i].z = geo[ind++]; +//#ifndef TBB2 + // std::cout << "L1Algo Input Magnetic field:" << z[i][0] << " " << B[i].x[0] << " " << B[i].y[0] << " " << B[i].z[0] + // << std::endl; +//#endif // TBB2 + // } + // //fVtxFieldRegion.Set(B[0], z[0], B[1], z[1], B[2], z[2]); + // //fVtxFieldValue = B[0]; + //} //vStations.clear(); - L1Station vStations[L1Parameters::kMaxNstations] _fvecalignment; - int NStations = static_cast<int>(geo[ind++]); - NMvdStations = static_cast<int>(geo[ind++]); // TODO: get rid of NMbdStations (S. Zh.) - NStsStations = static_cast<int>(geo[ind++]); // TODO: get rid of NStsStations (S. Zh.) + //L1Station vStations[L1Parameters::kMaxNstations] _fvecalignment; + //int NStations = static_cast<int>(geo[ind++]); + + //int NMvdStations = static_cast<int>(geo[ind++]); // TODO: get rid of NMbdStations (S. Zh.) + int nStationsSts = fInitManager.GetStationsNumber(static_cast<L1DetectorID>(1)); + fNstationsBeforePipe = fInitManager.GetStationsNumber(static_cast<L1DetectorID>(0)); + //int NStsStations = static_cast<int>(geo[ind++]); // TODO: get rid of NStsStations (S. Zh.) - fNfieldStations = NStsStations + NMvdStations; + fNfieldStations = nStationsSts + fNstationsBeforePipe; // TODO: Provide special getter for it (S.Zharko, 12.05.2022) if (fTrackingMode == kMcbm) { fNfieldStations = -1; } // cout << "N MVD & STS stations: " << NMvdStations << " " << NStations-NMvdStations << endl; -#ifndef TBB2 - std::cout << "L1Algo Input " << NStations << " Stations:" << std::endl; -#endif // TBB2 - for (int i = 0; i < NStations; i++) { - L1Station& st = vStations[i]; - st.type = geo[ind++]; - st.timeInfo = 1; - if (st.type == 1) st.timeInfo = 0; - st.z = geo[ind++]; - st.materialInfo.thick = geo[ind++]; - st.Rmin = geo[ind++]; - st.Rmax = geo[ind++]; - st.materialInfo.RL = geo[ind++]; - st.materialInfo.RadThick = st.materialInfo.thick / st.materialInfo.RL; - st.materialInfo.logRadThick = log(st.materialInfo.RadThick); - - double f_phi = geo[ind++]; - double f_sigma = geo[ind++]; - double b_phi = geo[ind++]; - double b_sigma = geo[ind++]; - double dt = geo[ind++]; //TODO: Add this field to L1BaseStationInfo and to ToString fcn (S.Zharko) - double c_f = cos(f_phi); - double s_f = sin(f_phi); - double c_b = cos(b_phi); - double s_b = sin(b_phi); - - st.frontInfo.cos_phi = c_f; - st.frontInfo.sin_phi = s_f; - st.frontInfo.sigma2 = f_sigma * f_sigma; - - st.backInfo.cos_phi = c_b; - st.backInfo.sin_phi = s_b; - st.backInfo.sigma2 = b_sigma * b_sigma; - st.dt = dt; - - if (fabs(b_phi - f_phi) < .1) { - double th = b_phi - f_phi; - double det = cos(th); - det *= det; - st.XYInfo.C00 = (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det; - st.XYInfo.C10 = -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det; - st.XYInfo.C11 = (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det; - //std::cout << "!!! det "<< det<<std::endl; - } - else { - double det = c_f * s_b - s_f * c_b; - det *= det; - st.XYInfo.C00 = (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det; - st.XYInfo.C10 = -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det; - st.XYInfo.C11 = (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det; - //std::cout << "!!! det "<< det<<std::endl; - } - //std::cout.precision(10); - //std::cout << "Station "<<i<<" " << st.XYInfo.C00[0]<<" "<<st.XYInfo.C11[0]<<" "<<st.XYInfo.C10[0]<<std::endl; - //std::cout << " "<< i<<" fsigma " << f_sigma<<" bsigma "<<b_sigma<<std::endl; - - // st.xInfo.cos_phi = c_f/(c_f*s_b - c_b*s_f); - // st.xInfo.sin_phi =-c_b/(c_f*s_b - c_b*s_f); - st.xInfo.cos_phi = -s_f / (c_f * s_b - c_b * s_f); - st.xInfo.sin_phi = s_b / (c_f * s_b - c_b * s_f); - st.xInfo.sigma2 = st.XYInfo.C00; - - st.yInfo.cos_phi = c_b / (c_b * s_f - c_f * s_b); - st.yInfo.sin_phi = -c_f / (c_b * s_f - c_f * s_b); - st.yInfo.sigma2 = st.XYInfo.C11; - //std::cout << "st.xInfo.cos_phi "<<st.xInfo.cos_phi<< " st.yInfo.cos_phi " << st.yInfo.cos_phi << std::endl; - //std::cout << "st.xInfo.sin_phi "<<st.xInfo.sin_phi<< " st.yInfo.sin_phi " << st.yInfo.sin_phi << std::endl; - - //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++]; - for (int iC = 0; iC < N; iC++) - st.fieldSlice.cy[iC] = geo[ind++]; - for (int iC = 0; iC < N; iC++) - st.fieldSlice.cz[iC] = geo[ind++]; -#ifndef TBB2 - std::cout << " " << st.z[0] << " " << st.materialInfo.thick[0] << " " << st.materialInfo.RL[0] << ", " << N - << " field coeff." << std::endl; - std::cout << " " << f_phi << " " << f_sigma << " " << b_phi << " " << b_sigma << std::endl; -#endif // TBB2 - } - - fTrackingLevel = static_cast<int>(geo[ind++]); - fMomentumCutOff = geo[ind++]; - fGhostSuppression = static_cast<int>(geo[ind++]); - - { - // fvec By0 = vStations[NStations-1].fieldSlice.cy[0]; - fvec z0 = vStations[NStations - 1].z; - fvec sy = 0., Sy = 0.; - for (int i = NStations - 1; i >= 0; i--) { - L1Station& st = vStations[i]; - fvec dz = st.z - z0; - fvec By = vStations[i].fieldSlice.cy[0]; - 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; - } - } +//#ifndef TBB2 + //std::cout << "L1Algo Input " << NStations << " Stations:" << std::endl; +//#endif // TBB2 + // for (int i = 0; i < NStations; i++) { + // L1Station& st = vStations[i]; + // st.type = geo[ind++]; + // st.timeInfo = 1; + // if (st.type == 1) st.timeInfo = 0; + // st.z = geo[ind++]; + // st.materialInfo.thick = geo[ind++]; + // st.Rmin = geo[ind++]; + // st.Rmax = geo[ind++]; + // st.materialInfo.RL = geo[ind++]; + // st.materialInfo.RadThick = st.materialInfo.thick / st.materialInfo.RL; + // st.materialInfo.logRadThick = log(st.materialInfo.RadThick); + + // double f_phi = geo[ind++]; + // double f_sigma = geo[ind++]; + // double b_phi = geo[ind++]; + // double b_sigma = geo[ind++]; + // double dt = geo[ind++]; //TODO: Add this field to L1BaseStationInfo and to ToString fcn (S.Zharko) + // double c_f = cos(f_phi); + // double s_f = sin(f_phi); + // double c_b = cos(b_phi); + // double s_b = sin(b_phi); + + // st.frontInfo.cos_phi = c_f; + // st.frontInfo.sin_phi = s_f; + // st.frontInfo.sigma2 = f_sigma * f_sigma; + + // st.backInfo.cos_phi = c_b; + // st.backInfo.sin_phi = s_b; + // st.backInfo.sigma2 = b_sigma * b_sigma; + // st.dt = dt; + + // if (fabs(b_phi - f_phi) < .1) { + // double th = b_phi - f_phi; + // double det = cos(th); + // det *= det; + // st.XYInfo.C00 = (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det; + // st.XYInfo.C10 = -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det; + // st.XYInfo.C11 = (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det; + // //std::cout << "!!! det "<< det<<std::endl; + // } + // else { + // double det = c_f * s_b - s_f * c_b; + // det *= det; + // st.XYInfo.C00 = (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det; + // st.XYInfo.C10 = -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det; + // st.XYInfo.C11 = (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det; + // //std::cout << "!!! det "<< det<<std::endl; + // } + // //std::cout.precision(10); + // //std::cout << "Station "<<i<<" " << st.XYInfo.C00[0]<<" "<<st.XYInfo.C11[0]<<" "<<st.XYInfo.C10[0]<<std::endl; + // //std::cout << " "<< i<<" fsigma " << f_sigma<<" bsigma "<<b_sigma<<std::endl; + + // // st.xInfo.cos_phi = c_f/(c_f*s_b - c_b*s_f); + // // st.xInfo.sin_phi =-c_b/(c_f*s_b - c_b*s_f); + // st.xInfo.cos_phi = -s_f / (c_f * s_b - c_b * s_f); + // st.xInfo.sin_phi = s_b / (c_f * s_b - c_b * s_f); + // st.xInfo.sigma2 = st.XYInfo.C00; + + // st.yInfo.cos_phi = c_b / (c_b * s_f - c_f * s_b); + // st.yInfo.sin_phi = -c_f / (c_b * s_f - c_f * s_b); + // st.yInfo.sigma2 = st.XYInfo.C11; + // //std::cout << "st.xInfo.cos_phi "<<st.xInfo.cos_phi<< " st.yInfo.cos_phi " << st.yInfo.cos_phi << std::endl; + // //std::cout << "st.xInfo.sin_phi "<<st.xInfo.sin_phi<< " st.yInfo.sin_phi " << st.yInfo.sin_phi << std::endl; + + // //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++]; + // for (int iC = 0; iC < N; iC++) + // st.fieldSlice.cy[iC] = geo[ind++]; + // for (int iC = 0; iC < N; iC++) + // st.fieldSlice.cz[iC] = geo[ind++]; +//#ifndef TBB2 +// std::cout << " " << st.z[0] << " " << st.materialInfo.thick[0] << " " << st.materialInfo.RL[0] << ", " << N +// << " field coeff." << std::endl; +// std::cout << " " << f_phi << " " << f_sigma << " " << b_phi << " " << b_sigma << std::endl; +//#endif // TBB2 + //} + + //fTrackingLevel = static_cast<int>(geo[ind++]); + //fMomentumCutOff = geo[ind++]; + //fGhostSuppression = static_cast<int>(geo[ind++]); + + //{ + // // fvec By0 = vStations[NStations-1].fieldSlice.cy[0]; + // fvec z0 = vStations[NStations - 1].z; + // fvec sy = 0., Sy = 0.; + // for (int i = NStations - 1; i >= 0; i--) { + // L1Station& st = vStations[i]; + // fvec dz = st.z - z0; + // fvec By = vStations[i].fieldSlice.cy[0]; + // 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; + // } + //} // for( int iS = 0; iS < NStations; ++iS ) { /// Grid is created for each station with the same step: xStep,yStep // L1Grid &grid = vGrid[iS]; // @@ -220,9 +225,9 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra // grid.Create(-1,1,-0.6,0.6,0.00317899,0.00105966); // } -#ifndef TBB2 - std::cout << "L1Algo initialized" << std::endl; -#endif // TBB2 +//#ifndef TBB2 + //std::cout << "L1Algo initialized" << std::endl; +//#endif // TBB2 // // NEW INITIALIZATION (BETA) @@ -236,52 +241,38 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra LOG(info) << "InitManager " << fInitManager.GetInitController().ToString(); // Get number of stations - int nStationsNew = fInitManager.GetStationsNumber(); + //int nStationsNew = fInitManager.GetStationsNumber(); // TODO: we must to get rid of station specification in the L1Algo (S.Zh.) - int nMvdStationsNew = fInitManager.GetStationsNumber(static_cast<L1DetectorID>(0)); - int nStsStationsNew = fInitManager.GetStationsNumber(static_cast<L1DetectorID>(1)); - int nFieldStationsNew = nMvdStationsNew + nStsStationsNew; + //int nMvdStationsNew = fInitManager.GetStationsNumber(static_cast<L1DetectorID>(0)); + //int nStsStationsNew = fInitManager.GetStationsNumber(static_cast<L1DetectorID>(1)); + //int nFieldStationsNew = nMvdStationsNew + nStsStationsNew; fNstations = fInitManager.GetStationsNumber(); // Get field near target - L1FieldValue vtxFieldValueNew = fInitManager.GetTargetFieldValue(); - L1FieldRegion vtxFieldRegionNew = fInitManager.GetTargetFieldRegion(); + fVtxFieldValue = fInitManager.GetTargetFieldValue(); + fVtxFieldRegion = fInitManager.GetTargetFieldRegion(); // Fill L1Station array fInitManager.TransferL1StationArray(fStations); - LOG(info) << "**********************************************************************"; - LOG(info) << "* New L1Algo initialization cross check (tmp log, to be removed!) *"; - LOG(info) << "**********************************************************************"; - LOG(info) << "** Number of stations (origial) **"; - LOG(info) << "\tTotal: " << NStations; - LOG(info) << "\tMVD: " << NMvdStations; - LOG(info) << "\tSTS: " << NStsStations; - LOG(info) << "\tField: " << fNfieldStations; + LOG(info) << " ***********************"; + LOG(info) << " * L1Algo parameters *"; + LOG(info) << " ***********************"; + LOG(info) << ""; + LOG(info) << "** Number of stations (new) **"; - LOG(info) << "\tTotal: " << fNstations; - LOG(info) << "\tMVD: " << nMvdStationsNew; - LOG(info) << "\tSTS: " << nStsStationsNew; - LOG(info) << "\tField: " << nFieldStationsNew; - - LOG(info) << "** Magnetic field near target (original)**"; - LOG(info) << "\tField Value: " << '\n' << vtxFieldValue.ToString(/*indent = */ 1) << '\n'; - LOG(info) << "\tField Region: " << '\n' << vtxFieldRegion.ToString(/*indent = */ 1) << '\n'; - - LOG(info) << "** Magnetic field near target (new)**"; - LOG(info) << "\tField Value: " << '\n' << vtxFieldValueNew.ToString(/*indent = */ 1) << '\n'; - LOG(info) << "\tField Region: " << '\n' << vtxFieldRegionNew.ToString(/*indent = */ 1) << '\n'; - - LOG(info) << "** Original L1Station array content **"; - int nStations = fInitManager.GetStationsNumber(); - for (int iSt = 0; iSt < NStations; ++iSt) { - LOG(info) << "Station Global No: " << iSt << ' ' << vStations[iSt].ToString(/*verbosity = */ 0); - } - LOG(info) << "** New L1Station array content **"; + LOG(info) << "\tTotal stations: " << fNstations; + LOG(info) << "\tStations before pipe: " << fNstationsBeforePipe; + LOG(info) << "\tStations in field region: " << fNfieldStations; + + LOG(info) << "** Magnetic field near target **"; + LOG(info) << "\tField Value: " << '\n' << fVtxFieldValue.ToString(/*indent = */ 1) << '\n'; + LOG(info) << "\tField Region: " << '\n' << fVtxFieldRegion.ToString(/*indent = */ 1) << '\n'; + + LOG(info) << "** L1 active stations list **"; for (int iSt = 0; iSt < fNstations; ++iSt) { - LOG(info) << "Station Global No: " << iSt; - LOG(info) << '\n' << fStations[iSt].ToString(/*verbosity = */ 0); + LOG(info) << " #" << iSt << " " << fStations[iSt].ToString(/*verbosity = */ 0); } // Print L1Parameters diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h index 5560a8695a..b62a2e4fcb 100644 --- a/reco/L1/L1Algo/L1Algo.h +++ b/reco/L1/L1Algo/L1Algo.h @@ -79,10 +79,13 @@ class L1AlgoEfficiencyPerformance; #endif typedef int Tindex; +using L1StationsArray_t = std::array<L1Station, L1Parameters::kMaxNstations>; + /// Central class of L1 tracking /// class L1Algo { public: + L1Algo(unsigned int nThreads = 1); L1Algo(const L1Algo&) = delete; @@ -220,24 +223,22 @@ public: void SetNThreads(unsigned int n); private: - int fNstations {0}; ///< number of all detector stations + int fNstations {0}; ///< number of all detector stations + int fNstationsBeforePipe {0}; ///< number of stations before pipe (MVD stations in CBM) + int fNfieldStations {0}; ///< number of stations in the field region + alignas(16) L1StationsArray_t fStations; ///< array of L1Station objects + public: - + /// Gets total number of stations used in tracking int GetNstations() const { return fNstations; } + /// Gets number of stations before the pipe (MVD stations in CBM) + int GetNstationsBeforePipe() const { return fNstationsBeforePipe; } + /// Gets number of stations situated in field region (MVD + STS in CBM) + int GetNfieldStations() const { return fNfieldStations;} + /// Gets reference to the stations array + const L1StationsArray_t& GetStations() const { return fStations; } - int NMvdStations {0}; ///< number of mvd stations - int NStsStations {0}; ///< number of sts stations - int fNfieldStations {0}; ///< number of stations in the field region - - - - // TODO: Replace _fvecalignment with C++11 alignas(16) attibute, see vStationsNew (S.Zh.) - //L1Station vStations[L1Parameters::kMaxNstations] _fvecalignment; // station info - -private: - alignas(16) std::array<L1Station, L1Parameters::kMaxNstations> fStations; public: - const std::array<L1Station, L1Parameters::kMaxNstations>& GetStations() const { return fStations; } L1Vector<L1Material> fRadThick {"fRadThick"}; // material for each station @@ -357,8 +358,8 @@ public: friend class CbmL1; - const L1FieldValue& GetVtxFieldValue() const { return vtxFieldValue; } - const L1FieldRegion& GetVtxFieldRegion() const { return vtxFieldRegion; } + const L1FieldValue& GetVtxFieldValue() const { return fVtxFieldValue; } + const L1FieldRegion& GetVtxFieldRegion() const { return fVtxFieldRegion; } /// ----- Hit-point-strips conversion routines ------ void GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, fscal& _z, const L1Station& sta); @@ -702,14 +703,14 @@ private: L1XYMeasurementInfo TargetXYInfo _fvecalignment {}; // target constraint [cm] - L1FieldRegion vtxFieldRegion _fvecalignment {}; // really doesn't used - L1FieldValue vtxFieldValue _fvecalignment {}; // field at teh vertex position. + L1FieldRegion fVtxFieldRegion _fvecalignment {}; // really doesn't used + L1FieldValue fVtxFieldValue _fvecalignment {}; // field at teh vertex position. // int TripNumThread; - int fTrackingLevel {0}; // currently not used - int fGhostSuppression {0}; // currently not used - float fMomentumCutOff {0}; // currently not used + int fTrackingLevel {2}; // currently not used + int fGhostSuppression {1}; // currently not used + float fMomentumCutOff {0.2}; // currently not used /// ----- Debug features ----- #ifdef PULLS diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx index 556ff6d87d..15bf8ec644 100644 --- a/reco/L1/L1Algo/L1CATrackFinder.cxx +++ b/reco/L1/L1Algo/L1CATrackFinder.cxx @@ -187,8 +187,8 @@ inline void L1Algo::f11( /// input 1st stage of singlet search // estimate field for singlet creation int istac = istal / 2; // "center" station - // if (istal >= NMvdStations) // to make it in the same way for both with and w\o mvd - // istac = (istal-NMvdStations)/2+NMvdStations; + // if (istal >= fNstationsBeforePipe) // to make it in the same way for both with and w\o mvd + // istac = (istal - fNstationsBeforePipe) / 2 + fNstationsBeforePipe; L1Station& stac = fStations[istac]; fvec zstac = stac.z; @@ -333,7 +333,7 @@ inline void L1Algo::f11( /// input 1st stage of singlet search else { fit.L1AddMaterial(T, fStations[ista].materialInfo, fMaxInvMom, 1); } - if (ista == NMvdStations - 1) fit.L1AddPipeMaterial(T, fMaxInvMom, 1); + if (ista == fNstationsBeforePipe - 1) fit.L1AddPipeMaterial(T, fMaxInvMom, 1); } // add left hit @@ -372,7 +372,7 @@ inline void L1Algo::f11( /// input 1st stage of singlet search else { fit.L1AddMaterial(T, stal.materialInfo, fMaxInvMom, 1); } - if ((istam >= NMvdStations) && (istal <= NMvdStations - 1)) { fit.L1AddPipeMaterial(T, fMaxInvMom, 1); } + if ((istam >= fNstationsBeforePipe) && (istal <= fNstationsBeforePipe - 1)) { fit.L1AddPipeMaterial(T, fMaxInvMom, 1); } fvec dz = zstam - zl; L1ExtrapolateTime(T, dz, stam.timeInfo); @@ -694,7 +694,7 @@ inline void L1Algo::f30( // input else { fit.L1AddMaterial(T2, stam.materialInfo, T2.qp, 1); } - if ((istar >= NMvdStations) && (istam <= NMvdStations - 1)) { fit.L1AddPipeMaterial(T2, T2.qp, 1); } + if ((istar >= fNstationsBeforePipe) && (istam <= fNstationsBeforePipe - 1)) { fit.L1AddPipeMaterial(T2, T2.qp, 1); } fvec dz2 = star.z - T2.z; L1ExtrapolateTime(T2, dz2, stam.timeInfo); @@ -1039,7 +1039,7 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction else { fit.L1AddMaterial(T, sta[ih].materialInfo, T.qp, 1); } - if (ista[ih] == NMvdStations - 1) fit.L1AddPipeMaterial(T, T.qp, 1); + if (ista[ih] == fNstationsBeforePipe - 1) fit.L1AddPipeMaterial(T, T.qp, 1); L1Filter(T, sta[ih].frontInfo, u[ih]); L1Filter(T, sta[ih].backInfo, v[ih]); @@ -1075,7 +1075,7 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction else { fit.L1AddMaterial(T, sta[ih].materialInfo, T.qp, 1); } - if (ista[ih] == NMvdStations) fit.L1AddPipeMaterial(T, T.qp, 1); + if (ista[ih] == fNstationsBeforePipe) fit.L1AddPipeMaterial(T, T.qp, 1); L1Filter(T, sta[ih].frontInfo, u[ih]); L1Filter(T, sta[ih].backInfo, v[ih]); @@ -1107,7 +1107,7 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction else { fit.L1AddMaterial(T, sta[ih].materialInfo, T.qp, 1); } - if (ista[ih] == NMvdStations + 1) fit.L1AddPipeMaterial(T, T.qp, 1); + if (ista[ih] == fNstationsBeforePipe + 1) fit.L1AddPipeMaterial(T, T.qp, 1); L1Filter(T, sta[ih].frontInfo, u[ih]); L1Filter(T, sta[ih].backInfo, v[ih]); @@ -1894,8 +1894,8 @@ void L1Algo::CATrackFinder() float SigmaTargetX = caIteration.GetTargetPosSigmaX(); float SigmaTargetY = caIteration.GetTargetPosSigmaY(); // target constraint [cm] - // Select magnetic field. For primary tracks - vtxFieldValue, for secondary tracks - st.fieldSlice - if (caIteration.IsPrimary()) { fTargB = vtxFieldValue; } + // Select magnetic field. For primary tracks - fVtxFieldValue, for secondary tracks - st.fieldSlice + if (caIteration.IsPrimary()) { fTargB = fVtxFieldValue; } else { fStations[0].fieldSlice.GetFieldValue(0, 0, fTargB); } // NOTE: calculates field fTargB in the center of 0th station @@ -1903,7 +1903,7 @@ void L1Algo::CATrackFinder() //if ((isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter) || (isec == kAllPrimIter) // || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter)) { // target - // fTargB = vtxFieldValue; + // fTargB = fVtxFieldValue; // if ((isec == kFastPrimIter) || (isec == kAllPrimIter) || (isec == kAllPrimEIter)) // SigmaTargetX = SigmaTargetY = 1; // target // else diff --git a/reco/L1/L1Algo/L1TrackExtender.cxx b/reco/L1/L1Algo/L1TrackExtender.cxx index 8b370d9955..8f0ac6d786 100644 --- a/reco/L1/L1Algo/L1TrackExtender.cxx +++ b/reco/L1/L1Algo/L1TrackExtender.cxx @@ -140,8 +140,8 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir, L1ExtrapolateTime(T, dz); fit.L1AddMaterial(T, sta.materialInfo, qp0, 1); - if ((step * ista <= step * (NMvdStations + (step + 1) / 2 - 1)) - && (step * ista_prev >= step * (NMvdStations + (step + 1) / 2 - 1 - step))) + if ((step * ista <= step * (fNstationsBeforePipe + (step + 1) / 2 - 1)) + && (step * ista_prev >= step * (fNstationsBeforePipe + (step + 1) / 2 - 1 - step))) fit.L1AddPipeMaterial(T, qp0, 1); fvec u = hit.u; diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx index 5304cc071c..2f9528a9dc 100644 --- a/reco/L1/L1Algo/L1TrackFitter.cxx +++ b/reco/L1/L1Algo/L1TrackFitter.cxx @@ -142,7 +142,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe. fit.L1AddMaterial(T, sta.materialInfo, qp0, ONE); } - // if (ista==NMvdStations-1) fit.L1AddPipeMaterial( T, qp0, 1); + // if (ista == fNstationsBeforePipe - 1) fit.L1AddPipeMaterial( T, qp0, 1); fvec u = hit.u; fvec v = hit.v; @@ -280,7 +280,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe. else { fit.L1AddMaterial(T, sta.materialInfo, qp0, ONE); } - // if (ista==NMvdStations) fit.L1AddPipeMaterial( T, qp0, 1); + // if (ista == fNstationsBeforePipe) fit.L1AddPipeMaterial( T, qp0, 1); L1Filter(T, sta.frontInfo, u); L1Filter(T, sta.backInfo, v); @@ -394,7 +394,7 @@ void L1Algo::L1KFTrackFitter() const int ista = (*fStripFlag)[hit.f] / 4; iSta[i] = ista; w[ista][iVec] = 1.; - if (ista > NMvdStations) w_time[ista][iVec] = 1.; + if (ista > fNstationsBeforePipe) w_time[ista][iVec] = 1.; u[ista][iVec] = hit.u; v[ista][iVec] = hit.v; @@ -512,7 +512,7 @@ void L1Algo::L1KFTrackFitter() L1Extrapolate(T, z[i], qp0, fld, &w1); - if (i == NMvdStations - 1) { + if (i == fNstationsBeforePipe - 1) { fit.L1AddPipeMaterial(T, qp0, wIn); fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(1.f), wIn); @@ -674,7 +674,7 @@ void L1Algo::L1KFTrackFitter() // T1.ExtrapolateLine( z[i]); - if (i == NMvdStations) { + if (i == fNstationsBeforePipe) { fit.L1AddPipeMaterial(T, qp0, wIn); fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn); @@ -952,7 +952,7 @@ void L1Algo::L1KFTrackFitterMuch() T1.Extrapolate(z[i], qp01, fld, &w1); - if (i == NMvdStations) { + if (i == fNstationsBeforePipe) { fit.L1AddPipeMaterial(T, qp0, wIn); fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn); diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx index f0927c0cab..93a28ce71a 100644 --- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx +++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx @@ -93,7 +93,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo) FairRootManager* fManger = FairRootManager::Instance(); TClonesArray* listStsHits = (TClonesArray*) fManger->GetObject("StsHit"); - int NMvdStations = CbmL1::Instance()->algo->NMvdStations; + int NMvdStations = CbmL1::Instance()->algo->GetNstationsBeforePipe(); TClonesArray* listMvdHits = 0; if (NMvdStations > 0.) listMvdHits = (TClonesArray*) fManger->GetObject("MvdHit"); @@ -406,7 +406,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<L1FieldRe CbmStsTrack* t[fvecLen]; int nStations = CbmL1::Instance()->algo->GetNstations(); - int NMvdStations = CbmL1::Instance()->algo->NMvdStations; + int NMvdStations = CbmL1::Instance()->algo->GetNstationsBeforePipe(); const L1Station* sta = CbmL1::Instance()->algo->GetStations().begin(); fvec* zSta = new fvec[nStations]; for (int iSta = 0; iSta < nStations; iSta++) @@ -572,7 +572,7 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<L1F FairRootManager* fManger = FairRootManager::Instance(); TClonesArray* listStsHits = (TClonesArray*) fManger->GetObject("StsHit"); TClonesArray* listMvdHits = 0; - int NMvdStations = CbmL1::Instance()->algo->NMvdStations; + int NMvdStations = CbmL1::Instance()->algo->GetNstationsBeforePipe(); if (NMvdStations > 0.) listMvdHits = (TClonesArray*) fManger->GetObject("MvdHit"); int nTracks_SIMD = fvecLen; @@ -644,7 +644,7 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks, FairRootManager* fManger = FairRootManager::Instance(); TClonesArray* listStsHits = (TClonesArray*) fManger->GetObject("StsHit"); TClonesArray* listMvdHits = 0; - int NMvdStations = CbmL1::Instance()->algo->NMvdStations; + int NMvdStations = CbmL1::Instance()->algo->GetNstationsBeforePipe(); if (NMvdStations > 0.) listMvdHits = (TClonesArray*) fManger->GetObject("MvdHit"); int nTracks_SIMD = fvecLen; -- GitLab