/** * @file * @author Christian Simon * @since 2017-08-08 */ #include "CbmTofMatchReco.h" #include "CbmMatch.h" #include "CbmTrackMatchNew.h" #include "CbmMCEventList.h" #include "CbmMCDataManager.h" #include "CbmMCDataArray.h" #include "CbmMCTrack.h" #include "CbmTofPoint.h" #include "CbmTofDigiExp.h" #include "CbmTofHit.h" #include "CbmTofTracklet.h" #include "CbmTofDef.h" #include "FairLogger.h" #include "FairRootManager.h" #include "TClonesArray.h" // --------------------------------------------------------------------------- CbmTofMatchReco::CbmTofMatchReco() : FairTask("TofMatchReco"), fTofPointsTB(NULL), fTofDigis(NULL), fTofHits(NULL), fTofHitDigiMatches(NULL), fTofHitPointMatches(NULL), fTofHitTrackMatches(NULL), fTofPointMatches(NULL), fTofTracks(NULL), fTofTrackMatches(NULL), fTofPoints(NULL), fMCTracks(NULL), fAccMCTracks(NULL), fTofAccMCTrackMatches(NULL), fTofAccMCTrackPointMatches(NULL), fMCEventList(NULL), fbPointsInTS(kFALSE), fbAlternativeBranchNames(kFALSE), fiFileIndex(0), fbMatchMCTracks(kFALSE), fbCreateHitMatches(kTRUE) { } // --------------------------------------------------------------------------- // --------------------------------------------------------------------------- CbmTofMatchReco::~CbmTofMatchReco() { } // --------------------------------------------------------------------------- // --------------------------------------------------------------------------- void CbmTofMatchReco::Exec(Option_t*) { if(fbCreateHitMatches) { fTofHitPointMatches->Delete(); if(fbPointsInTS) { fTofPointMatches->Delete(); } for (Int_t iHit = 0; iHit < fTofHits->GetEntriesFast(); iHit++) { CbmTofHit* tHit = dynamic_cast(fTofHits->At(iHit)); CbmMatch* tHitDigiMatch = dynamic_cast(fTofHitDigiMatches->At(iHit)); CbmMatch* tHitPointMatch = new((*fTofHitPointMatches)[iHit]) CbmMatch(); Int_t iNDigiLinks = tHitDigiMatch->GetNofLinks(); for(Int_t iDigiLink = 0; iDigiLink < iNDigiLinks; iDigiLink++) { const CbmLink& tDigiLink = tHitDigiMatch->GetLink(iDigiLink); CbmTofDigiExp* tDigi = dynamic_cast(fTofDigis->At(tDigiLink.GetIndex())); CbmMatch* tDigiMatch = tDigi->GetMatch(); Int_t iNMCLinks = tDigiMatch->GetNofLinks(); for(Int_t iMCLink = 0; iMCLink < iNMCLinks; iMCLink++) { CbmLink& tMCLink = tDigiMatch->GetLink(iMCLink); Int_t iOriginalFile = tMCLink.GetFile(); Int_t iOriginalEntry = tMCLink.GetEntry(); Int_t iOriginalIndex = tMCLink.GetIndex(); Float_t dOriginalWeight = tMCLink.GetWeight(); Int_t iOriginalID = tMCLink.GetUniqueID(); tMCLink.SetWeight(1.); // MC point link duplicates are avoided when using 'CbmMatch::AddLink' // which compares the link object to be inserted with all links already // attached to the 'CbmMatch' object w.r.t. fFile, fEntry and fIndex. // If a matching link is found its weight is incremented by the weight // of the link to be inserted. // The weight then (approximately) indicates how many digis of the cluster // the MC point has contributed to. Here, we do not distinguish between // primary MC signals, afterpulses or reflections. // Collect all contributions from detector dark rate which is not stored // in a TS due to 'fbPointsInTS == kFALSE'. The condition 'fFile == -1' // is only fulfilled for dark rate points which were not meant to be // stored in a TS. if(-1 == tMCLink.GetFile()) { // Make sure all transient dark rate contributions are combined in a // single 'CbmLink' object. tMCLink.SetFile(-tof::link_TransientDarkPoint); tMCLink.SetEntry(-tof::link_TransientDarkPoint); tMCLink.SetIndex(-tof::link_TransientDarkPoint); } // MC reference point is not available in the current TS. This applies // also to dark rate points which are linked to by 'fIndex == -1' if and // only if they could not be stored in a TS. In case they were not meant // to be stored in a TS in the first place, all corresponding link // references to them would show 'fIndex == 0'. // Add a default 'CbmLink' object with fUniqueID == tof::link_PointInDifferentTS // to the match object and count all such cases in fWeight. if(-1 == tMCLink.GetIndex()) { // Make sure all broken MC point references are combined in a // single 'CbmLink' object. tMCLink.SetFile(-tof::link_PointInDifferentTS); tMCLink.SetEntry(-tof::link_PointInDifferentTS); tMCLink.SetIndex(-tof::link_PointInDifferentTS); tMCLink.SetUniqueID(tof::link_PointInDifferentTS); } if(fbPointsInTS) { // While looping over the 'CbmTofHit' branch record the point to hit // forward references for each primary MC signal found in the entire // set of backlinks. Whenever such MC point information is found add // the hit that it is attached to to its dedicated 'CbmMatch' (same // index as the point in its branch). This way, some MC points in the // input branch are not looped over as they did not create any hit in // the detector. The 'TClonesArray' would return a NULL pointer if // one tried to access data at this index. The 'CbmLink::fWeight' // for each hit descending from the point hints at the cluster size. // // For QA purposes, also internal cluster points marked as // 'tof::point_Internal' in 'CbmTofClusterMC' are provided with // point-to-hit forward references here. if(tof::kiPrimarySignalRemainder == tMCLink.GetUniqueID()%tof::kiPrimarySignalDivisor && -1 < tMCLink.GetIndex()) { CbmMatch* tPointMatch(NULL); if(tMCLink.GetIndex() >= fTofPointMatches->GetEntriesFast()) { tPointMatch = new((*fTofPointMatches)[tMCLink.GetIndex()]) CbmMatch(); } else { tPointMatch = dynamic_cast(fTofPointMatches->At(tMCLink.GetIndex())); } if(!tPointMatch) { tPointMatch = new((*fTofPointMatches)[tMCLink.GetIndex()]) CbmMatch(); } tPointMatch->AddLink(1., iHit, FairRootManager::Instance()->GetEntryNr(), tMCLink.GetFile()); } } tHitPointMatch->AddLink(tMCLink); tMCLink.SetFile(iOriginalFile); tMCLink.SetEntry(iOriginalEntry); tMCLink.SetIndex(iOriginalIndex); tMCLink.SetWeight(dOriginalWeight); tMCLink.SetUniqueID(iOriginalID); } } } if(fbPointsInTS) { // Apparently, 'TClonesArray' buffers with "empty" slots cannot be filled // into the output tree. For this reason, all MC points that no hit match // object has been created for due to counter inefficiencies need to receive // one without any 'CbmLink' inserted. for(Int_t iPoint = 0; iPoint < fTofPointsTB->GetEntriesFast(); iPoint++) { if(!fTofPointMatches->At(iPoint)) { new((*fTofPointMatches)[iPoint]) CbmMatch(); } } } } if(fbMatchMCTracks) { fTofHitTrackMatches->Delete(); fAccMCTracks->Clear(); fTofAccMCTrackPointMatches->Delete(); std::map, Int_t> AccMCTrackMap; for(Int_t iHit = 0; iHit < fTofHits->GetEntriesFast(); iHit++) { CbmTofHit* tHit = dynamic_cast(fTofHits->At(iHit)); CbmMatch* tHitPointMatch = dynamic_cast(fTofHitPointMatches->At(iHit)); CbmMatch* tHitTrackMatch = new((*fTofHitTrackMatches)[iHit]) CbmMatch(); Int_t iNPointLinks = tHitPointMatch->GetNofLinks(); for(Int_t iPointLink = 0; iPointLink < iNPointLinks; iPointLink++) { const CbmLink& tPointLink = tHitPointMatch->GetLink(iPointLink); Int_t iPointLinkType = tPointLink.GetUniqueID(); Int_t iPointFileIndex = tPointLink.GetFile(); Int_t iPointEventIndex = tPointLink.GetEntry(); Int_t iPointArrayIndex = tPointLink.GetIndex(); // The MC point referred to by the hit-to-point backlink can be retrieved. if(-1 < iPointFileIndex && -1 < iPointEventIndex && -1 < iPointArrayIndex) { CbmTofPoint* tPoint(NULL); CbmLink tMCTrackLink = CbmLink(); Int_t iMCTrackLinkType(iPointLinkType); Int_t iMCTrackFileIndex(-1); Int_t iMCTrackEventIndex(-1); Int_t iMCTrackArrayIndex(-1); Float_t dMCTrackLinkWeight(1.); if(fbPointsInTS) { tPoint = dynamic_cast(fTofPointsTB->At(iPointArrayIndex)); iMCTrackFileIndex = tPoint->GetUniqueID(); iMCTrackEventIndex = tPoint->GetEventID(); iMCTrackArrayIndex = tPoint->GetTrackID(); } else { tPoint = dynamic_cast(fTofPoints->Get(iPointFileIndex, iPointEventIndex, iPointArrayIndex)); iMCTrackFileIndex = iPointFileIndex; iMCTrackEventIndex = iPointEventIndex; iMCTrackArrayIndex = tPoint->GetTrackID(); } // beam point or dark point (MC track not existing) if(-1 == iMCTrackArrayIndex) { // dark point (not assigned to any event) if(-1 == iMCTrackEventIndex) { // MC track link type is taken from the MC point link. iMCTrackFileIndex = -tof::link_TransientDarkPoint; iMCTrackEventIndex = -tof::link_TransientDarkPoint; iMCTrackArrayIndex = -tof::link_TransientDarkPoint; } // beam point else { iMCTrackLinkType = tof::link_BeamParticleTrack; iMCTrackFileIndex = -tof::link_BeamParticleTrack; iMCTrackEventIndex = -tof::link_BeamParticleTrack; iMCTrackArrayIndex = -tof::link_BeamParticleTrack; } tMCTrackLink.SetFile(iMCTrackFileIndex); tMCTrackLink.SetEntry(iMCTrackEventIndex); tMCTrackLink.SetIndex(iMCTrackArrayIndex); tMCTrackLink.SetWeight(dMCTrackLinkWeight); tMCTrackLink.SetUniqueID(iMCTrackLinkType); } // MC track can be retrieved. else { CbmMCTrack* tMCTrack = dynamic_cast(fMCTracks->Get(iMCTrackFileIndex, iMCTrackEventIndex, iMCTrackArrayIndex)); auto AccMCTrackID = std::make_tuple(iMCTrackFileIndex, iMCTrackEventIndex, iMCTrackArrayIndex); Int_t iAccMCTrackIndex(-1); // MC track has not been referred to before. if(AccMCTrackMap.find(AccMCTrackID) == AccMCTrackMap.end()) { iAccMCTrackIndex = fAccMCTracks->GetEntriesFast(); CbmMCTrack* tAccMCTrack = new((*fAccMCTracks)[iAccMCTrackIndex]) CbmMCTrack(*tMCTrack); AccMCTrackMap.emplace(AccMCTrackID, iAccMCTrackIndex); CbmMatch* tAccMCTrackPointMatch = new((*fTofAccMCTrackPointMatches)[iAccMCTrackIndex]) CbmMatch(); tAccMCTrackPointMatch->AddLink(1., iPointArrayIndex, iPointEventIndex, iPointFileIndex); // Convert the MC track start time relative to the time the beam // entered the target in the original event to a point in time // relative to the beginning of the entire simulation run. Due to // 'CbmMCTrack::fStartT' being a 'Double32_t' type this does // precision-wise not survive an I/O cycle! if(fMCEventList) { tAccMCTrack->SetStartT(tAccMCTrack->GetStartT() + fMCEventList->GetEventTime(iMCTrackEventIndex, iMCTrackFileIndex)); } } else { iAccMCTrackIndex = AccMCTrackMap.at(AccMCTrackID); // A single MC track can create MC points in several counters which // in turn refer to the same MC track. For this reason, a link is // added to the track-to-point forward match object even if the // MC track is already known, i.e. has been found in a hit-to-point // backlink. Owing to the clustering procedure, it is possible // that several hits hold references to the same MC point which // would then be looped over more than once here and appear in the // track-to-point forward match several times. This increases the // weight of the respective links but should not bother us. CbmMatch* tAccMCTrackPointMatch = dynamic_cast(fTofAccMCTrackPointMatches->At(iAccMCTrackIndex)); tAccMCTrackPointMatch->AddLink(1., iPointArrayIndex, iPointEventIndex, iPointFileIndex); } tMCTrackLink.SetFile(fiFileIndex); tMCTrackLink.SetEntry(FairRootManager::Instance()->GetEntryNr()); tMCTrackLink.SetIndex(iAccMCTrackIndex); tMCTrackLink.SetWeight(dMCTrackLinkWeight); tMCTrackLink.SetUniqueID(iMCTrackLinkType); } tHitTrackMatch->AddLink(tMCTrackLink); } else { tHitTrackMatch->AddLink(tPointLink); } } } if(fTofTracks) { fTofTrackMatches->Delete(); fTofAccMCTrackMatches->Delete(); for(Int_t iTracklet = 0; iTracklet < fTofTracks->GetEntriesFast(); iTracklet++) { CbmTofTracklet* tTracklet = dynamic_cast(fTofTracks->At(iTracklet)); CbmTrackMatchNew* tTrackMatch = new((*fTofTrackMatches)[iTracklet]) CbmTrackMatchNew(); for(Int_t iHit = 0; iHit < tTracklet->GetNofHits(); iHit++) { CbmMatch* tHitTrackMatch = dynamic_cast(fTofHitTrackMatches->At(tTracklet->GetHitIndex(iHit))); Int_t iNHitTrackLinks = tHitTrackMatch->GetNofLinks(); for(Int_t iHitTrackLink = 0; iHitTrackLink < iNHitTrackLinks; iHitTrackLink++) { const CbmLink& tHitTrackLink = tHitTrackMatch->GetLink(iHitTrackLink); tTrackMatch->AddLink(tHitTrackLink); Int_t iAccMCTrackIndex = tHitTrackLink.GetIndex(); // The link refers to a valid accepted MC track object. Establish a new // forward reference from that object to the current ToF tracklet. if(-1 < iAccMCTrackIndex) { CbmMatch* tAccMCTrackMatch(NULL); if(iAccMCTrackIndex >= fTofAccMCTrackMatches->GetEntriesFast()) { tAccMCTrackMatch = new((*fTofAccMCTrackMatches)[iAccMCTrackIndex]) CbmMatch(); } else { tAccMCTrackMatch = dynamic_cast(fTofAccMCTrackMatches->At(iAccMCTrackIndex)); } if(!tAccMCTrackMatch) { tAccMCTrackMatch = new((*fTofAccMCTrackMatches)[iAccMCTrackIndex]) CbmMatch(); } tAccMCTrackMatch->AddLink(1., iTracklet, FairRootManager::Instance()->GetEntryNr(), fiFileIndex); } } } const CbmLink& tTrackletRefLink = tTrackMatch->GetMatchedLink(); // The weight of a tracklet-to-track link indicates to how many hits // which the tracklet is composed of the referenced MC track contributed // to. Thus, the weight of the matched link (i.e. [one of] the link[s] with // the largest overall weight in the match object) gives the maximum number // of hits a single MC track referred to by the tracklet helped producing. // If this maximum weight equals the number of hits attached to the tracklet, // the tracklet is considered a "pure" reconstruction of an original MC track. Int_t iNTrueTrackletHits(static_cast(tTrackletRefLink.GetWeight())); // Make sure that a hit from the beam reference counter attached to // the tracklet is considered a "true" hit even though it does not // refer to the actual track particle. for(Int_t iHit = 0; iHit < tTracklet->GetNofHits(); iHit++) { CbmMatch* tHitTrackMatch = dynamic_cast(fTofHitTrackMatches->At(tTracklet->GetHitIndex(iHit))); Int_t iNHitTrackLinks = tHitTrackMatch->GetNofLinks(); for(Int_t iHitTrackLink = 0; iHitTrackLink < iNHitTrackLinks; iHitTrackLink++) { const CbmLink& tHitTrackLink = tHitTrackMatch->GetLink(iHitTrackLink); if(tof::link_BeamParticleTrack == tHitTrackLink.GetUniqueID()) { iNTrueTrackletHits++; break; } } } tTrackMatch->SetNofTrueHits(iNTrueTrackletHits); tTrackMatch->SetNofWrongHits(tTracklet->GetNofHits() - iNTrueTrackletHits); } } // So far, only MC tracks contributing to reconstructed ToF hits/tracklets have // been identified. To extract all accepted MC tracks within the reconstructed // event, however, it is necessary to scan the MC track affiliations of all // MC points available. If MC points are directly provided in an input // data array, one can find also those MC tracks which left no measurable // trace in the entire setup, i.e. for which no hit was created. They provide // the normalization of efficiency calculations for all geometrically accepted // tracks in the setup. If MC points are only accessible through hit-to-point // backlinks, determining the MC track reconstruction efficiency is thus not // possible. if(fbPointsInTS) { for(Int_t iPoint = 0; iPoint < fTofPointsTB->GetEntriesFast(); iPoint++) { CbmTofPoint* tPoint = dynamic_cast(fTofPointsTB->At(iPoint)); Int_t iMCTrackFileIndex = tPoint->GetUniqueID(); Int_t iMCTrackEventIndex = tPoint->GetEventID(); Int_t iMCTrackArrayIndex = tPoint->GetTrackID(); if(-1 < iMCTrackArrayIndex) { CbmMCTrack* tMCTrack = dynamic_cast(fMCTracks->Get(iMCTrackFileIndex, iMCTrackEventIndex, iMCTrackArrayIndex)); auto AccMCTrackID = std::make_tuple(iMCTrackFileIndex, iMCTrackEventIndex, iMCTrackArrayIndex); Int_t iAccMCTrackIndex(-1); // MC track has not been referred to before. if(AccMCTrackMap.find(AccMCTrackID) == AccMCTrackMap.end()) { iAccMCTrackIndex = fAccMCTracks->GetEntriesFast(); CbmMCTrack* tAccMCTrack = new((*fAccMCTracks)[iAccMCTrackIndex]) CbmMCTrack(*tMCTrack); AccMCTrackMap.emplace(AccMCTrackID, iAccMCTrackIndex); CbmMatch* tAccMCTrackPointMatch = new((*fTofAccMCTrackPointMatches)[iAccMCTrackIndex]) CbmMatch(); tAccMCTrackPointMatch->AddLink(1., iPoint, FairRootManager::Instance()->GetEntryNr(), fiFileIndex); // Convert the MC track start time relative to the time the beam // entered the target in the original event to a point in time // relative to the beginning of the entire simulation run. Due to // 'CbmMCTrack::fStartT' being a 'Double32_t' type this does // precision-wise not survive an I/O cycle! if(fMCEventList) { tAccMCTrack->SetStartT(tAccMCTrack->GetStartT() + fMCEventList->GetEventTime(iMCTrackEventIndex, iMCTrackFileIndex)); } } else { iAccMCTrackIndex = AccMCTrackMap.at(AccMCTrackID); CbmMatch* tAccMCTrackPointMatch = dynamic_cast(fTofAccMCTrackPointMatches->At(iAccMCTrackIndex)); tAccMCTrackPointMatch->AddLink(1., iPoint, FairRootManager::Instance()->GetEntryNr(), fiFileIndex); } } } // Create (empty) track-to-tracklet match objects also for accepted tracks // that could not be reconstructed. if(fTofTracks) { for(Int_t iTrack = 0; iTrack < fAccMCTracks->GetEntriesFast(); iTrack++) { if(!fTofAccMCTrackMatches->At(iTrack)) { new((*fTofAccMCTrackMatches)[iTrack]) CbmMatch(); } } } } } } // --------------------------------------------------------------------------- // --------------------------------------------------------------------------- InitStatus CbmTofMatchReco::Init() { // --- Get FairRootManager instance FairRootManager* tRootManager = FairRootManager::Instance(); if(!tRootManager) { LOG(ERROR)<<"FairRootManager not found."<(tRootManager->GetObject(tPointBranchName)); if(fTofPointsTB) { fbPointsInTS = kTRUE; } // --- Get input array (CbmTofDigiExp) // Ensure compatibility with realistically calibrated digis created by 'CbmTofTestBeamClusterizer'. fTofDigis = dynamic_cast(tRootManager->GetObject("TofCalDigi")); if(!fTofDigis) { fTofDigis = dynamic_cast(tRootManager->GetObject(tDigiBranchName)); if(!fTofDigis) { LOG(ERROR)<(tRootManager->GetObject(tHitBranchName)); if(!fTofHits) { LOG(ERROR)<(tRootManager->GetObject(tHitDigiMatchBranchName)); if(!fTofHitDigiMatches) { LOG(ERROR)<(tRootManager->GetObject(tHitPointMatchBranchName)); if(!fTofHitPointMatches) { fTofHitPointMatches = new TClonesArray("CbmMatch", 10000); tRootManager->Register(tHitPointMatchBranchName, "TOF", fTofHitPointMatches, IsOutputBranchPersistent(tHitPointMatchBranchName)); } else { fbCreateHitMatches = kFALSE; } if(fbPointsInTS) { // --- Register output array (CbmMatch) fTofPointMatches = dynamic_cast(tRootManager->GetObject(tPointHitMatchBranchName)); if(!fTofPointMatches) { if(kFALSE == fbCreateHitMatches) { LOG(ERROR)<Register(tPointHitMatchBranchName, "TOF", fTofPointMatches, IsOutputBranchPersistent(tPointHitMatchBranchName)); } } if(fbMatchMCTracks) { CbmMCDataManager* tMCManager = dynamic_cast(tRootManager->GetObject("MCDataManager")); if(!tMCManager) { LOG(ERROR)<<"Could not retrieve branch MCDataManager from FairRootManager."<InitBranch(tMCTrackBranchName.Data()); if(!fMCTracks) { LOG(ERROR)<InitBranch(tMCPointBranchName.Data()); if(!fTofPoints) { LOG(ERROR)<(tRootManager->GetObject("MCEventList.")); if(!fMCEventList) { // The input data arrays correspond to a reconstructed event. fMCEventList = dynamic_cast(tRootManager->GetObject("EventList.")); } fAccMCTracks = new TClonesArray("CbmMCTrack", 10000); tRootManager->Register(tAccMCTrackBranchName, "TOF", fAccMCTracks, IsOutputBranchPersistent(tAccMCTrackBranchName)); fTofAccMCTrackPointMatches = new TClonesArray("CbmMatch", 10000); tRootManager->Register(tAccMCTrackPointMatchBranchName, "TOF", fTofAccMCTrackPointMatches, IsOutputBranchPersistent(tAccMCTrackPointMatchBranchName)); fTofHitTrackMatches = new TClonesArray("CbmMatch", 10000); tRootManager->Register(tHitTrackMatchBranchName, "TOF", fTofHitTrackMatches, IsOutputBranchPersistent(tHitTrackMatchBranchName)); // --- Get input array (CbmTofTracklet) fTofTracks = dynamic_cast(tRootManager->GetObject(tTrackBranchName.Data())); if(fTofTracks) { fTofTrackMatches = new TClonesArray("CbmTrackMatchNew", 10000); tRootManager->Register(tTrackMatchBranchName, "TOF", fTofTrackMatches, IsOutputBranchPersistent(tTrackMatchBranchName)); fTofAccMCTrackMatches = new TClonesArray("CbmMatch", 10000); tRootManager->Register(tAccMCTrackMatchBranchName, "TOF", fTofAccMCTrackMatches, IsOutputBranchPersistent(tAccMCTrackMatchBranchName)); } } return kSUCCESS; } // --------------------------------------------------------------------------- ClassImp(CbmTofMatchReco)