From e1c026512cfbc8d14ab6d1353631cb6659dc8035 Mon Sep 17 00:00:00 2001 From: "s.zharko@gsi.de" <s.zharko@gsi.de> Date: Wed, 11 May 2022 13:42:32 +0200 Subject: [PATCH] L1Algo: minor documentation updates --- reco/L1/CbmL1.cxx | 22 +--- reco/L1/CbmL1.h | 40 ++++++- reco/L1/CbmL1ReadEvent.cxx | 172 ++++++++++++++++++++--------- reco/L1/L1Algo/L1Algo.h | 46 +++++--- reco/L1/L1Algo/L1CAIteration.cxx | 29 +++-- reco/L1/L1Algo/L1CAIteration.h | 74 +++++++------ reco/L1/L1Algo/L1CATrackFinder.cxx | 14 +-- reco/L1/L1Algo/L1Parameters.h | 22 +++- reco/L1/L1Algo/L1Track.h | 52 ++++++--- 9 files changed, 310 insertions(+), 161 deletions(-) diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index 5f27f4defa..f39bba9a8b 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -86,7 +86,7 @@ using std::ios; ClassImp(CbmL1) - static L1Algo algo_static _fvecalignment; + static L1Algo algo_static _fvecalignment; // TODO: gAlgo //L1AlgoInputData* fData_static _fvecalignment; @@ -1556,14 +1556,13 @@ 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(); L1Vector<std::pair<double, int>> SortStsHits("CbmL1::SortStsHits"); SortStsHits.reserve(listStsHits->GetEntriesFast()); - L1_SHOW(listStsHits->GetEntriesFast()); - float start_t = 10000000000; // TODO: move these values to CbmL1Parameters namespace (S.Zharko) @@ -1582,8 +1581,6 @@ void CbmL1::Reconstruct(CbmEvent* event) } TsStart = start_t; ///reco TS start time is set to smallest hit time - - L1_SHOW(TsStart); std::sort(SortStsHits.begin(), SortStsHits.end()); StsIndex.clear(); @@ -1593,18 +1590,13 @@ void CbmL1::Reconstruct(CbmEvent* event) StsIndex.push_back(j); }; - L1_SHOW(fLegacyEventMode); - L1_SHOW(TsStart); if (!fLegacyEventMode && fPerformance) { int nofEvents = fEventList->GetNofEvents(); - L1_SHOW(fEventList->GetNofEvents()); for (int iE = 0; iE < nofEvents; iE++) { int fileId = fEventList->GetFileIdByIndex(iE); int eventId = fEventList->GetEventIdByIndex(iE); - L1_SHOW(fileId); - L1_SHOW(eventId); vFileEvent.insert(DFSET::value_type(fileId, eventId)); } } @@ -1621,8 +1613,6 @@ void CbmL1::Reconstruct(CbmEvent* event) // repack data L1Vector<CbmL1Track> vRTracksCur("CbmL1::vRTracksCur"); // reconstructed tracks - L1_SHOW(vRTracksCur.size()); - L1_SHOW(vRTracksCur.capacity()); { int nHits = 0; int nSta = 1; @@ -1636,15 +1626,11 @@ void CbmL1::Reconstruct(CbmEvent* event) } vRTracksCur.reserve(10 + (2 * nHits) / nSta); } - L1_SHOW(vRTracksCur.size()); - L1_SHOW(vRTracksCur.capacity()); fTrackingTime = 0; while (areDataLeft) { - - L1_SHOW(areDataLeft); fData->Clear(); if (event) { @@ -1658,7 +1644,6 @@ void CbmL1::Reconstruct(CbmEvent* event) FstHitinTs = 0; } - L1_SHOW(fSTAPDataMode); if (fSTAPDataMode >= 2) { // 2,3 fData->ReadHitsFromFile(fSTAPDataDir.Data(), 1, fVerbose); @@ -1669,7 +1654,7 @@ void CbmL1::Reconstruct(CbmEvent* event) ReadEvent(fData, TsStart, TsLength, TsOverlap, FstHitinTs, areDataLeft, event); } - if (0) { // correct hits on MC // dbg + if constexpr (0) { // correct hits on MC // dbg TRandom3 random; L1Vector<int> strips("CbmL1::strips"); for (unsigned int iH = 0; iH < (*algo->vStsHits).size(); ++iH) { @@ -1717,7 +1702,6 @@ void CbmL1::Reconstruct(CbmEvent* event) } } - if (fPerformance) { HitMatch(); // calculate the max number of Hits\mcPoints on continuous(consecutive) stations diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h index 4944f6af50..745c94aa85 100644 --- a/reco/L1/CbmL1.h +++ b/reco/L1/CbmL1.h @@ -142,7 +142,7 @@ public: ~CbmL1(/*if (targetFieldSlice) delete;*/); - /// Gets a pointer to L1InitManager (for access in run_reco.C) + /// Gets a pointer to L1InitManager (for an access in run_reco.C) L1InitManager* GetInitManager() { return fpInitManager; } /// Gets a set of active detectors used in tracking // TODO: think about return (value, reference or const reference?) (S.Zh.) @@ -180,6 +180,8 @@ public: // void SetGhostSuppression( Bool_t b ){ fGhostSuppression= b; } // void SetDetectorEfficiency( Double_t eff ){ fDetectorEfficiency = eff; } + /// Reconstructs an event + /// \param event Pointer to current CbmEvent object void Reconstruct(CbmEvent* event = NULL); // bool ReadMCDataFromFile(const char work_dir[100], const int maxNEvent, const int iVerbose); @@ -201,17 +203,49 @@ private: /// Read information about hits, mcPoints and mcTracks into L1 classes + /// Repacks data from the external TClonesArray objects to the internal L1 arrays + /// \param fData_ Pointer to the target object containig L1Algo internal arrays of hits + /// \param TsStart Reference to the timeslice start time + /// \param TsLength Reference to the timeslice length + /// \param TsOverlap Reference to the timeslice overlap length (does not used at the moment) + /// \param FstHitinTs Index of the first hit in the time-lice + /// \param areDataLeft Flag: true - data were left after reading the sub-timeslice + /// \param event Pointer to the current CbmEvent object void ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, float& TsOverlap, int& FstHitinTs, bool& areDataLeft, CbmEvent* event = NULL); + /// Converts data from generic FairMCPoint based class to the CbmL1MCPoint (dummy method) + /// \param MC Pointer to a target CbmL1MCPoint object + /// \param iPoint Index of the point into the input MC points CbmMCDataArray object for the particular detector + /// \param MVD Index of the detector subsystem + /// \return flag: false - success, true - some errors occured bool ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int MVD); // help procedure + + /// Converts data from generic FairMCPoint based class to the CbmL1MCPoint + /// \param MC Pointer to a target CbmL1MCPoint object + /// \param iPoint Index of the point into the input MC points CbmMCDataArray object for the particular detector + /// \param file Index of the input file + /// \param event Index of the input event + /// \param MVD Index of the detector subsystem + /// \return flag: false - success, true - some errors occured + // TODO: Probably, we should replace input parameter MVD with the template (S.Zharko) bool ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int MVD); + // static bool compareZ(const int &a, const int &b ); // bool compareZ(const int &a, const int &b ); + + /// Fills the vMCTracks vector and the dFEI2vMCTracks set void Fill_vMCTracks(); - /// Input Performance - void HitMatch(); // Procedure for match hits and MCPoints. + /* + * Input Performance + */ + + /// Procedure for match hits and MCPoints. + /// Reads information about correspondence between hits and mcpoints and fill CbmL1MCPoint::hitIds and CbmL1Hit::mcPointIds arrays + /// should be called after fill of algo + void HitMatch(); + void FieldApproxCheck(); // Build histos with difference between Field map and approximated field void FieldIntegralCheck(); // Build 2D histo: dependence of the field integral on phi and theta void InputPerformance(); // Build histos about input data, like hit pulls, etc. diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx index 1f37ac3eac..176daae07a 100644 --- a/reco/L1/CbmL1ReadEvent.cxx +++ b/reco/L1/CbmL1ReadEvent.cxx @@ -67,28 +67,48 @@ static bool compareZ(const int& a, const int& b) return l1->Get_Z_vMCPoint(a) < l1->Get_Z_vMCPoint(b); } -struct TmpHit { // used for sort Hits before writing in the normal arrays - int iStripF, - iStripB; // indices of strips - int iStation; - int ExtIndex; // index of hit in the TClonesArray array ( negative for MVD ) - double u_front, u_back; // positions of strips - double x, y, z; // position of hit - double xmc, ymc, p; // mc-position of hit, momentum - double tx, ty; // mc-slope of mc point - double dx, dy, dxy; - double du, dv; - int iMC; // index of MCPoint in the vMCPoints array - double time = 0., dt = 1.e10; - int Det; +/// Local structure for a hit, containing both measured and MC information. +/// The structure is used to sort hits before writing them into L1 input arrays +/// +struct TmpHit { + int iStripF; ///< index of front strip + int iStripB; ///< index of back strip + int iStation; ///< index of station + int ExtIndex; ///< index of hit in the external TClonesArray array (NOTE: negative for MVD) + double u_front, u_back; ///< positions of strips + double x, y, z; ///< position of hit (Cortesian coordinates) + double xmc, ymc, p; ///< mc position of hit, total momentum + double tx, ty; ///< mc slopes of the mc point + double dx, dy, dxy; ///< hit position errors in Cortesian coordinates + double du, dv; ///< hit position errors in strips coordinates + int iMC; ///< index of MCPoint in the vMCPoints array + double time = 0.; ///< time of the hit + double dt = 1.e10; ///< time error of the hit + int Det; int id; int track; + /// Provides comparison of two hits. + /// If two hits belong to different stations, + /// the smallest hit belongs to the station with the smallest index. Otherwise, the smallest hit + /// has the smallest y coordinate + /// \param a Left hit + /// \param b Right hit + /// \return boolean: true - the left hit is smaller then the right one static bool Compare(const TmpHit& a, const TmpHit& b) { return (a.iStation < b.iStation) || ((a.iStation == b.iStation) && (a.y < b.y)); } + /// Creates a hit from the CbmL1MCPoint object + /// \param point constant reference to the input MC-point + /// \param det + /// \param nTmpHits + /// \param nStripF + /// \param ip + /// \param NStrips + /// \param st reference to the station info object + // TODO: Probably, L1Station& st parameter should be constant. Do we really want to modify a station here? (S.Zharko) void CreateHitFromPoint(const CbmL1MCPoint& point, int det, int nTmpHits, int nStripF, int ip, int& NStrips, L1Station& st) { @@ -130,6 +150,12 @@ struct TmpHit { // used for sort Hits before writing in the normal arrays iMC = ip; } + /// Sets randomized position and time of the hit + /// The positions are smeared within predefined errors dx, dy, dt; z coordinate + /// of the hit is known precisely + /// \param point constant reference to the input MC-point + /// \param st reference to the station info object + // TODO: Probably, L1Station& st parameter should be constant. Do we really want to modify a station here? (S.Zharko) void SetHitFromPoint(const CbmL1MCPoint& point, L1Station& st) { x = 0.5 * (point.xIn + point.xOut) + gRandom->Gaus(0, dx); @@ -141,7 +167,6 @@ struct TmpHit { // used for sort Hits before writing in the normal arrays } }; -/// Repack data from Clones Arrays to L1 arrays void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, float& /*TsOverlap*/, int& FstHitinTs, bool& areDataLeft, CbmEvent* event) { @@ -152,15 +177,15 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, fData_->Clear(); // clear arrays for next event - vMCPoints.clear(); - vMCPoints_in_Time_Slice.clear(); - vMCTracks.clear(); - vStsHits.clear(); - vRTracks.clear(); - vHitMCRef.clear(); - vHitStore.clear(); - dFEI2vMCPoints.clear(); - dFEI2vMCTracks.clear(); + vMCPoints.clear(); /* <CbmL1MCPoint> */ + vMCPoints_in_Time_Slice.clear(); /* <int> */ + vMCTracks.clear(); /* <CbmL1MCTrack> */ + vStsHits.clear(); /* <CbmL1Hit> */ + vRTracks.clear(); /* <CbmL1Track> */ + vHitMCRef.clear(); /* <int>: indexes of MC-points in vMCPoints (by index of algo->vStsHits) */ + vHitStore.clear(); /* <CbmL1HitStore> */ + dFEI2vMCPoints.clear(); /* dFEI vs MC-point index: dFEI = index * 10000 + fileID + eventID * 0.0001 */ + dFEI2vMCTracks.clear(); /* dFEI vs MC-track index: dFEI = index * 10000 + fileID + eventID * 0.0001 */ if (fVerbose >= 10) cout << "ReadEvent: clear is done." << endl; @@ -212,8 +237,18 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, int firstMuchPoint = 0; int firstTofPoint = 0; + /* + * MC hits and tracks gathering: START + * + * If performance is studied, i.e. fPerformanse > 0, MC-information is used for + * efficiencies calculation. Otherwise this + */ + if (fPerformance) { + // Fill vMCTracks and dFEI2vMCTracks Fill_vMCTracks(); + + // Fill vMCPoints, vMCPoints_in_Time_Slice and dFEI2vMCPoints vMCPoints.clear(); vMCPoints.reserve(5 * vMCTracks.size() * NStation); vMCPoints_in_Time_Slice.clear(); @@ -223,6 +258,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, Int_t iFile = set_it->first; Int_t iEvent = set_it->second; + // TODO: Probably to put this code in L1Algo interfaces (S.Zharko) if (fMvdPoints) { Int_t nMvdPointsInEvent = fMvdPoints->Size(iFile, iEvent); double maxDeviation = 0; @@ -230,7 +266,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, CbmL1MCPoint MC; if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 1)) { MC.iStation = -1; - L1Station* sta = algo->vStations; + L1Station* sta = algo->vStations; // TODO: Wrap it into interface algo->GetStations() (S.Zharko) double bestDist = 1.e20; for (Int_t iSt = 0; iSt < NMvdStations; iSt++) { // use z_in since z_out is sometimes very wrong @@ -447,10 +483,23 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, } //fPerformance + /* + * MC hits and tracks gathering: END + */ + + /* + * Measured hits gathering: START + * + * In this section the measured hits from different detector subsystems are reformatted according to TmpHit structure + * (NOTE: independent from particular detector design) and then pushed to the tmpHit vector. In the performance study mode + * matching with MC points is done + */ + int NStrips = 0; + // // get MVD hits - + // if (listMvdHits) { int firstDetStrip = NStrips; @@ -532,6 +581,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, } } + // + // Get STS hits + // if (listStsHits && (2 != fStsUseMcHit)) { Int_t nEntSts = (event ? event->GetNofData(ECbmDataType::kStsHit) : listStsHits->GetEntriesFast()); @@ -667,6 +719,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, } // for j } // if listStsHits + // + // Get MuCh hits + // if ((2 == fMuchUseMcHit) && fUseMUCH) { // SG!! create TRD hits from TRD points int firstDetStrip = NStrips; @@ -791,8 +846,10 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, nTrdHits++; } } - - + + // + // Get TRD hits + // if (fUseTRD && listTrdHits && (2 != fTrdUseMcHit)) { int firstDetStrip = NStrips; vector<bool> mcUsed(nTrdPoints, 0); @@ -912,6 +969,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, } // for j } // if listTrdHits + // + // Get ToF hits + // if ((2 == fTofUseMcHit) && fUseTOF) { // SG!! create TRD hits from TRD points int firstDetStrip = NStrips; @@ -1039,13 +1099,24 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, << nTrdHits << " tof " << nTofHits << endl; } - // sort hits + // total number of hits int nHits = nMvdHits + nStsHits + nMuchHits + nTrdHits + nTofHits; + /* + * Measured hits gathering: END + */ + /* + * Hits sorting + * + * Two hits are compared as follows. If the hits are measured with two different stations, the smallest hit has the smallest + * station ID. If the hits are measured within one station, the smallest hit has the smallest y position coordinate. + */ sort(tmpHits.begin(), tmpHits.end(), TmpHit::Compare); - // save strips in L1Algo + /* + * Save strips into L1Algo + */ fData_->NStsStrips = NStrips; fData_->fStripFlag.reset(NStrips, 0); int maxHitIndex = 0; @@ -1056,15 +1127,13 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, fData_->fStripFlag[th.iStripF] = flag; fData_->fStripFlag[th.iStripB] = flag; if (maxHitIndex < th.id) { maxHitIndex = th.id; } - } // ih if (fVerbose >= 10) { cout << "ReadEvent: strips are read." << endl; } /* - * Fill and then save vStsHits, vHitStore and vHitMCRef vectors as well as fData->vStsHits + * Fill and save vStsHits, vHitStore and vHitMCRef vectors as well as fData->vStsHits */ - int nEffHits = 0; SortedIndex.reset(maxHitIndex + 1); @@ -1160,20 +1229,19 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, fData_->GetStsHitsStopIndex()); if (fPerformance) { - if (fVerbose >= 10) cout << "HitMatch is done." << endl; if (fVerbose >= 10) cout << "MCPoints and MCTracks are saved." << endl; } - if (fVerbose >= 10) cout << "ReadEvent is done." << endl; - - } // void CbmL1::ReadEvent() - - +// +//-------------------------------------------------------------------------------------------------- +// void CbmL1::Fill_vMCTracks() { vMCTracks.clear(); + + // Count the total number of tracks in this event and reserve memory { Int_t nMCTracks = 0; for (DFSET::iterator set_it = vFileEvent.begin(); set_it != vFileEvent.end(); ++set_it) { @@ -1185,21 +1253,19 @@ void CbmL1::Fill_vMCTracks() } int fileEvent = 0; + /* Loop over MC event */ for (DFSET::iterator set_it = vFileEvent.begin(); set_it != vFileEvent.end(); ++set_it, ++fileEvent) { Int_t iFile = set_it->first; Int_t iEvent = set_it->second; auto header = dynamic_cast<FairMCEventHeader*>(fMcEventHeader->Get(iFile, iEvent)); assert(header); - if (fVerbose > 2) { - LOG(info) << "mc event vertex at " << header->GetX() << " " << header->GetY() << " " << header->GetZ(); - } Int_t nMCTrack = fMCTracks->Size(iFile, iEvent); - + /* Loop over MC tracks */ for (Int_t iMCTrack = 0; iMCTrack < nMCTrack; iMCTrack++) { CbmMCTrack* MCTrack = L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(iFile, iEvent, iMCTrack)); - if (!MCTrack) continue; + if (!MCTrack) { continue; } int mother_ID = MCTrack->GetMotherId(); @@ -1216,6 +1282,7 @@ void CbmL1::Fill_vMCTracks() mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); } + Int_t IND_Track = vMCTracks.size(); //or iMCTrack? CbmL1MCTrack T(mass, q, vr, vp, IND_Track, mother_ID, pdg); // CbmL1MCTrack T(mass, q, vr, vp, iMCTrack, mother_ID, pdg); @@ -1232,12 +1299,16 @@ void CbmL1::Fill_vMCTracks() } //Links } //Fill_vMCTracks - +// +//-------------------------------------------------------------------------------------------------- +// +// TODO: Probably, we should reduce code here, rewriting this funciton as a template from CbmMvdPoint (S.Zharko) bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int MVD) { TVector3 xyzI, PI, xyzO, PO; Int_t mcID = -1; Double_t time = 0.f; + if (MVD == 1) { CbmMvdPoint* pt = L1_DYNAMIC_CAST<CbmMvdPoint*>(fMvdPoints->Get(file, event, iPoint)); // file, event, object //CbmMvdPoint *pt = L1_DYNAMIC_CAST<CbmMvdPoint*> (Point); @@ -1388,12 +1459,13 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int M return 0; } - +// +//-------------------------------------------------------------------------------------------------- +// bool CbmL1::ReadMCPoint(CbmL1MCPoint* /*MC*/, int /*iPoint*/, int /*MVD*/) { return 0; } - -/// Procedure for match hits and MCPoints. -/// Read information about correspondence between hits and mcpoints and fill CbmL1MCPoint::hitIds and CbmL1Hit::mcPointIds arrays -/// should be called after fill of algo +// +//-------------------------------------------------------------------------------------------------- +// void CbmL1::HitMatch() { const int NHits = vStsHits.size(); diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h index 5fcbc99a02..d5037da614 100644 --- a/reco/L1/L1Algo/L1Algo.h +++ b/reco/L1/L1Algo/L1Algo.h @@ -19,14 +19,14 @@ class L1AlgoDraw; //#define XXX // time debug //#define COUNTERS // diff counters (hits, doublets, ... ) -#define MERGE_CLONES +// =====>dispatched<===== // #define MERGE_CLONES // #define TRACKS_FROM_TRIPLETS_ITERATION kAllPrimIter //#define HitErrors //#define GLOBAL //#define mCBM -#define LAST_ITERATION kAllSecIter +// =====>dispatched<===== // #define LAST_ITERATION kAllSecIter #define FIND_GAPED_TRACKS // use triplets with gaps // =====>dispatched<===== // #define USE_RL_TABLE #ifndef TRACKS_FROM_TRIPLETS @@ -34,7 +34,6 @@ class L1AlgoDraw; #endif //#define USE_EVENT_NUMBER //#endif -//#define MERGE_CLONES #include <array> @@ -152,21 +151,35 @@ public: Tindex fDupletPortionStopIndex[L1Parameters::kMaxNstations] {0}; // end of the duplet portions for the station L1Vector<Tindex> fDupletPortionSize {"L1Algo::fDupletPortionSize"}; // N duplets in a portion + /********************************************************************************************//** + * Temporary vectors used by the clone merger + * TODO: Probably, the subclass L1TrackMerger for clones merger would help to improve + * readability (S.Zharko) + ***********************************************************************************************/ + // - // Temporary vectors used by the clone merger + // Vectors that are parallel to fTracks // - // vectors that are parallel to fTracks - L1Vector<unsigned short> fMergerTrackFirstStation {"L1Algo::fMergerTrackFirstStation"}; // first station of a track - L1Vector<unsigned short> fMergerTrackLastStation {"L1Algo::fMergerTrackLastStation"}; // last station of a track - L1Vector<L1HitIndex_t> fMergerTrackFirstHit {"L1Algo::fMergerTrackFirstHit"}; // index of the first tracks hit - L1Vector<L1HitIndex_t> fMergerTrackLastHit {"L1Algo::fMergerTrackLastHit"}; // index of the last tracks hit - L1Vector<unsigned short> fMergerTrackNeighbour { - "L1Algo::fMergerTrackNeighbour"}; // track that can be merged with the given track - L1Vector<float> fMergerTrackChi2 {"L1Algo::fMergerTrackChi2"}; // chi2 of the merge - L1Vector<char> fMergerTrackIsStored {"L1Algo::fMergerTrackIsStored"}; // is the track already stored to the output - L1Vector<char> fMergerTrackIsDownstreamNeighbour { - "L1Algo::fMergerTrackIsDownstreamNeighbour"}; // is the track a downstream neighbor of another track - // other vectors + /// First station of a track + L1Vector<unsigned short> fMergerTrackFirstStation {"L1Algo::fMergerTrackFirstStation"}; + /// Last station of a track + L1Vector<unsigned short> fMergerTrackLastStation {"L1Algo::fMergerTrackLastStation"}; + /// Index of the first hit of a track + L1Vector<L1HitIndex_t> fMergerTrackFirstHit {"L1Algo::fMergerTrackFirstHit"}; + /// Index of the last hit of a track + L1Vector<L1HitIndex_t> fMergerTrackLastHit {"L1Algo::fMergerTrackLastHit"}; + /// Index (TODO:??) of a track that can be merge with the given track + L1Vector<unsigned short> fMergerTrackNeighbour { "L1Algo::fMergerTrackNeighbour"}; + /// Chi2 value of the track merging procedure + L1Vector<float> fMergerTrackChi2 {"L1Algo::fMergerTrackChi2"}; + /// Flag: is the given track already stored to the output + L1Vector<char> fMergerTrackIsStored {"L1Algo::fMergerTrackIsStored"}; + /// Flag: is the track a downstream neighbour of another track + L1Vector<char> fMergerTrackIsDownstreamNeighbour {"L1Algo::fMergerTrackIsDownstreamNeighbour"}; + // + // Utility vectors + // + /// Tracks after the merging procedure L1Vector<L1Track> fMergerTracksNew {"L1Algo::fMergerTracksNew"}; // vector of tracks after the merge L1Vector<L1HitIndex_t> fMergerRecoHitsNew {"L1Algo::fMergerRecoHitsNew"}; // vector of track hits after the merge @@ -426,6 +439,7 @@ private: void MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5]); void FilterTracks(fvec const r[5], fvec const C[15], fvec const m[5], fvec const V[15], fvec R[5], fvec W[15], fvec* chi2); + /// void CAMergeClones(); diff --git a/reco/L1/L1Algo/L1CAIteration.cxx b/reco/L1/L1Algo/L1CAIteration.cxx index 1d72350423..5dbf8e86ff 100644 --- a/reco/L1/L1Algo/L1CAIteration.cxx +++ b/reco/L1/L1Algo/L1CAIteration.cxx @@ -2,16 +2,23 @@ SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ +/***********************************************************************************************//** + * @file L1CAIteration.cxx + * @brief Definition of the L1CAIteration class methods + * @since 05.02.2022 + * @author S.Zharko <s.zharko@gsi.de> + ***************************************************************************************************/ + #include "L1CAIteration.h" #include <FairLogger.h> #include <sstream> - +// //---------------------------------------------------------------------------------------------------------------------- // L1CAIteration::L1CAIteration() noexcept { LOG(debug) << "L1CAIteration: Default constructor called for " << this; } - +// //---------------------------------------------------------------------------------------------------------------------- // L1CAIteration::L1CAIteration(const L1CAIteration& other) noexcept @@ -33,7 +40,7 @@ L1CAIteration::L1CAIteration(const L1CAIteration& other) noexcept { LOG(debug) << "L1CAIteration: Copy constructor called: " << &other << " was copied into " << this; } - +// //---------------------------------------------------------------------------------------------------------------------- // L1CAIteration::L1CAIteration(L1CAIteration&& other) noexcept @@ -41,18 +48,18 @@ L1CAIteration::L1CAIteration(L1CAIteration&& other) noexcept this->Swap(other); LOG(debug) << "L1CAIteration: Move constructor called: " << &other << " was moved into " << this; } - +// //---------------------------------------------------------------------------------------------------------------------- // L1CAIteration::L1CAIteration(const std::string& name) noexcept : fName(name) { LOG(debug) << "L1CAIteration: Constructor from name called for " << this; } - +// //---------------------------------------------------------------------------------------------------------------------- // L1CAIteration::~L1CAIteration() noexcept { LOG(debug) << "L1CAIteration: Destructor called for " << this; } - +// //---------------------------------------------------------------------------------------------------------------------- // L1CAIteration& L1CAIteration::operator=(const L1CAIteration& other) noexcept @@ -61,7 +68,7 @@ L1CAIteration& L1CAIteration::operator=(const L1CAIteration& other) noexcept LOG(debug) << "L1CAIteration: Copy operator= called: " << &other << " was copied into " << this; return *this; } - +// //---------------------------------------------------------------------------------------------------------------------- // L1CAIteration& L1CAIteration::operator=(L1CAIteration&& other) noexcept @@ -73,7 +80,7 @@ L1CAIteration& L1CAIteration::operator=(L1CAIteration&& other) noexcept LOG(debug) << "L1CAIteration: Move operator= called: " << &other << " was moved into " << this; return *this; } - +// //---------------------------------------------------------------------------------------------------------------------- // void L1CAIteration::Print(int verbosityLevel) const @@ -81,7 +88,7 @@ void L1CAIteration::Print(int verbosityLevel) const if (verbosityLevel == 0) { LOG(info) << " - " << fName; } if (verbosityLevel > 0) { LOG(info) << ToString(0); } } - +// //---------------------------------------------------------------------------------------------------------------------- // void L1CAIteration::SetTargetPosSigmaXY(float sigmaX, float sigmaY) @@ -89,7 +96,7 @@ void L1CAIteration::SetTargetPosSigmaXY(float sigmaX, float sigmaY) fTargetPosSigmaX = sigmaX; fTargetPosSigmaY = sigmaY; } - +// //---------------------------------------------------------------------------------------------------------------------- // void L1CAIteration::Swap(L1CAIteration& other) noexcept @@ -110,7 +117,7 @@ void L1CAIteration::Swap(L1CAIteration& other) noexcept std::swap(fTargetPosSigmaX, other.fTargetPosSigmaX); std::swap(fTargetPosSigmaY, other.fTargetPosSigmaY); } - +// //---------------------------------------------------------------------------------------------------------------------- // std::string L1CAIteration::ToString(int indentLevel) const diff --git a/reco/L1/L1Algo/L1CAIteration.h b/reco/L1/L1Algo/L1CAIteration.h index deed79ee74..8ff0f8a785 100644 --- a/reco/L1/L1Algo/L1CAIteration.h +++ b/reco/L1/L1Algo/L1CAIteration.h @@ -2,12 +2,12 @@ SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ - -/************************************************************************************************************ - * @file L1CAIteration.h - * @brief Class for L1 CA Track Finder Iteration - * @since 01.15.2022 - ***********************************************************************************************************/ +/***********************************************************************************************//** + * @file L1CAIteration.h + * @brief Declaration of the L1CAIteration class + * @since 05.02.2022 + * @author S.Zharko <s.zharko@gsi.de> + ***************************************************************************************************/ #ifndef L1CAIteration_h #define L1CAIteration_h 1 @@ -15,13 +15,18 @@ #include <bitset> #include <string> -/// Class L1CAIteration describes L1 Track finder iteration. Each iteration contains specific cuts and special -/// flags. +/// Class L1CAIteration describes L1 Track finder iteration. +/// Each iteration utilizes special physics cuts and run condition to find tracks of a particular +/// class (e.g., fast primary tracks or secondary electron tracks). Hits associated with tracks +/// reconstructed during current iteration are removed from the further iterations. /// class L1CAIteration { + /// Enumeration ControlFlag is used to keep flags, which controls behaviour of the track finder + /// iterations loop + /// enum class ControlFlag { - kePrimary, ///< true - track is primary, false - track is secondary (not primary) + kePrimary, ///> true - track is primary, false - track is secondary (not primary) keEnd }; using ControlFlags_t = std::bitset<static_cast<int>(ControlFlag::keEnd)>; @@ -37,26 +42,26 @@ public: L1CAIteration(const std::string& name) noexcept; /// Destructor ~L1CAIteration() noexcept; - /// Copy operator= + /// Copy assignment operator L1CAIteration& operator=(const L1CAIteration& other) noexcept; - /// Move operator= + /// Move assignment operator L1CAIteration& operator=(L1CAIteration&& other) noexcept; /// Gets doublet chi2 upper cut float GetDoubletChi2Cut() const { return fDoubletChi2Cut; } - /// + /// Gets correction for accounting overlaping and iff z float GetMaxDZ() const { return fMaxDZ; } /// Gets max considered q/p for tracks float GetMaxInvMom() const { return fMaxInvMom; } - /// Gets max slope (tx\ty) in 3d hit position of a triplet + /// Gets max slope (tx\ty) in 3D hit position of a triplet float GetMaxSlope() const { return fMaxSlope; } - /// Gets max slope (tx\ty) in prim vertex + /// Gets max slope (tx\ty) in primary vertex float GetMaxSlopePV() const { return fMaxSlopePV; } /// Gets the name of the iteration const std::string& GetName() const { return fName; } - /// + /// Gets size of region [TODO: units??] to attach new hits to the created track float GetPickGather() const { return fPickGather; } - /// + /// Gets min value of dp/dp_error, for which two tiplets are neighbours float GetPickNeighbour() const { return fPickNeighbour; } /// Gets sigma target position in X direction [cm] float GetTargetPosSigmaX() const { return fTargetPosSigmaX; } @@ -75,23 +80,25 @@ public: /// Sets doublet chi2 upper cut void SetDoubletChi2Cut(float input) { fDoubletChi2Cut = input; } - /// + /// Sets correction for accounting overlaping and iff z void SetMaxDZ(float input) { fMaxDZ = input; } /// Sets max considered q/p for tracks void SetMaxInvMom(float input) { fMaxInvMom = input; } - /// Sets max slope (tx\ty) in 3d hit position of a triplet + /// Sets max slope (tx\ty) in 3D hit position of a triplet void SetMaxSlope(float input) { fMaxSlope = input; } - /// Sets max slope (tx\ty) in prim vertex + /// Sets max slope (tx\ty) in primary vertex void SetMaxSlopePV(float input) { fMaxSlopePV = input; } /// Sets name of the iteration void SetName(const std::string& name) { fName = name; } - /// + /// Sets size of region [TODO: units??] to attach new hits to the created track void SetPickGather(float input) { fPickGather = input; } - /// + /// Sets min value of dp/dp_error, for which two tiplets are neighbours void SetPickNeighbour(float input) { fPickNeighbour = input; } /// Sets flag: primary tracks - true, secondary tracks - false void SetPrimary(bool flag) { fControlFlags[static_cast<int>(ControlFlag::kePrimary)] = flag; } - /// Sets sigma of target positions in XY plane (in cm) + /// Sets sigma of target positions in XY plane + /// \param sigmaX Sigma value in X direction [cm] + /// \param sigmaX Sigma value in Y direction [cm] void SetTargetPosSigmaXY(float sigmaX, float sigmaY); /// Sets track chi2 upper cut void SetTrackChi2Cut(float input) { fTrackChi2Cut = input; } @@ -101,27 +108,28 @@ public: /// Swap method void Swap(L1CAIteration& other) noexcept; /// String representation of the class contents + /// \param indentLevel Level of indentation for the text (in terms of \t symbols) std::string ToString(int indentLevel = 0) const; private: - //------------------------------------------------------------------------------------------------------------------- + /** Basic fields **/ std::string fName {""}; ///< Iteration name ControlFlags_t fControlFlags {}; ///< bitset flags to control iteration behaviour - //------------------------------------------------------------------------------------------------------------------- - // Track finder dependent cuts - //------------------------------------------------------------------------------------------------------------------- - + /** Track finder dependent cuts **/ // TODO: Iteratively change the literals to floats (S.Zharko) + // NOTE: For each new cut one should not forget to create a setter and a getter, insert the value + // initialization in the copy constructor and the Swap operator as well as a string repre- + // sentation to the ToString method (S.Zharko) float fTrackChi2Cut {10.f}; ///> track chi2 upper cut float fTripletChi2Cut {21.1075f}; ///> triplet chi2 upper cut float fDoubletChi2Cut {11.3449 * 2.f / 3.f}; ///> doublet chi2 upper cut - float fPickGather {3.0}; ///> size of region to attach new hits to the created track - float fPickNeighbour {5.0}; ///> min value of dp/dp_error, for which to tiplets are neighbours - float fMaxInvMom {1.0 / 0.5}; ///> max considered q/p for tracks - float fMaxSlopePV {1.1}; ///> max slope (tx\ty) in prim vertex - float fMaxSlope {2.748}; ///> max slope (tx\ty) in 3d hit position of a triplet - float fMaxDZ {0.f}; ///> Correction for accounting overlaping and iff z + float fPickGather {3.0}; ///> size of region to attach new hits to the created track [TODO: units??] + float fPickNeighbour {5.0}; ///> min value of dp/dp_error, for which two tiplets are neighbours + float fMaxInvMom {1.0 / 0.5}; ///> max considered q/p for tracks + float fMaxSlopePV {1.1}; ///> max slope (tx\ty) in primary vertex + float fMaxSlope {2.748}; ///> max slope (tx\ty) in 3D hit position of a triplet + float fMaxDZ {0.f}; ///> Correction for accounting overlaping and iff z [TODO: units??] float fTargetPosSigmaX {0}; ///> Constraint on target position in X direction [cm] float fTargetPosSigmaY {0}; ///> Constraint on target position in Y direction [cm] }; diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx index aecf7d6105..cac3f853b5 100644 --- a/reco/L1/L1Algo/L1CATrackFinder.cxx +++ b/reco/L1/L1Algo/L1CATrackFinder.cxx @@ -700,9 +700,9 @@ inline void L1Algo::f30( // input if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) // extrapolate to the right hit station { - if (istar <= fNfieldStations) L1Extrapolate(T2, star.z, T2.qp, f2); + if (istar <= fNfieldStations) L1Extrapolate(T2, star.z, T2.qp, f2); // Full extrapolation in the magnetic field else - L1ExtrapolateLine(T2, star.z); + L1ExtrapolateLine(T2, star.z); // Extrapolation with line () } else L1Extrapolate(T2, star.z, T2.qp, f2); @@ -715,7 +715,7 @@ inline void L1Algo::f30( // input if (fabs(T2.tx[i2_4]) > fMaxSlope) continue; if (fabs(T2.ty[i2_4]) > fMaxSlope) continue; - const fvec Pick_r22 = (fTripletChi2Cut - T2.chi2); + const fvec Pick_r22 = (fTripletChi2Cut - T2.chi2); const float& timeError = T2.C55[i2_4]; const float& time = T2.t[i2_4]; // find first possible hit @@ -746,7 +746,7 @@ inline void L1Algo::f30( // input //for (int irh = 0; irh < ( StsHitsUnusedStopIndex[istar] - StsHitsUnusedStartIndex[istar] ); irh++){ const L1HitPoint& hitr = vStsHits_r[irh]; -#ifdef USE_EVENT_NUMBER +#ifdef USE_EVENT_NUMBER // NOTE: if ((T2.n[i2_4] != hitr.n)) continue; #endif // USE_EVENT_NUMBER const fscal zr = hitr.Z(); @@ -2519,9 +2519,9 @@ void L1Algo::CATrackFinder() c_timerG.Start(); #endif -#ifdef MERGE_CLONES - CAMergeClones(); -#endif + if constexpr (L1Parameters::kIfMergeClones) { + CAMergeClones(); + } #ifdef XXX c_timerG.Stop(); diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h index c779f7436d..4ef0102e1b 100644 --- a/reco/L1/L1Algo/L1Parameters.h +++ b/reco/L1/L1Algo/L1Parameters.h @@ -20,9 +20,9 @@ class L1Parameters { public: - ///////////////////////////// - // COMPILE TIME CONSTANTS // - /////////////////////////// + /**********************//** + * COMPILE TIME CONSTANTS * + **************************/ // // Array sizes @@ -50,11 +50,25 @@ public: // Compile control flags // - /// Selector for the radiation length tables usage + /// Flag for the radiation length tables usage /// true - Radiational tables will be used, /// false - basic station material info is used static constexpr bool kIfUseRadLengthTable {true}; + /// Flag for calling the CAMergeClones procedure ... TODO + static constexpr bool kIfMergeClones {true}; + + /// Flag: debug mode for analyzing the doublets pergormance efficiencies + static constexpr bool kIfDebugDoubletsPerformance {false}; + /// Flag: debug mode for analyzing the tiplets pergormance efficiencies + static constexpr bool kIfDebugTripletsPerformance {false}; + /// Flag: debug mode for creating pools for triplets. + /// NOTE: this feature will work only if the L1Parameters::kIfDebugTipletsPerformace is true! + static constexpr bool kIfCreateTipletPulls {false}; + + + + public: /// Default constructor L1Parameters() = default; diff --git a/reco/L1/L1Algo/L1Track.h b/reco/L1/L1Algo/L1Track.h index 1cd1ea6216..ff53a42582 100644 --- a/reco/L1/L1Algo/L1Track.h +++ b/reco/L1/L1Algo/L1Track.h @@ -23,34 +23,50 @@ #include <limits> +/// Class representing a track in the L1 tracking algorithm +/// Track parameters vector: {x, y, tx, ty, q/p, z, t} +/// Covariation matrix: C[20] (C55) corresponds to the time variance +/// class L1Track { public: static constexpr float kNaN {std::numeric_limits<float>::signaling_NaN()}; - unsigned char NHits; - unsigned char n; - float Momentum {kNaN}, fTrackTime {kNaN}; - fscal TFirst[7] {kNaN}, CFirst[21] {kNaN}, TLast[7] {kNaN}, CLast[21] {kNaN}, Tpv[7] {kNaN}, Cpv[21] {kNaN}, - chi2 {kNaN}; - short int NDF; - - int FirstHitIndex, LastHitIndex; - int index; - int ista; + unsigned char NHits; ///> Number of hits in track + unsigned char n; ///> TODO: ?? + float Momentum {kNaN}; ///> TODO: ?? + float fTrackTime {kNaN}; ///> Track time + fscal TFirst[7] {kNaN}; ///> Track parameters on the first station + fscal CFirst[21] {kNaN}; ///> Track parameter covariation matrix elemenst on the first station + fscal TLast[7] {kNaN}; ///> Track parameters on the last station + fscal CLast[21] {kNaN}; ///> Track parameter covatiation matrix elements on the second station + fscal Tpv[7] {kNaN}; ///> Track parameters in the primary vertex + fscal Cpv[21] {kNaN}; ///> Track parameter covariation matrix elements in the primary vertex + fscal chi2 {kNaN}; ///> Track fit chi-square value + short int NDF; ///> Track fit NDF value + // TODO: shouldn't it be the L1HitIndex_t type instead of int? (S.Zharko) + int FirstHitIndex; ///> Track first hit index in the hits vector + int LastHitIndex; ///> Track last hit index in the hits vector + int index; ///> Index of the track + int ista; ///> TODO: ?? + /// Provides comparison of two L1Track objects + /// If two tracks have different number of hits, the smallest track has the largest number of hits. + /// If two tracks have the same numbers of hits and ... static bool compareCand(const L1Track& a, const L1Track& b) { - - if (a.NHits != b.NHits) return (a.NHits > b.NHits); - - if (a.ista != b.ista) return (a.ista < b.ista); - - else + if (a.NHits != b.NHits) { + return (a.NHits > b.NHits); + } + if (a.ista != b.ista) { + return (a.ista < b.ista); + } + else { return (a.chi2 < b.chi2); + } } - + /// Provides comparison for two tracks by the time variance static bool compare(const L1Track& a, const L1Track& b) { return (a.Cpv[20] <= b.Cpv[20]); } }; @@ -61,4 +77,4 @@ public: // else return false; // } -#endif +#endif // L1Track_H -- GitLab