diff --git a/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyTypes.h b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyTypes.h index 4766dc6787351..5fe2ab542e25c 100644 --- a/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyTypes.h +++ b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyTypes.h @@ -292,7 +292,8 @@ struct ClResTPC { struct RecPV { o2::dataformats::PrimaryVertex pv{}; o2::MCEventLabel mcEvLbl{}; - ClassDefNV(RecPV, 1); + int id = -1; + ClassDefNV(RecPV, 2); }; struct MCVertex { diff --git a/Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudy.cxx b/Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudy.cxx index 9dba400fe6edc..4834e7a268783 100644 --- a/Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudy.cxx +++ b/Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudy.cxx @@ -106,8 +106,10 @@ class TrackMCStudy : public Task bool addMCParticle(const MCTrack& mctr, const o2::MCCompLabel& lb, TParticlePDG* pPDG = nullptr); bool acceptMCCharged(const MCTrack& tr, const o2::MCCompLabel& lb, int followDec = -1); bool propagateToRefX(o2::track::TrackParCov& trcTPC, o2::track::TrackParCov& trcITS); - bool refitV0(int i, o2::dataformats::V0& v0, const o2::globaltracking::RecoContainer& recoData); + bool refitV0(int i, o2::dataformats::V0& v0, o2::dataformats::V0Index& v0id, const o2::globaltracking::RecoContainer& recoData); + bool fitV0(GTrackID prID0, GTrackID prID1, int idPV, o2::dataformats::V0& v0, o2::dataformats::V0Index& v0id, const o2::globaltracking::RecoContainer& recoData); void updateTimeDependentParams(ProcessingContext& pc); + o2::track::TrackParCov getCorrectedTPCTrack(const o2::tpc::TrackTPC& trTPC, float timeMUS); float getDCAYCut(float pt) const; const std::vector* mCurrMCTracks = nullptr; @@ -126,6 +128,7 @@ class TrackMCStudy : public Task bool mRecProcStage = false; //< flag that the MC particle was added only at the stage of reco tracks processing int mNTPCOccBinLength = 0; ///< TPC occ. histo bin length in TBs float mNTPCOccBinLengthInv = -1.f; + float mVDriftTB = -1.f; int mVerbose = 0; float mITSTimeBiasMUS = 0.f; float mITSROFrameLengthMUS = 0.f; ///< ITS RO frame in mus @@ -146,7 +149,7 @@ class TrackMCStudy : public Task int pdg = 0; int daughterFirst = -1; int daughterLast = -1; - int foundSVID = -1; + std::vector foundSVIdx; }; std::vector> mDecaysMaps; // for every parent particle to watch, store its label and entries of 1st/last decay product labels in mDecProdLblPool std::unordered_map mSelMCTracks; @@ -253,8 +256,6 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData) auto vtxRefs = recoData.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs auto prop = o2::base::Propagator::Instance(); int nv = vtxRefs.size(); - float vdriftTB = mTPCVDriftHelper.getVDriftObject().getVDrift() * o2::tpc::ParameterElectronics::Instance().ZbinWidth; // VDrift expressed in cm/TimeBin - float itsBias = 0.5 * mITSROFrameLengthMUS + o2::itsmft::DPLAlpideParam::Instance().roFrameBiasInBC * o2::constants::lhc::LHCBunchSpacingMUS; // ITS time is supplied in \mus as beginning of ROF prepareITSData(recoData); loadTPCOccMap(recoData); @@ -432,7 +433,7 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData) pvLbl = pvvecLbl[iv]; pvID = iv; if (pvLbl.isSet() && pvLbl.getEventID() < mMCVtVec.size()) { - mMCVtVec[pvLbl.getEventID()].recVtx.emplace_back(RecPV{pvvec[iv], pvLbl}); + mMCVtVec[pvLbl.getEventID()].recVtx.emplace_back(RecPV{pvvec[iv], pvLbl, iv}); } } const auto& vtref = vtxRefs[iv]; @@ -636,7 +637,7 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData) } return s; }; - for (int svID; svID < (int)v0s.size(); svID++) { + for (int svID = 0; svID < (int)v0s.size(); svID++) { const auto& v0idx = v0s[svID]; int nOKProngs = 0, realMCSVID = -1; int8_t decTypeID = -1; @@ -667,7 +668,7 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData) } if (nOKProngs == v0idx.getNProngs()) { // all prongs are from the decay of MC parent which deemed to be interesting, flag it LOGP(debug, "Decay {}/{} was found", decTypeID, realMCSVID); - mDecaysMaps[decTypeID][realMCSVID].foundSVID = svID; + mDecaysMaps[decTypeID][realMCSVID].foundSVIdx.push_back(svID); } } } @@ -675,6 +676,13 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData) // collect ITS/TPC cluster info for selected MC particles fillMCClusterInfo(recoData); + for (auto& mcVtx : mMCVtVec) { // sort rec.vertices in mult. order + std::sort(mcVtx.recVtx.begin(), mcVtx.recVtx.end(), [](const RecPV& lhs, const RecPV& rhs) { + return lhs.pv.getNContributors() > rhs.pv.getNContributors(); + }); + (*mDBGOut) << "mcVtxTree" << "mcVtx=" << mcVtx << "\n"; + } + // single tracks for (auto& entry : mSelMCTracks) { auto& trackFam = entry.second; @@ -682,7 +690,10 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData) } // decays + std::vector foundV0s, testV0s; + std::vector foundV0IDs, testV0IDs; std::vector decFam; + std::vector corrPVs; for (int id = 0; id < mNCheckDecays; id++) { std::string decTreeName = fmt::format("dec{}", params.decayPDG[id]); for (const auto& dec : mDecaysMaps[id]) { @@ -699,21 +710,51 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData) decFam.push_back(dtFamily); } if (!skip) { - o2::dataformats::V0 v0; - if (dec.foundSVID >= 0 && !refitV0(dec.foundSVID, v0, recoData)) { - v0.invalidate(); + corrPVs.clear(); + const auto& recPVs = mMCVtVec[dec.mother.getEventID()].recVtx; + for (const auto& rpv : recPVs) { + corrPVs.push_back(rpv.id); + } + foundV0s.clear(); + foundV0s.resize(dec.foundSVIdx.size()); + foundV0IDs.resize(dec.foundSVIdx.size()); + testV0s.clear(); + testV0IDs.clear(); + int corrSVPV = -1; + for (int ivf = 0; ivf < (int)dec.foundSVIdx.size(); ivf++) { + if (!refitV0(dec.foundSVIdx[ivf], foundV0s[ivf], foundV0IDs[ivf], recoData)) { + foundV0s[ivf].invalidate(); + } + if (recPVs.size() && recoData.getV0sIdx()[dec.foundSVIdx[ivf]].getVertexID() == recPVs[0].id) { + corrSVPV = ivf; + } + } + if (decFam.size() == 2 && recPVs.size()) { + + for (const auto& pr0 : decFam[0].recTracks) { + if (pr0.gid.getSource() == GTrackID::ITS) { // don't consider ITS only tracks + continue; + } + for (const auto& pr1 : decFam[1].recTracks) { + if (pr1.gid.getSource() == GTrackID::ITS) { // don't consider ITS only tracks + continue; + } + auto& v0t = testV0s.emplace_back(); + auto& v0tid = testV0IDs.emplace_back(); + if (!fitV0(pr0.gid, pr1.gid, recPVs[0].id, v0t, v0tid, recoData)) { + v0t.invalidate(); + } + } + } } - (*mDBGOut) << decTreeName.c_str() << "pdgPar=" << dec.pdg << "trPar=" << dec.parent << "prod=" << decFam << "found=" << dec.foundSVID << "sv=" << v0 << "\n"; + (*mDBGOut) << decTreeName.c_str() << "pdgPar=" << dec.pdg << "trPar=" << dec.parent << "prod=" << decFam << "foundSVIdx=" + << dec.foundSVIdx << "corrSVPV=" << corrSVPV << "corrPV=" << corrPVs + << "foundSVID=" << foundV0IDs << "foundSV=" << foundV0s + << "testSVID=" << testV0IDs << "testSV=" << testV0s + << "\n"; } } } - - for (auto& mcVtx : mMCVtVec) { // sort rec.vertices in mult. order - std::sort(mcVtx.recVtx.begin(), mcVtx.recVtx.end(), [](const RecPV& lhs, const RecPV& rhs) { - return lhs.pv.getNContributors() > rhs.pv.getNContributors(); - }); - (*mDBGOut) << "mcVtxTree" << "mcVtx=" << mcVtx << "\n"; - } } void TrackMCStudy::processTPCTrackRefs() @@ -973,6 +1014,7 @@ void TrackMCStudy::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) return; } if (mTPCVDriftHelper.accountCCDBInputs(matcher, obj)) { + mVDriftTB = mTPCVDriftHelper.getVDriftObject().getVDrift() * o2::tpc::ParameterElectronics::Instance().ZbinWidth; // VDrift expressed in cm/TimeBin return; } if (mTPCCorrMapsLoader.accountCCDBInputs(matcher, obj)) { @@ -1161,12 +1203,19 @@ bool TrackMCStudy::addMCParticle(const MCTrack& mcPart, const o2::MCCompLabel& l return true; } -bool TrackMCStudy::refitV0(int i, o2::dataformats::V0& v0, const o2::globaltracking::RecoContainer& recoData) +bool TrackMCStudy::refitV0(int i, o2::dataformats::V0& v0, o2::dataformats::V0Index& v0id, const o2::globaltracking::RecoContainer& recoData) { const auto& id = recoData.getV0sIdx()[i]; - auto seedP = recoData.getTrackParam(id.getProngID(0)); - auto seedN = recoData.getTrackParam(id.getProngID(1)); - bool isTPConly = (id.getProngID(0).getSource() == GTrackID::TPC) || (id.getProngID(1).getSource() == GTrackID::TPC); + return fitV0(id.getProngID(0), id.getProngID(1), id.getVertexID(), v0, v0id, recoData); +} + +bool TrackMCStudy::fitV0(GTrackID prID0, GTrackID prID1, int idPV, o2::dataformats::V0& v0, o2::dataformats::V0Index& v0id, const o2::globaltracking::RecoContainer& recoData) +{ + const auto& pv = recoData.getPrimaryVertex(idPV); + + auto seedP = recoData.getTrackParam(prID0); + auto seedN = recoData.getTrackParam(prID1); + bool isTPConlyP = (prID0.getSource() == GTrackID::TPC), isTPConlyN = (prID1.getSource() == GTrackID::TPC), isTPConly = isTPConlyP || isTPConlyN; const auto& svparam = o2::vertexing::SVertexerParams::Instance(); if (svparam.mTPCTrackPhotonTune && isTPConly) { mFitterV0.setMaxDZIni(svparam.mTPCTrackMaxDZIni); @@ -1174,7 +1223,8 @@ bool TrackMCStudy::refitV0(int i, o2::dataformats::V0& v0, const o2::globaltrack mFitterV0.setMaxChi2(svparam.mTPCTrackMaxChi2); mFitterV0.setCollinear(true); } - int nCand = mFitterV0.process(seedP, seedN); + int nCand = mFitterV0.process(isTPConlyP ? getCorrectedTPCTrack(recoData.getTPCTrack(prID0), pv.getTimeStamp().getTimeStamp()) : recoData.getTrackParam(prID0), + isTPConlyN ? getCorrectedTPCTrack(recoData.getTPCTrack(prID1), pv.getTimeStamp().getTimeStamp()) : recoData.getTrackParam(prID1)); if (svparam.mTPCTrackPhotonTune && isTPConly) { // restore // Reset immediately to the defaults mFitterV0.setMaxDZIni(svparam.maxDZIni); @@ -1182,6 +1232,9 @@ bool TrackMCStudy::refitV0(int i, o2::dataformats::V0& v0, const o2::globaltrack mFitterV0.setMaxChi2(svparam.maxChi2); mFitterV0.setCollinear(false); } + v0id.setVertexID(idPV); + v0id.setProngID(0, prID0); + v0id.setProngID(1, prID1); if (nCand == 0) { // discard this pair return false; } @@ -1196,7 +1249,6 @@ bool TrackMCStudy::refitV0(int i, o2::dataformats::V0& v0, const o2::globaltrack trNProp.getPxPyPzGlo(pN); std::array pV0 = {pP[0] + pN[0], pP[1] + pN[1], pP[2] + pN[2]}; auto p2V0 = pV0[0] * pV0[0] + pV0[1] * pV0[1] + pV0[2] * pV0[2]; - const auto& pv = recoData.getPrimaryVertex(id.getVertexID()); const auto v0XYZ = mFitterV0.getPCACandidatePos(cand); 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]; float cosPA = prodXYZv0 / std::sqrt((dx * dx + dy * dy + dz * dz) * p2V0); @@ -1239,6 +1291,19 @@ void TrackMCStudy::loadTPCOccMap(const o2::globaltracking::RecoContainer& recoDa } } +o2::track::TrackParCov TrackMCStudy::getCorrectedTPCTrack(const o2::tpc::TrackTPC& trTPC, float timeMUS) +{ + const float MUS2TPCBin = 1.f / (8 * o2::constants::lhc::LHCBunchSpacingMUS); + o2::track::TrackParCov tr{trTPC}; + if (trTPC.hasBothSidesClusters() || timeMUS < -1e6) { + return tr; + } + float tTB = timeMUS * MUS2TPCBin; + float dDrift = (tTB - trTPC.getTime0()) * mVDriftTB; + tr.setZ(trTPC.getZ() + (trTPC.hasASideClustersOnly() ? dDrift : -dDrift)); + return tr; +} + DataProcessorSpec getTrackMCStudySpec(GTrackID::mask_t srcTracks, GTrackID::mask_t srcClusters, const o2::tpc::CorrectionMapsLoaderGloOpts& sclOpts, bool checkSV) { std::vector outputs; diff --git a/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h b/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h index b933363bb352d..1df6e12b32a59 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h @@ -15,6 +15,8 @@ #ifndef O2_S_VERTEXER_H #define O2_S_VERTEXER_H +#define _DBGMC_SVERTEXER_ + #include "gsl/span" #include "DataFormatsCalibration/MeanVertexObject.h" #include "DataFormatsGlobalTracking/RecoContainer.h" @@ -197,6 +199,23 @@ class SVertexer std::vector> mFitterCasc; std::vector> mFitter3body; +#ifdef _DBGMC_SVERTEXER_ + bool mWatchHit = false; + bool mWatchHitVtxOK = false; + std::string mWatchHitRep{}; + o2::MCCompLabel mWatchLb; + std::vector mWatchLblVec; + bool checkLbl(o2::MCCompLabel& l) const + { + for (auto lb : mWatchLblVec) { + if (lb == l) { + return true; + } + } + return false; + } +#endif + PIDResponse mPIDresponse; int mNThreads = 1; diff --git a/Detectors/Vertexing/src/SVertexer.cxx b/Detectors/Vertexing/src/SVertexer.cxx index 1d48bcceb0097..6378d7b9d8b4c 100644 --- a/Detectors/Vertexing/src/SVertexer.cxx +++ b/Detectors/Vertexing/src/SVertexer.cxx @@ -24,6 +24,7 @@ #include "ReconstructionDataFormats/StrangeTrack.h" #include "CommonConstants/GeomConstants.h" #include "DataFormatsITSMFT/TrkClusRef.h" +#include "SimulationDataFormat/MCEventLabel.h" #ifdef WITH_OPENMP #include @@ -259,6 +260,22 @@ void SVertexer::produceOutput(o2::framework::ProcessingContext& pc) //__________________________________________________________________ void SVertexer::init() { +#ifdef _DBGMC_SVERTEXER_ + // must contain vector of {trackID, evID} pairs to watch, in the format + // std::vector> vTEID = {{trID0, evID0},{trID1, evID1} ...}. + // It is expected trID0 and trID1 are the labels of MC prongs + // #include "sv_labels_to_watch.inc" +#include "/data0/pp/LHC24af_550916MC1/tf8/dec22.txt" + // + // or just add locally something like + // std::vector> vTEID = {{930, 63}}; + + for (auto& p : vTEID) { + mWatchLblVec.emplace_back(p.first, p.second, 0); + mWatchLblVec.emplace_back(p.first + 1, p.second, 0); + } +#endif + // } //__________________________________________________________________ @@ -472,11 +489,31 @@ void SVertexer::buildT2V(const o2::globaltracking::RecoContainer& recoData) // a int it = vtref.getFirstEntry(), itLim = it + vtref.getEntries(); for (; it < itLim; it++) { auto tvid = trackIndex[it]; +#ifdef _DBGMC_SVERTEXER_ + mWatchHit = false; + mWatchHitVtxOK = false; + mWatchHitRep = ""; + if (mUseMC && mWatchLblVec.size()) { + mWatchLb = recoData.getTrackMCLabel(tvid); + mWatchHit = checkLbl(mWatchLb); + if (mWatchHit) { + int vtMCID = recoData.getPrimaryVertexMCLabel(iv).getEventID(); + mWatchHitVtxOK = vtMCID == mWatchLb.getEventID(); + mWatchHitRep = fmt::format("watched label {} for {}, MCVtx={}", mWatchLb.asString(), tvid.asString(), mWatchHitVtxOK); + LOGP(info, "Checking {} | PV#={}", mWatchHitRep, iv); + } + } +#endif if (!recoData.isTrackSourceLoaded(tvid.getSource())) { continue; } if (tvid.getSource() == GIndex::TPC) { if (mSVParams->mExcludeTPCtracks) { +#ifdef _DBGMC_SVERTEXER_ + if (mWatchHit) { + LOGP(info, "Reject {} due to the ExcludeTPCtracks", mWatchHitRep); + } +#endif continue; } // unconstrained TPC tracks require special treatment: there is no point in checking DCA to mean vertex since it is not precise, @@ -490,10 +527,20 @@ void SVertexer::buildT2V(const o2::globaltracking::RecoContainer& recoData) // a auto tref = tmap.find(tvid); if (tref != tmap.end()) { mTracksPool[tref->second.second][tref->second.first].vBracket.setMax(iv); // this track was already processed with other vertex, account the latter +#ifdef _DBGMC_SVERTEXER_ + if (mWatchHit) { + LOGP(info, "Skipping already processed ambiguous {}", mWatchHitRep); + } +#endif continue; } // was it already rejected? if (rejmap.find(tvid) != rejmap.end()) { +#ifdef _DBGMC_SVERTEXER_ + if (mWatchHit) { + LOGP(info, "Skipping already rejected ambiguous {}", mWatchHitRep); + } +#endif continue; } } @@ -541,6 +588,11 @@ void SVertexer::buildT2V(const o2::globaltracking::RecoContainer& recoData) // a if (tvid.isAmbiguous()) { rejmap[tvid] = true; } +#ifdef _DBGMC_SVERTEXER_ + if (mWatchHit) { + LOGP(info, "Reject {} in acceptTrack", mWatchHitRep); + } +#endif continue; } @@ -557,6 +609,11 @@ void SVertexer::buildT2V(const o2::globaltracking::RecoContainer& recoData) // a if (tvid.isAmbiguous()) { // track attached to >1 vertex, remember that it was already processed tmap[tvid] = {mTracksPool[posneg].size() - 1, posneg}; } +#ifdef _DBGMC_SVERTEXER_ + if (mWatchHit) { + LOGP(info, "Pooled {}", mWatchHitRep); + } +#endif } } // register 1st track of each charge for each vertex @@ -582,10 +639,42 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, auto& fitterV0 = mFitterV0[ithread]; // Fast rough cuts on pairs before feeding to DCAFitter, tracks are not in the same Frame or at same X bool isTPConly = (seedP.gid.getSource() == GIndex::TPC || seedN.gid.getSource() == GIndex::TPC); +#ifdef _DBGMC_SVERTEXER_ + bool watchHit = false, watchHitVtxOK = false; + o2::MCCompLabel watchLb0, watchLb1; + std::string watchPairRep{}; + if (mUseMC && mWatchLblVec.size()) { + watchLb0 = mRecoCont->getTrackMCLabel(seedP.gid); + watchLb1 = mRecoCont->getTrackMCLabel(seedN.gid); + watchHit = checkLbl(watchLb0) && checkLbl(watchLb1) && (watchLb0.getEventID() == watchLb1.getEventID() && std::abs(watchLb0.getTrackID() - watchLb1.getTrackID()) == 1); + if (watchHit) { + auto vlist = seedP.vBracket.getOverlap(seedN.vBracket); // indices of vertices shared by both seeds + std::string vlistStr{}; + for (int iv = vlist.getMin(); iv <= vlist.getMax(); iv++) { + int vtMCID = mRecoCont->getPrimaryVertexMCLabel(iv).getEventID(); + watchHitVtxOK = vtMCID == watchLb0.getEventID(); + vlistStr += fmt::format(" {}", iv); + if (watchHitVtxOK) { + vlistStr += "+"; + // break; + } + } + if (watchHitVtxOK) { + watchPairRep = fmt::format("watched labels {} of {} and {} of {}, MCVtx={} in {}", watchLb0.asString(), seedP.gid.asString(), watchLb1.asString(), seedN.gid.asString(), watchHitVtxOK, vlistStr); + LOGP(info, "Checking {}", watchPairRep); + } + } + } +#endif if (mSVParams->mTPCTrackPhotonTune && isTPConly) { // Check if Tgl is close enough if (std::abs(seedP.getTgl() - seedN.getTgl()) > mSVParams->maxV0TglAbsDiff) { LOG(debug) << "RejTgl"; +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {}, RegTgl:{}", watchPairRep, seedP.getTgl() - seedN.getTgl()); + } +#endif return false; } // Check in transverse plane @@ -599,8 +688,15 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, trkPosCircle.yC - trkEleCircle.yC); float r2r = trkPosCircle.rC + trkEleCircle.rC; float dcr = c2c - r2r; + if (std::abs(dcr) > mSVParams->mTPCTrackD2R) { LOG(debug) << "RejD2R " << c2c << " " << r2r << " " << dcr; +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} on dCR = {} (r2r={} c2c={}) (xyr: [{} {} {}] : [{} {} {}])", + watchPairRep, dcr, r2r, c2c, trkPosCircle.xC, trkPosCircle.yC, trkPosCircle.rC, trkEleCircle.xC, trkEleCircle.yC, trkEleCircle.rC); + } +#endif return false; } // Will the conversion point look somewhat reasonable @@ -609,6 +705,12 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, float dR = std::hypot(r2_r * trkPosCircle.xC + r1_r * trkEleCircle.xC, r2_r * trkPosCircle.yC + r1_r * trkEleCircle.yC); if (dR > mSVParams->mTPCTrackDR) { LOG(debug) << "RejDR" << dR; +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} on dR = {} (xyr: [{} {} {}] : [{} {} {}])", + watchPairRep, dR, trkPosCircle.xC, trkPosCircle.yC, trkPosCircle.rC, trkEleCircle.xC, trkEleCircle.yC, trkEleCircle.rC); + } +#endif return false; } @@ -631,6 +733,11 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, if (nCand == 0) { // discard this pair LOG(debug) << "RejDCAFitter"; +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {}: no candidates", watchPairRep); + } +#endif return false; } const auto& v0XYZ = fitterV0.getPCACandidate(); @@ -639,17 +746,32 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, float dxv0 = v0XYZ[0] - mMeanVertex.getX(), dyv0 = v0XYZ[1] - mMeanVertex.getY(), r2v0 = dxv0 * dxv0 + dyv0 * dyv0; if (r2v0 < mMinR2ToMeanVertex) { LOG(debug) << "RejMinR2ToMeanVertex"; +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} on RejMinR2ToMeanVertex check: r2v0={}", watchPairRep, r2v0); + } +#endif return false; } float rv0 = std::sqrt(r2v0), drv0P = rv0 - seedP.minR, drv0N = rv0 - seedN.minR; if (drv0P > mSVParams->causalityRTolerance || drv0P < -mSVParams->maxV0ToProngsRDiff || drv0N > mSVParams->causalityRTolerance || drv0N < -mSVParams->maxV0ToProngsRDiff) { LOG(debug) << "RejCausality " << drv0P << " " << drv0N; +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} on causality check: dr:{} {}", watchPairRep, drv0P, drv0N); + } +#endif return false; } const int cand = 0; if (!fitterV0.isPropagateTracksToVertexDone(cand) && !fitterV0.propagateTracksToVertex(cand)) { LOG(debug) << "RejProp failed"; +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} on propagateTracksToVertex", watchPairRep); + } +#endif return false; } const auto& trPProp = fitterV0.getTrack(0, cand); @@ -665,10 +787,20 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, float pt2V0 = pV0[0] * pV0[0] + pV0[1] * pV0[1], prodXYv0 = dxv0 * pV0[0] + dyv0 * pV0[1], tDCAXY = prodXYv0 / pt2V0; if (pt2V0 < mMinPt2V0) { // pt cut LOG(debug) << "RejPt2 " << pt2V0; +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} on MinPt2V0 check: {}", watchPairRep, pt2V0); + } +#endif return false; } if (pV0[2] * pV0[2] / pt2V0 > mMaxTgl2V0) { // tgLambda cut LOG(debug) << "RejTgL " << pV0[2] * pV0[2] / pt2V0; +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} on MaxTgl2V0 check: {}", watchPairRep, pV0[2] * pV0[2] / pt2V0); + } +#endif return false; } float p2V0 = pt2V0 + pV0[2] * pV0[2], ptV0 = std::sqrt(pt2V0); @@ -722,6 +854,11 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, if (!goodHyp && mSVParams->checkV0Hypothesis) { LOG(debug) << "RejHypo"; if (!checkFor3BodyDecays && !checkForCascade) { +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} as a bad hypothesis for >2 prong decays", watchPairRep); + } +#endif return false; } else { rejectAfter3BodyCheck = true; @@ -735,6 +872,11 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, if (dca2 > mMaxDCAXY2ToMeanVertexV0Casc || cosPAXY < mSVParams->minCosPAXYMeanVertexCascV0) { LOG(debug) << "Rej for cascade DCAXY2: " << dca2 << " << cosPAXY: " << cosPAXY; if (!checkFor3BodyDecays) { +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} as a bad cascade", watchPairRep); + } +#endif return false; } else { rejectAfter3BodyCheck = true; @@ -757,9 +899,22 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, if (mSVParams->mTPCTrackPhotonTune && isTPConly) { // Check for looser cut for tpc-only photons only if (dca2 > mSVParams->mTPCTrackMaxDCAXY2ToMeanVertex) { +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + auto tv0 = fitterV0.createParentTrackPar(); + LOGP(info, "Reject {} in TPCTrackMaxDCAXY2ToMeanVertex check: dca2={} : V0 XYZ:[{:.2f},{:.2f},{:.2f}], P:[{:.2f},{:.2f},{:.2f}] | {}", + watchPairRep, dca2, v0XYZ[0], v0XYZ[1], v0XYZ[2], pV0[0], pV0[1], pV0[2], tv0.asString()); + } +#endif return false; } } else { +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + auto tv0 = fitterV0.createParentTrackPar(); + LOGP(info, "Reject {} in mMaxDCAXY2ToMeanVertex ({}) and CosPAXYMeanVertex ({}) checks | {}", watchPairRep, dca2, cosPAXY, tv0.asString()); + } +#endif return false; } } @@ -793,6 +948,11 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, bestCosPA = cosPA; } if (!candFound) { +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} as no candidate PV is found", watchPairRep); + } +#endif return false; } if (bestCosPA < mSVParams->minCosPACascV0) { @@ -809,6 +969,11 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, n3bodyDecays += check3bodyDecays(v0Idxnew, v0new, rv0, pV0, p2V0, iP, POS, vlist, ithread); } if (rejectAfter3BodyCheck) { +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Reject {} after 3-body check", watchPairRep); + } +#endif return false; } @@ -851,7 +1016,11 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, mStrTracker->processV0(iv, v0new, v0Idxnew, ithread); } } - +#ifdef _DBGMC_SVERTEXER_ + if (watchHit && watchHitVtxOK) { + LOGP(info, "Watched {}, V0 status: {}", watchPairRep, mV0sIdxTmp[ithread].size() - nV0Ini != 0 ? "ACC" : "REJ"); + } +#endif return mV0sIdxTmp[ithread].size() - nV0Ini != 0; } @@ -1280,6 +1449,11 @@ void SVertexer::setNThreads(int n) bool SVertexer::processTPCTrack(const o2::tpc::TrackTPC& trTPC, GIndex gid, int vtxid) { if (mSVParams->mTPCTrackMaxX > 0. && trTPC.getX() > mSVParams->mTPCTrackMaxX) { +#ifdef _DBGMC_SVERTEXER_ + if (mWatchHit) { + LOGP(info, "Reject {} in processTPCTrack: TPCTrackMaxX check, {}", mWatchHitRep, trTPC.getX()); + } +#endif return true; } // if TPC trackis unconstrained, try to create in the tracks pool a clone constrained to vtxid vertex time. @@ -1306,6 +1480,11 @@ bool SVertexer::processTPCTrack(const o2::tpc::TrackTPC& trTPC, GIndex gid, int auto err = correctTPCTrack(trLoc, trTPC, twe.getTimeStamp(), twe.getTimeStampError()); if (err < 0) { mTracksPool[posneg].pop_back(); // discard +#ifdef _DBGMC_SVERTEXER_ + if (mWatchHit) { + LOGP(info, "Reject {} in processTPCTrack: correctTPCTrack fails", mWatchHitRep); + } +#endif return true; } @@ -1319,15 +1498,29 @@ bool SVertexer::processTPCTrack(const o2::tpc::TrackTPC& trTPC, GIndex gid, int o2::math_utils::CircleXYf_t trkCircle; trLoc.getCircleParams(mBz, trkCircle, sna, csa); float cR = std::hypot(trkCircle.xC, trkCircle.yC); - float drd2 = std::sqrt(cR * cR - trkCircle.rC * trkCircle.rC); - bool dRD2 = drd2 > mSVParams->mTPCTrackXY2Radius; + bool dRD2 = cR - trkCircle.rC > mSVParams->mTPCTrackXY2Radius; if (dCls || dDPV || dRD2) { mTracksPool[posneg].pop_back(); +#ifdef _DBGMC_SVERTEXER_ + if (mWatchHit) { + LOGP(info, "Reject {} in processTPCTrack: CheckClus:{} ({}) CheckDPv:{} ({}) CheckRD2:{} ({}) | {}", + mWatchHitRep, + dCls, + trTPC.getNClusters(), + dDPV, + trLoc.getX() * trLoc.getTgl() - trLoc.getZ() + vtx.getZ(), + dRD2, cR - trkCircle.rC > mSVParams->mTPCTrackXY2Radius, ((o2::track::TrackPar&)trLoc).asString()); + } +#endif return true; } } - +#ifdef _DBGMC_SVERTEXER_ + if (mWatchHit) { + LOGP(info, "Pooled TPC-only {}", mWatchHitRep); + } +#endif return true; }