diff --git a/algo/ca/core/data/CaHit.cxx b/algo/ca/core/data/CaHit.cxx index d719c97818afd6ebf7bf69cca247269d71f5fe18..fbc134729df6650245d3a0c9f0fe8928a092843f 100644 --- a/algo/ca/core/data/CaHit.cxx +++ b/algo/ca/core/data/CaHit.cxx @@ -14,6 +14,15 @@ using cbm::algo::ca::Hit; +std::string Hit::ToString() const +{ + std::stringstream ss; + ss << "Hit " << fId << " station " << fStation << " front key " << fFrontKey << " back key " << fBackKey << " X " + << fX << " Y " << fY << " Z " << fZ << " T " << fT << " dX2 " << fDx2 << " dY2 " << fDy2 << " dXY " << fDxy + << " dT2 " << fDt2 << " rangeX " << fRangeX << " rangeY " << fRangeY << " rangeT " << fRangeT; + return ss.str(); +} + std::string Hit::ToString(int verbose, bool bHeader) const { if (verbose < 1) { return ""; } diff --git a/algo/ca/core/data/CaHit.h b/algo/ca/core/data/CaHit.h index 218969ff1184db808d4ca9017098db5ac1a417cb..6a090f3c05ed5bade8d0c022a96ba9ad4d71d80b 100644 --- a/algo/ca/core/data/CaHit.h +++ b/algo/ca/core/data/CaHit.h @@ -10,6 +10,9 @@ #include <boost/serialization/access.hpp> +#include <sstream> +#include <string> + #include "CaSimd.h" namespace cbm::algo::ca @@ -122,10 +125,13 @@ namespace cbm::algo::ca /// Get the station index int Station() const { return fStation; } + /// \brief Simple string representation of the hit class + std::string ToString() const; + /// \brief String representation of the hit class /// \param verbose Verbosity level /// \param bHeader If true, prints the header - std::string ToString(int verbose = 2, bool bHeader = false) const; + std::string ToString(int verbose, bool bHeader) const; private: ///----------------------------------------------------------------------------- diff --git a/algo/ca/core/pars/CaMaterialMonitor.cxx b/algo/ca/core/pars/CaMaterialMonitor.cxx index 1586b983104f1f6e1f0c8de239bac4dc2916e309..c204f13f06927a240d4d5077bc4da5b7aa2e4c71 100644 --- a/algo/ca/core/pars/CaMaterialMonitor.cxx +++ b/algo/ca/core/pars/CaMaterialMonitor.cxx @@ -132,7 +132,7 @@ namespace cbm::algo::ca else { if (fNhitsTotal > 0) { ss << "No active material. "; } else { - ss << "No hits to identify active areas. "; + ss << "No hits yet to identify active areas. "; } } @@ -151,7 +151,7 @@ namespace cbm::algo::ca } } else { - ss << "No hit statistics. "; + ss << "No hit statistics yet. "; } return ss.str(); diff --git a/algo/ca/core/tracking/CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx index 27522753ddbdbe7cbb6de6de97bf56cce6a7bf90..90b3d382476e5010f22c12c92e290e9e7522f61e 100644 --- a/algo/ca/core/tracking/CaTrackFinder.cxx +++ b/algo/ca/core/tracking/CaTrackFinder.cxx @@ -69,11 +69,16 @@ void TrackFinder::FindTracks() } frAlgo.fvHitKeyFlags.reset(frAlgo.fInputData.GetNhitKeys(), 0); + frAlgo.fStripToTrack.reset(frAlgo.fInputData.GetNhitKeys(), -1); frAlgo.fHitTimeInfo.reset(frAlgo.fInputData.GetNhits()); + // TODO: move these values to Parameters namespace (S.Zharko) - const fscal tsLength = 10000; // length of sub-TS - const fscal tsOverlap = 15; // min length of overlap region + + // length of sub-TS + const fscal tsLength = (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) ? 500 : 10000; + + const fscal tsOverlap = 15; // min length of overlap region fscal minProtonMomentum = 0.1; fscal targX = frAlgo.GetParameters().GetTargetPositionX()[0]; @@ -82,6 +87,10 @@ void TrackFinder::FindTracks() fscal tsStart = std::numeric_limits<fscal>::max(); // starting time of sub-TS + fscal statTsStart = 0.; // end time of the TS + fscal statTsEnd = 0.; // end time of the TS + int statNhitsTotal = 0; + // calculate possible event time for the hits (frAlgo.fHitTimeInfo array) const int nDataStreams = frAlgo.fInputData.GetNdataStreams(); @@ -91,6 +100,7 @@ void TrackFinder::FindTracks() fscal maxTimeBeforeHit = std::numeric_limits<fscal>::lowest(); int nStreamHits = frAlgo.fInputData.GetStreamNhits(iStream); + statNhitsTotal += nStreamHits; for (int ih = 0; ih < nStreamHits; ++ih) { @@ -119,12 +129,21 @@ void TrackFinder::FindTracks() info.fEventTimeMax = std::numeric_limits<fscal>::max(); } + if (info.fEventTimeMin > 500 * 1.e6) { // cut hits with bogus start time > 500 ms + frAlgo.fvHitKeyFlags[h.FrontKey()] = 1; + frAlgo.fvHitKeyFlags[h.BackKey()] = 1; + LOG(warning) << "CA: skip bogus hit " << h.ToString(); + continue; + } + + if (maxTimeBeforeHit < info.fEventTimeMax) { maxTimeBeforeHit = info.fEventTimeMax; } info.fMaxTimeBeforeHit = maxTimeBeforeHit; if (tsStart > info.fEventTimeMax - 5) { tsStart = info.fEventTimeMax - 5; // 5 ns margin } + if (statTsEnd < info.fEventTimeMin) { statTsEnd = info.fEventTimeMin; } } fscal minTimeAfterHit = std::numeric_limits<fscal>::max(); @@ -138,9 +157,17 @@ void TrackFinder::FindTracks() } } + statTsStart = tsStart; + + LOG(info) << "CA tracker process time slice " << statTsStart / 1.e6 << " ms ... " << statTsEnd / 1.e6 << " ms"; + + if (statTsEnd > statTsStart + 500.e6) { // 500 ms maximal length of the TS + statTsEnd = statTsStart + 500.e6; + } + // cut data into sub-timeslices and process them one by one - bool areDataLeft = true; // is the whole TS processed + bool areUntouchedDataLeft = true; // is the whole TS processed ca::HitIndex_t sliceFirstHit[nDataStreams]; @@ -148,7 +175,12 @@ void TrackFinder::FindTracks() sliceFirstHit[iStream] = frAlgo.fInputData.GetStreamStartIndex(iStream); } - while (areDataLeft) { + int statNwindows = 0; + int statNhitsProcessed = 0; + int statLastLogTimeChunk = -1; + + while (areUntouchedDataLeft) { + frAlgo.fMonitor.IncrementCounter(ECounter::SubTS); // select the sub-slice hits @@ -156,10 +188,13 @@ void TrackFinder::FindTracks() frAlgo.fSliceHitIds[iS].clear(); } - areDataLeft = false; + areUntouchedDataLeft = false; // TODO: SG: skip empty regions and start the subslice with the earliest hit + statNwindows++; + int statNwindowHits = 0; + for (int iStream = 0; iStream < nDataStreams; ++iStream) { for (ca::HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < frAlgo.fInputData.GetStreamStopIndex(iStream); @@ -171,7 +206,7 @@ void TrackFinder::FindTracks() continue; } if (info.fEventTimeMin > tsStart + tsLength) { // the hit is too late for the sub slice - areDataLeft = true; + areUntouchedDataLeft = true; if (info.fMinTimeAfterHit > tsStart + tsLength) { // this hit and all later hits are out of the sub-slice break; @@ -180,6 +215,7 @@ void TrackFinder::FindTracks() else { if (tsStart <= info.fEventTimeMax) { // the hit belongs to the sub-slice frAlgo.fSliceHitIds[h.Station()].push_back(caHitId); + statNwindowHits++; if (info.fMaxTimeBeforeHit < tsStart + tsLength - tsOverlap) { // this hit and all hits before are before the overlap sliceFirstHit[iStream] = caHitId + 1; @@ -188,6 +224,27 @@ void TrackFinder::FindTracks() } // else the hit has been alread processed in previous sub-slices } } + statNhitsProcessed += statNwindowHits; + + // print the LOG for every 10 ms of data processed + { + int currentChunk = (int) ((tsStart - statTsStart) / 10.e6); + if (!areUntouchedDataLeft || currentChunk > statLastLogTimeChunk) { + statLastLogTimeChunk = currentChunk; + LOG(info) << "CA tracker process sliding window N " << statNwindows << ": time " << tsStart / 1.e6 << " ms + " + << tsLength / 1.e3 << " us) with " << statNwindowHits << " hits. " + << " Processed " << 100. * tsStart / (statTsEnd - statTsStart) << " % of the TS time and " + << 100. * statNhitsProcessed / statNhitsTotal << " % of TS hits"; + } + } + + if (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) { + // cut at 50 hits per station per 1 us. + int maxStationHits = (int) (50 * tsLength / 1.e3); + for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ista++) { + if ((int) frAlgo.fSliceHitIds[ista].size() > maxStationHits) { frAlgo.fSliceHitIds[ista].clear(); } + } + } frAlgo.fMonitorDataPerCall.StartTimer(ETimer::TrackFinder); frAlgo.fTrackFinderWindow.CaTrackFinderSlice(); @@ -218,7 +275,7 @@ void TrackFinder::FindTracks() } } - if (!areDataLeft) { // don't reject tracks in the overlap when no more data are left + if (!areUntouchedDataLeft) { // don't reject tracks in the overlap when no more data are left isTrackCompletelyInOverlap = 0; } diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.cxx b/algo/ca/core/tracking/CaTrackFinderWindow.cxx index e00da9e137bf0147f74294cabb5bbfae2a1f1273..4f0abbff39120d8bc07028b1aa967fa636d56914 100644 --- a/algo/ca/core/tracking/CaTrackFinderWindow.cxx +++ b/algo/ca/core/tracking/CaTrackFinderWindow.cxx @@ -144,9 +144,6 @@ void TrackFinderWindow::ReadWindowData() frAlgo.fHitFirstTriplet.reset(nHits, 0); frAlgo.fHitNtriplets.reset(nHits, 0); - // TODO: SG: reduce the array size to the number of keys in subslice - - frAlgo.fStripToTrack.reset(frAlgo.fInputData.GetNhitKeys(), -1); frAlgo.fSliceRecoTracks.clear(); frAlgo.fSliceRecoTracks.reserve(2 * nHits / frAlgo.GetParameters().GetNstationsActive()); @@ -417,7 +414,10 @@ void TrackFinderWindow::CaTrackFinderSlice() frAlgo.fTrackCandidates.clear(); - frAlgo.fStripToTrack.reset(frAlgo.fInputData.GetNhitKeys(), -1); + for (const auto& h : frAlgo.fWindowHits) { + frAlgo.fStripToTrack[h.FrontKey()] = -1; + frAlgo.fStripToTrack[h.BackKey()] = -1; + } //== Loop over triplets with the required level, find and store track candidates diff --git a/algo/ca/core/tracking/CaTripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx index 2b0e2433f7ccac02206ebab771df99bca9ea6b38..df63e6c1cd1a2932596428a112c9aa1fdfaf0f45 100644 --- a/algo/ca/core/tracking/CaTripletConstructor.cxx +++ b/algo/ca/core/tracking/CaTripletConstructor.cxx @@ -56,9 +56,11 @@ void TripletConstructor::InitStations(int istal, int istam, int istar) sta2 = fIstaM; sta1 = fIstaM - 1; if (sta0 >= sta1) { sta0 = sta1 - 1; } + if (sta0 < 0) { sta0 = 0; } + } + if (fAlgo->fParameters.GetNstationsActive() >= 3) { + assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < fAlgo->fParameters.GetNstationsActive()); } - - assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < fAlgo->fParameters.GetNstationsActive()); fFld1Sta[0] = &fAlgo->fParameters.GetStation(sta0); fFld1Sta[1] = &fAlgo->fParameters.GetStation(sta1); diff --git a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C index 618ab5e27234c3a6559f5dc7028e437173a48026..4a3b17afdfac67bbf3526cff758c74809184b64b 100644 --- a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C +++ b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C @@ -285,7 +285,7 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570, if (bEvB) { CbmTaskBuildRawEvents* evBuildRaw = new CbmTaskBuildRawEvents(); - + evBuildRaw->SetVerbose(3); //Choose between NoOverlap, MergeOverlap, AllowOverlap evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::AllowOverlap); //For time being it is needed. We will remove CbmMuchBeamTimeDigi. @@ -370,8 +370,11 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570, if (bSTS && geoSetup->IsActive(ECbmModuleId::kSts)) { CbmRecoSts* recoSts = new CbmRecoSts(); - recoSts->SetMode(ECbmRecoMode::EventByEvent); - + if (bEvB) { recoSts->SetMode(ECbmRecoMode::EventByEvent); } + else { + recoSts->SetMode(ECbmRecoMode::Timeslice); + } + recoSts->SetVerbose(3); recoSts->SetTimeCutDigisAbs(20.0); // cluster finder: time cut in ns recoSts->SetTimeCutClustersAbs(20.0); // hit finder: time cut in ns @@ -598,68 +601,9 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570, if (bL1) { run->AddTask(new CbmTrackingDetectorInterfaceInit()); - CbmKF* kalman = new CbmKF(); - run->AddTask(kalman); - - CbmL1* l1 = new CbmL1(); + CbmL1* l1 = new CbmL1("L1"); l1->SetMcbmMode(); - - // --- 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()); - } - - - // 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; + l1->SetVerbose(3); run->AddTask(l1); CbmL1GlobalTrackFinder* globalTrackFinder = new CbmL1GlobalTrackFinder(); diff --git a/macro/beamtime/mcbm2022/mcbm_reco.C b/macro/beamtime/mcbm2022/mcbm_reco.C index 88b58f490e108703a43224dec20e685c258eed8e..3a3f0cf3a43531dbb8ff6642de241d46a55fbcba 100644 --- a/macro/beamtime/mcbm2022/mcbm_reco.C +++ b/macro/beamtime/mcbm2022/mcbm_reco.C @@ -521,23 +521,16 @@ Bool_t mcbm_reco(UInt_t uRunId = 2391, run->AddTask(new CbmTrackingDetectorInterfaceInit()); - CbmKF* kalman = new CbmKF(); - run->AddTask(kalman); - CbmL1* l1 = new CbmL1("L1", 0); // <= Disable verbose mode l1->SetMcbmMode(); - - /// PAL, 08/02/2023: Is the following comment still valid with the current master? - // 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); + + CbmKF* kalman = new CbmKF(); + run->AddTask(kalman); } // ========================================================================= // === QA === diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index f3f7cffbab3c1b0429c556af06af5c76ec918f2e..5c141ee94737088391870439d77f37f76a27e5c9 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -214,7 +214,7 @@ InitStatus CbmL1::Init() fUseMUCH = 1; fUseTRD = 1; fUseTOF = 1; - fInitManager.DevSetIgnoreHitSearchAreas(true); // uncomment for debug + // fInitManager.DevSetIgnoreHitSearchAreas(true); // uncomment for debug } if (ca::Framework::TrackingMode::kGlobal == fTrackingMode) { @@ -502,8 +502,7 @@ InitStatus CbmL1::Init() std::vector<ca::StationInitializer> vStations; vStations.reserve(100); - bool bDisableTime = - (ca::Framework::TrackingMode::kMcbm == fTrackingMode) ? true : false; /// TMP, move to parameters! + bool bDisableTime = false; /// TMP, move to parameters! // *** MVD stations info *** if (fUseMVD) { @@ -534,7 +533,6 @@ InitStatus CbmL1::Init() auto stationInfo = ca::StationInitializer(ca::EDetectorID::kSts, iSt); // TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID stationInfo.SetStationType(0); // STS - stationInfo.SetTimeInfo(stsInterface->IsTimeInfoProvided(iSt)); stationInfo.SetTimeInfo(!bDisableTime && stsInterface->IsTimeInfoProvided(iSt)); stationInfo.SetFieldStatus(ca::Framework::TrackingMode::kMcbm == fTrackingMode ? 0 : 1); stationInfo.SetZref(stsInterface->GetZref(iSt)); @@ -757,7 +755,12 @@ void CbmL1::Reconstruct(CbmEvent* event) } } - if (fVerbose > 1) { cout << "\nCbmL1::Exec event " << fEventNo << " ...\n\n"; } + if (fVerbose > 1) { + if (event) { LOG(info) << "\nCbmL1::Exec event " << fEventNo << " ...\n"; } + else { + LOG(info) << "\nCbmL1::Exec time slice " << fEventNo << " ...\n"; + } + } // ----- Read data from branches and send data from IODataManager to ca::Framework ---------------------------------------- @@ -1102,9 +1105,7 @@ void CbmL1::DefineSTAPNames(const char* dirName) fSTAPDataPrefix = fSTAPDataPrefix.Strip(TString::EStripType::kBoth, '.'); TString sDirName = TString(dirName); - if (sDirName.Length() == 0) { - fSTAPDataDir = pathToRecoOutput.parent_path().string(); - } + if (sDirName.Length() == 0) { fSTAPDataDir = pathToRecoOutput.parent_path().string(); } else if (bfs::exists(sDirName.Data()) && bfs::is_directory(sDirName.Data())) { fSTAPDataDir = sDirName; } diff --git a/reco/L1/OffLineInterface/CbmL1GlobalFindTracksEvents.cxx b/reco/L1/OffLineInterface/CbmL1GlobalFindTracksEvents.cxx index 1aa86de3f1d03e2311f055899683f1a33926e87c..57fbf011511574cc84f08efef90e7f66d103d0a7 100644 --- a/reco/L1/OffLineInterface/CbmL1GlobalFindTracksEvents.cxx +++ b/reco/L1/OffLineInterface/CbmL1GlobalFindTracksEvents.cxx @@ -117,7 +117,7 @@ void CbmL1GlobalFindTracksEvents::Exec(Option_t* /*opt*/) } //? event branch present else { // Timeslice reconstruction without events - ProcessEvent(nullptr); + result = ProcessEvent(nullptr); nHits = result.first; nTracks = result.second; }