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)
97 void process(const
o2::globaltracking::RecoContainer& recoData);
100 void processTPCTrackRefs();
101 void loadTPCOccMap(const
o2::globaltracking::RecoContainer& recoData);
102 void fillMCClusterInfo(const
o2::globaltracking::RecoContainer& recoData);
103 void prepareITSData(const
o2::globaltracking::RecoContainer& recoData);
104 bool processMCParticle(
int src,
int ev,
int trid);
105 bool addMCParticle(const
MCTrack& mctr, const
o2::
MCCompLabel& lb, TParticlePDG* pPDG =
nullptr);
108 bool refitV0(
int i,
o2::dataformats::
V0&
v0, const
o2::globaltracking::RecoContainer& recoData);
110 float getDCAYCut(
float pt) const;
112 const
std::vector<
o2::
MCTrack>* mCurrMCTracks =
nullptr;
113 TVector3 mCurrMCVertex;
114 o2::tpc::VDriftHelper mTPCVDriftHelper{};
116 std::shared_ptr<DataRequest> mDataRequest;
117 std::shared_ptr<o2::base::GRPGeomRequest> mGGCCDBRequest;
118 std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOut;
119 std::vector<float> mTBinClOcc;
120 std::vector<float> mTBinClOccHist;
121 std::vector<long> mIntBC;
122 std::vector<float> mTPCOcc;
123 std::vector<int> mITSOcc;
124 bool mCheckSV =
false;
125 bool mRecProcStage =
false;
126 int mNTPCOccBinLength = 0;
127 float mNTPCOccBinLengthInv = -1.f;
129 float mITSTimeBiasMUS = 0.f;
130 float mITSROFrameLengthMUS = 0.f;
131 float mTPCTBinMUS = 0.;
133 int mNCheckDecays = 0;
137 std::vector<int> mITSROF;
138 std::vector<TBracket> mITSROFBracket;
139 std::vector<o2::MCCompLabel> mDecProdLblPool;
140 std::vector<MCVertex> mMCVtVec{};
146 int daughterFirst = -1;
147 int daughterLast = -1;
150 std::vector<std::vector<DecayRef>> mDecaysMaps;
151 std::unordered_map<o2::MCCompLabel, TrackFamily> mSelMCTracks;
152 std::unordered_map<o2::MCCompLabel, std::pair<int, int>> mSelTRefIdx;
153 std::vector<o2::track::TrackPar> mSelTRefs;
155 static constexpr float MaxSnp = 0.9;
163 mDBGOut = std::make_unique<o2::utils::TreeStreamRedirector>(
"trackMCStudy.root",
"recreate");
164 mVerbose = ic.
options().
get<
int>(
"device-verbosity");
167 for (
int id = 0;
id <
sizeof(
params.decayPDG) /
sizeof(
int);
id++) {
168 if (
params.decayPDG[
id] < 0) {
173 mDecaysMaps.resize(mNCheckDecays);
174 mTPCCorrMapsLoader.
init(ic);
180 for (
int i = 0;
i < mNCheckDecays;
i++) {
181 mDecaysMaps[
i].clear();
183 mDecProdLblPool.clear();
185 mCurrMCTracks =
nullptr;
188 updateTimeDependentParams(pc);
189 mRecProcStage =
false;
198 bool updateMaps =
false;
204 LOGP(info,
"Updating TPC fast transform map with new VDrift factor of {} wrt reference {} and DriftTimeOffset correction {} wrt {} from source {}",
214 static bool initOnceDone =
false;
219 LOGP(info,
"VertexTrackMatcher ITSROFrameLengthMUS:{}", mITSROFrameLengthMUS);
222 mTPCTBinMUS = elParam.ZbinWidth;
227 mFitterV0.setUseAbsDCA(svparam.useAbsDCA);
228 mFitterV0.setPropagateToPCA(
false);
229 mFitterV0.setMaxR(svparam.maxRIni);
230 mFitterV0.setMinParamChange(svparam.minParamChange);
231 mFitterV0.setMinRelChi2Change(svparam.minRelChi2Change);
232 mFitterV0.setMaxDZIni(svparam.maxDZIni);
233 mFitterV0.setMaxDXYIni(svparam.maxDXYIni);
234 mFitterV0.setMaxChi2(svparam.maxChi2);
236 mFitterV0.setUsePropagator(svparam.usePropagator);
237 mFitterV0.setRefitWithMatCorr(svparam.refitWithMatCorr);
238 mFitterV0.setMaxStep(svparam.maxStep);
239 mFitterV0.setMaxSnp(svparam.maxSnp);
240 mFitterV0.setMinXSeed(svparam.minXSeed);
247 constexpr float SQRT12Inv = 0.288675f;
254 int nv = vtxRefs.size();
258 prepareITSData(recoData);
259 loadTPCOccMap(recoData);
260 auto getITSPatt = [&](
GTrackID gid, uint8_t& ncl) {
264 ncl = itsTrf.getNClusters();
265 for (
int il = 0; il < 7; il++) {
266 if (itsTrf.hasHitOnLayer(il)) {
273 for (
int il = 0; il < 7; il++) {
274 if (itsTr.hasHitOnLayer(il)) {
285 uint8_t clSect = 0, clRow = 0;
289 trc.getClusterReference(clRefs, trc.getNClusterReferences() - 1, clSect, clRow, clIdx);
290 const auto& clus = tpcClusAcc.clusters[clSect][clRow][clIdx];
291 int padFromEdge =
int(clus.getPad()), npads = o2::gpu::GPUTPCGeometry::NPads(clRow);
292 if (padFromEdge > npads / 2) {
293 padFromEdge = npads - 1 - padFromEdge;
295 tref.padFromEdge = uint8_t(padFromEdge);
296 tref.lowestPadRow = clRow;
303 const auto* TPCClMClab = recoData.
inputsTPCclusters->clusterIndex.clustersMCTruth;
305 for (
int ic = 0; ic < trc.getNClusterReferences(); ic++) {
306 uint8_t clSect = 0, clRow = 0;
308 trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
309 auto labels = TPCClMClab->getLabels(clIdx + TPCClusterIdxStruct.clusterOffset[clSect][clRow]);
310 for (
auto& lbl : labels) {
325 unsigned int rofCount = 0;
327 for (
const auto& mcIR : mcEvRecords) {
328 long tbc = mcIR.differenceInBC(recoData.
startIR);
329 auto& mcVtx = mMCVtVec.emplace_back();
331 mcVtx.ID = mIntBC.size();
332 mIntBC.push_back(tbc);
333 int occBin = tbc / 8 * mNTPCOccBinLengthInv;
334 mTPCOcc.push_back(occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]));
336 long gbc = mcIR.toLong();
337 while (rofCount < ITSClusROFRec.size()) {
338 long rofbcMin = ITSClusROFRec[rofCount].getBCData().toLong() + ITSTimeBias, rofbcMax = rofbcMin + ITSROFLen;
339 if (gbc < rofbcMin) {
340 mITSOcc.push_back(0);
341 }
else if (gbc < rofbcMax) {
342 mITSOcc.push_back(ITSClusROFRec[rofCount].getNEntries());
349 if (mNTPCOccBinLengthInv > 0.f) {
350 mcVtx.occTPCV.resize(
params.nOccBinsDrift);
351 int grp = TMath::Max(1, TMath::Nint(
params.nTBPerOccBin * mNTPCOccBinLengthInv));
352 for (
int ib = 0; ib <
params.nOccBinsDrift; ib++) {
354 int tbs = occBin + TMath::Nint(ib *
params.nTBPerOccBin * mNTPCOccBinLengthInv);
355 for (
int ig = 0; ig < grp; ig++) {
356 if (tbs >= 0 && tbs <
int(mTBinClOccHist.size())) {
357 smb += mTBinClOccHist[tbs];
361 mcVtx.occTPCV[ib] = smb;
364 if (rofCount >= ITSClusROFRec.size()) {
365 mITSOcc.push_back(0);
370 int curSrcMC = 0, curEvMC = 0;
371 for (curSrcMC = 0; curSrcMC < (
int)mcReader.
getNSources(); curSrcMC++) {
373 LOGP(info,
"Source {}", curSrcMC);
376 bool okAccVtx =
true;
377 if (nev != (
int)mMCVtVec.size()) {
378 LOGP(
debug,
"source {} has {} events while {} MC vertices were booked", curSrcMC, nev, mMCVtVec.size());
380 if (nev > (
int)mMCVtVec.size()) {
384 for (curEvMC = 0; curEvMC < nev; curEvMC++) {
386 LOGP(info,
"Event {}", curEvMC);
388 mCurrMCTracks = &mcReader.
getTracks(curSrcMC, curEvMC);
391 auto&
pos = mMCVtVec[curEvMC].pos;
393 pos[0] = mCurrMCVertex.X();
394 pos[1] = mCurrMCVertex.Y();
395 pos[2] = mCurrMCVertex.Z();
398 for (
int itr = 0; itr < mCurrMCTracks->size(); itr++) {
399 processMCParticle(curSrcMC, curEvMC, itr);
404 for (
int id = 0;
id < mNCheckDecays;
id++) {
405 LOGP(info,
"Decay PDG={} : {} entries",
params.decayPDG[
id], mDecaysMaps[
id].size());
410 mRecProcStage =
true;
411 for (
int iv = 0; iv < nv; iv++) {
413 LOGP(info,
"processing PV {} of {}", iv, nv);
417 if (iv < (
int)pvvecLbl.size()) {
418 pvLbl = pvvecLbl[iv];
421 mMCVtVec[pvLbl.
getEventID()].recVtx.emplace_back(
RecPV{pvvec[iv], pvLbl});
424 const auto& vtref = vtxRefs[iv];
430 int idMin = vtref.getFirstEntryOfSource(is), idMax = idMin + vtref.getEntriesOfSource(is);
431 for (
int i = idMin;
i < idMax;
i++) {
432 auto vid = trackIndex[
i];
434 if (trc.getPt() <
params.minPt || std::abs(trc.getTgl()) >
params.maxTgl) {
440 auto entry = mSelMCTracks.find(lbl);
441 if (
entry == mSelMCTracks.end()) {
442 if (lbl.getSourceID() != curSrcMC || lbl.getEventID() != curEvMC) {
443 curSrcMC = lbl.getSourceID();
444 curEvMC = lbl.getEventID();
445 mCurrMCTracks = &mcReader.
getTracks(curSrcMC, curEvMC);
448 if (!acceptMCCharged((*mCurrMCTracks)[lbl.getTrackID()], lbl)) {
451 entry = mSelMCTracks.find(lbl);
453 auto& trackFamily =
entry->second;
454 if (vid.isAmbiguous()) {
455 if (trackFamily.contains(vid)) {
459 auto& trf = trackFamily.recTracks.emplace_back();
469 if (lblITS == trackFamily.mcTrackInfo.label) {
473 if (trcITSF.getPt() <
params.minPt || std::abs(trcITSF.getTgl()) >
params.maxTgl) {
476 auto entryOfFake = mSelMCTracks.find(lblITS);
477 if (entryOfFake == mSelMCTracks.end()) {
480 auto& trackFamilyOfFake = entryOfFake->second;
481 auto& trfOfFake = trackFamilyOfFake.recTracks.emplace_back();
486 LOGP(info,
"Matched rec track {} to MC track {}", vid.asString(),
entry->first.asString());
495 LOGP(info,
"collected {} MC tracks", mSelMCTracks.size());
496 if (
params.minTPCRefsToExtractClRes > 0) {
497 processTPCTrackRefs();
501 for (
auto&
entry : mSelMCTracks) {
502 auto& trackFam =
entry.second;
503 auto& tracks = trackFam.recTracks;
505 if (tracks.empty()) {
509 LOGP(info,
"Processing MC track#{} {} -> {} reconstructed tracks", mcnt - 1,
entry.first.asString(), tracks.size());
512 std::sort(tracks.begin(), tracks.end(), [](
const RecTrack& lhs,
const RecTrack& rhs) {
513 const auto mskL = lhs.gid.getSourceDetectorsMask();
514 const auto mskR = rhs.gid.getSourceDetectorsMask();
515 bool itstpcL = mskL[DetID::ITS] && mskL[DetID::TPC], itstpcR = mskR[DetID::ITS] && mskR[DetID::TPC];
516 if (itstpcL && !itstpcR) {
519 return lhs.gid.getSource() > rhs.gid.getSource();
523 for (
auto& tref : tracks) {
524 if (tref.gid.isSourceSet()) {
530 auto msk = tref.gid.getSourceDetectorsMask();
533 tref.pattITS = getITSPatt(gidSet[
GTrackID::ITS], tref.nClITS);
534 if (trackFam.entITS < 0) {
535 trackFam.entITS = tcnt;
538 if (lblITS.isFake()) {
541 if (lblITS == trackFam.mcTrackInfo.label) {
542 trackFam.entITSFound = tcnt;
550 if (msk[
DetID::TPC] && trackFam.entITSTPC < 0) {
551 trackFam.entITSTPC = tcnt;
559 tref.nClTPC = trtpc.getNClusters();
560 getLowestPadrow(trtpc, tref);
561 flagTPCClusters(trtpc,
entry.first);
562 if (trackFam.entTPC < 0) {
563 trackFam.entTPC = tcnt;
564 trackFam.tpcT0 = trtpc.getTime0();
570 float ts = 0, terr = 0;
575 const auto& itsBra = mITSROFBracket[mITSROF[tref.gid.getIndex()]];
576 tref.ts =
timeEst{itsBra.mean(), itsBra.delta() * SQRT12Inv};
579 LOGP(info,
"Invalid entry {} of {} getTrackMCLabel {}", tcnt, tracks.size(), tref.gid.asString());
583 if (trackFam.entITS > -1 && trackFam.entTPC > -1) {
584 auto vidITS = tracks[trackFam.entITS].gid;
585 auto vidTPC = tracks[trackFam.entTPC].gid;
588 if (propagateToRefX(trcTPC, trcITS)) {
589 trackFam.trackITSProp = trcITS;
590 trackFam.trackTPCProp = trcTPC;
592 trackFam.trackITSProp.invalidate();
593 trackFam.trackTPCProp.invalidate();
596 trackFam.trackITSProp.invalidate();
597 trackFam.trackTPCProp.invalidate();
603 auto v0s = recoData.getV0sIdx();
606 s += fmt::format(
" par {} Ntpccl={} Nitscl={} ",
f.mcTrackInfo.pdgParent,
f.mcTrackInfo.nTPCCl,
f.mcTrackInfo.nITSCl);
607 for (
auto& t :
f.recTracks) {
608 s += t.gid.asString();
613 for (
int svID; svID < (
int)v0s.size(); svID++) {
614 const auto& v0idx = v0s[svID];
615 int nOKProngs = 0, realMCSVID = -1;
616 int8_t decTypeID = -1;
617 for (
int ipr = 0; ipr < v0idx.getNProngs(); ipr++) {
618 auto mcl = recoData.getTrackMCLabel(v0idx.getProngID(ipr));
619 auto itl = mSelMCTracks.find(mcl);
620 if (itl == mSelMCTracks.end()) {
624 auto& trackFamily = itl->second;
625 int decayParentIndex = trackFamily.mcTrackInfo.parentEntry;
626 if (decayParentIndex < 0) {
630 realMCSVID = decayParentIndex;
631 decTypeID = trackFamily.mcTrackInfo.parentDecID;
633 LOGP(
debug,
"Prong{} {} comes from {}/{}", ipr, prpr(trackFamily), decTypeID, realMCSVID);
636 if (realMCSVID != decayParentIndex || decTypeID != trackFamily.mcTrackInfo.parentDecID) {
639 LOGP(
debug,
"Prong{} {} comes from {}/{}", ipr, prpr(trackFamily), decTypeID, realMCSVID);
642 if (nOKProngs == v0idx.getNProngs()) {
643 LOGP(
debug,
"Decay {}/{} was found", decTypeID, realMCSVID);
644 mDecaysMaps[decTypeID][realMCSVID].foundSVID = svID;
650 fillMCClusterInfo(recoData);
653 for (
auto&
entry : mSelMCTracks) {
654 auto& trackFam =
entry.second;
655 (*mDBGOut) <<
"tracks" <<
"tr=" << trackFam <<
"\n";
659 std::vector<TrackFamily> decFam;
660 for (
int id = 0;
id < mNCheckDecays;
id++) {
661 std::string decTreeName = fmt::format(
"dec{}",
params.decayPDG[
id]);
662 for (
const auto& dec : mDecaysMaps[
id]) {
665 for (
int idd = dec.daughterFirst; idd <= dec.daughterLast; idd++) {
666 auto dtLbl = mDecProdLblPool[idd];
667 const auto& dtFamily = mSelMCTracks[dtLbl];
668 if (dtFamily.mcTrackInfo.pdgParent != dec.pdg) {
669 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);
673 decFam.push_back(dtFamily);
677 if (dec.foundSVID >= 0 && !refitV0(dec.foundSVID,
v0, recoData)) {
680 (*mDBGOut) << decTreeName.c_str() <<
"pdgPar=" << dec.pdg <<
"trPar=" << dec.parent <<
"prod=" << decFam <<
"found=" << dec.foundSVID <<
"sv=" <<
v0 <<
"\n";
685 for (
auto& mcVtx : mMCVtVec) {
686 std::sort(mcVtx.recVtx.begin(), mcVtx.recVtx.end(), [](
const RecPV& lhs,
const RecPV& rhs) {
687 return lhs.pv.getNContributors() > rhs.pv.getNContributors();
689 (*mDBGOut) <<
"mcVtxTree" <<
"mcVtx=" << mcVtx <<
"\n";
693void TrackMCStudy::processTPCTrackRefs()
695 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};
696 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};
697 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};
699 for (
auto&
entry : mSelMCTracks) {
700 auto lb =
entry.first;
701 auto trspan = mcReader.getTrackRefs(lb.getSourceID(), lb.getEventID(), lb.getTrackID());
702 int q =
entry.second.mcTrackInfo.track.getCharge();
706 int ref0entry = mSelTRefs.size(), nrefsSel = 0;
707 for (
const auto& trf : trspan) {
708 if (trf.getDetectorId() != 1) {
711 float pT = std::sqrt(trf.Px() * trf.Px() + trf.Py() * trf.Py());
715 float secX, secY,
phi = std::atan2(trf.Y(), trf.X());
718 float phiPt = std::atan2(trf.Py(), trf.Px());
719 o2::math_utils::bringTo02Pi(phiPt);
720 auto dphiPt = phiPt - alpsec[sector];
728 float tgL = trf.Pz() / pT;
729 std::array<float, 5> pars = {secY, trf.Z(), std::sin(dphiPt), tgL, q / pT};
730 auto& refTrack = mSelTRefs.emplace_back(secX, alpsec[sector], pars);
731 refTrack.setUserField(uint16_t(sector));
734 if (nrefsSel <
params.minTPCRefsToExtractClRes) {
735 mSelTRefs.resize(ref0entry);
738 mSelTRefIdx[lb] = std::make_pair(ref0entry, ref0entry + nrefsSel);
747 const auto* TPCClMClab = recoData.
inputsTPCclusters->clusterIndex.clustersMCTruth;
751 for (uint8_t sector = 0; sector < 36; sector++) {
752 for (uint8_t
row = 0;
row < 152;
row++) {
753 unsigned int offs = TPCClusterIdxStruct.clusterOffset[sector][
row];
754 for (
unsigned int icl0 = 0; icl0 < TPCClusterIdxStruct.nClusters[sector][
row]; icl0++) {
755 const auto labels = TPCClMClab->getLabels(icl0 + offs);
757 for (
const auto& lbl : labels) {
758 if (!lbl.isValid()) {
763 const auto& clus = TPCClusterIdxStruct.clusters[sector][
row][icl0];
764 int tbinH =
int(clus.getTime() * mNTPCOccBinLengthInv);
765 clRes.contTracks.clear();
766 bool doClusRes = (
params.minTPCRefsToExtractClRes > 0) && (
params.rejectClustersResStat <= 0. || gRandom->Rndm() <
params.rejectClustersResStat);
767 for (
auto lbl : labels) {
768 bool corrAttach = lbl.isFake();
769 lbl.setFakeFlag(
false);
770 auto entry = mSelMCTracks.find(lbl);
771 if (
entry == mSelMCTracks.end()) {
774 auto& mctr =
entry->second.mcTrackInfo;
776 if (
row > mctr.maxTPCRow) {
777 mctr.maxTPCRow =
row;
778 mctr.maxTPCRowSect = sector;
780 }
else if (
row == 0 && mctr.nUsedPadRows == 0) {
783 if (
row < mctr.minTPCRow) {
784 mctr.minTPCRow =
row;
785 mctr.minTPCRowSect = sector;
787 if (mctr.minTPCRowSect == sector &&
row > mctr.maxTPCRowInner) {
788 mctr.maxTPCRowInner =
row;
795 auto entTRefIDsIt = mSelTRefIdx.find(lbl);
796 if (entTRefIDsIt == mSelTRefIdx.end()) {
800 mTPCCorrMapsLoader.Transform(sector,
row, clus.getPad(), clus.getTime(), xc, yc, zc, mctr.bcInTF / 8.);
802 const auto& entTRefIDs = entTRefIDsIt->second;
804 int entIDBelow = -1, entIDAbove = -1;
805 float xBelow = -1e6, xAbove = 1e6;
807 for (
int entID = entTRefIDs.first; entID < entTRefIDs.second; entID++) {
808 const auto& refTr = mSelTRefs[entID];
809 if (refTr.getUserField() != sector % 18) {
812 if ((refTr.getX() < xc) && (refTr.getX() > xBelow) && (refTr.getX() > xc -
params.maxTPCRefExtrap)) {
813 xBelow = refTr.getX();
816 if ((refTr.getX() > xc) && (refTr.getX() < xAbove) && (refTr.getX() < xc +
params.maxTPCRefExtrap)) {
817 xAbove = refTr.getX();
821 if ((entIDBelow < 0 && entIDAbove < 0) || (
params.requireTopBottomRefs && (entIDBelow < 0 || entIDAbove < 0))) {
826 bool okBelow = entIDBelow >= 0 && prop->PropagateToXBxByBz((tparBelow = mSelTRefs[entIDBelow]), xc, 0.99, 2.);
827 bool okAbove = entIDAbove >= 0 && prop->PropagateToXBxByBz((tparAbove = mSelTRefs[entIDAbove]), xc, 0.99, 2.);
828 if ((!okBelow && !okAbove) || (
params.requireTopBottomRefs && (!okBelow || !okAbove))) {
833 auto& clCont = clRes.contTracks.emplace_back();
834 clCont.corrAttach = corrAttach;
836 clCont.below = {mSelTRefs[entIDBelow].getX(), tparBelow.getY(), tparBelow.getZ()};
837 clCont.snp += tparBelow.getSnp();
838 clCont.tgl += tparBelow.getTgl();
839 clCont.q2pt += tparBelow.getQ2Pt();
843 clCont.above = {mSelTRefs[entIDAbove].getX(), tparAbove.getY(), tparAbove.getZ()};
844 clCont.snp += tparAbove.getSnp();
845 clCont.tgl += tparAbove.getTgl();
846 clCont.q2pt += tparAbove.getQ2Pt();
850 if (clRes.contTracks.size() == 1) {
851 int occBin = mctr.bcInTF / 8 * mNTPCOccBinLengthInv;
852 clRes.occ = occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]);
854 clCont.xyz = {xc, yc, zc};
861 clRes.contTracks.pop_back();
865 if (clRes.getNCont()) {
868 clRes.qtot = clus.getQtot();
869 clRes.qmax = clus.getQmax();
870 clRes.flags = clus.getFlags();
871 clRes.sigmaTimePacked = clus.sigmaTimePacked;
872 clRes.sigmaPadPacked = clus.sigmaPadPacked;
873 clRes.ncont = ncontLb;
878 }
else if (tbinH >=
int(mTBinClOccHist.size())) {
879 tbinH = (
int)mTBinClOccHist.size() - 1;
881 clRes.occBin = mTBinClOccHist[tbinH];
883 (*mDBGOut) <<
"clres" <<
"clr=" << clRes <<
"\n";
891 for (
unsigned int icl = 0; icl <
ITSClusters.size(); icl++) {
892 const auto labels = mcITSClusters->getLabels(icl);
893 for (
const auto& lbl : labels) {
894 auto entry = mSelMCTracks.find(lbl);
895 if (
entry == mSelMCTracks.end()) {
898 auto& mctr =
entry->second.mcTrackInfo;
907 bool refReached =
false;
908 constexpr float TgHalfSector = 0.17632698f;
916 if (fabs(trcTPC.getY()) <
par.XMatchingRef * TgHalfSector) {
924 if (!trcTPC.rotate(alphaNew) != 0) {
932 float alp = trcTPC.getAlpha();
949 if (mTPCVDriftHelper.accountCCDBInputs(matcher, obj)) {
952 if (mTPCCorrMapsLoader.accountCCDBInputs(matcher, obj)) {
956 LOG(info) <<
"ITS Alpide param updated";
958 par.printKeyValues();
970 int nROFs = ITSTrackROFRec.size();
972 mITSROFBracket.clear();
973 mITSROF.reserve(ITSTracksArray.size());
974 mITSROFBracket.reserve(ITSTracksArray.size());
975 for (
int irof = 0; irof < nROFs; irof++) {
976 const auto& rofRec = ITSTrackROFRec[irof];
977 long nBC = rofRec.getBCData().differenceInBC(recoData.
startIR);
979 float tMax = tMin + mITSROFrameLengthMUS;
980 mITSROFBracket.emplace_back(tMin, tMax);
981 for (
int it = 0; it < rofRec.getNEntries(); it++) {
982 mITSROF.push_back(irof);
994bool TrackMCStudy::processMCParticle(
int src,
int ev,
int trid)
996 const auto& mcPart = (*mCurrMCTracks)[trid];
997 int pdg = mcPart.GetPdgCode();
1003 if (mcPart.T() <
params.decayMotherMaxT) {
1004 for (
int id = 0;
id < mNCheckDecays;
id++) {
1005 if (
params.decayPDG[
id] == std::abs(pdg)) {
1011 auto& decayPool = mDecaysMaps[decay];
1012 int idd0 = mcPart.getFirstDaughterTrackId(), idd1 = mcPart.getLastDaughterTrackId();
1013 int dtStart = mDecProdLblPool.size(), dtEnd = -1;
1017 for (
int idd = idd0; idd <= idd1; idd++) {
1018 const auto& product = (*mCurrMCTracks)[idd];
1020 if (!acceptMCCharged(product, lbld, decay)) {
1022 mDecProdLblPool.resize(dtStart);
1025 mDecProdLblPool.push_back(lbld);
1029 dtEnd = mDecProdLblPool.size();
1030 for (
int dtid = dtStart; dtid < dtEnd; dtid++) {
1031 mSelMCTracks[mDecProdLblPool[dtid]].mcTrackInfo.parentEntry = decayPool.size();
1032 mSelMCTracks[mDecProdLblPool[dtid]].mcTrackInfo.parentDecID = int8_t(decay);
1035 std::array<float, 3> xyz{(float)mcPart.GetStartVertexCoordinatesX(), (float)mcPart.GetStartVertexCoordinatesY(), (float)mcPart.GetStartVertexCoordinatesZ()};
1036 std::array<float, 3> pxyz{(float)mcPart.GetStartVertexMomentumX(), (float)mcPart.GetStartVertexMomentumY(), (float)mcPart.GetStartVertexMomentumZ()};
1037 decayPool.emplace_back(DecayRef{lbl,
1039 mcPart.GetPdgCode(), dtStart, dtEnd});
1041 LOGP(info,
"Adding MC parent pdg={} {}, with prongs in {}:{} range", pdg, lbl.asString(), dtStart, dtEnd);
1049 if (mSelMCTracks.find(lbl) == mSelMCTracks.end()) {
1050 res = acceptMCCharged(mcPart, lbl);
1063 if (mVerbose > 1 && followDecay > -1) {
1064 LOGP(info,
"rejecting decay {} prong : pdg={}, pT={}, tgL={}, r={}", followDecay, tr.
GetPdgCode(), tr.
GetPt(), tr.
GetTgl(), std::sqrt(tr.
R2()));
1069 float r2 = dx * dx + dy * dy;
1070 float posTgl2 = r2 > 1 && std::abs(dz) < 20 ? dz * dz / r2 : 0;
1071 if (posTgl2 >
params.maxPosTglMC *
params.maxPosTglMC) {
1072 if (mVerbose > 1 && followDecay > -1) {
1073 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));
1077 if (
params.requireITSorTPCTrackRefs) {
1080 for (
const auto& trf : trspan) {
1095 if (pPDG->Charge() == 0.) {
1098 return addMCParticle(tr, lb, pPDG);
1109 auto& mcEntry = mSelMCTracks[lb];
1110 mcEntry.mcTrackInfo.pdg = mcPart.
GetPdgCode();
1111 mcEntry.mcTrackInfo.track =
o2::track::TrackPar(xyz, pxyz, TMath::Nint(pPDG->Charge() / 3),
true);
1112 mcEntry.mcTrackInfo.label = lb;
1113 mcEntry.mcTrackInfo.bcInTF = mIntBC[lb.
getEventID()];
1114 mcEntry.mcTrackInfo.occTPC = mTPCOcc[lb.
getEventID()];
1115 mcEntry.mcTrackInfo.occITS = mITSOcc[lb.
getEventID()];
1116 mcEntry.mcTrackInfo.occTPCV = mMCVtVec[lb.
getEventID()].occTPCV;
1117 if (mRecProcStage) {
1118 mcEntry.mcTrackInfo.setAddedAtRecStage();
1121 mcEntry.mcTrackInfo.setPrimary();
1126 const auto& mcPartPar = (*mCurrMCTracks)[moth];
1127 mcEntry.mcTrackInfo.pdgParent = mcPartPar.GetPdgCode();
1145 if (svparam.mTPCTrackPhotonTune && isTPConly) {
1146 mFitterV0.setMaxDZIni(svparam.mTPCTrackMaxDZIni);
1147 mFitterV0.setMaxDXYIni(svparam.mTPCTrackMaxDXYIni);
1148 mFitterV0.setMaxChi2(svparam.mTPCTrackMaxChi2);
1149 mFitterV0.setCollinear(
true);
1151 int nCand = mFitterV0.process(seedP, seedN);
1152 if (svparam.mTPCTrackPhotonTune && isTPConly) {
1154 mFitterV0.setMaxDZIni(svparam.maxDZIni);
1155 mFitterV0.setMaxDXYIni(svparam.maxDXYIni);
1156 mFitterV0.setMaxChi2(svparam.maxChi2);
1157 mFitterV0.setCollinear(
false);
1163 if (!mFitterV0.isPropagateTracksToVertexDone(cand) && !mFitterV0.propagateTracksToVertex(cand)) {
1166 const auto& trPProp = mFitterV0.getTrack(0, cand);
1167 const auto& trNProp = mFitterV0.getTrack(1, cand);
1168 std::array<float, 3> pP{}, pN{};
1169 trPProp.getPxPyPzGlo(pP);
1170 trNProp.getPxPyPzGlo(pN);
1171 std::array<float, 3> pV0 = {pP[0] + pN[0], pP[1] + pN[1], pP[2] + pN[2]};
1172 auto p2V0 = pV0[0] * pV0[0] + pV0[1] * pV0[1] + pV0[2] * pV0[2];
1174 const auto v0XYZ = mFitterV0.getPCACandidatePos(cand);
1175 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];
1176 float cosPA = prodXYZv0 / std::sqrt((dx * dx + dy * dy + dz * dz) * p2V0);
1178 v0.
setDCA(mFitterV0.getChi2AtPCACandidate(cand));
1188 auto TPCRefitter = std::make_unique<o2::gpu::GPUO2InterfaceRefit>(&recoData.
inputsTPCclusters->clusterIndex, &mTPCCorrMapsLoader, prop->getNominalBz(),
1190 mNTPCOccBinLength = TPCRefitter->getParam()->rec.tpc.occupancyMapTimeBins;
1192 if (mNTPCOccBinLength > 1 && TPCOccMap.size()) {
1193 mNTPCOccBinLengthInv = 1. / mNTPCOccBinLength;
1196 mTBinClOcc.resize(nTPCOccBins);
1197 mTBinClOccHist.resize(nTPCOccBins);
1198 float sm = 0., tb = 0.5 * mNTPCOccBinLength;
1199 for (
int i = 0;
i < nTPCOccBins;
i++) {
1200 mTBinClOccHist[
i] = TPCRefitter->getParam()->GetUnscaledMult(tb);
1201 tb += mNTPCOccBinLength;
1203 for (
int i = nTPCOccBins;
i--;) {
1204 sm += mTBinClOccHist[
i];
1205 if (
i + sumBins < nTPCOccBins) {
1206 sm -= mTBinClOccHist[
i + sumBins];
1211 mTBinClOcc.resize(1);
1212 mTBinClOccHist.resize(1);
1218 std::vector<OutputSpec> outputs;
1220 {
"device-verbosity", VariantType::Int, 0, {
"Verbosity level"}},
1221 {
"dcay-vs-pt", VariantType::String,
"0.0105 + 0.0350 / pow(x, 1.1)", {
"Formula for global tracks DCAy vs pT cut"}},
1222 {
"min-tpc-clusters", VariantType::Int, 60, {
"Cut on TPC clusters"}},
1223 {
"max-tpc-dcay", VariantType::Float, 2.f, {
"Cut on TPC dcaY"}},
1224 {
"max-tpc-dcaz", VariantType::Float, 2.f, {
"Cut on TPC dcaZ"}},
1225 {
"min-x-prop", VariantType::Float, 6.f, {
"track should be propagated to this X at least"}}};
1226 auto dataRequest = std::make_shared<DataRequest>();
1228 dataRequest->requestTracks(srcTracks, useMC);
1229 dataRequest->requestClusters(srcClusters, useMC);
1230 dataRequest->requestPrimaryVertices(useMC);
1232 dataRequest->requestSecondaryVertices(useMC);
1236 auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(
false,
1242 dataRequest->inputs,
1249 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)
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)
recalculate inverse correction
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
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
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
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"