diff --git a/reco/KF/CbmKFParticleFinderQa.cxx b/reco/KF/CbmKFParticleFinderQa.cxx index 418d76c4ed39ddca0a2c95eb8fd02caf53b21687..851c894d3e3e34f623b170f62a25b53e768cc498 100644 --- a/reco/KF/CbmKFParticleFinderQa.cxx +++ b/reco/KF/CbmKFParticleFinderQa.cxx @@ -45,30 +45,7 @@ using std::vector; CbmKFParticleFinderQa::CbmKFParticleFinderQa(const char* name, Int_t iVerbose, const KFParticleTopoReconstructor* tr, TString outFileName) : FairTask(name, iVerbose) - , fMCTracksBranchName("MCTrack") - , fTrackMatchBranchName("StsTrackMatch") - , fMCTrackArray(nullptr) - , fMCTrackArrayEvent(nullptr) - , fEventList(nullptr) - , fTrackMatchArray(nullptr) - , fRecParticles(nullptr) - , fMCParticles(nullptr) - , fMatchParticles(nullptr) - , fSaveParticles(false) - , fSaveMCParticles(false) - , fTimeSliceMode(false) , fOutFileName(outFileName) - , fOutFile(nullptr) - , fEfffileName("Efficiency.txt") - , fTopoPerformance(nullptr) - , fPrintFrequency(100) - , fNEvents(0) - , fTime() - , fSuperEventAnalysis(false) - , fReferenceResults("./") - , fDecayToAnalyse(-1) - , fCheckDecayQA(false) - , fTestOk(false) { for (Int_t i = 0; i < 5; i++) fTime[i] = 0; @@ -80,8 +57,9 @@ CbmKFParticleFinderQa::CbmKFParticleFinderQa(const char* name, Int_t iVerbose, c TDirectory* curDirectory = gDirectory; if (!(fOutFileName == "")) fOutFile = new TFile(fOutFileName.Data(), "RECREATE"); - else + else { fOutFile = gFile; + } fTopoPerformance->CreateHistos("KFTopoReconstructor", fOutFile, tr->GetKFParticleFinder()->GetReconstructionList()); gFile = curFile; @@ -109,12 +87,12 @@ InitStatus CbmKFParticleFinderQa::Init() return kERROR; } - //check the mode - fTimeSliceMode = 0; - if (ioman->CheckBranch("CbmEvent")) fTimeSliceMode = 1; + // check the mode + fLegacyEventMode = 0; + if (!ioman->CheckBranch("CbmEvent") && ioman->CheckBranch("MCTrack")) { fLegacyEventMode = 1; } - //MC Tracks - if (fTimeSliceMode) { + // MC Tracks + if (!fLegacyEventMode) { FairRootManager* fManger = FairRootManager::Instance(); CbmMCDataManager* mcManager = (CbmMCDataManager*) fManger->GetObject("MCDataManager"); if (mcManager == 0) Error("CbmKFParticleFinderQa::Init", "MC Data Manager not found!"); @@ -126,18 +104,21 @@ InitStatus CbmKFParticleFinderQa::Init() return kERROR; } - fEventList = (CbmMCEventList*) ioman->GetObject("MCEventList."); - if (fEventList == 0) { + fMcEventList = (CbmMCEventList*) ioman->GetObject("MCEventList."); + if (fMcEventList == 0) { Error("CbmKFParticleFinderQa::Init", "MC Event List not found!"); return kERROR; } + + fRecoEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent")); + if (nullptr == fRecoEvents) { LOG(warn) << GetName() << ": No event array! Will process entire tree."; } } else { fMCTrackArrayEvent = (TClonesArray*) ioman->GetObject("MCTrack"); } - //Track match - fTrackMatchArray = (TClonesArray*) ioman->GetObject(fTrackMatchBranchName); + // Track match + fTrackMatchArray = (TClonesArray*) ioman->GetObject("StsTrackMatch"); if (fTrackMatchArray == 0) { Error("CbmKFParticleFinderQa::Init", "track match array not found!"); return kERROR; @@ -158,181 +139,210 @@ InitStatus CbmKFParticleFinderQa::Init() fMatchParticles = new TClonesArray("KFParticleMatch", 100); ioman->Register("KFParticleMatch", "KFParticle", fMatchParticles, IsOutputBranchPersistent("KFParticleMatch")); } + return kSUCCESS; } void CbmKFParticleFinderQa::Exec(Option_t* /*opt*/) { - if (!fSuperEventAnalysis) { - if (fSaveParticles) fRecParticles->Delete(); - if (fSaveMCParticles) { - fMCParticles->Delete(); - fMatchParticles->Delete(); + if (fSuperEventAnalysis) { + LOG(error) << GetName() << " SuperEventAnalysis option currently doesn't work"; + return; + } + + if (fSaveParticles) fRecParticles->Delete(); + if (fSaveMCParticles) { + fMCParticles->Delete(); + fMatchParticles->Delete(); + } + + // make sure that the number of events in time slice is 1 + + int nRecoEvents = 1; + int nMCEvents = 1; + + if (!fLegacyEventMode) { + if (fRecoEvents) { nRecoEvents = fRecoEvents->GetEntriesFast(); } + nMCEvents = fMcEventList->GetNofEvents(); + } + + if (nMCEvents != 1 || nRecoEvents != 1) { + LOG(warning) << GetName() << " the task doesn't properly work with more than one event in the time slice"; + } + + vector<KFMCTrack> mcTracks; + vector<int> mcToKFmcMap[nMCEvents]; + + int mcIndexOffset = 0; + + for (int iMCEvent = 0; iMCEvent < nMCEvents; iMCEvent++) { + + CbmLink mcEventLink = fMcEventList->GetEventLinkByIndex(iMCEvent); + int nMCTracks = 0; + if (!fLegacyEventMode) { + nMCTracks = fMCTrackArray->Size(mcEventLink); + mcToKFmcMap[iMCEvent].resize(nMCTracks, -1); + } + else { + nMCTracks = fMCTrackArrayEvent->GetEntriesFast(); } - int nMCEvents = 1; - if (fTimeSliceMode) nMCEvents = fEventList->GetNofEvents(); + for (Int_t iMC = 0; iMC < nMCTracks; iMC++) { - vector<KFMCTrack> mcTracks; - vector<vector<vector<int>>> indexMap(1); - indexMap[0].resize(nMCEvents); + CbmMCTrack* cbmMCTrack; + CbmLink mcTrackLink = mcEventLink; + mcTrackLink.SetIndex(iMC); + if (!fLegacyEventMode) { cbmMCTrack = (CbmMCTrack*) fMCTrackArray->Get(mcTrackLink); } + else { + cbmMCTrack = (CbmMCTrack*) fMCTrackArrayEvent->At(iMC); + } - int mcIndexOffset = 0; + KFMCTrack kfMCTrack; + kfMCTrack.SetX(cbmMCTrack->GetStartX()); + kfMCTrack.SetY(cbmMCTrack->GetStartY()); + kfMCTrack.SetZ(cbmMCTrack->GetStartZ()); + kfMCTrack.SetPx(cbmMCTrack->GetPx()); + kfMCTrack.SetPy(cbmMCTrack->GetPy()); + kfMCTrack.SetPz(cbmMCTrack->GetPz()); + + Int_t pdg = cbmMCTrack->GetPdgCode(); + Double_t q = 1; + if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) + q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0; + else if (pdg == 1000010020) + q = 1; + else if (pdg == -1000010020) + q = -1; + else if (pdg == 1000010030) + q = 1; + else if (pdg == -1000010030) + q = -1; + else if (pdg == 1000020030) + q = 2; + else if (pdg == -1000020030) + q = -2; + else if (pdg == 1000020040) + q = 2; + else if (pdg == -1000020040) + q = -2; + else + q = 0; + Double_t p = cbmMCTrack->GetP(); - for (int iMCEvent = 0; iMCEvent < nMCEvents; iMCEvent++) { - int nMCTracks = 0; - if (fTimeSliceMode) nMCTracks = fMCTrackArray->Size(0, iMCEvent); + if (cbmMCTrack->GetMotherId() >= 0) kfMCTrack.SetMotherId(cbmMCTrack->GetMotherId() + mcIndexOffset); else - nMCTracks = fMCTrackArrayEvent->GetEntriesFast(); - - if (fTimeSliceMode) indexMap[0][iMCEvent].resize(nMCTracks); - - for (Int_t iMC = 0; iMC < nMCTracks; iMC++) { - CbmMCTrack* cbmMCTrack; - if (fTimeSliceMode) cbmMCTrack = (CbmMCTrack*) fMCTrackArray->Get(0, iMCEvent, iMC); - else - cbmMCTrack = (CbmMCTrack*) fMCTrackArrayEvent->At(iMC); - - KFMCTrack kfMCTrack; - kfMCTrack.SetX(cbmMCTrack->GetStartX()); - kfMCTrack.SetY(cbmMCTrack->GetStartY()); - kfMCTrack.SetZ(cbmMCTrack->GetStartZ()); - kfMCTrack.SetPx(cbmMCTrack->GetPx()); - kfMCTrack.SetPy(cbmMCTrack->GetPy()); - kfMCTrack.SetPz(cbmMCTrack->GetPz()); - - Int_t pdg = cbmMCTrack->GetPdgCode(); - Double_t q = 1; - if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) - q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0; - else if (pdg == 1000010020) - q = 1; - else if (pdg == -1000010020) - q = -1; - else if (pdg == 1000010030) - q = 1; - else if (pdg == -1000010030) - q = -1; - else if (pdg == 1000020030) - q = 2; - else if (pdg == -1000020030) - q = -2; - else if (pdg == 1000020040) - q = 2; - else if (pdg == -1000020040) - q = -2; - else - q = 0; - Double_t p = cbmMCTrack->GetP(); - - if (cbmMCTrack->GetMotherId() >= 0) kfMCTrack.SetMotherId(cbmMCTrack->GetMotherId() + mcIndexOffset); - else - kfMCTrack.SetMotherId(-iMCEvent - 1); - kfMCTrack.SetQP(q / p); - kfMCTrack.SetPDG(pdg); - kfMCTrack.SetNMCPoints(0); - - if (fTimeSliceMode) indexMap[0][iMCEvent][iMC] = mcTracks.size(); - mcTracks.push_back(kfMCTrack); - } + kfMCTrack.SetMotherId(-iMCEvent - 1); + kfMCTrack.SetQP(q / p); + kfMCTrack.SetPDG(pdg); + kfMCTrack.SetNMCPoints(0); - mcIndexOffset += nMCTracks; + if (!fLegacyEventMode) { mcToKFmcMap[iMCEvent][iMC] = mcTracks.size(); } + mcTracks.push_back(kfMCTrack); } - Int_t ntrackMatches = fTrackMatchArray->GetEntriesFast(); + mcIndexOffset += nMCTracks; + } - vector<int> trackMatch(ntrackMatches, -1); + Int_t ntrackMatches = fTrackMatchArray->GetEntriesFast(); + vector<int> trackMatch(ntrackMatches, -1); + + for (int iTr = 0; iTr < ntrackMatches; iTr++) { + + CbmTrackMatchNew* stsTrackMatch = (CbmTrackMatchNew*) fTrackMatchArray->At(iTr); - for (int iTr = 0; iTr < ntrackMatches; iTr++) { - CbmTrackMatchNew* stsTrackMatch = (CbmTrackMatchNew*) fTrackMatchArray->At(iTr); - - if (stsTrackMatch->GetNofLinks() == 0) continue; - Float_t bestWeight = 0.f; - Float_t totalWeight = 0.f; - Int_t mcTrackId = -1; - CbmLink link; - for (int iLink = 0; iLink < stsTrackMatch->GetNofLinks(); iLink++) { - totalWeight += stsTrackMatch->GetLink(iLink).GetWeight(); - if (stsTrackMatch->GetLink(iLink).GetWeight() > bestWeight) { - bestWeight = stsTrackMatch->GetLink(iLink).GetWeight(); - int iMCTrack = stsTrackMatch->GetLink(iLink).GetIndex(); - link = stsTrackMatch->GetLink(iLink); - - if (fTimeSliceMode) mcTrackId = indexMap[link.GetFile()][link.GetEntry()][iMCTrack]; - else - mcTrackId = stsTrackMatch->GetLink(iLink).GetIndex(); + if (stsTrackMatch->GetNofLinks() == 0) continue; + Float_t bestWeight = 0.f; + Float_t totalWeight = 0.f; + Int_t mcTrackId = -1; + CbmLink link; + + for (int iLink = 0; iLink < stsTrackMatch->GetNofLinks(); iLink++) { + totalWeight += stsTrackMatch->GetLink(iLink).GetWeight(); + if (stsTrackMatch->GetLink(iLink).GetWeight() > bestWeight) { + bestWeight = stsTrackMatch->GetLink(iLink).GetWeight(); + int iMCTrack = stsTrackMatch->GetLink(iLink).GetIndex(); + link = stsTrackMatch->GetLink(iLink); + if (!fLegacyEventMode) { + int eventIndex = fMcEventList->GetEventIndex(link); + if (eventIndex >= 0) { mcTrackId = mcToKFmcMap[eventIndex][iMCTrack]; } + } + else { + mcTrackId = stsTrackMatch->GetLink(iLink).GetIndex(); } } - if (bestWeight / totalWeight < 0.7) continue; - // if(mcTrackId >= nMCTracks || mcTrackId < 0) - // { - // std::cout << "Sts Matching is wrong! StsTackId = " << mcTrackId << " N mc tracks = " << nMCTracks << std::endl; - // continue; - // } - - if (TMath::Abs(mcTracks[mcTrackId].PDG()) > 4000 - && !(TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000010020 - || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000010030 - || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000020030 - || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000020040)) - continue; - mcTracks[mcTrackId].SetReconstructed(); - trackMatch[iTr] = mcTrackId; } - fTopoPerformance->SetMCTracks(mcTracks); - fTopoPerformance->SetTrackMatch(trackMatch); + if (mcTrackId < 0) continue; - fTopoPerformance->CheckMCTracks(); - fTopoPerformance->MatchTracks(); - fTopoPerformance->FillHistos(); + if (bestWeight / totalWeight < 0.7) continue; + // if(mcTrackId >= nMCTracks || mcTrackId < 0) + // { + // std::cout << "Sts Matching is wrong! StsTackId = " << mcTrackId << " N mc tracks = " << nMCTracks << std::endl; + // continue; + // } - fNEvents++; - fTime[4] += fTopoPerformance->GetTopoReconstructor()->Time(); - for (int iT = 0; iT < 4; iT++) - fTime[iT] += fTopoPerformance->GetTopoReconstructor()->StatTime(iT); - if (fNEvents % fPrintFrequency == 0) { - std::cout << "Topo reconstruction time" - << " Real = " << std::setw(10) << fTime[4] / fNEvents * 1.e3 << " ms" << std::endl; - std::cout << " Init " << fTime[0] / fNEvents * 1.e3 << " ms" << std::endl; - std::cout << " PV Finder " << fTime[1] / fNEvents * 1.e3 << " ms" << std::endl; - std::cout << " Sort Tracks " << fTime[2] / fNEvents * 1.e3 << " ms" << std::endl; - std::cout << " KF Particle Finder " << fTime[3] / fNEvents * 1.e3 << " ms" << std::endl; - } + if (TMath::Abs(mcTracks[mcTrackId].PDG()) > 4000 + && !(TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000010020 || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000010030 + || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000020030 + || TMath::Abs(mcTracks[mcTrackId].PDG()) == 1000020040)) + continue; - // save particles to a ROOT file - if (fSaveParticles) { - for (unsigned int iP = 0; iP < fTopoPerformance->GetTopoReconstructor()->GetParticles().size(); iP++) { - new ((*fRecParticles)[iP]) KFParticle(fTopoPerformance->GetTopoReconstructor()->GetParticles()[iP]); - } + mcTracks[mcTrackId].SetReconstructed(); + trackMatch[iTr] = mcTrackId; + } + + fTopoPerformance->SetMCTracks(mcTracks); + fTopoPerformance->SetTrackMatch(trackMatch); + + fTopoPerformance->CheckMCTracks(); + fTopoPerformance->MatchTracks(); + fTopoPerformance->FillHistos(); + + fNEvents++; + fTime[4] += fTopoPerformance->GetTopoReconstructor()->Time(); + for (int iT = 0; iT < 4; iT++) + fTime[iT] += fTopoPerformance->GetTopoReconstructor()->StatTime(iT); + if (fNEvents % fPrintFrequency == 0) { + std::cout << "Topo reconstruction time" + << " Real = " << std::setw(10) << fTime[4] / fNEvents * 1.e3 << " ms" << std::endl; + std::cout << " Init " << fTime[0] / fNEvents * 1.e3 << " ms" << std::endl; + std::cout << " PV Finder " << fTime[1] / fNEvents * 1.e3 << " ms" << std::endl; + std::cout << " Sort Tracks " << fTime[2] / fNEvents * 1.e3 << " ms" << std::endl; + std::cout << " KF Particle Finder " << fTime[3] / fNEvents * 1.e3 << " ms" << std::endl; + } + + // save particles to a ROOT file + if (fSaveParticles) { + for (unsigned int iP = 0; iP < fTopoPerformance->GetTopoReconstructor()->GetParticles().size(); iP++) { + new ((*fRecParticles)[iP]) KFParticle(fTopoPerformance->GetTopoReconstructor()->GetParticles()[iP]); } + } - if (fSaveMCParticles) { - for (unsigned int iP = 0; iP < fTopoPerformance->GetTopoReconstructor()->GetParticles().size(); iP++) { - new ((*fMatchParticles)[iP]) KFParticleMatch(); - KFParticleMatch* p = (KFParticleMatch*) (fMatchParticles->At(iP)); - - Short_t matchType = 0; - int iMCPart = -1; - if (!(fTopoPerformance->ParticlesMatch()[iP].IsMatchedWithPdg())) //background - { - if (fTopoPerformance->ParticlesMatch()[iP].IsMatched()) { - iMCPart = fTopoPerformance->ParticlesMatch()[iP].GetBestMatchWithPdg(); - matchType = 1; - } - } - else { + if (fSaveMCParticles) { + for (unsigned int iP = 0; iP < fTopoPerformance->GetTopoReconstructor()->GetParticles().size(); iP++) { + new ((*fMatchParticles)[iP]) KFParticleMatch(); + KFParticleMatch* p = (KFParticleMatch*) (fMatchParticles->At(iP)); + + Short_t matchType = 0; + int iMCPart = -1; + if (!(fTopoPerformance->ParticlesMatch()[iP].IsMatchedWithPdg())) //background + { + if (fTopoPerformance->ParticlesMatch()[iP].IsMatched()) { iMCPart = fTopoPerformance->ParticlesMatch()[iP].GetBestMatchWithPdg(); - matchType = 2; + matchType = 1; } - - p->SetMatch(iMCPart); - p->SetMatchType(matchType); } - - for (unsigned int iP = 0; iP < fTopoPerformance->MCParticles().size(); iP++) { - new ((*fMCParticles)[iP]) KFMCParticle(fTopoPerformance->MCParticles()[iP]); + else { + iMCPart = fTopoPerformance->ParticlesMatch()[iP].GetBestMatchWithPdg(); + matchType = 2; } + + p->SetMatch(iMCPart); + p->SetMatchType(matchType); + } + + for (unsigned int iP = 0; iP < fTopoPerformance->MCParticles().size(); iP++) { + new ((*fMCParticles)[iP]) KFMCParticle(fTopoPerformance->MCParticles()[iP]); } } } diff --git a/reco/KF/CbmKFParticleFinderQa.h b/reco/KF/CbmKFParticleFinderQa.h index b66c511718ccce22e09fb71857d3c13f152a5ca6..10cd4e95023cd6aeab8e6a05be3e0685e53de108 100644 --- a/reco/KF/CbmKFParticleFinderQa.h +++ b/reco/KF/CbmKFParticleFinderQa.h @@ -30,8 +30,6 @@ public: ~CbmKFParticleFinderQa(); void SetEffFileName(const TString& name) { fEfffileName = name; } - void SetMCTrackBranchName(const TString& name) { fMCTracksBranchName = name; } - void SetTrackMatchBranchName(const TString& name) { fTrackMatchBranchName = name; } virtual InitStatus Init(); virtual void Exec(Option_t* opt); @@ -59,45 +57,43 @@ private: void FitDecayQAHistograms(float sigma[14], const bool saveReferenceResults = false) const; void CheckDecayQA(); - //names of input branches - TString fMCTracksBranchName; //! Name of the input TCA with MC tracks - TString fTrackMatchBranchName; //! Name of the input TCA with track match - //input branches - CbmMCDataArray* fMCTrackArray; //mc tracks - TClonesArray* fMCTrackArrayEvent; - CbmMCEventList* fEventList; //mc event list in timeslice - TClonesArray* fTrackMatchArray; //track match + TClonesArray* fRecoEvents {nullptr}; //! Array of CbmEvent objects + CbmMCDataArray* fMCTrackArray {nullptr}; //mc tracks + TClonesArray* fMCTrackArrayEvent {nullptr}; + CbmMCEventList* fMcEventList {nullptr}; //mc event list in timeslice + TClonesArray* fTrackMatchArray {nullptr}; //track match // output arrays of particles - TClonesArray* fRecParticles; // output array of KF Particles - TClonesArray* fMCParticles; // output array of MC Particles - TClonesArray* fMatchParticles; // output array of match objects + TClonesArray* fRecParticles {nullptr}; // output array of KF Particles + TClonesArray* fMCParticles {nullptr}; // output array of MC Particles + TClonesArray* fMatchParticles {nullptr}; // output array of match objects - Bool_t fSaveParticles; - Bool_t fSaveMCParticles; + Bool_t fSaveParticles {false}; + Bool_t fSaveMCParticles {false}; - bool fTimeSliceMode; + bool fLegacyEventMode {false}; // event-by-event mode where data is stored in an old-fasion way //output file with histograms - TString fOutFileName; - TFile* fOutFile; - TString fEfffileName; + TString fOutFileName {"CbmKFParticleFinderQa.root"}; + TFile* fOutFile {nullptr}; + TString fEfffileName {"Efficiency.txt"}; + //KF Particle QA - KFTopoPerformance* fTopoPerformance; + KFTopoPerformance* fTopoPerformance {nullptr}; - Int_t fPrintFrequency; - Int_t fNEvents; + Int_t fPrintFrequency {100}; + Int_t fNEvents {0}; Double_t fTime[5]; //for super event analysis - bool fSuperEventAnalysis; + bool fSuperEventAnalysis {false}; //for tests - TString fReferenceResults; - int fDecayToAnalyse; - bool fCheckDecayQA; - bool fTestOk; + TString fReferenceResults {"./"}; + int fDecayToAnalyse {-1}; + bool fCheckDecayQA {false}; + bool fTestOk {false}; ClassDef(CbmKFParticleFinderQa, 1); };