#include "CbmMuchMatchTracks.h" #include "CbmMuchCluster.h" #include "CbmMuchDigiMatch.h" #include "CbmMuchPixelHit.h" #include "CbmMuchTrack.h" #include "CbmTrackMatch.h" #include "FairMCPoint.h" #include "FairRootManager.h" #include "TClonesArray.h" #include <iomanip> #include <iostream> #include <map> CbmMuchMatchTracks::CbmMuchMatchTracks() : FairTask("CbmMuchMatchTracks") , fTracks(NULL) , fPoints(NULL) , fPixelHits(NULL) , fMatches(NULL) , fPixelDigiMatches(NULL) , fClusters(NULL) , fNofHits(0) , fNofTrueHits(0) , fNofWrongHits(0) , fNofFakeHits(0) , fNEvents(0) {} CbmMuchMatchTracks::~CbmMuchMatchTracks() {} InitStatus CbmMuchMatchTracks::Init() { FairRootManager* ioman = FairRootManager::Instance(); if (ioman == NULL) Fatal("CbmMuchMatchTracks::Init", "RootManager not instantised!"); fPixelHits = (TClonesArray*) ioman->GetObject("MuchPixelHit"); if (fPixelHits == NULL) Fatal("CbmMuchMatchTracks::Init", "No MuchPixelHit array!"); fTracks = (TClonesArray*) ioman->GetObject("MuchTrack"); if (fTracks == NULL) Fatal("CbmMuchMatchTracks::Init", "No MuchTrack array!"); fPoints = (TClonesArray*) ioman->GetObject("MuchPoint"); if (fPoints == NULL) Fatal("CbmMuchMatchTracks::Init", "No MuchPoint array!"); fPixelDigiMatches = (TClonesArray*) ioman->GetObject("MuchDigiMatch"); if (fPixelDigiMatches == NULL) Fatal("CbmMuchMatchTracks::Init", "No MuchDigiMatch array!"); fClusters = (TClonesArray*) ioman->GetObject("MuchCluster"); if (fClusters == NULL) Info("CbmMuchMatchTracks::Init", "No cluster array -- simple hit to digi matching will be used"); fMatches = new TClonesArray("CbmTrackMatch", 100); ioman->Register("MuchTrackMatch", "MUCH", fMatches, IsOutputBranchPersistent("MuchTrackMatch")); return kSUCCESS; } void CbmMuchMatchTracks::Exec(Option_t*) { fMatches->Clear(); Int_t nofTracks = fTracks->GetEntriesFast(); for (Int_t iTrack = 0; iTrack < nofTracks; iTrack++) { // Loop over tracks // std::map stores MC track id to number of contributions of this MC track std::map<Int_t, Int_t> matchMap; CbmMuchTrack* pTrack = static_cast<CbmMuchTrack*>(fTracks->At(iTrack)); if (pTrack == NULL) continue; Int_t nofHits = pTrack->GetNofHits(); for (Int_t iHit = 0; iHit < nofHits; iHit++) { // Loop over hits HitType hitType = pTrack->GetHitType(iHit); if (hitType == kMUCHPIXELHIT) { ExecPixel(matchMap, pTrack->GetHitIndex(iHit)); } else { TObject::Fatal("CbmMuchMatchTracks::Exec", "Hit type not supported!"); } } // Loop over hits Int_t nofTrue = 0; Int_t bestMcTrackId = -1; Int_t nPoints = 0; for (std::map<Int_t, Int_t>::iterator it = matchMap.begin(); it != matchMap.end(); it++) { if (it->first != -1 && it->second >= nofTrue) { bestMcTrackId = it->first; nofTrue = it->second; } nPoints += it->second; } Int_t nofFake = 0; Int_t nofWrong = nofHits - nofTrue - nofFake; Int_t nofMcTracks = matchMap.size() - 1; new ((*fMatches)[iTrack]) CbmTrackMatch(bestMcTrackId, nofTrue, nofWrong, nofFake, nofMcTracks); fNofHits += nofHits; fNofTrueHits += nofTrue; fNofWrongHits += nofWrong; fNofFakeHits += nofFake; if (fVerbose > 1) std::cout << "iTrack=" << iTrack << " mcTrack=" << bestMcTrackId << " nPoints=" << nPoints << " nofTrue=" << nofTrue << " nofWrong=" << nofWrong << " nofFake=" << nofFake << " nofMcTracks=" << nofMcTracks << std::endl; } // Loop over tracks fNEvents++; } void CbmMuchMatchTracks::Finish() { Double_t trueHits = 100. * Double_t(fNofTrueHits) / Double_t(fNofHits); Double_t wrongHits = 100. * Double_t(fNofWrongHits) / Double_t(fNofHits); Double_t fakeHits = 100. * Double_t(fNofFakeHits) / Double_t(fNofHits); std::cout << "=================================================" << std::endl; std::cout << "===== " << GetName() << ": Run summary " << std::endl; std::cout << "True hits: " << trueHits << "%" << std::endl; std::cout << "Wrong hits: " << wrongHits << "%" << std::endl; std::cout << "Fake hits: " << fakeHits << "%" << std::endl; std::cout << "=================================================" << std::endl; } void CbmMuchMatchTracks::ExecPixel(std::map<Int_t, Int_t>& matchMap, Int_t index) { // std::set stores MC track indices contributed to a certain hit std::set<Int_t> mcIdHit; CbmMuchPixelHit* hit = static_cast<CbmMuchPixelHit*>(fPixelHits->At(index)); if (hit == NULL) return; Int_t clusterId = hit->GetRefId(); CbmMuchCluster* cluster = static_cast<CbmMuchCluster*>(fClusters->At(clusterId)); if (cluster == NULL) return; for (Int_t iDigi = 0; iDigi < cluster->GetNofDigis(); iDigi++) { Int_t digiId = cluster->GetDigi(iDigi); CbmMuchDigiMatch* digiMatch = static_cast<CbmMuchDigiMatch*>(fPixelDigiMatches->At(digiId)); if (digiMatch == NULL) continue; for (Int_t iPoint = 0; iPoint < digiMatch->GetNofLinks(); iPoint++) { Int_t pointIndex = digiMatch->GetLink(iPoint).GetIndex(); if (pointIndex < 0) { // Fake or background hit mcIdHit.insert(-1); continue; } FairMCPoint* point = static_cast<FairMCPoint*>(fPoints->At(pointIndex)); if (point == NULL) continue; mcIdHit.insert(point->GetTrackID()); } } // loop over digis for (std::set<Int_t>::iterator it = mcIdHit.begin(); it != mcIdHit.end(); it++) { matchMap[*it]++; } } ClassImp(CbmMuchMatchTracks);