diff --git a/macro/beamtime/mcbm2022/mcbm_event_reco.C b/macro/beamtime/mcbm2022/mcbm_event_reco.C index 734dc982f981bc2425b69cb092f040e1eb46711c..46c5eeeba201408b69e00d65d6a9b43738eb2981 100644 --- a/macro/beamtime/mcbm2022/mcbm_event_reco.C +++ b/macro/beamtime/mcbm2022/mcbm_event_reco.C @@ -852,77 +852,70 @@ Bool_t mcbm_event_reco(UInt_t uRunId = 2391, // === L1 === // ========================================================================= - // CbmKF* kalman = new CbmKF(); - // run->AddTask(kalman); - // CbmL1* l1 = new CbmL1(); - // l1->SetLegacyEventMode(1); - // l1->SetMcbmMode(); - // l1->SetUseHitErrors(1); - // if (strcmp(geoSetupTag.data(), "mcbm_beam_2021_07_surveyed") == 0) l1->SetMissingHits(1); - // - // // --- Material budget file names - // TString mvdGeoTag; - // if (geoSetup->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 (geoSetup->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 (geoSetup->GetGeoTag(ECbmModuleId::kMuch, muchGeoTag)) { - // - // // --- Parameter file name - // TString geoTag; - // geoSetup->GetGeoTag(ECbmModuleId::kMuch, geoTag); - // Int_t muchFlag = 0; - // if (geoTag.Contains("mcbm")) muchFlag = 1; - // - // 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 (geoSetup->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 (geoSetup->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()); - // } - // - - // Workaround to get it running: - // 1) Change fUseGlobal in line 129 of CbmStsParSetModule.h to - // Bool_t fUseGlobal = kTRUE; - // 2) Change fUseGlobal in line 114 of CbmStsParSetSensor.h to - // Bool_t fUseGlobal = kTRUE; - //run->AddTask(l1); - - // CbmL1GlobalTrackFinder* globalTrackFinder = new CbmL1GlobalTrackFinder(); - // FairTask* globalFindTracks = new CbmL1GlobalFindTracksEvents(globalTrackFinder); - //run->AddTask(globalFindTracks); + run->AddTask(new CbmTrackingDetectorInterfaceInit()); + + + 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 (geoSetup->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 (geoSetup->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 (geoSetup->GetGeoTag(ECbmModuleId::kMuch, muchGeoTag)) { + + // --- Parameter file name + TString geoTag; + geoSetup->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 trdGeoTag; + if (geoSetup->GetGeoTag(ECbmModuleId::kTrd, trdGeoTag)) { + TString parFile = gSystem->Getenv("VMCWORKDIR"); + parFile = parFile + "/parameters/trd/trd_matbudget_v22a_mcbm.root"; //trd_matbudget_" + trdGeoTag + ".root "; + std::cout << "Using material budget file " << parFile << std::endl; + l1->SetTrdMaterialBudgetFileName(parFile.Data()); + } + + TString tofGeoTag; + if (geoSetup->GetGeoTag(ECbmModuleId::kTof, tofGeoTag)) { + TString parFile = gSystem->Getenv("VMCWORKDIR"); + // parFile = parFile + "/parameters/tof/tof_matbudget_" + tofGeoTag + ".root "; + parFile = parFile + "/parameters/tof/tof_matbudget_v21d_mcbm.root"; + std::cout << "Using material budget file " << parFile << std::endl; + l1->SetTofMaterialBudgetFileName(parFile.Data()); + } + run->AddTask(l1); + + CbmL1GlobalTrackFinder* globalTrackFinder = new CbmL1GlobalTrackFinder(); + FairTask* globalFindTracks = new CbmL1GlobalFindTracksEvents(globalTrackFinder); + run->AddTask(globalFindTracks); // ========================================================================= // === QA === diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index f80d034c45dac1eccc1b680c23ef04c99d959217..7f6c60d73e8491834221eb9fa1c060c3740ce8e0 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -184,7 +184,7 @@ InitStatus CbmL1::Init() // Init L1 algo core - fpAlgo = &gAlgo; + fpAlgo = &gAlgo; fpAlgo->Init(fUseHitErrors, fTrackingMode, fMissingHits); fHistoDir = gROOT->mkdir("L1"); @@ -523,6 +523,8 @@ InitStatus CbmL1::Init() auto stationInfo = L1BaseStationInfo(L1DetectorID::kSts, iSt); stationInfo.SetStationType(0); // STS stationInfo.SetTimeInfo(stsInterface->IsTimeInfoProvided(iSt)); + stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) + : false); stationInfo.SetTimeResolution(stsInterface->GetTimeResolution(iSt)); stationInfo.SetFieldStatus(L1Algo::TrackingMode::kMcbm == fTrackingMode ? 0 : 1); stationInfo.SetZ(stsInterface->GetZ(iSt)); @@ -548,6 +550,8 @@ InitStatus CbmL1::Init() auto stationInfo = L1BaseStationInfo(L1DetectorID::kMuch, iSt); stationInfo.SetStationType(2); // MuCh stationInfo.SetTimeInfo(muchInterface->IsTimeInfoProvided(iSt)); + stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) + : false); stationInfo.SetTimeResolution(muchInterface->GetTimeResolution(iSt)); stationInfo.SetFieldStatus(0); stationInfo.SetZ(muchInterface->GetZ(iSt)); @@ -573,6 +577,8 @@ InitStatus CbmL1::Init() auto stationInfo = L1BaseStationInfo(L1DetectorID::kTrd, iSt); stationInfo.SetStationType((iSt == 1 || iSt == 3) ? 6 : 3); // MuCh stationInfo.SetTimeInfo(trdInterface->IsTimeInfoProvided(iSt)); + stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) + : false); stationInfo.SetTimeResolution(trdInterface->GetTimeResolution(iSt)); stationInfo.SetFieldStatus(0); stationInfo.SetZ(trdInterface->GetZ(iSt)); @@ -608,6 +614,8 @@ InitStatus CbmL1::Init() auto stationInfo = L1BaseStationInfo(L1DetectorID::kTof, iSt); stationInfo.SetStationType(4); stationInfo.SetTimeInfo(tofInterface->IsTimeInfoProvided(iSt)); + stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) + : false); stationInfo.SetTimeResolution(tofInterface->GetTimeResolution(iSt)); stationInfo.SetFieldStatus(0); stationInfo.SetZ(tofInterface->GetZ(iSt)); diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx index a015b73cf0b565bc5ab3fcbb122d91ae3ea73254..d6d081f62dce67713685d550e68c951e33ed786b 100644 --- a/reco/L1/CbmL1ReadEvent.cxx +++ b/reco/L1/CbmL1ReadEvent.cxx @@ -48,6 +48,7 @@ #include "TDatabasePDG.h" #include "TRandom.h" + #include <iostream> using std::cout; @@ -186,8 +187,8 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int fvMCTracks.clear(); /* <CbmL1MCTrack> */ fvExternalHits.clear(); /* <CbmL1Hit> */ fvRecoTracks.clear(); /* <CbmL1Track> */ // TODO: Move from this function to more suitable place (S.Zharko) - fvHitPointIndexes.clear(); /* <int>: indexes of MC-points in fvMCPoints (by index of fpAlgo->vHits) */ - fvHitStore.clear(); /* <CbmL1HitStore> */ + fvHitPointIndexes.clear(); /* <int>: indexes of MC-points in fvMCPoints (by index of fpAlgo->vHits) */ + fvHitStore.clear(); /* <CbmL1HitStore> */ fmMCPointsLinksMap.clear(); fmMCTracksLinksMap.clear(); @@ -551,9 +552,7 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int if (fpMvdHitMatches) { CbmMatch* mvdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fpMvdHitMatches->At(j)); if (mvdHitMatch->GetNofLinks() > 0) - if (mvdHitMatch->GetLink(0).GetIndex() < nMvdPoints) { - th.iMC = mvdHitMatch->GetLink(0).GetIndex(); - } + if (mvdHitMatch->GetLink(0).GetIndex() < nMvdPoints) { th.iMC = mvdHitMatch->GetLink(0).GetIndex(); } } } //if( h.MC_Point >=0 ) // DEBUG !!!! @@ -639,10 +638,7 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int /// stop if reco TS ends and many hits left if (!event) - if ((th.time > (TsStart + TsLength)) && ((nEntSts - hitIndex) > 300)) { - - break; - } + if ((th.time > (TsStart + TsLength)) && ((nEntSts - hitIndex) > 300)) { break; } TVector3 pos, err; mh->Position(pos); @@ -804,17 +800,22 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int vector<bool> mcUsed(nTrdPoints, 0); - for (int iHit = 0; iHit < fpTrdHits->GetEntriesFast(); iHit++) { + Int_t nEntTrd = (event ? event->GetNofData(ECbmDataType::kTrdHit) : fpTrdHits->GetEntriesFast()); + + for (int iHit = 0; iHit < nEntTrd; iHit++) { TmpHit th; - CbmTrdHit* mh = L1_DYNAMIC_CAST<CbmTrdHit*>(fpTrdHits->At(iHit)); + Int_t hitIndex = 0; + hitIndex = (event ? event->GetIndex(ECbmDataType::kTrdHit, iHit) : iHit); + + CbmTrdHit* mh = L1_DYNAMIC_CAST<CbmTrdHit*>(fpTrdHits->At(hitIndex)); if ((L1Algo::TrackingMode::kGlobal == fTrackingMode) && (int) mh->GetClassType() != 1) { // SGtrd2d!! skip TRD-1D hit continue; } - th.ExtIndex = iHit; + th.ExtIndex = hitIndex; th.Det = 3; th.id = tmpHits.size(); @@ -927,13 +928,20 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int int firstDetStrip = NStrips; - for (int j = 0; j < fpTofHits->GetEntriesFast(); j++) { + Int_t nEntTof = (event ? event->GetNofData(ECbmDataType::kTofHit) : fpTofHits->GetEntriesFast()); + + for (int j = 0; j < nEntTof; j++) { + + + Int_t hitIndex = 0; + hitIndex = (event ? event->GetIndex(ECbmDataType::kTofHit, j) : j); + TmpHit th; - CbmTofHit* mh = L1_DYNAMIC_CAST<CbmTofHit*>(fpTofHits->At(j)); + CbmTofHit* mh = L1_DYNAMIC_CAST<CbmTofHit*>(fpTofHits->At(hitIndex)); - th.ExtIndex = j; + th.ExtIndex = hitIndex; th.Det = 4; th.id = tmpHits.size(); @@ -1035,7 +1043,7 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int int maxHitIndex = 0; for (int ih = 0; ih < nHits; ih++) { - TmpHit& th = tmpHits[ih]; + TmpHit& th = tmpHits[ih]; if (maxHitIndex < th.id) { maxHitIndex = th.id; } } // ih @@ -1264,8 +1272,8 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int M if (!pt) return 1; if (!fLegacyEventMode) { - Double_t StartTime = fTimeSlice->GetStartTime(); // not used - Double_t EndTime = fTimeSlice->GetEndTime(); // not used + Double_t StartTime = fTimeSlice->GetStartTime(); // not used + Double_t EndTime = fTimeSlice->GetEndTime(); // not used Double_t Time_MC_point = pt->GetTime() + fpEventList->GetEventTime(event, file); // not used if (StartTime > 0) @@ -1388,7 +1396,7 @@ void CbmL1::HitMatch() } } - Float_t bestWeight = 0.f; + Float_t bestWeight = 0.f; for (Int_t iLink = 0; iLink < stsHitMatch.GetNofLinks(); iLink++) { const CbmLink& link = stsHitMatch.GetLink(iLink); Int_t iFile = link.GetFile();