diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index bbb0e86c29cf3228afcfe45548c32b52a69bfc9d..52aa3bb0775addc455bcd47cc01e175a23d426f0 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -71,7 +71,8 @@ set(SRCS L1Algo/utils/L1AlgoPulls.cxx catools/CaToolsMCData.cxx catools/CaToolsMCPoint.cxx - catools/CaToolsPerformance.cxx + catools/CaToolsMCTrack.cxx + catools/CaToolsQa.cxx catools/CaToolsDebugger.cxx catools/CaToolsWindowFinder.cxx catools/CaToolsWFExpression.cxx diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx index d2e0e038764c4cfd04b566b9e1059a99b3d7cc0e..e76e9e51a8432f01316e63b95c2fabcbb25708d6 100644 --- a/reco/L1/CbmCaMCModule.cxx +++ b/reco/L1/CbmCaMCModule.cxx @@ -11,6 +11,7 @@ #include "CbmEvent.h" #include "CbmL1.h" // for L1DetectorID +#include "CbmL1Constants.h" #include "CbmL1Hit.h" #include "CbmMCDataManager.h" #include "CbmMCDataObject.h" @@ -154,6 +155,55 @@ void CbmCaMCModule::InitEvent(CbmEvent* /*pEvent*/) fMCData.Clear(); this->ReadMCTracks(); this->ReadMCPoints(); + + // ------ Prepare tracks: set point indexes and redefine indexes from external to internal containers + for (auto& aTrk : fMCData.GetTrackContainer()) { + aTrk.SortPointIndexes( + [&](const int& iPl, const int& iPr) { return fMCData.GetPoint(iPl).GetZ() < fMCData.GetPoint(iPr).GetZ(); }); + } +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaMCModule::InitTrackInfo(const L1Vector<CbmL1Hit>& vHits) +{ + // ----- Initialize stations arrangement and hit indexes + fMCData.InitTrackInfo(vHits); + + // ----- Define reconstructable and additional flags + for (auto& aTrk : fMCData.GetTrackContainer()) { + bool isRec = true; // is track reconstructable + + // Cut on momentum + isRec &= aTrk.GetP() > CbmL1Constants::MinRecoMom; + + // Cut on max number of points on station + isRec &= aTrk.GetMaxNofPointsOnStation() <= 4; + + bool isAdd = isRec; // is track additional + + // Cut on number of stations + switch (fPerformanceMode) { + case 1: isRec &= aTrk.GetNofConsStationsWithHit() >= CbmL1Constants::MinNStations; break; + case 2: isRec &= aTrk.GetTotNofStationsWithHit() >= CbmL1Constants::MinNStations; break; + case 3: isRec &= aTrk.GetNofConsStationsWithPoint() >= CbmL1Constants::MinNStations; break; + default: LOG(fatal) << "CA MC Module: invalid performance mode " << fPerformanceMode; + } + + if (aTrk.GetNofPoints() > 0) { + isAdd &= aTrk.GetNofConsStationsWithHit() == aTrk.GetTotNofStationsWithHit(); + isAdd &= aTrk.GetNofConsStationsWithPoint() == aTrk.GetTotNofStationsWithHit(); + isAdd &= aTrk.GetTotNofStationsWithHit() == aTrk.GetTotNofStationsWithPoint(); + isAdd &= aTrk.GetNofConsStationsWithHit() >= 3; + isAdd &= fMCData.GetPoint(aTrk.GetPointIndexes()[0]).GetStationId() == 0; + isAdd &= !isRec; + } + else { + isAdd = false; + } + aTrk.SetFlagReconstructable(isRec); + aTrk.SetFlagAdditional(isAdd); + } } // --------------------------------------------------------------------------------------------------------------------- @@ -161,17 +211,22 @@ void CbmCaMCModule::InitEvent(CbmEvent* /*pEvent*/) void CbmCaMCModule::Finish() { std::cout << "\033[1;32mFinishing performance\033[0m\n"; } -// ***************************************** -// ** Hits and MC-points matching ** -// ***************************************** +// ********************************** +// ** Reco and MC matching ** +// ********************************** // --------------------------------------------------------------------------------------------------------------------- // -void CbmCaMCModule::MatchPointsWithHits(L1Vector<CbmL1Hit>& vHits) +void CbmCaMCModule::MatchPointsWithHits(L1Vector<CbmL1Hit>& vInputExtHits) { - int nHits = vHits.size(); + // FIXME: Cleaning up makes it safer, but probably is not needed after old code removal + for (auto& hit : vInputExtHits) { + hit.mcPointIds.clear(); + } + + int nHits = vInputExtHits.size(); for (int iH = 0; iH < nHits; ++iH) { - auto& hit = vHits[iH]; + auto& hit = vInputExtHits[iH]; int iP = fMCData.GetPointIndexOfHit(iH); if (iP > -1) { hit.mcPointIds.push_back_no_warning(iP); @@ -245,7 +300,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kSts>(int iHit) const if (link.GetWeight() > bestWeight) { bestWeight = link.GetWeight(); - iPoint = fMCData.FindGlobalPointIndex(CalcGlobMCPointIndex(iIndex, L1DetectorID::kSts), iEvent, iFile); + iPoint = fMCData.FindInternalPointIndex(CalcGlobMCPointIndex(iIndex, L1DetectorID::kSts), iEvent, iFile); assert(iPoint != -1); } } @@ -270,7 +325,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMuch>(int iHit) const } int iFile = hitMatchMuch->GetLink(iLink).GetFile(); int iEvent = hitMatchMuch->GetLink(iLink).GetEntry(); - iPoint = fMCData.FindGlobalPointIndex(CalcGlobMCPointIndex(iMc, L1DetectorID::kMuch), iEvent, iFile); + iPoint = fMCData.FindInternalPointIndex(CalcGlobMCPointIndex(iMc, L1DetectorID::kMuch), iEvent, iFile); assert(iPoint != -1); } } @@ -296,7 +351,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHit) const int iFile = hitMatch->GetLink(0).GetFile(); int iEvent = hitMatch->GetLink(0).GetEntry(); - iPoint = fMCData.FindGlobalPointIndex(CalcGlobMCPointIndex(iMc, L1DetectorID::kTrd), iEvent, iFile); + iPoint = fMCData.FindInternalPointIndex(CalcGlobMCPointIndex(iMc, L1DetectorID::kTrd), iEvent, iFile); assert(iPoint != -1); } } @@ -318,7 +373,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHit) const if (iMc < 0) return iPoint; assert(iMc >= 0 && iMc < fpTofPoints->Size(iFile, iEvent)); - int iPointPrelim = fMCData.FindGlobalPointIndex(CalcGlobMCPointIndex(iMc, L1DetectorID::kTof), iEvent, iFile); + int iPointPrelim = fMCData.FindInternalPointIndex(CalcGlobMCPointIndex(iMc, L1DetectorID::kTof), iEvent, iFile); if (iPointPrelim == -1) { continue; } iPoint = iPointPrelim; } @@ -326,6 +381,63 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHit) const return iPoint; } +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaMCModule::MatchRecoAndMCTracks(L1Vector<CbmL1Track>& vRecoTracks, L1Vector<CbmL1Hit>& vInputExtHits) +{ + // NOTE: Potentially there is a case, when + + for (int iTre = 0; iTre < int(vRecoTracks.size()); ++iTre) { + auto& trkRe = vRecoTracks[iTre]; + + // ----- Count number of hits from each MC track belonging to this reconstructed track + auto& mNofHitsVsMCTrkID = trkRe.hitMap; + mNofHitsVsMCTrkID.clear(); + for (int iH : trkRe.Hits) { + for (int iP : vInputExtHits[iH].mcPointIds) { + assert(iP > -1); // Should never be triggered + int iTmc = fMCData.GetPoint(iP).GetTrackId(); + if (mNofHitsVsMCTrkID.find(iTmc) == mNofHitsVsMCTrkID.end()) { mNofHitsVsMCTrkID[iTmc] = 1; } + else { + mNofHitsVsMCTrkID[iTmc] += 1; + } + + } // loop over global point ids stored for hit: end + } // loop over hit ids stored for a reconstructed track: end + + + // ----- Calculate track max purity + // NOTE: Maximal purity is a maximum fraction of hits from a single MC track. A reconstructed track can be matched + // to several MC tracks, because it uses hits from different particle. Purity shows, how fully the true track + // was reconstructed. + int nHitsTrkRe = trkRe.Hits.size(); // number of hits in a given reconstructed track + double maxPurity = 0.; // [0, 1] + for (auto& item : mNofHitsVsMCTrkID) { + int iTmc = item.first; + int nHitsTrkMc = item.second; + + if (iTmc < 0) { continue; } + + if (double(nHitsTrkMc) > double(nHitsTrkRe) * maxPurity) { maxPurity = double(nHitsTrkMc) / double(nHitsTrkRe); } + + ca::tools::MCTrack& trkMc = fMCData.GetTrack(iTmc); + + // Match reconstructed and MC tracks, if purity with this MC track passes the threshold + if (double(nHitsTrkMc) >= CbmL1Constants::MinPurity * double(nHitsTrkRe)) { + trkMc.AddRecoTrackIndex(iTre); + trkRe.AddMCTrackIndex(iTmc); + } + // If purity is low, but some hits of a given MC track are presented in the reconstructed track + else { + trkMc.AddTouchTrackIndex(iTre); + } + + // Update max purity of the reconstructed track + trkRe.SetMaxPurity(maxPurity); + } // loop over hit map: end + } // loop over reconstructed tracks: end +} + // *********************** // ** Accessors ** @@ -465,7 +577,7 @@ void CbmCaMCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray* } } - int iTrack = fMCData.FindGlobalTrackIndex(pExtPoint->GetTrackID(), iEvent, iFile); + int iTrack = fMCData.FindInternalTrackIndex(pExtPoint->GetTrackID(), iEvent, iFile); assert(iTrack > -1); if (iStSelected != -1) { int iStActive = fpParameters->GetStationIndexActive(iStSelected, L1DetectorID::kTof); @@ -484,9 +596,11 @@ void CbmCaMCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray* if (vTofPointToTrack[iStLocGeo][iTrk] < 0) { continue; } ca::tools::MCPoint aPoint; if (FillMCPoint<L1DetectorID::kTof>(vTofPointToTrack[iStLocGeo][iTrk], iEvent, iFile, aPoint)) { - aPoint.SetStationID(iStActive); - int iGlobP = this->CalcGlobMCPointIndex(vTofPointToTrack[iStLocGeo][iTrk], L1DetectorID::kTof); - fMCData.AddPoint(aPoint, iGlobP, iEvent, iFile); + aPoint.SetStationId(iStActive); + aPoint.SetExternalId(CalcGlobMCPointIndex(vTofPointToTrack[iStLocGeo][iTrk], L1DetectorID::kTof)); + int iPInt = fMCData.GetNofPoints(); // assign an index of point in internal container + if (aPoint.GetTrackId() > -1) { fMCData.GetTrack(aPoint.GetTrackId()).AddPointIndex(iPInt); } + fMCData.AddPoint(aPoint); ++fvNofPointsUsed[int(L1DetectorID::kTof)]; } } // loop over tracks: end @@ -550,31 +664,51 @@ void CbmCaMCModule::ReadMCTracks() LOG(warn) << "CbmCaMCModule: track with (iF, iE, iT) = " << iFile << ' ' << iEvent << ' ' << iTrkExt << " not found"; } - TVector3 vtx; - TLorentzVector mom4; - pExtMCTrk->GetStartVertex(vtx); - pExtMCTrk->Get4Momentum(mom4); - int pdg = pExtMCTrk->GetPdgCode(); - unsigned procID = pExtMCTrk->GetGeantProcessId(); - int motherID = pExtMCTrk->GetMotherId(); - double charge = 0.; // in [e] - double mass = 0.; // in [GeV/c2] - - // Initialize charge and mass from PDG - auto pPdgBase = TDatabasePDG::Instance()->GetParticle(pdg); - if (pPdgBase) { - charge = pPdgBase->Charge() / 3.; - mass = pPdgBase->Mass(); + // Create a CbmL1MCTrack + ca::tools::MCTrack aTrk; + + aTrk.SetId(fMCData.GetNofTracks()); // assign current number of tracks read so far as an ID of a new track + aTrk.SetExternalId(iTrkExt); // external index of track is its index from CbmMCTrack objects container + aTrk.SetEventId(iEvent); + aTrk.SetFileId(iFile); + + aTrk.SetStartT(pExtMCTrk->GetStartT()); + aTrk.SetStartX(pExtMCTrk->GetStartX()); + aTrk.SetStartY(pExtMCTrk->GetStartY()); + aTrk.SetStartZ(pExtMCTrk->GetStartZ()); + + aTrk.SetPx(pExtMCTrk->GetPx()); + aTrk.SetPy(pExtMCTrk->GetPy()); + aTrk.SetPz(pExtMCTrk->GetPz()); + + aTrk.SetFlagSignal(aTrk.IsPrimary() && (aTrk.GetStartZ() > (pEvtHeader->GetZ() + 1.e-10))); + + // In CbmMCTrack mass is defined from ROOT PDG data base. If track is left from ion, and its pdg is not registered + // in the data base, its mass is calculated as A times proton mass. + aTrk.SetMass(pExtMCTrk->GetMass()); + + // The charge in CbmMCTrack is similarly to mass defined from ROOT PDG data base. The value of charge there is + // given in the units of 1/3e (absolute value of d-quark charge). In ca::tools::MCTrack we recalculate this + // value to the units of e. + aTrk.SetCharge(pExtMCTrk->GetCharge() / 3.); + + // Set index of mother track. We assume, that mother track was recorded into the internal array before + int extMotherId = pExtMCTrk->GetMotherId(); + if (extMotherId < 0) { + // This is a primary track, and it's mother ID equals -1 or -2. This value is taken also as internal mother ID + aTrk.SetMotherId(extMotherId); + } + else { + // This is a secondary track, mother ID should be recalculated for the internal array + int motherId = fMCData.FindInternalTrackIndex(extMotherId, iEvent, iFile); + assert(motherId > -1); + aTrk.SetMotherId(motherId); } - // Create a CbmL1MCTrack - int iTrkInt = fMCData.GetNofTracks(); // index of track in the MCData container - CbmL1MCTrack aTrk(mass, charge, vtx, mom4, iTrkInt, motherID, pdg, procID); - aTrk.time = pExtMCTrk->GetStartT(); - aTrk.iFile = iFile; - aTrk.iEvent = iEvent; - aTrk.isSignal = aTrk.IsPrimary() && (aTrk.z > pEvtHeader->GetZ() + 1.e-10); - fMCData.AddTrack(aTrk, iTrkExt, iEvent, iFile); + aTrk.SetProcessId(pExtMCTrk->GetGeantProcessId()); + aTrk.SetPdgCode(pExtMCTrk->GetPdgCode()); + + fMCData.AddTrack(aTrk); } // Loop over MC tracks: end } // Loop over MC events: end } diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index 66926522b6242d1bcc763c5b5b58c11b1cb092b8..2a66c06dafb86fe78eaee50be50ab2ce2cc5d5be 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -1009,6 +1009,7 @@ void CbmL1::Reconstruct(CbmEvent* event) if (fPerformance) { if (fVerbose > 1) { cout << "Performance..." << endl; } TrackMatch(); + fpMCModule->MatchRecoAndMCTracks(fvRecoTracks, fvExternalHits); EfficienciesPerformance(); HistoPerformance(); TrackFitPerformance(); diff --git a/reco/L1/CbmL1Hit.h b/reco/L1/CbmL1Hit.h index 94f3facaa9bef5bbe01cb464d2d15fbfd37599c0..02007d3c923bd2badfa392802f84ab5b95a0069b 100644 --- a/reco/L1/CbmL1Hit.h +++ b/reco/L1/CbmL1Hit.h @@ -25,16 +25,16 @@ public: /// class CbmL1HitDebugInfo { // TODO: SZh 21.09.2022: Replace instances of this class with L1Hit public: - int ExtIndex; - int iStation; - double x; - double y; - double time; - double dx; - double dy; - double dt; - double dxy; - int Det; + int ExtIndex; ///< index of hit in the external branch + int iStation; ///< index of station in active stations array + double x; ///< x coordinate of position [cm] + double y; ///< y coordinate of position [cm] + double time; ///< hit time [ns] + double dx; ///< x coordinate error [cm] + double dy; ///< y coordinate error [cm] + double dt; ///< time error [ns] + double dxy; ///< covariance between x and y [cm2] + int Det; ///< detector subsystem ID L1Vector<int> mcPointIds {"CbmL1HitDebugInfo::mcPointIds"}; ///< indices of CbmL1MCPoint in L1->fvMCPoints array }; diff --git a/reco/L1/CbmL1MCTrack.h b/reco/L1/CbmL1MCTrack.h index cd11625bd325f91a6920507e38f9f51d84c005ed..530fa73d2f10ac0c8e0cde4a0f5db4509e62db90 100644 --- a/reco/L1/CbmL1MCTrack.h +++ b/reco/L1/CbmL1MCTrack.h @@ -53,14 +53,18 @@ public: void Init(); void AddRecoTrack(CbmL1Track* rTr) { rTracks.push_back_no_warning(rTr); } + void AddRecoTrackIndex(int iT) { rTrackIndexes.push_back_no_warning(iT); } L1Vector<CbmL1Track*>& GetRecoTracks() { return rTracks; } int GetNClones() const { return rTracks.size() - 1; } bool IsReconstructed() const { return rTracks.size(); } void AddTouchTrack(CbmL1Track* tTr) { tTracks.push_back_no_warning(tTr); } + void AddTouchTrackIndex(int iT) { tTrackIndexes.push_back_no_warning(iT); } bool IsDisturbed() const { return tTracks.size(); } void SetIsReconstructable(bool v) { isReconstructable = v; } + const auto& GetRecoTrackIndexes() const { return rTrackIndexes; } + const auto& GetTouchTrackIndexes() const { return tTrackIndexes; } friend class CbmL1; @@ -108,10 +112,14 @@ private: bool isAdditional = false; // is not reconstructable, but stil interesting // next members filled and used in Performance - L1Vector<CbmL1Track*> rTracks {"CbmL1MCTrack::rTracks"}; // array of assosiated recoTracks + L1Vector<CbmL1Track*> rTracks {"CbmL1MCTrack::rTracks"}; // array of associated recoTracks L1Vector<CbmL1Track*> tTracks {"CbmL1MCTrack::tTracks"}; // array of recoTracks - // wich aren't assosiated with this mcTrack, + // which aren't associated with this mcTrack, // but use some hits from it. + + // NOTE: SZh 14.12.2022: on the replacement from rTracks and tTracks + L1Vector<int> rTrackIndexes = {"CbmL1MCTrack::rTrackIndexes"}; // array of associated recoTrack indexes + L1Vector<int> tTrackIndexes = {"CbmL1MCTrack::tTrackIndexes"}; // ..... }; diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx index 777b344ae2c1f9568051eec25be342b1f7bfb1e8..38432f6307aadf3c7163a0963736958eea81e99e 100644 --- a/reco/L1/CbmL1Performance.cxx +++ b/reco/L1/CbmL1Performance.cxx @@ -72,6 +72,13 @@ void CbmL1::TrackMatch() map<int, CbmL1MCTrack*> pMCTrackMap; pMCTrackMap.clear(); + // ***** DEBUG: Start + for (int iTrk = 0; iTrk < int(fvMCTracks.size()); ++iTrk) { + assert(iTrk == fvMCTracks[iTrk].ID); + } + + // ***** DEBUG: End + // fill pMCTrackMap for (vector<CbmL1MCTrack>::iterator i = fvMCTracks.begin(); i != fvMCTracks.end(); ++i) { CbmL1MCTrack& MC = *i; @@ -128,15 +135,17 @@ void CbmL1::TrackMatch() if (double(posIt->second) > max_percent * double(hitsum)) max_percent = double(posIt->second) / double(hitsum); // set relation to the mcTrack - if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) continue; + if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) { continue; } CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first]; if (double(posIt->second) >= CbmL1Constants::MinPurity * double(hitsum)) { // found correspondent MCTrack pmtra->AddRecoTrack(prtra); + pmtra->AddRecoTrackIndex(iR); prtra->AddMCTrack(pmtra); } else { pmtra->AddTouchTrack(prtra); + pmtra->AddTouchTrackIndex(iR); } } // for hitmap prtra->SetMaxPurity(max_percent); diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx index 38103534caa06be1530dfaed2f2aa4ff97c5b239..5b7b083f9c017db2345418d7d05103ab6c40c55c 100644 --- a/reco/L1/CbmL1ReadEvent.cxx +++ b/reco/L1/CbmL1ReadEvent.cxx @@ -414,7 +414,6 @@ void CbmL1::ReadEvent(CbmEvent* event) MC.ID = itTrack->second; fvMCTracks[MC.ID].Points.push_back_no_warning(fvMCPoints.size()); - fmMCPointsLinksMap[CbmL1LinkKey(iMC, iEvent, iFile)] = fvMCPoints.size(); fvMCPoints.push_back(MC); fNpointsMvd++; @@ -609,6 +608,7 @@ void CbmL1::ReadEvent(CbmEvent* event) } } //time_slice + for (unsigned int iTr = 0; iTr < fvMCTracks.size(); iTr++) { sort(fvMCTracks[iTr].Points.begin(), fvMCTracks[iTr].Points.end(), compareMcPointZ); @@ -1269,12 +1269,14 @@ void CbmL1::Fill_vMCTracks() Int_t pdg = MCTrack->GetPdgCode(); unsigned int processID = MCTrack->GetGeantProcessId(); Double_t q = 0, mass = 0.; - if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) { - q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0; - mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); - } + //if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) { + // q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0; + // mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); + //} // TODO: Add light nuclei (d, t, He3, He4): they are common in tracking but not accounted in TDatabasePDG (S.Zharko) + q = MCTrack->GetCharge() / 3.; + mass = MCTrack->GetMass(); Int_t iTrack = fvMCTracks.size(); //or iMCTrack? CbmL1MCTrack T(mass, q, vr, vp, iTrack, mother_ID, pdg, processID); @@ -1452,6 +1454,8 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i MC->q = particlePDG->Charge() / 3.0; MC->mass = particlePDG->Mass(); } + + // TODO: Add light nuclei (d, t, He3, He4): they are common in tracking but not accounted in TDatabasePDG (S.Zharko) return 0; @@ -1463,6 +1467,16 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i // TODO: Duplicated code (from CbmL1::ReadEvent) void CbmL1::HitMatch() { + // Clear contents + for (auto& hit : fvExternalHits) { + hit.mcPointIds.clear(); + } + + for (auto& point : fvMCPoints) { + point.hitIds.clear(); + } + + // Fill new contents const int NHits = fvExternalHits.size(); for (int iH = 0; iH < NHits; iH++) { CbmL1HitDebugInfo& hit = fvHitDebugInfo[iH]; diff --git a/reco/L1/CbmL1Track.h b/reco/L1/CbmL1Track.h index 8ba6383d035b103239486d66ecf5fdd35af19317..b7434e01bc957f38e83a2ca008a18335ef2545df 100644 --- a/reco/L1/CbmL1Track.h +++ b/reco/L1/CbmL1Track.h @@ -40,14 +40,35 @@ public: int GetNOfHits() { return Hits.size(); } void AddMCTrack(CbmL1MCTrack* mcTr) { mcTracks.push_back(mcTr); } + + /// Adds an index of MC track index + void AddMCTrackIndex(int iMT) { fvMcTrackIndexes.push_back_no_warning(iMT); } + + /// Clears the contents of matched MC track indexes (and pointers) + void ClearMatchedMCTracks() + { + mcTracks.clear(); + fvMcTrackIndexes.clear(); + } + vector<CbmL1MCTrack*>& GetMCTracks() { return mcTracks; } CbmL1MCTrack* GetMCTrack() { return mcTracks[0]; } int GetNMCTracks() { return mcTracks.size(); } bool IsGhost() { return (maxPurity < CbmL1Constants::MinPurity); } - void SetMaxPurity(double maxPurity_) { maxPurity = maxPurity_; } + + /// Gets a reference to MC track indexes + const auto& GetMCTrackIndexes() const { return fvMcTrackIndexes; } + + /// Gets max purity double GetMaxPurity() { return maxPurity; } + /// Sets max purity + /// NOTE: max purity is calculated as a ratio of max number of hits left by an actual track and the + /// total number of hits in the track + void SetMaxPurity(double maxPurity_) { maxPurity = maxPurity_; } + + static bool compareChi2(const CbmL1Track& a, const CbmL1Track& b) { return (a.chi2 < b.chi2); } static bool comparePChi2(const CbmL1Track* a, const CbmL1Track* b) { return (a->chi2 < b->chi2); } @@ -62,10 +83,15 @@ public: double fTrackTime; - map<int, int> hitMap; // how many hits from each mcTrack belong to current recoTrack + map<int, int> hitMap; // N hits (second) from each mcTrack (first is a MC track ID) belong to current recoTrack + // FIXME: SZh 14.12.2022: map => unordered_map private: // next members filled and used in Performance - vector<CbmL1MCTrack*> mcTracks; // array of assosiated mc Tracks. Should be only one. + vector<CbmL1MCTrack*> mcTracks; // array of assosiated recoTracks. Should be only one. + + L1Vector<int> fvMcTrackIndexes = {"CbmL1Track::fvMcTrackIndexes"}; // global indexes of MC tracks + // NOTE: mcTracks should be replaced with fvMcTrackIndexes + double maxPurity; // maximum persent of hits, which belong to one mcTrack. }; diff --git a/reco/L1/catools/CaToolsLinkKey.h b/reco/L1/catools/CaToolsLinkKey.h index 78c91dea0531b9029ce6ff6c117d9f388b572901..7f9ce0560f3e411e0e9abbc244fd9a399c11e754 100644 --- a/reco/L1/catools/CaToolsLinkKey.h +++ b/reco/L1/catools/CaToolsLinkKey.h @@ -29,9 +29,9 @@ namespace ca::tools return lhs.fFile == rhs.fFile && lhs.fEvent == rhs.fEvent && lhs.fIndex == rhs.fIndex; } - int fIndex = undef::kI32; - int fEvent = undef::kI32; - int fFile = undef::kI32; + int fIndex = undef::kI32; ///< Index of MC point/track in external data structures + int fEvent = undef::kI32; ///< Index of MC event + int fFile = undef::kI32; ///< Index of MC file }; } // namespace ca::tools diff --git a/reco/L1/catools/CaToolsMCData.cxx b/reco/L1/catools/CaToolsMCData.cxx index 85c1fb0e27635948c635a8b36aea1a938bc55cdc..64573973b6c1ceb535a777068a8ea985e8cc1139 100644 --- a/reco/L1/catools/CaToolsMCData.cxx +++ b/reco/L1/catools/CaToolsMCData.cxx @@ -9,6 +9,8 @@ #include "CaToolsMCData.h" +#include "CbmL1Hit.h" + #include <iomanip> #include <sstream> #include <utility> // for std::move @@ -66,17 +68,17 @@ void MCData::Swap(MCData& other) noexcept // --------------------------------------------------------------------------------------------------------------------- // -void MCData::AddPoint(const MCPoint& point, int index, int event, int file) +void MCData::AddPoint(const MCPoint& point) { - fmPointLinkMap[LinkKey(index, event, file)] = static_cast<int>(fvPoints.size()); + fmPointLinkMap[point.GetLinkKey()] = point.GetId(); fvPoints.push_back(point); } // --------------------------------------------------------------------------------------------------------------------- // -void MCData::AddTrack(const CbmL1MCTrack& track, int index, int event, int file) +void MCData::AddTrack(const MCTrack& track) { - fmTrackLinkMap[LinkKey(index, event, file)] = static_cast<int>(fvTracks.size()); + fmTrackLinkMap[track.GetLinkKey()] = track.GetId(); fvTracks.push_back(track); } @@ -91,6 +93,26 @@ void MCData::Clear() fmTrackLinkMap.clear(); } +// --------------------------------------------------------------------------------------------------------------------- +// +void MCData::InitTrackInfo(const L1Vector<CbmL1Hit>& vHits) +{ + for (auto& aTrk : fvTracks) { + // Assign hits to tracks + aTrk.ClearHitIndexes(); + auto& vHitIds = aTrk.GetHitIndexes(); + for (int iP : aTrk.GetPointIndexes()) { + const auto& point = fvPoints[iP]; + for (int iH : point.GetHitIndexes()) { + if (std::find(vHitIds.begin(), vHitIds.end(), iH) == vHitIds.end()) { aTrk.AddHitIndex(iH); } + } + } + // Initialize arrangements of points and hits within stations + aTrk.InitPointsInfo(fvPoints); + aTrk.InitHitsInfo(vHits); + } +} + // --------------------------------------------------------------------------------------------------------------------- // std::string MCData::ToString(int verbose) const @@ -109,10 +131,13 @@ std::string MCData::ToString(int verbose) const << "px [GeV/c]" << ' ' << setw(10) << setfill(' ') << "py [GeV/c]" << ' ' << setw(10) << setfill(' ') << "pz [GeV/c]" << '\n'; for (int i = 0; i < kNofTracksToPrint; ++i) { - msg << setw(10) << setfill(' ') << fvTracks[i].ID << ' ' << setw(10) << setfill(' ') << fvTracks[i].mother_ID - << ' ' << setw(10) << setfill(' ') << fvTracks[i].pdg << ' ' << setw(10) << setfill(' ') << fvTracks[i].z - << ' ' << setw(10) << setfill(' ') << fvTracks[i].px << ' ' << setw(10) << setfill(' ') << fvTracks[i].py - << ' ' << setw(10) << setfill(' ') << fvTracks[i].pz << '\n'; + msg << setw(10) << setfill(' ') << fvTracks[i].GetId() << ' '; + msg << setw(10) << setfill(' ') << fvTracks[i].GetMotherId() << ' '; + msg << setw(10) << setfill(' ') << fvTracks[i].GetPdgCode() << ' '; + msg << setw(10) << setfill(' ') << fvTracks[i].GetStartZ() << ' '; + msg << setw(10) << setfill(' ') << fvTracks[i].GetPx() << ' '; + msg << setw(10) << setfill(' ') << fvTracks[i].GetPy() << ' '; + msg << setw(10) << setfill(' ') << fvTracks[i].GetPz() << '\n'; } } return msg.str(); diff --git a/reco/L1/catools/CaToolsMCData.h b/reco/L1/catools/CaToolsMCData.h index 8f4584864b48e6c2dfef90b2a8cc2ce509024fac..1fab4a312506b38cea6bec6bc1796f53ef35738f 100644 --- a/reco/L1/catools/CaToolsMCData.h +++ b/reco/L1/catools/CaToolsMCData.h @@ -10,13 +10,12 @@ #ifndef CaToolsMCData_h #define CaToolsMCData_h 1 -#include "CbmL1MCTrack.h" - #include <numeric> #include <string> #include "CaToolsLinkKey.h" #include "CaToolsMCPoint.h" +#include "CaToolsMCTrack.h" #include "L1Constants.h" #include "L1Vector.h" @@ -56,47 +55,45 @@ namespace ca::tools /// \param iPoint Index of MC point /// \param - /// Adds an MC point + /// Adds an MC point to points container and a corresponding link key to the point index map /// \param point MC point object - /// \param index Index of point in a link - /// \param event Event index of a link - /// \param file File index of a link - void AddPoint(const MCPoint& point, int index, int event, int file); + void AddPoint(const MCPoint& point); - /// Adds an MC track + /// Adds an MC track to tracks container and a corresponding link key to the link index map /// \param track MC track object - /// \param index Index of point in a link - /// \param event Event index of a link - /// \param file File index of a link - void AddTrack(const CbmL1MCTrack& track, int index, int event, int file); + void AddTrack(const MCTrack& track); /// Clears contents void Clear(); - /// Finds an MC point global index by given local index, event and file of a link - /// \param index Local index of MC point - /// \param event Event of link - /// \param file File of link - int FindGlobalPointIndex(int index, int event, int file) const + /// Finds an index of MC point in internal point container + /// \param index Index of MC point in external point container + /// \param event Index of MC event + /// \param file Index of MC file + /// \return Index of MC point in internal point container within event/TS + /// If the point is not found the function returns -1 + int FindInternalPointIndex(int index, int event, int file) const { auto it = fmPointLinkMap.find(LinkKey(index, event, file)); return (it != fmPointLinkMap.cend()) ? it->second : -1; } - /// Finds an MC track global index by given local index, event and file of a link - /// \param index Local index of track - /// \param event Index of event - /// \param file Index of file - int FindGlobalTrackIndex(int index, int event, int file) const + /// Finds an index of MC track in internal track container + /// \param index Index of MC track in external track container + /// \param event Index of MC event + /// \param file Index of MC file + /// \return Index of MC track in internal track container within event/TS + /// If the track is not found, the function returns -1 + int FindInternalTrackIndex(int index, int event, int file) const { auto it = fmTrackLinkMap.find(LinkKey(index, event, file)); return (it != fmTrackLinkMap.cend()) ? it->second : -1; } - /// Gets number of tracks + /// Gets number of tracks in this event/TS int GetNofTracks() const { return fvTracks.size(); } - /// Gets number of points + /// Gets number of points in this event/TS int GetNofPoints() const { return fvPoints.size(); } /// Gets a reference to MC point by its index @@ -114,12 +111,25 @@ namespace ca::tools /// \param iHit Index of hit auto GetPointIndexOfHit(int iHit) const { return fvPointIndexOfHit[iHit]; } - /// Gets a reference to MC track by its index + /// Gets a reference to MC track by its internal index const auto& GetTrack(int idx) const { return fvTracks[idx]; } + /// Gets a mutual reference to MC track by its internal index + auto& GetTrack(int idx) { return fvTracks[idx]; } + /// Gets a reference to the vector of tracks const auto& GetTrackContainer() const { return fvTracks; } + /// Gets a mutual reference to the vector of tracks + auto& GetTrackContainer() { return fvTracks; } + + /// \brief Initialize information about points and hits association with MC track + /// \param vHits Vector of hit objects + /// Initialize tracks: defines indexes of hits and points related to the track, calculates max number of points and + /// hits on a station, number of consecutive stations containing a hit or point and number of stations and points + /// with hits. + void InitTrackInfo(const L1Vector<CbmL1Hit>& vHits); + /// Registers index of point for a given index of hit /// \param iHit Index of hit /// \param iPoint Index of point @@ -146,12 +156,12 @@ namespace ca::tools std::string ToString(int verbose = 1) const; private: - // ********************** - // ** Member variables ** - // ********************** + // ****************************** + // ** Member variables ** + // ****************************** - L1Vector<MCPoint> fvPoints = {"ca::tools::MCData::fvMCPoints"}; ///< Container of points - L1Vector<CbmL1MCTrack> fvTracks = {"ca::tools::MCData::fvMCTracks"}; ///< Container of tracks + L1Vector<MCPoint> fvPoints = {"ca::tools::MCData::fvPoints"}; ///< Container of points + L1Vector<MCTrack> fvTracks = {"ca::tools::MCData::fvTracks"}; ///< Container of tracks /// Correspondence of MC point index to the global hit index L1Vector<int> fvPointIndexOfHit = {"ca::tools::MCData::fvPointIndexOfHit"}; diff --git a/reco/L1/catools/CaToolsMCPoint.cxx b/reco/L1/catools/CaToolsMCPoint.cxx index 8c53161671075a21b8c65ad282154b7f92f1a1e0..3df6aaea5d13ed8cbd7ce203cd3ff577ac7272f6 100644 --- a/reco/L1/catools/CaToolsMCPoint.cxx +++ b/reco/L1/catools/CaToolsMCPoint.cxx @@ -15,16 +15,6 @@ using namespace ca::tools; -// --------------------------------------------------------------------------------------------------------------------- -// -void MCPoint::SetLinkAddress(int pointID, int eventID, int fileID) -{ - // TODO: SZh 21.11.2022: Test uniqueness of the pointID among different subsystems in link definition - fPointID = pointID; - fEventID = eventID; - fFileID = fileID; -} - // --------------------------------------------------------------------------------------------------------------------- // std::string MCPoint::ToString(int verbose, bool printHeader) const @@ -49,6 +39,7 @@ std::string MCPoint::ToString(int verbose, bool printHeader) const msg << std::setw(14) << std::setfill(' ') << "zOut [cm]" << '|'; msg << std::setw(14) << std::setfill(' ') << "p [GeV/c]" << '|'; msg << std::setw(10) << std::setfill(' ') << "point ID" << ' '; + msg << std::setw(10) << std::setfill(' ') << "point ID (ext)" << ' '; msg << std::setw(10) << std::setfill(' ') << "event ID" << ' '; msg << std::setw(10) << std::setfill(' ') << "file ID" << ' '; } @@ -56,9 +47,9 @@ std::string MCPoint::ToString(int verbose, bool printHeader) const } else { if (verbose > 0) { - msg << std::setw(10) << std::setfill(' ') << fTrackID << ' '; - msg << std::setw(10) << std::setfill(' ') << fMotherID << '|'; - msg << std::setw(10) << std::setfill(' ') << fStationID << '|'; + msg << std::setw(10) << std::setfill(' ') << fTrackId << ' '; + msg << std::setw(10) << std::setfill(' ') << fMotherId << '|'; + msg << std::setw(10) << std::setfill(' ') << fStationId << '|'; msg << std::setw(10) << std::setfill(' ') << fPdgCode << ' '; msg << std::setw(10) << std::setfill(' ') << fMass << ' '; msg << std::setw(10) << std::setfill(' ') << fCharge << '|'; @@ -70,9 +61,10 @@ std::string MCPoint::ToString(int verbose, bool printHeader) const msg << std::setw(14) << std::setfill(' ') << fPosIn[2] << ' '; msg << std::setw(14) << std::setfill(' ') << fPosOut[2] << '|'; msg << std::setw(14) << std::setfill(' ') << this->GetP() << '|'; - msg << std::setw(10) << std::setfill(' ') << fPointID << ' '; - msg << std::setw(10) << std::setfill(' ') << fEventID << ' '; - msg << std::setw(10) << std::setfill(' ') << fFileID << ' '; + msg << std::setw(10) << std::setfill(' ') << fId << ' '; + msg << std::setw(10) << std::setfill(' ') << fLinkKey.fIndex << ' '; + msg << std::setw(10) << std::setfill(' ') << fLinkKey.fEvent << ' '; + msg << std::setw(10) << std::setfill(' ') << fLinkKey.fFile << ' '; } } } diff --git a/reco/L1/catools/CaToolsMCPoint.h b/reco/L1/catools/CaToolsMCPoint.h index 465b8569c6864d82e77f43b3fcddab6d1ee17d26..83aec7054bfa68cdd14eb2d93cbbf960d7808b3d 100644 --- a/reco/L1/catools/CaToolsMCPoint.h +++ b/reco/L1/catools/CaToolsMCPoint.h @@ -12,9 +12,11 @@ #include <string> +#include "CaToolsLinkKey.h" #include "L1Undef.h" #include "L1Vector.h" + enum class L1DetectorID; /// Class describes a Monte-Carlo point used in CA tracking QA analysis @@ -41,9 +43,9 @@ namespace ca::tools MCPoint& operator=(MCPoint&&) = default; - /// Adds hit ID - /// \param hitID A hit index in the external hits array - void AddHitID(int hitID) { fHitIDs.push_back_no_warning(hitID); } + /// Adds index of hits from the container of hits of event/TS + /// \param iH A hit index in the external hits container of event/TS + void AddHitID(int iH) { fvHitIndexes.push_back_no_warning(iH); } // ********************* @@ -54,19 +56,28 @@ namespace ca::tools double GetCharge() const { return fCharge; } /// Gets detector ID - L1DetectorID GetDetectorID() const { return fDetectorID; } + L1DetectorID GetDetectorId() const { return fDetectorId; } + + /// Gets MC event ID + int GetEventId() const { return fLinkKey.fEvent; } + + /// Gets MC file ID + int GetFileId() const { return fLinkKey.fFile; } + + /// Gets index of this point in internal CA container + int GetId() const { return fId; } - /// Gets event ID - int GetEventID() const { return fEventID; } + /// Gets container of matched hit indexes + const auto& GetHitIndexes() const { return fvHitIndexes; } - /// Gets file ID - int GetFileID() const { return fFileID; } + /// Gets link key + LinkKey GetLinkKey() const { return fLinkKey; } /// Gets mass of the particle [GeV/c2] double GetMass() const { return fMass; } /// Gets mother ID of the track - int GetMotherID() const { return fMotherID; } + int GetMotherId() const { return fMotherId; } /// Gets track momentum absolute value at reference z of station [GeV/c] double GetP() const { return std::sqrt(fMom[0] * fMom[0] + fMom[1] * fMom[1] + fMom[2] * fMom[2]); } @@ -77,9 +88,6 @@ namespace ca::tools /// Gets track momentum absolute value at entrance to station [GeV/c] double GetPIn() const { return std::sqrt(fMomIn[0] * fMomIn[0] + fMomIn[1] * fMomIn[1] + fMomIn[2] * fMomIn[2]); } - /// Gets point ID - int GetPointID() const { return fPointID; } - /// Gets track momentum absolute value at exit of station [GeV/c] double GetPOut() const { @@ -114,13 +122,13 @@ namespace ca::tools double GetPzOut() const { return fMomOut[2]; } /// Gets global ID of the active tracking station - int GetStationID() const { return fStationID; } + int GetStationId() const { return fStationId; } /// Gets time [ns] double GetTime() const { return fTime; } - /// Gets ID of track - int GetTrackID() const { return fTrackID; } + /// Gets ID of track from the internal CA MC track container (within event/TS) + int GetTrackId() const { return fTrackId; } /// Gets x coordinate at reference z of station [cm] double GetX() const { return fPos[0]; } @@ -149,11 +157,6 @@ namespace ca::tools /// Gets z coordinate at exit of station [cm] double GetZOut() const { return fPosOut[2]; } - - /// Gets Hit ID - int MapHitID(int iHitID) { return fHitIDs[iHitID]; } - - // ********************* // ** Setters ** // ********************* @@ -162,19 +165,25 @@ namespace ca::tools void SetCharge(double charge) { fCharge = charge; } /// Sets detector ID - void SetDetectorID(L1DetectorID detID) { fDetectorID = detID; } + void SetDetectorId(L1DetectorID detId) { fDetectorId = detId; } + + /// Sets index of MC event containing this point + void SetEventId(int eventId) { fLinkKey.fEvent = eventId; } + + /// Sets index of this point in external data structures + void SetExternalId(int id) { fLinkKey.fIndex = id; } + + /// Sets index of MC file containing this point + void SetFileId(int fileId) { fLinkKey.fFile = fileId; } + + /// Sets index of this point in the CA internal structure + void SetId(int id) { fId = id; } /// Sets particle mass [GeV/c2] void SetMass(double mass) { fMass = mass; } - /// Sets mother ID - void SetMotherID(int motherID) { fMotherID = motherID; } - - /// Sets link address - /// \param pointID Index of - /// \param eventID - /// \param fileID - void SetLinkAddress(int pointID, int eventID, int fileID); + /// Sets index of mother track in the internal CA data structures + void SetMotherId(int motherId) { fMotherId = motherId; } /// Sets PDG code void SetPdgCode(int pdg) { fPdgCode = pdg; } @@ -207,13 +216,13 @@ namespace ca::tools void SetPzOut(double pz) { fMomOut[2] = pz; } /// Sets global index of active station - void SetStationID(int stationID) { fStationID = stationID; } + void SetStationId(int stationId) { fStationId = stationId; } /// Sets time [ns] void SetTime(double time) { fTime = time; } - /// Sets track ID - void SetTrackID(int trackID) { fTrackID = trackID; } + /// Sets track ID in the CA internal track container (within event/TS) + void SetTrackId(int trackId) { fTrackId = trackId; } /// Sets x coordinate at reference z of station [cm] void SetX(double x) { fPos[0] = x; } @@ -262,22 +271,21 @@ namespace ca::tools std::array<double, 3> fMomIn = {undef::kD, undef::kD, undef::kD}; ///< Momentum at entrance to station [cm] std::array<double, 3> fMomOut = {undef::kD, undef::kD, undef::kD}; ///< Momentum at exit of station [cm] - int fPdgCode = undef::kI32; ///< Particle PDG code - double fMass = undef::kD; ///< Particle mass [GeV/c2] - double fCharge = undef::kD; ///< Particle charge [e] - double fTime = undef::kD; ///< Point time [ns] + double fMass = undef::kD; ///< Particle mass [GeV/c2] + double fCharge = undef::kD; ///< Particle charge [e] + double fTime = undef::kD; ///< Point time [ns] - int fTrackID = undef::kI32; ///< Track ID - int fMotherID = undef::kI32; ///< Mother particle ID (NOTE: probably is redundant) - int fStationID = undef::kI32; ///< Global index of active tracking station + LinkKey fLinkKey = {undef::kI32, undef::kI32, undef::kI32}; ///< Link key of point - L1DetectorID fDetectorID; ///< Detector ID of MC point + int fPdgCode = undef::kI32; ///< Particle PDG code + int fId = undef::kI32; ///< Index of MC point in the external MC point container + int fTrackId = undef::kI32; ///< Index of associated MC track in CA internal track container within TS/event + int fMotherId = undef::kI32; ///< Index of mother track in CA internal data structures (within event/TS) + int fStationId = undef::kI32; ///< Global index of active tracking station - int fPointID = undef::kI32; ///< TODO - int fEventID = undef::kI32; - int fFileID = undef::kI32; + L1DetectorID fDetectorId; ///< Detector ID of MC point - L1Vector<int> fHitIDs {"ca::tools::MCPoint::fHitIDs"}; + L1Vector<int> fvHitIndexes {"ca::tools::MCPoint::fvHitIndexes"}; }; } // namespace ca::tools diff --git a/reco/L1/catools/CaToolsMCTrack.cxx b/reco/L1/catools/CaToolsMCTrack.cxx new file mode 100644 index 0000000000000000000000000000000000000000..1a07b81dee6b97a713680082b81ea4fd0cf30375 --- /dev/null +++ b/reco/L1/catools/CaToolsMCTrack.cxx @@ -0,0 +1,150 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov, Sergei Zharko [committer] */ + +/// \file CaToolsMCTrack.h +/// \brief Class represents a MC track for CA tracking QA and performance (header) +/// \since 22.11.2022 +/// \author S.Zharko <s.zharko@gsi.de> (based on CbmL1MCTrack class by I.Kisel and S.Gorbunov) + +#include "CaToolsMCTrack.h" + +#include "CbmL1Hit.h" + +#include "CaToolsMCPoint.h" + +using namespace ca::tools; + + +// --------------------------------------------------------------------------------------------------------------------- +// +void MCTrack::Clear() +{ + MCTrack other(std::move(*this)); + ClearPointIndexes(); + ClearHitIndexes(); + ClearRecoTrackIndexes(); + ClearTouchTrackIndexes(); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void MCTrack::InitHitsInfo(const L1Vector<CbmL1Hit>& vHits) +{ + // NOTE: vHits must be sorted over stations! + fMaxNofHitsOnStation = 0; + fTotNofStationsWithHit = 0; + fNofConsStationsWithHit = 0; + + int iStPrev = -1; // Station index of previous hit + int currNofConsStationsWithHit = 0; // Current number of consecutive stations with a point in MC track + int currMaxNofHitsOnStation = 0; // Current max number of hits on station + + // Loop over hit indexes assigned to this track + for (int iH : fvHitIndexes) { + int iSt = vHits[iH].iStation; // current index of active station + // Check if the track is going to the backward direction and is not reconstructable + if (iSt < iStPrev) { + fMaxNofHitsOnStation = 0; + fNofConsStationsWithHit = 0; + //fTotNofStationsWithHit = 0; // TODO: Check, if this is correct + return; + } + + // Check if this hit is on the same station as the previous one and update max number of hits within a station a + // station + if (iSt == iStPrev) { currMaxNofHitsOnStation++; } // the same station + else { // the next station (reset hits counter and update number of stations) + if (currMaxNofHitsOnStation > fMaxNofHitsOnStation) { fMaxNofHitsOnStation = currMaxNofHitsOnStation; } + currMaxNofHitsOnStation = 1; + fTotNofStationsWithHit++; + } + + // Check if this hit is on the next station comparing with the previous hit and update the number of consecutive + // stations + if (iSt - iStPrev == 1) { currNofConsStationsWithHit++; } + else if (iSt - iStPrev > 1) { + if (currNofConsStationsWithHit > fNofConsStationsWithHit) { + fNofConsStationsWithHit = currNofConsStationsWithHit; + } + currNofConsStationsWithHit = 1; + } + + iStPrev = iSt; + } // Loop over hit indexes assigned to this track: end + if (currMaxNofHitsOnStation > fMaxNofHitsOnStation) { fMaxNofHitsOnStation = currMaxNofHitsOnStation; } + if (currNofConsStationsWithHit > fNofConsStationsWithHit) { fNofConsStationsWithHit = currNofConsStationsWithHit; } +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void MCTrack::InitPointsInfo(const L1Vector<MCPoint>& vPoints) +{ + fMaxNofPointsOnStation = 0; + fMaxNofPointsOnSensor = 0; + fTotNofStationsWithPoint = 0; + + int currMaxNofPointsOnStation = 0; // current max number of points within station + int currMaxNofPointsOnSensor = 0; // current max number of points within sensor + int iStPrev = -1; // index of station for previous point + float zPrev = 0; // Z position of previous point + + // Loop over point indexes assigned to this track + for (int iP : fvPointIndexes) { + const auto& point = vPoints[iP]; // current MC point object + int iSt = point.GetStationId(); // current index of active station + + // Check if this point is on the same station as the previous one and update max number of points within a station + if (iSt == iStPrev) { currMaxNofPointsOnStation++; } // the same station + else { // the next station (reset points counter and update number of stations) + if (currMaxNofPointsOnStation > fMaxNofPointsOnStation) { fMaxNofPointsOnStation = currMaxNofPointsOnStation; } + currMaxNofPointsOnStation = 1; + fTotNofStationsWithPoint++; + } + + // Check if this point is on the same sensor as the previous one and update max number of points within a sensor + // NOTE: points sometimes have different z positions within a sensor, so the max number of points within a sensor + // might be calculated incorrectly. (TODO: remove this variable?) + if (point.GetZ() == zPrev) { currMaxNofPointsOnSensor++; } + else { + if (currMaxNofPointsOnSensor > fMaxNofPointsOnSensor) { fMaxNofPointsOnSensor = currMaxNofPointsOnSensor; } + currMaxNofPointsOnSensor = 1; + } + + iStPrev = iSt; + zPrev = point.GetZ(); + } // Loop over point indexes assigned to this track + + if (currMaxNofPointsOnStation > fMaxNofPointsOnStation) { fMaxNofPointsOnStation = currMaxNofPointsOnStation; } + if (currMaxNofPointsOnSensor > fMaxNofPointsOnSensor) { fMaxNofPointsOnSensor = currMaxNofPointsOnSensor; } + + fNofConsStationsWithPoint = 0; + int currNofConsStationsWithPoint = 0; // current number of consecutive stations with points + iStPrev = -1; + // Loop over point indexes assigned to this track + for (int iP : fvPointIndexes) { + int iSt = vPoints[iP].GetStationId(); + // Check if this point is on the next station comparing with the previous point and update the number of consecutive + // stations + if (iSt - iStPrev == 1) { currNofConsStationsWithPoint++; } + else if (iSt - iStPrev > 1) { + if (currNofConsStationsWithPoint > fNofConsStationsWithPoint) { + fNofConsStationsWithPoint = currNofConsStationsWithPoint; + } + currNofConsStationsWithPoint = 1; + } + + if (iSt <= iStPrev) { continue; } // Tracks propagating in backward direction + iStPrev = iSt; + } + if (currNofConsStationsWithPoint > fNofConsStationsWithPoint) { + fNofConsStationsWithPoint = currNofConsStationsWithPoint; + } +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void MCTrack::SortPointIndexes(const std::function<bool(const int& lhs, const int& rhs)>& cmpFn) +{ + std::sort(fvPointIndexes.begin(), fvPointIndexes.end(), cmpFn); +} diff --git a/reco/L1/catools/CaToolsMCTrack.h b/reco/L1/catools/CaToolsMCTrack.h new file mode 100644 index 0000000000000000000000000000000000000000..da2a1bcef887725f1ecc43b371d2407b48af7dce --- /dev/null +++ b/reco/L1/catools/CaToolsMCTrack.h @@ -0,0 +1,365 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov, Sergei Zharko [committer] */ + +/// \file CaToolsMCTrack.h +/// \brief Class represents a MC track for CA tracking QA and performance (implementation) +/// \since 22.11.2022 +/// \author S.Zharko <s.zharko@gsi.de> (based on CbmL1MCTrack class by I.Kisel and S.Gorbunov) + +#ifndef CaToolsMCTrack_h +#define CaToolsMCTrack_h 1 + +#include <functional> + +#include "CaToolsLinkKey.h" +#include "L1Undef.h" +#include "L1Vector.h" + +class CbmL1Hit; + +namespace ca::tools +{ + class MCPoint; + class MCTrack { + public: + /// Default constructor + MCTrack() = default; + + /// Destructor + ~MCTrack() = default; + + /// Copy constructor + MCTrack(const MCTrack&) = default; + + /// Move constructor + MCTrack(MCTrack&&) = default; + + /// Copy assignment operator + MCTrack& operator=(const MCTrack&) = default; + + /// Move assignment operator + MCTrack& operator=(MCTrack&&) = default; + + /// Adds index of point, associated with this MC track + /// \param iP Index of MC point in the container of internal MC points for this event/TS + void AddPointIndex(int iP) { fvPointIndexes.push_back_no_warning(iP); } + + /// Adds index of hit, associated with this MC track + /// \param iH Index of hit in the container of internal hits for this event/TS + void AddHitIndex(int iH) { fvHitIndexes.push_back_no_warning(iH); } + + /// Adds index of reconstructed track, associated with this MC track + /// \param iTre Index of track in the array of reconstructed tracks of the TS + void AddRecoTrackIndex(int iTre) { fvRecoTrackIndexes.push_back_no_warning(iTre); } + + /// Adds index of reconstructed track, which is not associated with this MC track but contains some hits, produced + /// by this MC track + /// \param iTre Index of track in the array of reconstructed tracks of the TS + void AddTouchTrackIndex(int iTre) { fvTouchTrackIndexes.push_back_no_warning(iTre); } + + /// Clears contents + void Clear(); + + /// Clears container of point indexes + void ClearPointIndexes() { fvPointIndexes.clear(); } + + /// Clears container of hit indexes + void ClearHitIndexes() { fvHitIndexes.clear(); } + + /// Clears container of reconstructed track indexes + void ClearRecoTrackIndexes() { fvRecoTrackIndexes.clear(); } + + /// Clears container of indexes of reconstructed tracks, which are not assigned to this MC track, but contain some + /// of its hits + void ClearTouchTrackIndexes() { fvTouchTrackIndexes.clear(); } + + // ********************* + // ** Getters ** + // ********************* + + /// Gets charge [e] + double GetCharge() const { return fCharge; } + + /// Gets total energy [GeV] + double GetE() const { return std::sqrt(GetP() * GetP() + fMass * fMass); } + + /// Gets kinetic energy [GeV] + double GetEkin() const { return GetE() - fMass; } + + /// Gets index of MC event containing this track in external data structures + int GetEventId() const { return fLinkKey.fEvent; } + + /// Gets index of the track in the external data structures + int GetExternalId() const { return fLinkKey.fIndex; } + + /// Gets index of MC file containing this track in external data structures + int GetFileId() const { return fLinkKey.fFile; } + + /// Gets a reference to associated hit indexes + const auto& GetHitIndexes() const { return fvHitIndexes; } + + /// Gets index of track + int GetId() const { return fId; } + + /// Gets link key + LinkKey GetLinkKey() const { return fLinkKey; } + + /// Gets particle mass [GeV/c2] + double GetMass() const { return fMass; } + + /// Gets max number of hits within a station + int GetMaxNofHitsOnStation() const { return fMaxNofHitsOnStation; } + + /// Gets max number of points within a station + int GetMaxNofPointsOnStation() const { return fMaxNofPointsOnStation; } + + /// Gets index of mother track in CA internal data structures + int GetMotherId() const { return fMotherId; } + + /// Gets number of consecutive stations with hits + int GetNofConsStationsWithHit() const { return fNofConsStationsWithHit; } + + /// Gets number of consecutive stations with MC points + int GetNofConsStationsWithPoint() const { return fNofConsStationsWithPoint; } + + /// Gets number of hits + int GetNofHits() const { return fvHitIndexes.size(); } + + /// Gets number of points + int GetNofPoints() const { return fvPointIndexes.size(); } + + /// Gets number of assigned reconstructed tracks + int GetNofRecoTracks() const { return fvRecoTrackIndexes.size(); } + + /// Gets number of reconstructed tracks, which contain hits from this MC track + int GetNofTouchTracks() const { return fvTouchTrackIndexes.size(); } + + /// Gets absolute momentum [GeV/c] + double GetP() const { return std::sqrt(fMom[0] * fMom[0] + fMom[1] * fMom[1] + fMom[2] * fMom[2]); } + + /// Gets PDG encoding + int GetPdgCode() const { return fPdgCode; } + + /// Gets azimuthal angle [rad] + double GetPhi() const; + + /// Gets a reference to associated point indexes + const auto& GetPointIndexes() const { return fvPointIndexes; } + + /// Gets process ID + unsigned GetProcessId() const { return fProcId; } + + /// Gets transverse momentum [GeV/c] + double GetPt() const { return std::sqrt(fMom[0] * fMom[0] + fMom[1] * fMom[1]); } + + /// Gets x component of momentum [GeV/c] + double GetPx() const { return fMom[0]; } + + /// Gets x component of momentum [GeV/c] + double GetPy() const { return fMom[1]; } + + /// Gets x component of momentum [GeV/c] + double GetPz() const { return fMom[2]; } + + /// Gets rapidity + double GetRapidity() const { return 0.5 * std::log((GetE() + fMom[2]) / (GetE() - fMom[2])); } + + /// Gets a reference to vector associated reconstructed track indexes + const auto& GetRecoTrackIndexes() const { return fvRecoTrackIndexes; } + + /// Gets particle speed [c] + double GetSpeed() const { return GetP() / GetE(); } + + /// Gets time of the track vertex [ns] + double GetStartT() const { return fTime; } + + /// Gets x component of the track vertex [cm] + double GetStartX() const { return fPos[0]; } + + /// Gets y component of the track vertex [cm] + double GetStartY() const { return fPos[1]; } + + /// Gets z component of the track vertex [cm] + double GetStartZ() const { return fPos[2]; } + + /// Gets total number of stations with hits + int GetTotNofStationsWithHit() const { return fTotNofStationsWithHit; } + + /// Gets total number of stations with MC points + int GetTotNofStationsWithPoint() const { return fTotNofStationsWithPoint; } + + /// Gets a reference to vector of reconstructed track indexes, not associated with this MC track but containing some + /// hits, produced by this MC track + const auto& GetTouchTrackIndexes() const { return fvTouchTrackIndexes; } + + /// \brief Initializes information about MC track hits arrangement within stations + /// Defines: + /// #1) Number of stations with hits + /// #2) Maximal number of hits within one station + /// #3) Number of consecutive stations with a hit in MC track + /// \param vHits Vector of hits for a given TS + void InitHitsInfo(const L1Vector<CbmL1Hit>& vHits); + + /// \brief Initializes information about MC track points arrangement within stations + /// Defines: + /// #1) Number of stations with points + /// #2) Maximal number of points within one station + /// #3) Maximal number of points within one sensor (with same z-position) + /// #4) Number of consecutive stations with a point in MC track + /// \param vPoints Vector of points for a given TS + void InitPointsInfo(const L1Vector<MCPoint>& vPoints); + + + // ******************* + // ** Flags ** + // ******************* + + /// Returns .... TODO + bool IsAdditional() const { return fIsAdditional; } + + /// Returns true, if this MC track + bool IsDisturbed() const { return fvTouchTrackIndexes.size(); } + + /// Returns flag, if the track is primary (process ID is 0 either mother ID is -1) + bool IsPrimary() const { return fProcId == 0; } + + /// Returns flag, if track is reconstructable. The definition of particle reconstructability depends on the + /// performance mode and should be determined in the experiment-specific MC module + bool IsReconstructable() const { return fIsReconstructable; } + + bool IsReconstructed() const { return fvRecoTrackIndexes.size(); } + + /// Returns flag, if the track comes from a real particle (true), or from generated noise (false) + bool IsSignal() const { return fIsSignal; } + + + // ********************* + // ** Setters ** + // ********************* + + /// Sets charge [e] + void SetCharge(double q) { fCharge = q; } + + /// Sets index of MC event containing this track in external data structures + void SetEventId(int iEvent) { fLinkKey.fEvent = iEvent; } + + /// Sets index of track in external data structures + /// \note This index should be used to navigate through links + void SetExternalId(int id) { fLinkKey.fIndex = id; } + + /// Sets index of MC file containing this track in external data structures + void SetFileId(int iFile) { fLinkKey.fFile = iFile; } + + /// Sets index of track in the CA internal data structure (within event/TS) + void SetId(int id) { fId = id; } + + /// Sets index of mother track + /// \note Mother ID should refer to the internal index of track + void SetMotherId(int motherId) { fMotherId = motherId; } + + /// Sets particle mass [GeV/c2] + void SetMass(double mass) { fMass = mass; } + + /// Sets PDG encoding + void SetPdgCode(int pdg) { fPdgCode = pdg; } + + /// Sets process ID + void SetProcessId(unsigned procId) { fProcId = procId; } + + /// Sets x component of momentum [GeV/c] + void SetPx(double px) { fMom[0] = px; } + + /// Sets x component of momentum [GeV/c] + void SetPy(double py) { fMom[1] = py; } + + /// Sets x component of momentum [GeV/c] + void SetPz(double pz) { fMom[2] = pz; } + + /// Sets flag, if the track is additional (not reconstructable but still interesting) + void SetFlagAdditional(bool isAdditional) { fIsAdditional = isAdditional; } + + /// Sets flag, if the track is reconstructable + void SetFlagReconstructable(bool isReconstructable) { fIsReconstructable = isReconstructable; } + + /// Sets flag, if the track comes from signal + void SetFlagSignal(bool isSignal) { fIsSignal = isSignal; } + + /// Sets time of the track vertex [ns] + void SetStartT(double t) { fTime = t; } + + /// Sets x component of the track vertex [cm] + void SetStartX(double x) { fPos[0] = x; } + + /// Sets y component of the track vertex [cm] + void SetStartY(double y) { fPos[1] = y; } + + /// Sets z component of the track vertex [cm] + void SetStartZ(double z) { fPos[2] = z; } + + /// Sorts points inside track by a provided hit comparison function function + /// \param cmpFn Functional object to compare mcPoints + void SortPointIndexes(const std::function<bool(const int& lhs, const int& rhs)>& cmpFn); + + private: + // **************************** + // ** Data variables ** + // **************************** + + double fMass = undef::kD; ///< Particle mass [GeV/c2] + double fCharge = undef::kD; ///< Particle charge [e] + double fTime = undef::kD; ///< Time of track [cm] + + std::array<double, 3> fPos = {undef::kD, undef::kD, undef::kD}; ///< Track vertex components [cm] + std::array<double, 3> fMom = {undef::kD, undef::kD, undef::kD}; ///< Momentum components [GeV/c] + + int fPdgCode = undef::kI32; ///< PDG encoding + unsigned fProcId = undef::kU32; ///< Process ID (from ROOT::TProcessID) + + // Track address + int fId = undef::kI32; ///< Index of MC track in internal container for TS/event + int fMotherId = undef::kI32; ///< Index of mother MC track in the external tracks container + + LinkKey fLinkKey = {undef::kI32, undef::kI32, + undef::kI32}; ///< A link key of this track in the external data structures + + bool fIsSignal = false; ///< If the track comes from signal + bool fIsReconstructable = false; ///< If track is reconstructable + bool fIsAdditional = false; ///< If track is not reconstructable, but still interesting + + // Arrangement of hits and points within the stations + int fNofConsStationsWithHit = undef::kI32; ///< Number of consecutive stations with hits + int fNofConsStationsWithPoint = undef::kI32; ///< Number of consecutive stations with points + int fTotNofStationsWithHit = undef::kI32; ///< Total number of stations with hits + int fTotNofStationsWithPoint = undef::kI32; ///< Total number of stations with MC points + int fMaxNofPointsOnStation = undef::kI32; ///< Max number of MC points on a station + int fMaxNofPointsOnSensor = undef::kI32; ///< Max number of MC points with same Z (means on same sensor) + int fMaxNofHitsOnStation = undef::kI32; ///< Max number of hits on a station + + L1Vector<int> fvPointIndexes = {"ca::tools::fvPointIndexes"}; ///< Indexes of MC points in ext.container + L1Vector<int> fvHitIndexes = {"ca::tools::fvHitIndexes"}; ///< Indexes of hits in int.container + L1Vector<int> fvRecoTrackIndexes = {"ca::tools::fvRecoTrackIndexes"}; ///< Indexes of associated reco tracks + L1Vector<int> fvTouchTrackIndexes = {"ca::tools::fvTouchTrackIndexes"}; ///< Pointers to non-associated tracks, + ///< which use hits from this track + }; + +} // namespace ca::tools + + +// ************************************ +// ** Inline implementations ** +// ************************************ + +inline double ca::tools::MCTrack::GetPhi() const +{ + // NOTE: Implementation was taken from ROOT TMath::Atan2 + constexpr double kRootPi = 3.14159265358979323846; // Value of Pi, used in ROOT TMath + if (fMom[0] != 0) { return std::atan2(-fMom[1], -fMom[0]); } + if (fMom[1] == 0) { return 0; } + if (fMom[1] < 0) { return kRootPi / 2; } + else { + return -kRootPi / 2; + } +} + +#endif // CaToolsMCTrack_h diff --git a/reco/L1/catools/CaToolsPerformance.h b/reco/L1/catools/CaToolsPerformance.h deleted file mode 100644 index 4f2cb15858cfbd4034459d1ea92b29b1807c483f..0000000000000000000000000000000000000000 --- a/reco/L1/catools/CaToolsPerformance.h +++ /dev/null @@ -1,64 +0,0 @@ -/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergey Gorbunov, Sergei Zharko [committer] */ - -/// \file CaToolsPerformance.h -/// \brief Tracking performance class (header) -/// \since 23.09.2022 -/// \author S.Zharko <s.zharko@gsi.de> - -#ifndef CaToolsPerformance_h -#define CaToolsPerformance_h 1 - -#include "CaToolsMCData.h" - -class L1Parameters; - -namespace ca -{ - namespace tools - { - /// Class ca::tools::Performance defines CA tracking internal performance measurement - /// - class Performance { - public: - // ***************************************** - // ** Constructors and destructor ** - // ***************************************** - - - /// Default constructor - Performance() = default; - - /// Destructor - ~Performance() = default; - - /// Copy constructor - Performance(const Performance& other) = delete; - - /// Move constructor - Performance(Performance&& other) = delete; - - /// Copy assignment operator - Performance& operator=(const Performance& other) = delete; - - /// Move assignment operator - Performance& operator=(Performance&& other) = delete; - - - // *********************** - // ** Accessors ** - // *********************** - - /// Receives MC data object from MCDataManager - void ReceiveMCData(MCData&& mcData); - - - private: - // const L1Parameters* fpParameters = nullptr; ///< Instance of the tracking core class parameters - MCData fMCData = {}; ///< Object of MC data - }; - } // namespace tools -} // namespace ca - -#endif // CaToolsPerformance_h diff --git a/reco/L1/catools/CaToolsPerformance.cxx b/reco/L1/catools/CaToolsQa.cxx similarity index 70% rename from reco/L1/catools/CaToolsPerformance.cxx rename to reco/L1/catools/CaToolsQa.cxx index f814df8d3ec5ee83be98842d235e711d677ee146..1f693a8995d2b39251ccb3f4701999ee4ed8f41b 100644 --- a/reco/L1/catools/CaToolsPerformance.cxx +++ b/reco/L1/catools/CaToolsQa.cxx @@ -2,12 +2,12 @@ SPDX-License-Identifier: GPL-3.0-only Authors: Sergey Gorbunov, Sergei Zharko [committer] */ -/// \file CaToolsPerformance.h +/// \file CaToolsQa.h /// \brief Tracking performance class (implementation) /// \since 23.09.2022 /// \author S.Zharko <s.zharko@gsi.de> -#include "CaToolsPerformance.h" +#include "CaToolsQa.h" #include <utility> // for std::move @@ -15,7 +15,12 @@ using namespace ca::tools; // --------------------------------------------------------------------------------------------------------------------- // -void Performance::ReceiveMCData(MCData&& mcData) { fMCData = std::move(mcData); } +void Qa::InitHistograms() {} // --------------------------------------------------------------------------------------------------------------------- // +void Qa::FillHistograms() {} + +// --------------------------------------------------------------------------------------------------------------------- +// +void Qa::SaveHistograms() {} diff --git a/reco/L1/catools/CaToolsQa.h b/reco/L1/catools/CaToolsQa.h new file mode 100644 index 0000000000000000000000000000000000000000..87a58b31b96f14db37ca05894056e0b225cc293b --- /dev/null +++ b/reco/L1/catools/CaToolsQa.h @@ -0,0 +1,92 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov, Sergei Zharko [committer] */ + +/// \file CaToolsQa.h +/// \brief Tracking performance class (header) +/// \since 23.09.2022 +/// \author S.Zharko <s.zharko@gsi.de> + +#ifndef CaToolsQa_h +#define CaToolsQa_h 1 + +#include "L1Vector.h" + +class L1Parameters; +class TDirectory; +class TH1; +class TH1F; +class TH2F; +class TProfile; + +namespace ca::tools +{ + class MCData; + + /// Class ca::tools::Qa defines CA tracking internal performance measurement + /// + class Qa { + public: + // ***************************************** + // ** Constructors and destructor ** + // ***************************************** + + /// Default constructor + Qa() = default; + + /// Destructor + ~Qa() = default; + + /// Copy constructor + Qa(const Qa& other) = delete; + + /// Move constructor + Qa(Qa&& other) = delete; + + /// Copy assignment operator + Qa& operator=(const Qa& other) = delete; + + /// Move assignment operator + Qa& operator=(Qa&& other) = delete; + + /// Initializes histograms + /// \note should be called in the beginning of the run + void InitHistograms(); + + /// Fills histograms + /// \note should be called in the end of the run + void FillHistograms(); + + /// Saves histograms to the file + /// \note Should be called in the end of the run + void SaveHistograms(); + + // *********************** + // ** Accessors ** + // *********************** + + /// Sets pointer to MC data object + void SetMCData(const MCData* pMCData) { fpMCData = pMCData; } + + /// Sets pointer to tracking parameters object + void SetParameters(const L1Parameters* pParameters) { fpParameters = pParameters; } + + + private: + const L1Parameters* fpParameters = nullptr; ///< Instance of the tracking core class parameters + const MCData* fpMCData = nullptr; ///< Container for MC information on the event + + // *********************** + // ** Histograms ** + // *********************** + + TDirectory* fpHistoDirectory = nullptr; ///< directory object to store performance histograms + + TH1F* fphRegMcMom = nullptr; + TH1F* fphAccMCMom = nullptr; + + L1Vector<TH1*> fvpHistoList = {"ca::tools::Qa::fvpHistoList"}; ///< List of histograms for iteration over them + }; +} // namespace ca::tools + +#endif // CaToolsQa_h