13#include <TStopwatch.h>
55#include "GPUParam.inc"
59#include <unordered_map>
77using VTIndexV = std::pair<int, o2::dataformats::VtxTrackIndex>;
87 : mDataRequest(dr), mGGCCDBRequest(gr), mTracksSrc(
src), mCheckSV(checkSV)
98 void process(const
o2::globaltracking::RecoContainer& recoData);
101 void processTPCTrackRefs();
102 void loadTPCOccMap(const
o2::globaltracking::RecoContainer& recoData);
103 void fillMCClusterInfo(const
o2::globaltracking::RecoContainer& recoData);
104 void prepareITSData(const
o2::globaltracking::RecoContainer& recoData);
105 bool processMCParticle(
int src,
int ev,
int trid);
106 bool addMCParticle(const
MCTrack& mctr, const
o2::
MCCompLabel& lb, TParticlePDG* pPDG =
nullptr);
109 bool refitV0(
int i,
o2::dataformats::
V0&
v0, const
o2::globaltracking::RecoContainer& recoData);
111 float getDCAYCut(
float pt) const;
113 const
std::vector<
o2::
MCTrack>* mCurrMCTracks =
nullptr;
114 TVector3 mCurrMCVertex;
115 o2::tpc::VDriftHelper mTPCVDriftHelper{};
117 std::shared_ptr<DataRequest> mDataRequest;
118 std::shared_ptr<o2::base::GRPGeomRequest> mGGCCDBRequest;
119 std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOut;
120 std::vector<float> mTBinClOcc;
121 std::vector<float> mTBinClOccHist;
122 std::vector<long> mIntBC;
123 std::vector<float> mTPCOcc;
124 std::vector<int> mITSOcc;
125 bool mCheckSV =
false;
126 bool mRecProcStage =
false;
127 int mNTPCOccBinLength = 0;
128 float mNTPCOccBinLengthInv = -1.f;
130 float mITSTimeBiasMUS = 0.f;
131 float mITSROFrameLengthMUS = 0.f;
132 float mTPCTBinMUS = 0.;
134 int mNCheckDecays = 0;
138 std::vector<int> mITSROF;
139 std::vector<TBracket> mITSROFBracket;
140 std::vector<o2::MCCompLabel> mDecProdLblPool;
141 std::vector<MCVertex> mMCVtVec{};
147 int daughterFirst = -1;
148 int daughterLast = -1;
151 std::vector<std::vector<DecayRef>> mDecaysMaps;
152 std::unordered_map<o2::MCCompLabel, TrackFamily> mSelMCTracks;
153 std::unordered_map<o2::MCCompLabel, std::pair<int, int>> mSelTRefIdx;
154 std::vector<o2::track::TrackPar> mSelTRefs;
156 static constexpr float MaxSnp = 0.9;
164 mDBGOut = std::make_unique<o2::utils::TreeStreamRedirector>(
"trackMCStudy.root",
"recreate");
165 mVerbose = ic.
options().
get<
int>(
"device-verbosity");
168 for (
int id = 0;
id <
sizeof(
params.decayPDG) /
sizeof(
int);
id++) {
169 if (
params.decayPDG[
id] < 0) {
174 mDecaysMaps.resize(mNCheckDecays);
175 mTPCCorrMapsLoader.
init(ic);
181 for (
int i = 0;
i < mNCheckDecays;
i++) {
182 mDecaysMaps[
i].clear();
184 mDecProdLblPool.clear();
186 mCurrMCTracks =
nullptr;
189 updateTimeDependentParams(pc);
190 mRecProcStage =
false;
199 bool updateMaps =
false;
205 LOGP(info,
"Updating TPC fast transform map with new VDrift factor of {} wrt reference {} and DriftTimeOffset correction {} wrt {} from source {}",
215 static bool initOnceDone =
false;
220 LOGP(info,
"VertexTrackMatcher ITSROFrameLengthMUS:{}", mITSROFrameLengthMUS);
223 mTPCTBinMUS = elParam.ZbinWidth;
228 mFitterV0.setUseAbsDCA(svparam.useAbsDCA);
229 mFitterV0.setPropagateToPCA(
false);
230 mFitterV0.setMaxR(svparam.maxRIni);
231 mFitterV0.setMinParamChange(svparam.minParamChange);
232 mFitterV0.setMinRelChi2Change(svparam.minRelChi2Change);
233 mFitterV0.setMaxDZIni(svparam.maxDZIni);
234 mFitterV0.setMaxDXYIni(svparam.maxDXYIni);
235 mFitterV0.setMaxChi2(svparam.maxChi2);
237 mFitterV0.setUsePropagator(svparam.usePropagator);
238 mFitterV0.setRefitWithMatCorr(svparam.refitWithMatCorr);
239 mFitterV0.setMaxStep(svparam.maxStep);
240 mFitterV0.setMaxSnp(svparam.maxSnp);
241 mFitterV0.setMinXSeed(svparam.minXSeed);
248 constexpr float SQRT12Inv = 0.288675f;
255 int nv = vtxRefs.size();
259 prepareITSData(recoData);
260 loadTPCOccMap(recoData);
261 auto getITSPatt = [&](
GTrackID gid, uint8_t& ncl) {
265 ncl = itsTrf.getNClusters();
266 for (
int il = 0; il < 7; il++) {
267 if (itsTrf.hasHitOnLayer(il)) {
274 for (
int il = 0; il < 7; il++) {
275 if (itsTr.hasHitOnLayer(il)) {
286 uint8_t clSect = 0, clRow = 0, lowestR = -1;
291 for (
int ic = 0; ic < trc.getNClusterReferences(); ic++) {
292 trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
293 if (clRow < lowestR) {
297 unsigned int absoluteIndex = tpcClusAcc.clusterOffset[clSect][clRow] + clIdx;
302 tref.lowestPadRow = lowestR;
303 const auto& clus = tpcClusAcc.clusters[clSect][clRow][clIdx];
304 int padFromEdge =
int(clus.getPad()), npads = o2::gpu::GPUTPCGeometry::NPads(clRow);
305 if (padFromEdge > npads / 2) {
306 padFromEdge = npads - 1 - padFromEdge;
308 tref.padFromEdge = uint8_t(padFromEdge);
309 trc.getClusterReference(clRefs, 0, clSect, clRow, clIdx);
310 tref.rowMaxTPC = clRow;
317 const auto* TPCClMClab = recoData.
inputsTPCclusters->clusterIndex.clustersMCTruth;
319 for (
int ic = 0; ic < trc.getNClusterReferences(); ic++) {
320 uint8_t clSect = 0, clRow = 0;
322 trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
323 auto labels = TPCClMClab->getLabels(clIdx + TPCClusterIdxStruct.clusterOffset[clSect][clRow]);
324 for (
auto& lbl : labels) {
339 unsigned int rofCount = 0;
341 for (
const auto& mcIR : mcEvRecords) {
342 long tbc = mcIR.differenceInBC(recoData.
startIR);
343 auto& mcVtx = mMCVtVec.emplace_back();
345 mcVtx.ID = mIntBC.size();
346 mIntBC.push_back(tbc);
347 int occBin = tbc / 8 * mNTPCOccBinLengthInv;
348 mTPCOcc.push_back(occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]));
350 long gbc = mcIR.toLong();
351 while (rofCount < ITSClusROFRec.size()) {
352 long rofbcMin = ITSClusROFRec[rofCount].getBCData().toLong() + ITSTimeBias, rofbcMax = rofbcMin + ITSROFLen;
353 if (gbc < rofbcMin) {
354 mITSOcc.push_back(0);
355 }
else if (gbc < rofbcMax) {
356 mITSOcc.push_back(ITSClusROFRec[rofCount].getNEntries());
363 if (mNTPCOccBinLengthInv > 0.f) {
364 mcVtx.occTPCV.resize(
params.nOccBinsDrift);
365 int grp = TMath::Max(1, TMath::Nint(
params.nTBPerOccBin * mNTPCOccBinLengthInv));
366 for (
int ib = 0; ib <
params.nOccBinsDrift; ib++) {
368 int tbs = occBin + TMath::Nint(ib *
params.nTBPerOccBin * mNTPCOccBinLengthInv);
369 for (
int ig = 0; ig < grp; ig++) {
370 if (tbs >= 0 && tbs <
int(mTBinClOccHist.size())) {
371 smb += mTBinClOccHist[tbs];
375 mcVtx.occTPCV[ib] = smb;
378 if (rofCount >= ITSClusROFRec.size()) {
379 mITSOcc.push_back(0);
384 int curSrcMC = 0, curEvMC = 0;
385 for (curSrcMC = 0; curSrcMC < (
int)mcReader.
getNSources(); curSrcMC++) {
387 LOGP(info,
"Source {}", curSrcMC);
390 bool okAccVtx =
true;
391 if (nev != (
int)mMCVtVec.size()) {
392 LOGP(
debug,
"source {} has {} events while {} MC vertices were booked", curSrcMC, nev, mMCVtVec.size());
394 if (nev > (
int)mMCVtVec.size()) {
398 for (curEvMC = 0; curEvMC < nev; curEvMC++) {
400 LOGP(info,
"Event {}", curEvMC);
402 mCurrMCTracks = &mcReader.
getTracks(curSrcMC, curEvMC);
405 auto&
pos = mMCVtVec[curEvMC].pos;
407 pos[0] = mCurrMCVertex.X();
408 pos[1] = mCurrMCVertex.Y();
409 pos[2] = mCurrMCVertex.Z();
412 for (
int itr = 0; itr < mCurrMCTracks->size(); itr++) {
413 processMCParticle(curSrcMC, curEvMC, itr);
418 for (
int id = 0;
id < mNCheckDecays;
id++) {
419 LOGP(info,
"Decay PDG={} : {} entries",
params.decayPDG[
id], mDecaysMaps[
id].size());
424 mRecProcStage =
true;
425 for (
int iv = 0; iv < nv; iv++) {
427 LOGP(info,
"processing PV {} of {}", iv, nv);
431 if (iv < (
int)pvvecLbl.size()) {
432 pvLbl = pvvecLbl[iv];
435 mMCVtVec[pvLbl.
getEventID()].recVtx.emplace_back(
RecPV{pvvec[iv], pvLbl});
438 const auto& vtref = vtxRefs[iv];
444 int idMin = vtref.getFirstEntryOfSource(is), idMax = idMin + vtref.getEntriesOfSource(is);
445 for (
int i = idMin;
i < idMax;
i++) {
446 auto vid = trackIndex[
i];
448 if (trc.getPt() <
params.minPt || std::abs(trc.getTgl()) >
params.maxTgl) {
454 auto entry = mSelMCTracks.find(lbl);
455 if (
entry == mSelMCTracks.end()) {
456 if (lbl.getSourceID() != curSrcMC || lbl.getEventID() != curEvMC) {
457 curSrcMC = lbl.getSourceID();
458 curEvMC = lbl.getEventID();
459 mCurrMCTracks = &mcReader.
getTracks(curSrcMC, curEvMC);
462 if (!acceptMCCharged((*mCurrMCTracks)[lbl.getTrackID()], lbl)) {
465 entry = mSelMCTracks.find(lbl);
467 auto& trackFamily =
entry->second;
468 if (vid.isAmbiguous()) {
469 if (trackFamily.contains(vid)) {
473 auto& trf = trackFamily.recTracks.emplace_back();
483 if (lblITS == trackFamily.mcTrackInfo.label) {
487 if (trcITSF.getPt() <
params.minPt || std::abs(trcITSF.getTgl()) >
params.maxTgl) {
490 auto entryOfFake = mSelMCTracks.find(lblITS);
491 if (entryOfFake == mSelMCTracks.end()) {
494 auto& trackFamilyOfFake = entryOfFake->second;
495 auto& trfOfFake = trackFamilyOfFake.recTracks.emplace_back();
500 LOGP(info,
"Matched rec track {} to MC track {}", vid.asString(),
entry->first.asString());
509 LOGP(info,
"collected {} MC tracks", mSelMCTracks.size());
510 if (
params.minTPCRefsToExtractClRes > 0 ||
params.storeTPCTrackRefs) {
511 processTPCTrackRefs();
515 for (
auto&
entry : mSelMCTracks) {
516 auto& trackFam =
entry.second;
517 auto& tracks = trackFam.recTracks;
519 if (tracks.empty()) {
523 LOGP(info,
"Processing MC track#{} {} -> {} reconstructed tracks", mcnt - 1,
entry.first.asString(), tracks.size());
526 std::sort(tracks.begin(), tracks.end(), [](
const RecTrack& lhs,
const RecTrack& rhs) {
527 const auto mskL = lhs.gid.getSourceDetectorsMask();
528 const auto mskR = rhs.gid.getSourceDetectorsMask();
529 bool itstpcL = mskL[DetID::ITS] && mskL[DetID::TPC], itstpcR = mskR[DetID::ITS] && mskR[DetID::TPC];
530 if (itstpcL && !itstpcR) {
533 return lhs.gid.getSource() > rhs.gid.getSource();
535 if (
params.storeTPCTrackRefs) {
536 auto rft = mSelTRefIdx.find(
entry.first);
537 if (rft != mSelTRefIdx.end()) {
538 auto rfent = rft->second;
539 for (
int irf = rfent.first; irf < rfent.second; irf++) {
540 trackFam.mcTrackInfo.trackRefsTPC.push_back(mSelTRefs[irf]);
546 for (
auto& tref : tracks) {
547 if (tref.gid.isSourceSet()) {
553 auto msk = tref.gid.getSourceDetectorsMask();
556 tref.pattITS = getITSPatt(gidSet[
GTrackID::ITS], tref.nClITS);
557 if (trackFam.entITS < 0) {
558 trackFam.entITS = tcnt;
561 if (lblITS.isFake()) {
564 if (lblITS == trackFam.mcTrackInfo.label) {
565 trackFam.entITSFound = tcnt;
573 if (msk[
DetID::TPC] && trackFam.entITSTPC < 0) {
574 trackFam.entITSTPC = tcnt;
582 tref.nClTPC = trtpc.getNClusters();
583 if (trtpc.hasBothSidesClusters()) {
586 fillTPCClusterInfo(trtpc, tref);
587 flagTPCClusters(trtpc,
entry.first);
588 if (trackFam.entTPC < 0) {
589 trackFam.entTPC = tcnt;
590 trackFam.tpcT0 = trtpc.getTime0();
596 float ts = 0, terr = 0;
601 const auto& itsBra = mITSROFBracket[mITSROF[tref.gid.getIndex()]];
602 tref.ts =
timeEst{itsBra.mean(), itsBra.delta() * SQRT12Inv};
605 LOGP(info,
"Invalid entry {} of {} getTrackMCLabel {}", tcnt, tracks.size(), tref.gid.asString());
609 if (trackFam.entITS > -1 && trackFam.entTPC > -1) {
614 if (propagateToRefX(trcTPC, trcITS)) {
615 trackFam.trackITSProp = trcITS;
616 trackFam.trackTPCProp = trcTPC;
618 trackFam.trackITSProp.invalidate();
619 trackFam.trackTPCProp.invalidate();
622 trackFam.trackITSProp.invalidate();
623 trackFam.trackTPCProp.invalidate();
629 auto v0s = recoData.getV0sIdx();
632 s += fmt::format(
" par {} Ntpccl={} Nitscl={} ",
f.mcTrackInfo.pdgParent,
f.mcTrackInfo.nTPCCl,
f.mcTrackInfo.nITSCl);
633 for (
auto& t :
f.recTracks) {
634 s += t.gid.asString();
639 for (
int svID; svID < (
int)v0s.size(); svID++) {
640 const auto& v0idx = v0s[svID];
641 int nOKProngs = 0, realMCSVID = -1;
642 int8_t decTypeID = -1;
643 for (
int ipr = 0; ipr < v0idx.getNProngs(); ipr++) {
644 auto mcl = recoData.getTrackMCLabel(v0idx.getProngID(ipr));
645 auto itl = mSelMCTracks.find(mcl);
646 if (itl == mSelMCTracks.end()) {
650 auto& trackFamily = itl->second;
651 int decayParentIndex = trackFamily.mcTrackInfo.parentEntry;
652 if (decayParentIndex < 0) {
656 realMCSVID = decayParentIndex;
657 decTypeID = trackFamily.mcTrackInfo.parentDecID;
659 LOGP(
debug,
"Prong{} {} comes from {}/{}", ipr, prpr(trackFamily), decTypeID, realMCSVID);
662 if (realMCSVID != decayParentIndex || decTypeID != trackFamily.mcTrackInfo.parentDecID) {
665 LOGP(
debug,
"Prong{} {} comes from {}/{}", ipr, prpr(trackFamily), decTypeID, realMCSVID);
668 if (nOKProngs == v0idx.getNProngs()) {
669 LOGP(
debug,
"Decay {}/{} was found", decTypeID, realMCSVID);
670 mDecaysMaps[decTypeID][realMCSVID].foundSVID = svID;
676 fillMCClusterInfo(recoData);
679 for (
auto&
entry : mSelMCTracks) {
680 auto& trackFam =
entry.second;
681 (*mDBGOut) <<
"tracks" <<
"tr=" << trackFam <<
"\n";
685 std::vector<TrackFamily> decFam;
686 for (
int id = 0;
id < mNCheckDecays;
id++) {
687 std::string decTreeName = fmt::format(
"dec{}",
params.decayPDG[
id]);
688 for (
const auto& dec : mDecaysMaps[
id]) {
691 for (
int idd = dec.daughterFirst; idd <= dec.daughterLast; idd++) {
692 auto dtLbl = mDecProdLblPool[idd];
693 const auto& dtFamily = mSelMCTracks[dtLbl];
694 if (dtFamily.mcTrackInfo.pdgParent != dec.pdg) {
695 LOGP(error,
"{}-th decay (pdg={}): {} in {}:{} range refers to MC track with pdgParent = {}",
id,
params.decayPDG[
id], idd, dec.daughterFirst, dec.daughterLast, dtFamily.mcTrackInfo.pdgParent);
699 decFam.push_back(dtFamily);
703 if (dec.foundSVID >= 0 && !refitV0(dec.foundSVID,
v0, recoData)) {
706 (*mDBGOut) << decTreeName.c_str() <<
"pdgPar=" << dec.pdg <<
"trPar=" << dec.parent <<
"prod=" << decFam <<
"found=" << dec.foundSVID <<
"sv=" <<
v0 <<
"\n";
711 for (
auto& mcVtx : mMCVtVec) {
712 std::sort(mcVtx.recVtx.begin(), mcVtx.recVtx.end(), [](
const RecPV& lhs,
const RecPV& rhs) {
713 return lhs.pv.getNContributors() > rhs.pv.getNContributors();
715 (*mDBGOut) <<
"mcVtxTree" <<
"mcVtx=" << mcVtx <<
"\n";
719void TrackMCStudy::processTPCTrackRefs()
721 constexpr float alpsec[18] = {0.174533, 0.523599, 0.872665, 1.221730, 1.570796, 1.919862, 2.268928, 2.617994, 2.967060, 3.316126, 3.665191, 4.014257, 4.363323, 4.712389, 5.061455, 5.410521, 5.759587, 6.108652};
722 constexpr float sinAlp[18] = {0.173648, 0.500000, 0.766044, 0.939693, 1.000000, 0.939693, 0.766044, 0.500000, 0.173648, -0.173648, -0.500000, -0.766044, -0.939693, -1.000000, -0.939693, -0.766044, -0.500000, -0.173648};
723 constexpr float cosAlp[18] = {0.984808, 0.866025, 0.642788, 0.342020, 0.000000, -0.342020, -0.642788, -0.866025, -0.984808, -0.984808, -0.866025, -0.642788, -0.342020, -0.000000, 0.342020, 0.642788, 0.866025, 0.984808};
725 for (
auto&
entry : mSelMCTracks) {
726 auto lb =
entry.first;
727 auto trspan = mcReader.getTrackRefs(lb.getSourceID(), lb.getEventID(), lb.getTrackID());
728 int q =
entry.second.mcTrackInfo.track.getCharge();
732 int ref0entry = mSelTRefs.size(), nrefsSel = 0;
733 for (
const auto& trf : trspan) {
734 if (trf.getDetectorId() != 1) {
737 float pT = std::sqrt(trf.Px() * trf.Px() + trf.Py() * trf.Py());
741 float secX, secY,
phi = std::atan2(trf.Y(), trf.X());
744 float phiPt = std::atan2(trf.Py(), trf.Px());
745 o2::math_utils::bringTo02Pi(phiPt);
746 auto dphiPt = phiPt - alpsec[sector];
754 float tgL = trf.Pz() / pT;
755 std::array<float, 5> pars = {secY, trf.Z(), std::sin(dphiPt), tgL, q / pT};
756 auto& refTrack = mSelTRefs.emplace_back(secX, alpsec[sector], pars);
757 refTrack.setUserField(uint16_t(sector));
760 if (nrefsSel <
params.minTPCRefsToExtractClRes) {
761 mSelTRefs.resize(ref0entry);
764 mSelTRefIdx[lb] = std::make_pair(ref0entry, ref0entry + nrefsSel);
773 const auto* TPCClMClab = recoData.
inputsTPCclusters->clusterIndex.clustersMCTruth;
777 for (uint8_t
row = 0;
row < 152;
row++) {
778 for (uint8_t sector = 0; sector < 36; sector++) {
779 unsigned int offs = TPCClusterIdxStruct.clusterOffset[sector][
row];
780 for (
unsigned int icl0 = 0; icl0 < TPCClusterIdxStruct.nClusters[sector][
row]; icl0++) {
781 const auto labels = TPCClMClab->getLabels(icl0 + offs);
783 for (
const auto& lbl : labels) {
784 if (!lbl.isValid()) {
789 const auto& clus = TPCClusterIdxStruct.clusters[sector][
row][icl0];
790 int tbinH =
int(clus.getTime() * mNTPCOccBinLengthInv);
791 clRes.contTracks.clear();
792 bool doClusRes = (
params.minTPCRefsToExtractClRes > 0) && (
params.rejectClustersResStat <= 0. || gRandom->Rndm() <
params.rejectClustersResStat);
793 for (
auto lbl : labels) {
794 bool corrAttach = lbl.isFake();
795 lbl.setFakeFlag(
false);
796 auto entry = mSelMCTracks.find(lbl);
797 if (
entry == mSelMCTracks.end()) {
800 auto& mctr =
entry->second.mcTrackInfo;
802 if (
row > mctr.maxTPCRow) {
803 mctr.maxTPCRow =
row;
804 mctr.maxTPCRowSect = sector;
806 }
else if (
row == 0 && mctr.nUsedPadRows == 0) {
809 if (
row < mctr.minTPCRow) {
810 mctr.minTPCRow =
row;
811 mctr.minTPCRowSect = sector;
813 if (mctr.minTPCRowSect == sector &&
row > mctr.maxTPCRowInner) {
814 mctr.maxTPCRowInner =
row;
821 auto entTRefIDsIt = mSelTRefIdx.find(lbl);
822 if (entTRefIDsIt == mSelTRefIdx.end()) {
826 mTPCCorrMapsLoader.Transform(sector,
row, clus.getPad(), clus.getTime(), xc, yc, zc, mctr.bcInTF / 8.);
828 const auto& entTRefIDs = entTRefIDsIt->second;
830 int entIDBelow = -1, entIDAbove = -1;
831 float xBelow = -1e6, xAbove = 1e6;
833 for (
int entID = entTRefIDs.first; entID < entTRefIDs.second; entID++) {
834 const auto& refTr = mSelTRefs[entID];
835 if (refTr.getUserField() != sector % 18) {
838 if ((refTr.getX() < xc) && (refTr.getX() > xBelow) && (refTr.getX() > xc -
params.maxTPCRefExtrap)) {
839 xBelow = refTr.getX();
842 if ((refTr.getX() > xc) && (refTr.getX() < xAbove) && (refTr.getX() < xc +
params.maxTPCRefExtrap)) {
843 xAbove = refTr.getX();
847 if ((entIDBelow < 0 && entIDAbove < 0) || (
params.requireTopBottomRefs && (entIDBelow < 0 || entIDAbove < 0))) {
852 bool okBelow = entIDBelow >= 0 && prop->PropagateToXBxByBz((tparBelow = mSelTRefs[entIDBelow]), xc, 0.99, 2.);
853 bool okAbove = entIDAbove >= 0 && prop->PropagateToXBxByBz((tparAbove = mSelTRefs[entIDAbove]), xc, 0.99, 2.);
854 if ((!okBelow && !okAbove) || (
params.requireTopBottomRefs && (!okBelow || !okAbove))) {
859 auto& clCont = clRes.contTracks.emplace_back();
860 clCont.corrAttach = corrAttach;
862 clCont.below = {mSelTRefs[entIDBelow].getX(), tparBelow.getY(), tparBelow.getZ()};
863 clCont.snp += tparBelow.getSnp();
864 clCont.tgl += tparBelow.getTgl();
865 clCont.q2pt += tparBelow.getQ2Pt();
869 clCont.above = {mSelTRefs[entIDAbove].getX(), tparAbove.getY(), tparAbove.getZ()};
870 clCont.snp += tparAbove.getSnp();
871 clCont.tgl += tparAbove.getTgl();
872 clCont.q2pt += tparAbove.getQ2Pt();
876 if (clRes.contTracks.size() == 1) {
877 int occBin = mctr.bcInTF / 8 * mNTPCOccBinLengthInv;
878 clRes.occ = occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]);
880 clCont.xyz = {xc, yc, zc};
887 clRes.contTracks.pop_back();
891 if (clRes.getNCont()) {
894 clRes.qtot = clus.getQtot();
895 clRes.qmax = clus.getQmax();
896 clRes.flags = clus.getFlags();
897 clRes.sigmaTimePacked = clus.sigmaTimePacked;
898 clRes.sigmaPadPacked = clus.sigmaPadPacked;
899 clRes.ncont = ncontLb;
904 }
else if (tbinH >=
int(mTBinClOccHist.size())) {
905 tbinH = (
int)mTBinClOccHist.size() - 1;
907 clRes.occBin = mTBinClOccHist[tbinH];
909 (*mDBGOut) <<
"clres" <<
"clr=" << clRes <<
"\n";
917 for (
unsigned int icl = 0; icl <
ITSClusters.size(); icl++) {
918 const auto labels = mcITSClusters->getLabels(icl);
919 for (
const auto& lbl : labels) {
920 auto entry = mSelMCTracks.find(lbl);
921 if (
entry == mSelMCTracks.end()) {
924 auto& mctr =
entry->second.mcTrackInfo;
933 bool refReached =
false;
934 constexpr float TgHalfSector = 0.17632698f;
942 if (fabs(trcTPC.getY()) <
par.XMatchingRef * TgHalfSector) {
950 if (!trcTPC.rotate(alphaNew) != 0) {
958 float alp = trcTPC.getAlpha();
975 if (mTPCVDriftHelper.accountCCDBInputs(matcher, obj)) {
978 if (mTPCCorrMapsLoader.accountCCDBInputs(matcher, obj)) {
982 LOG(info) <<
"ITS Alpide param updated";
984 par.printKeyValues();
996 int nROFs = ITSTrackROFRec.size();
998 mITSROFBracket.clear();
999 mITSROF.reserve(ITSTracksArray.size());
1000 mITSROFBracket.reserve(ITSTracksArray.size());
1001 for (
int irof = 0; irof < nROFs; irof++) {
1002 const auto& rofRec = ITSTrackROFRec[irof];
1003 long nBC = rofRec.getBCData().differenceInBC(recoData.
startIR);
1005 float tMax = tMin + mITSROFrameLengthMUS;
1006 mITSROFBracket.emplace_back(tMin, tMax);
1007 for (
int it = 0; it < rofRec.getNEntries(); it++) {
1008 mITSROF.push_back(irof);
1020bool TrackMCStudy::processMCParticle(
int src,
int ev,
int trid)
1022 const auto& mcPart = (*mCurrMCTracks)[trid];
1023 int pdg = mcPart.GetPdgCode();
1029 if (mcPart.T() <
params.decayMotherMaxT) {
1030 for (
int id = 0;
id < mNCheckDecays;
id++) {
1031 if (
params.decayPDG[
id] == std::abs(pdg)) {
1037 auto& decayPool = mDecaysMaps[decay];
1038 int idd0 = mcPart.getFirstDaughterTrackId(), idd1 = mcPart.getLastDaughterTrackId();
1039 int dtStart = mDecProdLblPool.size(), dtEnd = -1;
1043 for (
int idd = idd0; idd <= idd1; idd++) {
1044 const auto& product = (*mCurrMCTracks)[idd];
1046 if (!acceptMCCharged(product, lbld, decay)) {
1048 mDecProdLblPool.resize(dtStart);
1051 mDecProdLblPool.push_back(lbld);
1055 dtEnd = mDecProdLblPool.size();
1056 for (
int dtid = dtStart; dtid < dtEnd; dtid++) {
1057 mSelMCTracks[mDecProdLblPool[dtid]].mcTrackInfo.parentEntry = decayPool.size();
1058 mSelMCTracks[mDecProdLblPool[dtid]].mcTrackInfo.parentDecID = int8_t(decay);
1061 std::array<float, 3> xyz{(float)mcPart.GetStartVertexCoordinatesX(), (float)mcPart.GetStartVertexCoordinatesY(), (float)mcPart.GetStartVertexCoordinatesZ()};
1062 std::array<float, 3> pxyz{(float)mcPart.GetStartVertexMomentumX(), (float)mcPart.GetStartVertexMomentumY(), (float)mcPart.GetStartVertexMomentumZ()};
1063 decayPool.emplace_back(DecayRef{lbl,
1065 mcPart.GetPdgCode(), dtStart, dtEnd});
1067 LOGP(info,
"Adding MC parent pdg={} {}, with prongs in {}:{} range", pdg, lbl.asString(), dtStart, dtEnd);
1075 if (mSelMCTracks.find(lbl) == mSelMCTracks.end()) {
1076 res = acceptMCCharged(mcPart, lbl);
1089 if (mVerbose > 1 && followDecay > -1) {
1090 LOGP(info,
"rejecting decay {} prong : pdg={}, pT={}, tgL={}, r={}", followDecay, tr.
GetPdgCode(), tr.
GetPt(), tr.
GetTgl(), std::sqrt(tr.
R2()));
1095 float r2 = dx * dx + dy * dy;
1096 float posTgl2 = r2 > 1 && std::abs(dz) < 20 ? dz * dz / r2 : 0;
1097 if (posTgl2 >
params.maxPosTglMC *
params.maxPosTglMC) {
1098 if (mVerbose > 1 && followDecay > -1) {
1099 LOGP(info,
"rejecting decay {} prong : pdg={}, pT={}, tgL={}, dr={}, dz={} r={}, z={}, posTgl={}", followDecay, tr.
GetPdgCode(), tr.
GetPt(), tr.
GetTgl(), std::sqrt(r2), dz, std::sqrt(tr.
R2()), tr.
GetStartVertexCoordinatesZ(), std::sqrt(posTgl2));
1103 if (
params.requireITSorTPCTrackRefs) {
1106 for (
const auto& trf : trspan) {
1121 if (pPDG->Charge() == 0.) {
1124 return addMCParticle(tr, lb, pPDG);
1135 auto& mcEntry = mSelMCTracks[lb];
1136 mcEntry.mcTrackInfo.pdg = mcPart.
GetPdgCode();
1137 mcEntry.mcTrackInfo.track =
o2::track::TrackPar(xyz, pxyz, TMath::Nint(pPDG->Charge() / 3),
true);
1138 mcEntry.mcTrackInfo.label = lb;
1139 mcEntry.mcTrackInfo.bcInTF = mIntBC[lb.
getEventID()];
1140 mcEntry.mcTrackInfo.occTPC = mTPCOcc[lb.
getEventID()];
1141 mcEntry.mcTrackInfo.occITS = mITSOcc[lb.
getEventID()];
1142 mcEntry.mcTrackInfo.occTPCV = mMCVtVec[lb.
getEventID()].occTPCV;
1143 if (mRecProcStage) {
1144 mcEntry.mcTrackInfo.setAddedAtRecStage();
1147 mcEntry.mcTrackInfo.setPrimary();
1152 const auto& mcPartPar = (*mCurrMCTracks)[moth];
1153 mcEntry.mcTrackInfo.pdgParent = mcPartPar.GetPdgCode();
1171 if (svparam.mTPCTrackPhotonTune && isTPConly) {
1172 mFitterV0.setMaxDZIni(svparam.mTPCTrackMaxDZIni);
1173 mFitterV0.setMaxDXYIni(svparam.mTPCTrackMaxDXYIni);
1174 mFitterV0.setMaxChi2(svparam.mTPCTrackMaxChi2);
1175 mFitterV0.setCollinear(
true);
1177 int nCand = mFitterV0.process(seedP, seedN);
1178 if (svparam.mTPCTrackPhotonTune && isTPConly) {
1180 mFitterV0.setMaxDZIni(svparam.maxDZIni);
1181 mFitterV0.setMaxDXYIni(svparam.maxDXYIni);
1182 mFitterV0.setMaxChi2(svparam.maxChi2);
1183 mFitterV0.setCollinear(
false);
1189 if (!mFitterV0.isPropagateTracksToVertexDone(cand) && !mFitterV0.propagateTracksToVertex(cand)) {
1192 const auto& trPProp = mFitterV0.getTrack(0, cand);
1193 const auto& trNProp = mFitterV0.getTrack(1, cand);
1194 std::array<float, 3> pP{}, pN{};
1195 trPProp.getPxPyPzGlo(pP);
1196 trNProp.getPxPyPzGlo(pN);
1197 std::array<float, 3> pV0 = {pP[0] + pN[0], pP[1] + pN[1], pP[2] + pN[2]};
1198 auto p2V0 = pV0[0] * pV0[0] + pV0[1] * pV0[1] + pV0[2] * pV0[2];
1200 const auto v0XYZ = mFitterV0.getPCACandidatePos(cand);
1201 float dx = v0XYZ[0] - pv.getX(), dy = v0XYZ[1] - pv.getY(), dz = v0XYZ[2] - pv.getZ(), prodXYZv0 = dx * pV0[0] + dy * pV0[1] + dz * pV0[2];
1202 float cosPA = prodXYZv0 / std::sqrt((dx * dx + dy * dy + dz * dz) * p2V0);
1204 v0.
setDCA(mFitterV0.getChi2AtPCACandidate(cand));
1214 auto TPCRefitter = std::make_unique<o2::gpu::GPUO2InterfaceRefit>(&recoData.
inputsTPCclusters->clusterIndex, &mTPCCorrMapsLoader, prop->getNominalBz(),
1216 mNTPCOccBinLength = TPCRefitter->getParam()->rec.tpc.occupancyMapTimeBins;
1218 if (mNTPCOccBinLength > 1 && TPCOccMap.size()) {
1219 mNTPCOccBinLengthInv = 1. / mNTPCOccBinLength;
1222 mTBinClOcc.resize(nTPCOccBins);
1223 mTBinClOccHist.resize(nTPCOccBins);
1224 float sm = 0., tb = 0.5 * mNTPCOccBinLength;
1225 for (
int i = 0;
i < nTPCOccBins;
i++) {
1226 mTBinClOccHist[
i] = TPCRefitter->getParam()->GetUnscaledMult(tb);
1227 tb += mNTPCOccBinLength;
1229 for (
int i = nTPCOccBins;
i--;) {
1230 sm += mTBinClOccHist[
i];
1231 if (
i + sumBins < nTPCOccBins) {
1232 sm -= mTBinClOccHist[
i + sumBins];
1237 mTBinClOcc.resize(1);
1238 mTBinClOccHist.resize(1);
1244 std::vector<OutputSpec> outputs;
1246 {
"device-verbosity", VariantType::Int, 0, {
"Verbosity level"}},
1247 {
"dcay-vs-pt", VariantType::String,
"0.0105 + 0.0350 / pow(x, 1.1)", {
"Formula for global tracks DCAy vs pT cut"}},
1248 {
"min-tpc-clusters", VariantType::Int, 60, {
"Cut on TPC clusters"}},
1249 {
"max-tpc-dcay", VariantType::Float, 2.f, {
"Cut on TPC dcaY"}},
1250 {
"max-tpc-dcaz", VariantType::Float, 2.f, {
"Cut on TPC dcaZ"}},
1251 {
"min-x-prop", VariantType::Float, 6.f, {
"track should be propagated to this X at least"}}};
1252 auto dataRequest = std::make_shared<DataRequest>();
1254 dataRequest->requestTracks(srcTracks, useMC);
1255 dataRequest->requestClusters(srcClusters, useMC);
1256 dataRequest->requestPrimaryVertices(useMC);
1258 dataRequest->requestSecondaryVertices(useMC);
1262 auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(
false,
1268 dataRequest->inputs,
1275 AlgorithmSpec{adaptFromTask<TrackMCStudy>(dataRequest, ggRequest, srcTracks, sclOpts, checkSV)},
Helper class to access load maps from CCDB.
Defintions for N-prongs secondary vertex fit.
Definition of the GeometryManager class.
Helper for geometry and GRP related CCDB requests.
Global index for barrel track: provides provenance (detectors combination), index in respective array...
Utility functions for MC particles.
Configurable params for TPC ITS matching.
Definition of the Names Generator class.
Definition of the parameter class for the detector electronics.
Wrapper container for different reconstructed object types.
o2::track::TrackParCov TrackParCov
Configurable params for secondary vertexer.
Result of refitting TPC-ITS matched track.
Reference on ITS/MFT clusters set.
Helper class to extract VDrift from different sources.
Referenc on track indices contributing to the vertex, with possibility chose tracks from specific sou...
void setFakeFlag(bool v=true)
std::string asString() const
Double_t GetStartVertexMomentumZ() const
Double_t GetStartVertexMomentumX() const
Double_t GetStartVertexCoordinatesY() const
Double_t GetStartVertexCoordinatesZ() const
Double_t R2() const
production radius squared
Double_t GetStartVertexMomentumY() const
Double_t GetStartVertexCoordinatesX() const
Int_t GetPdgCode() const
Accessors.
Int_t getMotherTrackId() const
static TDatabasePDG * Instance()
void checkUpdates(o2::framework::ProcessingContext &pc)
static GRPGeomHelper & instance()
void setRequest(std::shared_ptr< GRPGeomRequest > req)
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
static const TrackMCStudyConfig & Instance()
Static class with identifiers, bitmasks and names for ALICE detectors.
T get(const char *key) const
ConfigParamRegistry const & options()
void setLumiScaleType(int32_t v)
void setLumiScaleMode(int32_t v)
void setCheckCTPIDCConsistency(bool v)
static constexpr int getLayer(int chipSW)
static bool isPhysicalPrimary(o2::MCTrack const &p, std::vector< o2::MCTrack > const &pcontainer)
std::vector< o2::InteractionTimeRecord > & getEventRecords(bool withQED=false)
bool initFromDigitContext(std::string_view filename)
DigitizationContext const * getDigitizationContext() const
size_t getNEvents(int source) const
Get number of events.
o2::dataformats::MCEventHeader const & getMCEventHeader(int source, int event) const
retrieves the MCEventHeader for a given eventID and sourceID
size_t getNSources() const
Get number of sources.
std::vector< MCTrack > const & getTracks(int source, int event) const
variant returning all tracks for source and event at once
void extractCCDBInputs(o2::framework::ProcessingContext &pc)
void updateVDrift(float vdriftCorr, float vdrifRef, float driftTimeOffset=0)
static void requestCCDBInputs(std::vector< o2::framework::InputSpec > &inputs, std::vector< o2::framework::ConfigParamSpec > &options, const CorrectionMapsLoaderGloOpts &gloOpts)
void init(o2::framework::InitContext &ic)
static void requestCCDBInputs(std::vector< o2::framework::InputSpec > &inputs, bool laser=true, bool itstpcTgl=true)
void extractCCDBInputs(o2::framework::ProcessingContext &pc, bool laser=true, bool itstpcTgl=true)
const VDriftCorrFact & getVDriftObject() const
static std::string_view getSourceName(Source s)
TrackMCStudy(std::shared_ptr< DataRequest > dr, std::shared_ptr< o2::base::GRPGeomRequest > gr, GTrackID::mask_t src, const o2::tpc::CorrectionMapsLoaderGloOpts &sclOpts, bool checkSV)
void endOfStream(EndOfStreamContext &ec) final
This is invoked whenever we have an EndOfStream event.
void init(InitContext &ic) final
void finaliseCCDB(ConcreteDataMatcher &matcher, void *obj) final
void run(ProcessingContext &pc) final
void process(const o2::globaltracking::RecoContainer &recoData)
~TrackMCStudy() final=default
GLenum const GLfloat * params
constexpr double LHCBunchSpacingMUS
constexpr int LHCMaxBunches
constexpr double LHCBunchSpacingNS
Node par(int index)
Parameters.
Defining PrimaryVertex explicitly as messageable.
std::vector< ConfigParamSpec > Options
detail::Bracket< float > Bracketf_t
float angle2Alpha(float phi)
int angle2Sector(float phi)
std::tuple< float, float > rotateZInv(float xG, float yG, float snAlp, float csAlp)
std::pair< int, o2::dataformats::VtxTrackIndex > VTIndexV
o2::dataformats::VtxTrackRef V2TRef
o2::framework::DataProcessorSpec getTrackMCStudySpec(o2::dataformats::GlobalTrackID::mask_t srcTracks, o2::dataformats::GlobalTrackID::mask_t srcClus, const o2::tpc::CorrectionMapsLoaderGloOpts &sclOpts, bool checkSV)
create a processor spec
o2::dataformats::VtxTrackIndex VTIndex
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Defining DataPointCompositeObject explicitly as copiable.
auto getITSTracks() const
GTrackID getITSContributorGID(GTrackID source) const
bool isTrackSourceLoaded(int src) const
o2::InteractionRecord startIR
GlobalIDSet getSingleDetectorRefs(GTrackID gidx) const
const o2::tpc::TrackTPC & getTPCTrack(GTrackID id) const
auto getITSClustersMCLabels() const
auto getPrimaryVertices() const
auto getPrimaryVertexMatchedTracks() const
const o2::tpc::ClusterNativeAccess & getTPCClusters() const
auto getITSABRefs() const
auto getTPCTracksClusterRefs() const
auto getPrimaryVertexMatchedTrackRefs() const
o2::MCCompLabel getTrackMCLabel(GTrackID id) const
GTrackID getTPCContributorGID(GTrackID source) const
auto getITSClustersROFRecords() const
const o2::track::TrackParCov & getTrackParam(GTrackID gidx) const
gsl::span< const unsigned char > clusterShMapTPC
externally set TPC clusters sharing map
void collectData(o2::framework::ProcessingContext &pc, const DataRequest &request)
const o2::track::TrackParCov & getTrackParamOut(GTrackID gidx) const
auto getPrimaryVertexMCLabels() const
const o2::dataformats::PrimaryVertex & getPrimaryVertex(int i) const
const o2::its::TrackITS & getITSTrack(GTrackID gid) const
std::unique_ptr< o2::tpc::internal::getWorkflowTPCInput_ret > inputsTPCclusters
auto getITSTracksROFRecords() const
void getTrackTime(GTrackID gid, float &t, float &tErr) const
gsl::span< const unsigned int > occupancyMapTPC
externally set TPC clusters occupancy map
auto getITSClusters() const
int lumiType
what estimator to used for corrections scaling: 0: no scaling, 1: CTP, 2: IDC
bool checkCTPIDCconsistency
int lumiMode
what corrections method to use: 0: classical scaling, 1: Using of the derivative map,...
float refTimeOffset
additive time offset reference (\mus)
float refVDrift
reference vdrift for which factor was extracted
float getTimeOffset() const
float timeOffsetCorr
additive time offset correction (\mus)
float corrFact
drift velocity correction factor (multiplicative)
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"