24 #include "FairRuntimeDb.h"
30 map<pair<ECbmModuleId, Int_t>, pair<UInt_t, set<Int_t>>>
54 int canSkipHits = 0.3 * fMinTrackLength;
55 fMinTrackLength -= canSkipHits;
57 FairRootManager* ioman = FairRootManager::Instance();
59 if (0 == ioman) LOG(fatal) <<
"No FairRootManager";
64 if (0 == mcManager) LOG(fatal) <<
"No MC data manager";
68 if (0 == fMCTracks) LOG(fatal) <<
"No MC tracks in the input file";
70 for (
int i = 0; fMCTracks->Size(0,
i) >= 0; ++
i) {
71 gMCTracks.push_back(vector<MCTrackDesc>());
73 auto nofMcTracks = fMCTracks->Size(0,
i);
75 eventTracks.resize(nofMcTracks);
77 for (Int_t j = 0; j < nofMcTracks; ++j) {
80 static_cast<const CbmMCTrack*
>(fMCTracks->Get(0,
i, j));
86 fGlobalTracks =
static_cast<TClonesArray*
>(ioman->GetObject(
"GlobalTrack"));
88 if (0 == fGlobalTracks) LOG(fatal) <<
"No global tracks in the input file";
91 fStsTracks =
static_cast<TClonesArray*
>(ioman->GetObject(
"StsTrack"));
93 if (0 == fStsTracks) LOG(fatal) <<
"No sts tracks in the input file";
95 fStsHits =
static_cast<TClonesArray*
>(ioman->GetObject(
"StsHit"));
97 if (0 == fStsHits) LOG(fatal) <<
"No sts hits in the input file";
99 fStsClusters =
static_cast<TClonesArray*
>(ioman->GetObject(
"StsCluster"));
101 if (0 == fStsClusters) LOG(fatal) <<
"No sts clusters in the input file";
103 fStsDigis =
static_cast<TClonesArray*
>(ioman->GetObject(
"StsDigi"));
105 if (0 == fStsDigis) LOG(fatal) <<
"No sts digis in the input file";
108 static_cast<TClonesArray*
>(ioman->GetObject(
"StsDigiMatch"));
110 if (0 == fStsDigiMatches)
111 LOG(fatal) <<
"No sts digi matches in the input file";
113 fStsPoints = mcManager->
InitBranch(
"StsPoint");
115 if (0 == fStsPoints) LOG(fatal) <<
"No sts MC points in the input file";
117 for (Int_t
i = 0; fStsPoints->Size(0,
i) >= 0; ++
i) {
118 Int_t nofPoints = fStsPoints->
Size(0,
i);
121 for (Int_t j = 0; j < nofPoints; ++j) {
123 static_cast<const CbmStsPoint*
>(fStsPoints->Get(0,
i, j));
124 Int_t trackId = stsPoint->GetTrackID();
125 Int_t stationNumber =
127 auto& trackDesk =
tracks[trackId];
135 fMuchTracks =
static_cast<TClonesArray*
>(ioman->GetObject(
"MuchTrack"));
137 if (0 == fMuchTracks) LOG(fatal) <<
"No much tracks in the input file";
139 fMuchHits =
static_cast<TClonesArray*
>(ioman->GetObject(
"MuchPixelHit"));
141 if (0 == fMuchHits) LOG(fatal) <<
"No much hits in the input file";
143 fMuchClusters =
static_cast<TClonesArray*
>(ioman->GetObject(
"MuchCluster"));
145 if (0 == fMuchClusters) LOG(fatal) <<
"No much clusters in the input file";
147 fMuchDigis =
static_cast<TClonesArray*
>(ioman->GetObject(
"MuchDigi"));
149 if (0 == fMuchDigis) LOG(fatal) <<
"No much digis in the input file";
152 static_cast<TClonesArray*
>(ioman->GetObject(
"MuchDigiMatch"));
154 if (0 == fMuchDigiMatches)
155 LOG(fatal) <<
"No much digi matches in the input file";
157 fMuchPoints = mcManager->
InitBranch(
"MuchPoint");
159 if (0 == fMuchPoints) LOG(fatal) <<
"No much MC points in the input file";
161 for (Int_t
i = 0; fMuchPoints->Size(0,
i) >= 0; ++
i) {
162 Int_t nofPoints = fMuchPoints->
Size(0,
i);
165 for (Int_t j = 0; j < nofPoints; ++j) {
168 Int_t trackId = muchPoint->GetTrackID();
169 int muchStationNumber =
173 int stationNumber = muchStationNumber * 3 + layerNumber;
174 auto& trackDesk =
tracks[trackId];
176 stationNumber, set<Int_t>()};
182 fTrdTracks =
static_cast<TClonesArray*
>(ioman->GetObject(
"TrdTrack"));
184 if (0 == fTrdTracks) LOG(fatal) <<
"No trd tracks in the input file";
186 fTrdHits =
static_cast<TClonesArray*
>(ioman->GetObject(
"TrdHit"));
188 if (0 == fTrdHits) LOG(fatal) <<
"No trd hits in the input file";
190 fTrdClusters =
static_cast<TClonesArray*
>(ioman->GetObject(
"TrdCluster"));
192 if (0 == fTrdClusters) LOG(fatal) <<
"No global tracks in the input file";
195 static_cast<TClonesArray*
>(ioman->GetObject(
"TrdDigiMatch"));
197 if (0 == fTrdDigiMatches)
198 LOG(fatal) <<
"No trd hit to digi matches in the input file";
200 fTrdPoints = mcManager->
InitBranch(
"TrdPoint");
202 if (0 == fTrdPoints) LOG(fatal) <<
"No trd MC points in the input file";
204 for (Int_t
i = 0; fTrdPoints->Size(0,
i) >= 0; ++
i) {
205 Int_t nofPoints = fTrdPoints->
Size(0,
i);
208 for (Int_t j = 0; j < nofPoints; ++j) {
210 static_cast<const CbmTrdPoint*
>(fTrdPoints->Get(0,
i, j));
211 Int_t trackId = trdPoint->GetTrackID();
214 auto& trackDesk =
tracks[trackId];
222 fTofHits =
static_cast<TClonesArray*
>(ioman->GetObject(
"TofHit"));
224 if (0 == fTofHits) LOG(fatal) <<
"No tof hits in the input file";
227 static_cast<TClonesArray*
>(ioman->GetObject(
"TofDigiMatch"));
229 if (0 == fTofHitDigiMatches)
230 LOG(fatal) <<
"No tof hit to digi matches in the input file";
232 fTofDigiPointMatches =
233 static_cast<TClonesArray*
>(ioman->GetObject(
"TofDigiMatchPoints"));
235 if (0 == fTofDigiPointMatches)
236 LOG(fatal) <<
"No tof digi to point matches in the input file";
238 fTofDigis =
static_cast<TClonesArray*
>(ioman->GetObject(
"TofDigi"));
240 if (0 == fTofDigis) LOG(fatal) <<
"No tof digis in the input file";
242 fTofPoints = mcManager->
InitBranch(
"TofPoint");
244 if (0 == fTofPoints) LOG(fatal) <<
"No tof MC points in the input file";
246 for (Int_t
i = 0; fTofPoints->Size(0,
i) >= 0; ++
i) {
247 Int_t nofPoints = fTofPoints->
Size(0,
i);
250 for (Int_t j = 0; j < nofPoints; ++j) {
252 static_cast<const CbmTofPoint*
>(fTofPoints->Get(0,
i, j));
253 Int_t trackId = tofPoint->GetTrackID();
254 auto& trackDesk =
tracks[trackId];
271 Int_t nofGlobalTracks = fGlobalTracks->GetEntriesFast();
273 for (Int_t
i = 0;
i < nofGlobalTracks; ++
i) {
281 map<void*, set<pair<ECbmModuleId, Int_t>>> mcTrackPtrs;
282 Int_t nofTrackHits = 0;
286 static_cast<const CbmStsTrack*
>(fStsTracks->At(stsTrackIndex));
288 nofTrackHits += nofStsHits;
290 for (Int_t j = 0; j < nofStsHits; ++j) {
293 static_cast<const CbmStsHit*
>(fStsHits->At(stsHitInd));
298 static_cast<const CbmStsCluster*
>(fStsClusters->At(frontClusterInd));
301 for (Int_t k = 0; k < nofFrontDigis; ++k) {
302 Int_t stsDigiInd = frontCluster->
GetDigi(k);
305 static_cast<const CbmMatch*
>(fStsDigiMatches->At(stsDigiInd));
308 for (Int_t
m = 0;
m < nofLinks; ++
m) {
313 fStsPoints->Get(0, eventId, mcPointId));
314 Int_t trackId = stsPoint->GetTrackID();
315 auto& mcTrack = eventMCTracks[trackId];
319 if (recoIter != mcTrack.mcPointMap.end())
320 recoIter->second.second.insert(
i);
327 static_cast<const CbmStsCluster*
>(fStsClusters->At(backClusterInd));
330 for (Int_t k = 0; k < nofBackDigis; ++k) {
331 Int_t stsDigiInd = backCluster->
GetDigi(k);
334 static_cast<const CbmMatch*
>(fStsDigiMatches->At(stsDigiInd));
337 for (Int_t
m = 0;
m < nofLinks; ++
m) {
342 fStsPoints->Get(0, eventId, mcPointId));
343 Int_t trackId = stsPoint->GetTrackID();
344 auto& mcTrack = eventMCTracks[trackId];
348 if (recoIter != mcTrack.mcPointMap.end())
349 recoIter->second.second.insert(
i);
359 static_cast<const CbmMuchTrack*
>(fMuchTracks->At(muchTrackIndex));
361 nofTrackHits += nofMuchHits;
363 for (Int_t j = 0; j < nofMuchHits; ++j) {
370 Int_t clusterId = muchHit->
GetRefId();
375 for (Int_t k = 0; k < nofDigis; ++k) {
376 Int_t digiId = cluster->
GetDigi(k);
379 static_cast<const CbmMatch*
>(fMuchDigiMatches->At(digiId));
382 for (Int_t
m = 0;
m < nofLinks; ++
m) {
387 fMuchPoints->Get(0, eventId, mcPointId));
388 Int_t trackId = muchPoint->GetTrackID();
389 auto& mcTrack = eventMCTracks[trackId];
390 auto recoIter = mcTrack.mcPointMap.find(
393 if (recoIter != mcTrack.mcPointMap.end())
394 recoIter->second.second.insert(
i);
404 static_cast<const CbmTrdTrack*
>(fTrdTracks->At(trdTrackIndex));
406 nofTrackHits += nofTrdHits;
408 for (Int_t j = 0; j < nofTrdHits; ++j) {
411 static_cast<const CbmTrdHit*
>(fTrdHits->At(trdHitInd));
412 Int_t clusterId = trdHit->
GetRefId();
414 static_cast<const CbmTrdCluster*
>(fTrdClusters->At(clusterId));
417 for (Int_t k = 0; k < nofDigis; ++k) {
418 Int_t digiId = cluster->
GetDigi(k);
420 static_cast<const CbmMatch*
>(fTrdDigiMatches->At(digiId));
423 for (Int_t
m = 0;
m < nofLinks; ++
m) {
428 fTrdPoints->Get(0, eventId, mcPointId));
429 Int_t trackId = trdPoint->GetTrackID();
430 auto& mcTrack = eventMCTracks[trackId];
434 if (recoIter != mcTrack.mcPointMap.end())
435 recoIter->second.second.insert(
i);
446 static_cast<const CbmMatch*
>(fTofHitDigiMatches->At(tofHitIndex));
449 for (Int_t j = 0; j < nofTofDigis; ++j) {
451 Int_t digiInd = digiLink.
GetIndex();
453 static_cast<const CbmMatch*
>(fTofDigiPointMatches->At(digiInd));
456 for (Int_t k = 0; k < nofPoints; ++k) {
458 Int_t mcEventId = pointLink.
GetEntry();
459 Int_t mcPointId = pointLink.
GetIndex();
461 fTofPoints->Get(0, mcEventId, mcPointId));
462 Int_t trackId = tofPoint->GetTrackID();
463 auto& mcTrack = eventMCTracks[trackId];
467 if (recoIter != mcTrack.mcPointMap.end())
468 recoIter->second.second.insert(
i);
477 for (
const auto& j : mcTrackPtrs) {
478 if (j.second.size() >= size_t(0.7 * nofTrackHits)) {
489 int nofRefMCTRacks = 0;
490 int nofMatchedRefMCTracks = 0;
491 uint maxTrackLength = 0;
494 for (
const auto& j :
i) {
495 set<pair<ECbmModuleId, UInt_t>> stationNofs;
497 for (
const auto& k : j.mcPointMap)
498 stationNofs.insert(make_pair(k.first.first, k.second.first));
500 if (stationNofs.size() >
static_cast<size_t>(maxTrackLength))
501 maxTrackLength = stationNofs.size();
503 if (stationNofs.size() <
static_cast<size_t>(fMinTrackLength))
continue;
505 map<Int_t, set<pair<ECbmModuleId, UInt_t>>> recoToStationMatches;
507 for (
const auto& k : j.mcPointMap) {
508 for (
const auto&
m : k.second.second)
509 recoToStationMatches[
m].insert(
510 make_pair(k.first.first, k.second.first));
513 int nofMatchedReco = 0;
515 for (
const auto& k : recoToStationMatches) {
516 if (k.second.size() >=
size_t(stationNofs.size() * 0.7))
520 if (nofMatchedReco > 1)
gNofClones += nofMatchedReco - 1;
525 if (fPrimaryParticlePdg >= 0
526 && (j.parentInd < 0 ||
i[j.parentInd].pdg != fPrimaryParticlePdg))
531 if (nofMatchedReco > 0) ++nofMatchedRefMCTracks;
535 cout <<
"Maximum MC track length: " << maxTrackLength << endl;
537 double res = 100 * nofMatchedRefMCTracks;
538 res /= nofRefMCTRacks;
539 cout <<
"The reconstruction efficiency: " << res <<
" % = 100 * "
540 << nofMatchedRefMCTracks <<
" / " << nofRefMCTRacks << endl;
544 cout <<
"The share of ghosts: " << res <<
" % = 100 * (" <<
gNofRecoTracks
549 cout <<
"The share of clones: " << res <<
" % = 100 * " <<
gNofClones <<
" / "
555 FairRun::Instance()->GetRuntimeDb()->getContainer(
"CbmBinnedSettings"));