diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h index 2da65db31d0eaa14f926c792cf6b465138794e16..c981b36076f98d068353f1414f98dcb4faffe171 100644 --- a/reco/L1/CbmL1.h +++ b/reco/L1/CbmL1.h @@ -356,14 +356,11 @@ private: void ReadEvent(CbmEvent* event = NULL); /// 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 iDet Index of the detector subsystem /// \param file Index of the input file /// \param event Index of the input event - /// \param iDet Index of the detector subsystem - /// \return flag: false - success, true - some errors occurred - // TODO: Probably, we should replace input parameter MVD with the template (S.Zharko) - bool ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int iDet); + /// \param iPoint Index of the point into the input MC points CbmMCDataArray object for the particular detector + void ReadMCPoint(L1DetectorID det, int iPoint, int file, int event); // static bool compareZ(const int &a, const int &b ); // bool compareZ(const int &a, const int &b ); diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx index bd3cf8d243b91393a29167924d6fb8785d6c5ec3..4568b0721258f533ad9464d027a23b137be7b2c0 100644 --- a/reco/L1/CbmL1Performance.cxx +++ b/reco/L1/CbmL1Performance.cxx @@ -1523,8 +1523,9 @@ void CbmL1::TrackFitPerformance() h_fit_chi2->Fill(it->chi2 / it->NDF); // last TRD point + /* do { - break; // only produce these plots in debug mode + break; // only produce these plots in debug mode if (it->GetNMCTracks() < 1) { break; } if (!fpTrdPoints) break; @@ -1575,7 +1576,7 @@ void CbmL1::TrackFitPerformance() tr.copyFromArrays(it->TLast, it->CLast); FillFitHistos(tr, mcP, isTimeFitted, h_fitTof); } while (0); // tof point - +*/ } // reco track } // void CbmL1::TrackFitPerformance() diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx index 6688a3095144d58c2e73f2788f44ce9f174923ff..b0f83de98b441c6d699db62f6bf2bee386951b8f 100644 --- a/reco/L1/CbmL1ReadEvent.cxx +++ b/reco/L1/CbmL1ReadEvent.cxx @@ -209,6 +209,7 @@ void CbmL1::ReadEvent(CbmEvent* event) fNpointsTrdAll = 0; fNpointsTofAll = 0; + // count total amount of mc points to set up the global indexing 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; @@ -222,37 +223,33 @@ void CbmL1::ReadEvent(CbmEvent* event) 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; + if (fUseMVD && fpMvdPoints) { - Int_t fNpointsMvdInEvent = fpMvdPoints->Size(iFile, iEvent); - for (Int_t iMC = 0; iMC < fNpointsMvdInEvent; iMC++) { - CbmL1MCPoint MC; - if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 0)) { fvMCPoints.push_back(MC); } + Int_t nMC = fpMvdPoints->Size(iFile, iEvent); + for (Int_t iMC = 0; iMC < nMC; iMC++) { + ReadMCPoint(L1DetectorID::kMvd, iFile, iEvent, iMC); } - } // fpMvdPoints + } if (fUseSTS && fpStsPoints) { - - Int_t nMC = fpStsPoints->Size(iFile, iEvent); + Int_t nMC = fpStsPoints->Size(iFile, iEvent); for (Int_t iMC = 0; iMC < nMC; iMC++) { - CbmL1MCPoint MC; - if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 1)) { fvMCPoints.push_back(MC); } + ReadMCPoint(L1DetectorID::kSts, iFile, iEvent, iMC); } - } // fpStsPoints + } if (fUseMUCH && fpMuchPoints) { Int_t nMC = fpMuchPoints->Size(iFile, iEvent); for (Int_t iMC = 0; iMC < nMC; iMC++) { - CbmL1MCPoint MC; - if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 2)) { fvMCPoints.push_back(MC); } + ReadMCPoint(L1DetectorID::kMuch, iFile, iEvent, iMC); } - } // fpMuchPoints + } if (fUseTRD && fpTrdPoints) { for (Int_t iMC = 0; iMC < fpTrdPoints->Size(iFile, iEvent); iMC++) { - CbmL1MCPoint MC; - if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 3)) { fvMCPoints.push_back(MC); } + ReadMCPoint(L1DetectorID::kTrd, iFile, iEvent, iMC); } - } // fpTrdPoints + } if (fUseTOF && fpTofPoints) { @@ -304,10 +301,7 @@ void CbmL1::ReadEvent(CbmEvent* event) int rpcAddr = pExtPoint->GetDetectorID() << kNofBitsRpcAddress; // New point if (rpcAddrCurr != rpcAddr || iTmcCurr != iTmc) { - if (iPointSelected != -1) { - CbmL1MCPoint MC; - if (!ReadMCPoint(&MC, iPointSelected, iFile, iEvent, 4)) { fvMCPoints.push_back(MC); } - } + if (iPointSelected != -1) { ReadMCPoint(L1DetectorID::kTof, iFile, iEvent, iPointSelected); } iTmcCurr = iTmc; rpcAddrCurr = rpcAddr; auto key = std::make_pair(iTmc, rpcAddr); @@ -332,15 +326,11 @@ void CbmL1::ReadEvent(CbmEvent* event) } } // Add the last point - if (iPointSelected != -1) { - CbmL1MCPoint MC; - if (!ReadMCPoint(&MC, iPointSelected, iFile, iEvent, 4)) { fvMCPoints.push_back(MC); } - } + if (iPointSelected != -1) { ReadMCPoint(L1DetectorID::kTof, iFile, iEvent, iPointSelected); } } // iP } } //time_slice - for (unsigned int iTr = 0; iTr < fvMCTracks.size(); iTr++) { sort(fvMCTracks[iTr].Points.begin(), fvMCTracks[iTr].Points.end(), compareMcPointZ); @@ -887,198 +877,146 @@ void CbmL1::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 iDet) +void CbmL1::ReadMCPoint(L1DetectorID iDet, int file, int event, int iPoint) { - TVector3 xyzI, PI, xyzO, PO; - Int_t mcID = -1; - Double_t time = 0.f; + FairMCPoint* fairPoint = nullptr; + TVector3 xyzO, PO; - double eventTime = fpMcEventList->GetEventTime(event, file); - int iStLoc = -1; - if (iDet == 0) { + int iStLoc = -1; + if (L1DetectorID::kMvd == iDet) { CbmMvdPoint* pt = L1_DYNAMIC_CAST<CbmMvdPoint*>(fpMvdPoints->Get(file, event, iPoint)); - - if (!pt) return 1; - pt->Position(xyzI); - pt->Momentum(PI); + assert(pt); + fairPoint = pt; pt->PositionOut(xyzO); pt->MomentumOut(PO); - mcID = pt->GetTrackID(); - time = eventTime + pt->GetTime(); iStLoc = CbmMvdTrackingInterface::Instance()->GetTrackingStationIndex(pt); } - if (iDet == 1) { + if (L1DetectorID::kSts == iDet) { CbmStsPoint* pt = L1_DYNAMIC_CAST<CbmStsPoint*>(fpStsPoints->Get(file, event, iPoint)); - if (!pt) return 1; - { - 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; - - if (EndTime > 0) - if (Time_MC_point > EndTime) return 1; - } - - pt->Position(xyzI); - pt->Momentum(PI); + assert(pt); + fairPoint = pt; pt->PositionOut(xyzO); pt->MomentumOut(PO); - mcID = pt->GetTrackID(); - time = eventTime + pt->GetTime(); - iStLoc = CbmStsTrackingInterface::Instance()->GetTrackingStationIndex(pt->GetDetectorID()); + iStLoc = CbmStsTrackingInterface::Instance()->GetTrackingStationIndex(pt); } - if (iDet == 2) { + if (L1DetectorID::kMuch == iDet) { CbmMuchPoint* pt = L1_DYNAMIC_CAST<CbmMuchPoint*>(fpMuchPoints->Get(file, event, iPoint)); - if (!pt) return 1; - { - 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; - - if (EndTime > 0) - if (Time_MC_point > EndTime) return 1; - } - - pt->Position(xyzI); - pt->Momentum(PI); + assert(pt); + fairPoint = pt; pt->PositionOut(xyzO); - pt->Momentum(PO); - mcID = pt->GetTrackID(); - time = eventTime + pt->GetTime(); - iStLoc = CbmMuchTrackingInterface::Instance()->GetTrackingStationIndex(pt->GetDetectorID()); + pt->MomentumOut(PO); + iStLoc = CbmMuchTrackingInterface::Instance()->GetTrackingStationIndex(pt); } - if (iDet == 3) { + if (L1DetectorID::kTrd == iDet) { CbmTrdPoint* pt = L1_DYNAMIC_CAST<CbmTrdPoint*>(fpTrdPoints->Get(file, event, iPoint)); - - if (!pt) return 1; - { - 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; - - if (EndTime > 0) - if (Time_MC_point > EndTime) return 1; - } - - pt->Position(xyzI); - pt->Momentum(PI); + assert(pt); + fairPoint = pt; pt->PositionOut(xyzO); pt->MomentumOut(PO); - mcID = pt->GetTrackID(); - - time = eventTime + pt->GetTime(); - iStLoc = CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(pt->GetDetectorID()); + iStLoc = CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(pt); } - if (iDet == 4) { + if (L1DetectorID::kTof == iDet) { CbmTofPoint* pt = L1_DYNAMIC_CAST<CbmTofPoint*>(fpTofPoints->Get(file, event, iPoint)); - if (!pt) return 1; - { - 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; - - if (EndTime > 0) - if (Time_MC_point > EndTime) return 1; - } - - pt->Position(xyzI); - pt->Momentum(PI); + assert(pt); + fairPoint = pt; pt->Position(xyzO); pt->Momentum(PO); - mcID = pt->GetTrackID(); - time = eventTime + pt->GetTime(); - iStLoc = CbmTofTrackingInterface::Instance()->GetTrackingStationIndex(pt->GetDetectorID()); + iStLoc = CbmTofTrackingInterface::Instance()->GetTrackingStationIndex(pt); } - MC->iStation = fpAlgo->GetParameters()->GetStationIndexActive(iStLoc, static_cast<L1DetectorID>(iDet)); - if (MC->iStation < 0) { return 1; } + CbmL1MCPoint MC; + + MC.iStation = fpAlgo->GetParameters()->GetStationIndexActive(iStLoc, iDet); + if (MC.iStation < 0) { return; } + + TVector3 xyzI, PI; + fairPoint->Position(xyzI); + fairPoint->Momentum(PI); TVector3 xyz = .5 * (xyzI + xyzO); TVector3 P = .5 * (PI + PO); - MC->x = xyz.X(); - MC->y = xyz.Y(); - MC->z = xyz.Z(); - MC->px = P.X(); - MC->py = P.Y(); - MC->pz = P.Z(); - MC->xIn = xyzI.X(); - MC->yIn = xyzI.Y(); - MC->zIn = xyzI.Z(); - MC->pxIn = PI.X(); - MC->pyIn = PI.Y(); - MC->pzIn = PI.Z(); - MC->xOut = xyzO.X(); - MC->yOut = xyzO.Y(); - MC->zOut = xyzO.Z(); - MC->pxOut = PO.X(); - MC->pyOut = PO.Y(); - MC->pzOut = PO.Z(); - MC->p = sqrt(fabs(MC->px * MC->px + MC->py * MC->py + MC->pz * MC->pz)); - auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(mcID, event, file)); + + MC.x = xyz.X(); + MC.y = xyz.Y(); + MC.z = xyz.Z(); + MC.px = P.X(); + MC.py = P.Y(); + MC.pz = P.Z(); + MC.xIn = xyzI.X(); + MC.yIn = xyzI.Y(); + MC.zIn = xyzI.Z(); + MC.pxIn = PI.X(); + MC.pyIn = PI.Y(); + MC.pzIn = PI.Z(); + MC.xOut = xyzO.X(); + MC.yOut = xyzO.Y(); + MC.zOut = xyzO.Z(); + MC.pxOut = PO.X(); + MC.pyOut = PO.Y(); + MC.pzOut = PO.Z(); + MC.p = sqrt(fabs(MC.px * MC.px + MC.py * MC.py + MC.pz * MC.pz)); + + CbmMCTrack* MCTrack = L1_DYNAMIC_CAST<CbmMCTrack*>(fpMcTracks->Get(file, event, fairPoint->GetTrackID())); + + auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(fairPoint->GetTrackID(), event, file)); assert(itTrack != fmMCTracksLinksMap.cend()); - MC->ID = itTrack->second; + MC.ID = itTrack->second; + + MC.pointId = iPoint; + MC.file = file; + MC.event = event; - MC->pointId = iPoint; - MC->file = file; - MC->event = event; + MC.time = fpMcEventList->GetEventTime(event, file) + fairPoint->GetTime(); - MC->time = time; + assert(MC.ID >= 0); + if (MC.ID < 0) return; - if (MC->ID < 0) return 1; - CbmMCTrack* MCTrack = L1_DYNAMIC_CAST<CbmMCTrack*>(fpMcTracks->Get(file, event, MC->ID)); - if (!MCTrack) return 1; - MC->pdg = MCTrack->GetPdgCode(); - MC->mother_ID = MCTrack->GetMotherId(); - TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(MC->pdg); + if (!MCTrack) { return; } + MC.pdg = MCTrack->GetPdgCode(); + MC.mother_ID = MCTrack->GetMotherId(); + + TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(MC.pdg); if (particlePDG) { - MC->q = particlePDG->Charge() / 3.0; - MC->mass = particlePDG->Mass(); + MC.q = particlePDG->Charge() / 3.0; + MC.mass = particlePDG->Mass(); } // END OF CUTS. Please, don't put any cuts below! - fvMCTracks[MC->ID].Points.push_back_no_warning(fvMCPoints.size()); int iPointInt = iPoint; switch (iDet) { - case 0: ++fNpointsMvd; break; - case 1: + case L1DetectorID::kMvd: { + ++fNpointsMvd; + break; + } + case L1DetectorID::kSts: iPointInt += fNpointsMvdAll; ++fNpointsSts; break; - case 2: + case L1DetectorID::kMuch: iPointInt += fNpointsMvdAll + fNpointsStsAll; ++fNpointsMuch; break; - case 3: + case L1DetectorID::kTrd: iPointInt += fNpointsMvdAll + fNpointsStsAll + fNpointsMuchAll; ++fNpointsTrd; break; - case 4: - iPointInt += fNpointsMvdAll + +fNpointsStsAll + fNpointsMuchAll + fNpointsTrdAll; + case L1DetectorID::kTof: + iPointInt += fNpointsMvdAll + fNpointsStsAll + fNpointsMuchAll + fNpointsTrdAll; ++fNpointsTof; break; + case L1DetectorID::kEND: break; } - fmMCPointsLinksMap[CbmL1LinkKey(iPointInt, event, file)] = fvMCPoints.size(); - return 0; + fmMCPointsLinksMap[CbmL1LinkKey(iPointInt, event, file)] = fvMCPoints.size(); + fvMCTracks[MC.ID].Points.push_back_no_warning(fvMCPoints.size()); + fvMCPoints.push_back(MC); }