Skip to content

GitLab

Commit e01908fe authored by Valentina Akishina's avatar Valentina Akishina
Browse files

L1: add broken triplets for mCBM

parent 33e953a2
Hide whitespace changes
Inline Side-by-side
Showing with 102 additions and 93 deletions
macro/mcbm/mcbm_reco_event.C
View file @ e01908fe
......@@ -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 --------------------------------------------
......
reco/L1/CbmL1.cxx
View file @ e01908fe
......@@ -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";
......
reco/L1/CbmL1.h
View file @ e01908fe
......@@ -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
......
reco/L1/CbmL1ReadEvent.cxx
View file @ e01908fe
......@@ -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];
......
reco/L1/L1Algo/L1Algo.cxx
View file @ e01908fe
......@@ -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++){
......
reco/L1/L1Algo/L1Algo.h
View file @ e01908fe
......@@ -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}};
......
reco/L1/L1Algo/L1CATrackFinder.cxx
View file @ e01908fe
......@@ -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]
);
);
}
} //
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment