diff --git a/reco/L1/CbmL1Constants.h b/reco/L1/CbmL1Constants.h index 6ef8bf3e8800eaab0a2753459583d9864b8580d1..cb17176fdaa22eabe896e9944056b9f5bcfebe0e 100644 --- a/reco/L1/CbmL1Constants.h +++ b/reco/L1/CbmL1Constants.h @@ -8,10 +8,10 @@ namespace CbmL1Constants { /// Performance constants - const double MinRecoMom = 0.1; // Extra set of tracks = (MinRecoMom, MinRefMom) - const double MinFastMom = 1.; // Primary set of tracks = (MinRefMom, +inf) //All reco tracks = (MinRecoMom, +inf) - const double MinPurity = 0.7; // min NStationInRecoTrack/NStationInMCTrack - const double MinNStations = 4.; // min number of stations in track to be reconstructable + const double MinRecoMom = 0.1; // Extra set of tracks = (MinRecoMom, MinRefMom) + const double MinFastMom = 1.; // Primary set of tracks = (MinRefMom, +inf) //All reco tracks = (MinRecoMom, +inf) + const double MinPurity = 0.7; // min NStationInRecoTrack/NStationInMCTrack + const int MinNStations = 4; // min number of stations on track to be reconstructable } // namespace CbmL1Constants #endif diff --git a/reco/L1/CbmL1MCTrack.cxx b/reco/L1/CbmL1MCTrack.cxx index d4df61048e0f883e8b41442653a8238b67d21304..6e52c5c939b322100184e09ecd1dc6b28542a01b 100644 --- a/reco/L1/CbmL1MCTrack.cxx +++ b/reco/L1/CbmL1MCTrack.cxx @@ -74,9 +74,8 @@ void CbmL1MCTrack::Init() } CalculateMCCont(); - CalculateHitCont(); + CountHitStations(); CalculateMaxNStaMC(); - CalculateMaxNStaHits(); CalculateIsReconstructable(); } @@ -113,66 +112,38 @@ void CbmL1MCTrack::CalculateMCCont() istaold = ista; } if (nMCContStations < ncont) nMCContStations = ncont; -}; // void CbmL1MCTrack::CalculateMCCont() - -void CbmL1MCTrack::CalculateHitCont() -{ - CbmL1* L1 = CbmL1::Instance(); - L1Algo* algo = L1->fpAlgo; - - int nhits = Hits.size(); - nHitContStations = 0; - int istaold = -1, ncont = 0; - { - for (int ih = 0; ih < nhits; ih++) { - int jh = Hits[ih]; - const L1Hit& h = algo->GetInputData().GetHit(jh); - int ista = h.iSt; - - if (ista - istaold == 1) ncont++; - else if (ista - istaold > 1) { - if (nHitContStations < ncont) nHitContStations = ncont; - ncont = 1; - } - - if (!(ista >= istaold)) { // tracks going in backward direction are not reconstructable - nHitContStations = 0; - return; - } - if (ista == istaold) continue; // backward direction - istaold = ista; - } - } - if (nHitContStations < ncont) nHitContStations = ncont; - +} -}; // void CbmL1MCTrack::CalculateHitCont() -void CbmL1MCTrack::CalculateMaxNStaHits() +void CbmL1MCTrack::CountHitStations() { CbmL1* L1 = CbmL1::Instance(); - maxNStaHits = 0; - nStations = 0; - int lastSta = -1; - int cur_maxNStaHits = 0; + int stationNhits[L1Constants::size::kMaxNstations] {0}; + for (unsigned int iH = 0; iH < Hits.size(); iH++) { CbmL1HitDebugInfo& sh = L1->fvHitDebugInfo[Hits[iH]]; - if (sh.iStation == lastSta) { cur_maxNStaHits++; } - else { // new station - if (!(sh.iStation > lastSta)) { // tracks going in backward direction are not reconstructable - maxNStaHits = 0; - return; - } - if (cur_maxNStaHits > maxNStaHits) maxNStaHits = cur_maxNStaHits; - cur_maxNStaHits = 1; - lastSta = sh.iStation; + stationNhits[sh.iStation]++; + } + + maxNStaHits = 0; + nStations = 0; + nHitContStations = 0; + int ncont = 0; + + for (int ista = 0; ista < L1->fpAlgo->GetParameters()->GetNstationsActive(); ista++) { + if (maxNStaHits < stationNhits[ista]) { maxNStaHits = stationNhits[ista]; } + if (stationNhits[ista] > 0) { nStations++; + ncont++; } - }; - if (cur_maxNStaHits > maxNStaHits) maxNStaHits = cur_maxNStaHits; - // cout << pdg << " " << p << " " << Hits.size() << " > " << maxNStaHits << endl; -}; // void CbmL1MCTrack::CalculateHitCont() + else { + if (nHitContStations < ncont) { nHitContStations = ncont; } + ncont = 0; + } + } + if (nHitContStations < ncont) nHitContStations = ncont; +} void CbmL1MCTrack::CalculateMaxNStaMC() { diff --git a/reco/L1/CbmL1MCTrack.h b/reco/L1/CbmL1MCTrack.h index 6dcac4bd72199ad97f644c2d630bcfe00deea5cf..330a4a20a6a10baeb8f03afa3c780b10be133b8e 100644 --- a/reco/L1/CbmL1MCTrack.h +++ b/reco/L1/CbmL1MCTrack.h @@ -69,8 +69,7 @@ public: private: void CalculateMCCont(); - void CalculateHitCont(); - void CalculateMaxNStaHits(); + void CountHitStations(); void CalculateMaxNStaMC(); void CalculateIsReconstructable(); @@ -89,6 +88,7 @@ public: int iFile = -1; int iEvent = -1; int mother_ID = -1; + int chainID = -1; // ID of the first particle in the decay chain int pdg = -1; unsigned int process_ID = (unsigned int) -1; bool isSignal {0}; diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx index 62e199a02c91af733f34891e36c8cb33d0fec9a3..777b344ae2c1f9568051eec25be342b1f7bfb1e8 100644 --- a/reco/L1/CbmL1Performance.cxx +++ b/reco/L1/CbmL1Performance.cxx @@ -54,6 +54,7 @@ #include <map> #include <vector> +#include "CaToolsDebugger.h" #include "L1Algo/L1Algo.h" #include "L1Algo/L1Def.h" #include "L1Algo/L1Fit.h" @@ -74,12 +75,12 @@ void CbmL1::TrackMatch() // fill pMCTrackMap for (vector<CbmL1MCTrack>::iterator i = fvMCTracks.begin(); i != fvMCTracks.end(); ++i) { CbmL1MCTrack& MC = *i; - if (pMCTrackMap.find(MC.ID) == pMCTrackMap.end()) { pMCTrackMap.insert(pair<int, CbmL1MCTrack*>(MC.ID, &MC)); } else { cout << "*** L1: Track ID " << MC.ID << " encountered twice! ***" << endl; } } + // -- prepare information about reconstructed tracks const int nRTracks = fvRecoTracks.size(); for (int iR = 0; iR < nRTracks; iR++) { @@ -91,6 +92,7 @@ void CbmL1::TrackMatch() // count how many hits from each mcTrack belong to current recoTrack map<int, int>& hitmap = prtra->hitMap; // how many hits from each mcTrack belong to current recoTrack + map<int, int> hitmapChain; for (vector<int>::iterator ih = (prtra->Hits).begin(); ih != (prtra->Hits).end(); ++ih) { const int nMCPoints = fvHitDebugInfo[*ih].mcPointIds.size(); @@ -98,12 +100,20 @@ void CbmL1::TrackMatch() int iMC = fvHitDebugInfo[*ih].mcPointIds[iP]; // cout<<iMC<<" iMC"<<endl; - int ID = -1; - if (iMC >= 0) ID = fvMCPoints[iMC].ID; + int ID = -1; + int chainID = -1; + if (iMC >= 0) { + ID = fvMCPoints[iMC].ID; + chainID = fvMCTracks[ID].chainID; + } if (hitmap.find(ID) == hitmap.end()) hitmap[ID] = 1; else { hitmap[ID] += 1; } + if (hitmapChain.find(chainID) == hitmapChain.end()) hitmapChain[chainID] = 1; + else { + hitmapChain[chainID] += 1; + } } // for iPoint } // for iHit @@ -130,6 +140,21 @@ void CbmL1::TrackMatch() } } // for hitmap prtra->SetMaxPurity(max_percent); + + if (0) { // debug + + double max_percentChain = 0.0; // [%]. maximum persent of hits, which belong to one mcTrack. + for (map<int, int>::iterator posIt = hitmapChain.begin(); posIt != hitmapChain.end(); + posIt++) { // loop over all touched MCTracks + if (posIt->first < 0) continue; // not a MC track - based on fake hits + // count max-purity + if (double(posIt->second) > max_percentChain * double(hitsum)) { + max_percentChain = double(posIt->second) / double(hitsum); + } + } // for hitmap + prtra->SetMaxPurity(max_percentChain); + } + } // for reco tracks } // @@ -291,10 +316,16 @@ void CbmL1::EfficienciesPerformance() ntra.fOutDir = fTableDir; // Setup a pointer for output directory L1_NTRA.fOutDir = fTableDir; // Save average efficiencies + ca::tools::Debugger::Instance().AddNtuple("ghost", "it:ih:p:x:y:z:t:dx:dy"); + static int statNghost = 0; + for (vector<CbmL1Track>::iterator rtraIt = fvRecoTracks.begin(); rtraIt != fvRecoTracks.end(); ++rtraIt) { ntra.ghosts += rtraIt->IsGhost(); if (0 && rtraIt->IsGhost()) { // Debug. - cout << " L1: ghost track: nhits " << rtraIt->GetNOfHits() << " p " << 1. / rtraIt->T[5] << " purity " + L1TrackPar tr; + tr.copyFromArrays(rtraIt->T, rtraIt->C); + + cout << " L1: ghost track: nhits " << rtraIt->GetNOfHits() << " p " << 1. / rtraIt->T[4] << " purity " << rtraIt->GetMaxPurity() << " | hits "; for (map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++) { cout << " (" << posIt->second << " " << posIt->first << ") "; @@ -302,9 +333,18 @@ void CbmL1::EfficienciesPerformance() cout << endl; for (map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++) { CbmL1MCTrack& t = fvMCTracks[posIt->first]; - cout << "mc " << posIt->first << " pdg " << t.pdg << " mother: " << t.mother_ID; + cout << "mc " << posIt->first << " pdg " << t.pdg << " mother: " << t.mother_ID << " chain " << t.chainID; cout << " n mc stations: " << t.NMCStations() << endl; } + for (unsigned int i = 0; i < rtraIt->Hits.size(); i++) { + const L1Hit& h = fpAlgo->fInputData.GetHit(rtraIt->Hits[i]); + const CbmL1HitDebugInfo& s = fvHitDebugInfo[rtraIt->Hits[i]]; + cout << " x y z t " << s.x << " " << s.y << " " << h.z << " dx " << s.dx << " dy " << s.dy << std::endl; + ca::tools::Debugger::Instance().FillNtuple("ghost", statNghost, i, fabs(1. / tr.qp[0]), s.x, s.y, h.z, h.t, + s.dx, s.dy); + } + tr.Print(0); + statNghost++; } } @@ -798,7 +838,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa h_ghost_chi2->Fill(prtra->chi2 / prtra->NDF); h_ghost_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF)); } - else { + else if (prtra->GetNMCTracks() > 0) { if (prtra->GetMCTrack()[0].IsReconstructable()) { h_reco_chi2->Fill(prtra->chi2 / prtra->NDF); h_reco_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF)); @@ -1321,7 +1361,10 @@ void CbmL1::TrackFitPerformance() FillFitHistos(tr, mcP, isTimeFitted, h_fitL); } - { // vertex + do { // vertex + + if (it->GetNMCTracks() < 1) { break; } + CbmL1MCTrack mc = *(it->GetMCTracks()[0]); fit.SetTrack(it->T, it->C); @@ -1452,13 +1495,15 @@ void CbmL1::TrackFitPerformance() h_fitPV[16]->Fill(tr.vi[0] * L1TrackPar::kClightNs); } } - } + } while (0); h_fit_chi2->Fill(it->chi2 / it->NDF); // last TRD point do { break; // only produce these plots in debug mode + if (it->GetNMCTracks() < 1) { break; } + if (!fpTrdPoints) break; const CbmL1MCTrack& mcTrack = *(it->GetMCTracks()[0]); int nTrdPoints = fpTrdPoints->Size(mcTrack.iFile, mcTrack.iEvent); @@ -1484,7 +1529,9 @@ void CbmL1::TrackFitPerformance() // last TOF point do { break; // only produce these plots in debug mode - // only produce these plots in debug mode + + if (it->GetNMCTracks() < 1) { break; } + if (!fpTofPoints) break; const CbmL1MCTrack& mcTrack = *(it->GetMCTracks()[0]); int nTofPoints = fpTofPoints->Size(mcTrack.iFile, mcTrack.iEvent); @@ -2078,7 +2125,6 @@ void CbmL1::InputPerformance() const CbmTrdHit* sh = L1_DYNAMIC_CAST<CbmTrdHit*>(fpTrdHits->At(iH)); const CbmMatch* hm = L1_DYNAMIC_CAST<CbmMatch*>(fpTrdHitMatches->At(iH)); - if (hm->GetNofLinks() == 0) continue; if (hm->GetNofLinks() != 1) continue; // only check single-linked hits diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx index 62aa75bd8c0ece89f4a0f4134f9d9cd9a0c10f1a..a6e37619fc3a31a38866f67687c23e4f834dea8c 100644 --- a/reco/L1/CbmL1ReadEvent.cxx +++ b/reco/L1/CbmL1ReadEvent.cxx @@ -695,7 +695,7 @@ void CbmL1::ReadEvent(CbmEvent* event) th.v = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; } th.Det = 0; - th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMvd>(th.ExtIndex) : -1; + th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMvd>(hitIndex) : -1; if (1 == fMvdUseMcHit && th.iMC >= 0) { th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation)); } @@ -984,9 +984,9 @@ void CbmL1::ReadEvent(CbmEvent* event) std::tie(th.u, th.v) = st.ConvXYtoUV<double>(th.x, th.y); - th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kTrd>(th.ExtIndex) : -1; + th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kTrd>(hitIndex) : -1; - if (1 == fTrdUseMcHit) { // replace hit by MC points + if (1 == fTrdUseMcHit) { // replace hit by MC point assert(fPerformance && fpTrdHitMatches); @@ -995,11 +995,11 @@ void CbmL1::ReadEvent(CbmEvent* event) assert(iMcTrd >= 0 && iMcTrd < fNpointsTrd); if (mcUsed[iMcTrd]) continue; // one hit per MC point bool doSmear = true; - if (0) { // SGtrd2d!! debug: artificial errors - th.dx = .1; - th.du = .1; - th.dy = .1; - th.dv = .1; + if (1) { // SGtrd2d!! set artificial errors when the input errors are unrealistically large + if (th.dx > 0.3) { th.dx = 0.3; } + if (th.dy > 0.2) { th.dy = 0.2; } + th.du = th.dx; + th.dv = th.dy; } th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation), doSmear); mcUsed[iMcTrd] = 1; @@ -1099,7 +1099,7 @@ void CbmL1::ReadEvent(CbmEvent* event) th.u = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; th.v = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; - th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kTof>(th.ExtIndex) : -1; + th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kTof>(hitIndex) : -1; if (1 == fTofUseMcHit && th.iMC > -1) { th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation)); @@ -1235,9 +1235,8 @@ void CbmL1::Fill_vMCTracks() fvMCTracks.reserve(nMCTracks); } - int fileEvent = 0; /* Loop over MC event */ - for (DFSET::iterator set_it = fvSelectedMcEvents.begin(); set_it != fvSelectedMcEvents.end(); ++set_it, ++fileEvent) { + for (DFSET::iterator set_it = fvSelectedMcEvents.begin(); set_it != fvSelectedMcEvents.end(); ++set_it) { Int_t iFile = set_it->first; Int_t iEvent = set_it->second; @@ -1254,6 +1253,14 @@ void CbmL1::Fill_vMCTracks() int mother_ID = MCTrack->GetMotherId(); + int chainID = iMCTrack; + int chainParent = mother_ID; + while (chainParent >= 0) { + chainID = chainParent; + CbmMCTrack* parent = L1_DYNAMIC_CAST<CbmMCTrack*>(fpMcTracks->Get(iFile, iEvent, chainParent)); + chainParent = parent->GetMotherId(); + } + if (mother_ID < 0 && mother_ID != -2) mother_ID = -iEvent - 1; TVector3 vr; TLorentzVector vp; @@ -1272,6 +1279,8 @@ void CbmL1::Fill_vMCTracks() Int_t iTrack = fvMCTracks.size(); //or iMCTrack? CbmL1MCTrack T(mass, q, vr, vp, iTrack, mother_ID, pdg, processID); + T.chainID = chainID; + // CbmL1MCTrack T(mass, q, vr, vp, iMCTrack, mother_ID, pdg); T.time = eventTime + MCTrack->GetStartT(); T.iFile = iFile; @@ -1283,8 +1292,8 @@ void CbmL1::Fill_vMCTracks() fmMCTracksLinksMap[CbmL1LinkKey(iMCTrack, iEvent, iFile)] = fvMCTracks.size() - 1; } //iTracks } //Links - } //Fill_vMCTracks + // //-------------------------------------------------------------------------------------------------- // @@ -1298,8 +1307,7 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i double eventTime = fpMcEventList->GetEventTime(event, file); if (iDet == 0) { - CbmMvdPoint* pt = L1_DYNAMIC_CAST<CbmMvdPoint*>(fpMvdPoints->Get(file, event, iPoint)); // file, event, object - //CbmMvdPoint *pt = L1_DYNAMIC_CAST<CbmMvdPoint*> (Point); + CbmMvdPoint* pt = L1_DYNAMIC_CAST<CbmMvdPoint*>(fpMvdPoints->Get(file, event, iPoint)); if (!pt) return 1; pt->Position(xyzI); @@ -1311,7 +1319,7 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i } if (iDet == 1) { - CbmStsPoint* pt = L1_DYNAMIC_CAST<CbmStsPoint*>(fpStsPoints->Get(file, event, iPoint)); // file, event, object + CbmStsPoint* pt = L1_DYNAMIC_CAST<CbmStsPoint*>(fpStsPoints->Get(file, event, iPoint)); if (!pt) return 1; { Double_t StartTime = fTimeSlice->GetStartTime(); @@ -1335,7 +1343,7 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i if (iDet == 2) { - CbmMuchPoint* pt = L1_DYNAMIC_CAST<CbmMuchPoint*>(fpMuchPoints->Get(file, event, iPoint)); // file, event, object + CbmMuchPoint* pt = L1_DYNAMIC_CAST<CbmMuchPoint*>(fpMuchPoints->Get(file, event, iPoint)); if (!pt) return 1; { Double_t StartTime = fTimeSlice->GetStartTime(); @@ -1359,13 +1367,13 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i if (iDet == 3) { - CbmTrdPoint* pt = L1_DYNAMIC_CAST<CbmTrdPoint*>(fpTrdPoints->Get(file, event, iPoint)); // file, event, object + CbmTrdPoint* pt = L1_DYNAMIC_CAST<CbmTrdPoint*>(fpTrdPoints->Get(file, event, iPoint)); if (!pt) return 1; { - Double_t StartTime = fTimeSlice->GetStartTime(); // not used - Double_t EndTime = fTimeSlice->GetEndTime(); // not used - Double_t Time_MC_point = eventTime + pt->GetTime(); // not used + Double_t StartTime = fTimeSlice->GetStartTime(); + Double_t EndTime = fTimeSlice->GetEndTime(); + Double_t Time_MC_point = eventTime + pt->GetTime(); if (StartTime > 0) if (Time_MC_point < StartTime) return 1; @@ -1384,7 +1392,7 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i } if (iDet == 4) { - CbmTofPoint* pt = L1_DYNAMIC_CAST<CbmTofPoint*>(fpTofPoints->Get(file, event, iPoint)); // file, event, object + CbmTofPoint* pt = L1_DYNAMIC_CAST<CbmTofPoint*>(fpTofPoints->Get(file, event, iPoint)); if (!pt) return 1; { Double_t StartTime = fTimeSlice->GetStartTime(); @@ -1464,5 +1472,5 @@ void CbmL1::HitMatch() hit.mcPointIds.push_back_no_warning(iP); fvMCPoints[iP].hitIds.push_back_no_warning(iH); } - } // for hits + } } diff --git a/reco/L1/CbmL1Track.h b/reco/L1/CbmL1Track.h index e5275831aa65903f7d64253ff4ca9231959db4ea..8ba6383d035b103239486d66ecf5fdd35af19317 100644 --- a/reco/L1/CbmL1Track.h +++ b/reco/L1/CbmL1Track.h @@ -22,6 +22,7 @@ #ifndef CbmL1Track_H #define CbmL1Track_H +#include "CbmL1Constants.h" #include "CbmL1MCTrack.h" #include "CbmL1TrackPar.h" @@ -42,7 +43,7 @@ public: vector<CbmL1MCTrack*>& GetMCTracks() { return mcTracks; } CbmL1MCTrack* GetMCTrack() { return mcTracks[0]; } int GetNMCTracks() { return mcTracks.size(); } - bool IsGhost() { return !(mcTracks.size()); } + bool IsGhost() { return (maxPurity < CbmL1Constants::MinPurity); } void SetMaxPurity(double maxPurity_) { maxPurity = maxPurity_; } double GetMaxPurity() { return maxPurity; } @@ -64,7 +65,7 @@ public: map<int, int> hitMap; // how many hits from each mcTrack belong to current recoTrack private: // next members filled and used in Performance - vector<CbmL1MCTrack*> mcTracks; // array of assosiated recoTracks. Should be only one. + vector<CbmL1MCTrack*> mcTracks; // array of assosiated mc Tracks. Should be only one. double maxPurity; // maximum persent of hits, which belong to one mcTrack. };