diff --git a/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h b/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h index b933363bb352d..904ea6441d0b4 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h @@ -38,6 +38,14 @@ #include "TPCFastTransform.h" #include "DataFormatsTPC/PIDResponse.h" +//add for 3body check +#include "SimulationDataFormat/MCEventLabel.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "CommonUtils/TreeStream.h" +#include "CommonUtils/TreeStreamRedirector.h" +#include +#include + namespace o2 { namespace tpc @@ -153,9 +161,11 @@ class SVertexer private: template void extractPVReferences(const TVI& v0s, TR& vtx2V0Refs, const TCI& cascades, TR& vtx2CascRefs, const T3I& vtxs3, TR& vtx2body3Refs); - bool checkV0(const TrackCand& seed0, const TrackCand& seed1, int iP, int iN, int ithread); + //bool checkV0(const TrackCand& seed0, const TrackCand& seed1, int iP, int iN, int ithread); + bool checkV0(const o2::globaltracking::RecoContainer& recoData, const TrackCand& seed0, const TrackCand& seed1, int iP, int iN, int ithread); // for 3body debug int checkCascades(const V0Index& v0Idx, const V0& v0, float rv0, std::array pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread); - int check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, std::array pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread); + //int check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, std::array pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread); + int check3bodyDecays(const o2::globaltracking::RecoContainer& recoData, const V0Index& v0Idx, const V0& v0, float rv0, std::array pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread); // for 3body debug void setupThreads(); void buildT2V(const o2::globaltracking::RecoContainer& recoTracks); void updateTimeDependentParams(); @@ -224,6 +234,10 @@ class SVertexer bool mEnableCascades = true; bool mEnable3BodyDecays = false; bool mUseMC = false; + + //For Debug + std::unique_ptr svDebug; + bool debugFor3BodyDecays = true; }; } // namespace vertexing diff --git a/Detectors/Vertexing/src/SVertexer.cxx b/Detectors/Vertexing/src/SVertexer.cxx index 1d48bcceb0097..50e20d4a3476b 100644 --- a/Detectors/Vertexing/src/SVertexer.cxx +++ b/Detectors/Vertexing/src/SVertexer.cxx @@ -76,10 +76,27 @@ void SVertexer::process(const o2::globaltracking::RecoContainer& recoData, o2::f #else int iThread = 0; #endif - checkV0(seedP, seedN, itp, itn, iThread); + // checkV0(seedP, seedN, itp, itn, iThread); + checkV0(recoData, seedP, seedN, itp, itn, iThread); // for 3body debug } } + // Debug Check for the TrackPool + for (int itp = 0; itp < ntrP; itp++) { + auto& seedP = mTracksPool[POS][itp]; + auto mclabelP = recoData.getTrackMCLabel(seedP.gid); + (*svDebug) << "PosTrackPool" + << "PosTrackMCGID=" << mclabelP + << "\n"; + } + for (int itn = 0; itn < ntrN; itn++) { + auto& seedN = mTracksPool[NEG][itn]; + auto mclabelN = recoData.getTrackMCLabel(seedN.gid); + (*svDebug) << "NegTrackPool" + << "NegTrackMCGID=" << mclabelN + << "\n"; + } + produceOutput(pc); } @@ -254,11 +271,15 @@ void SVertexer::produceOutput(o2::framework::ProcessingContext& pc) } extractPVReferences(v0sIdx, v0Refs, cascsIdx, cascRefs, body3Idx, vtx3bodyRefs); + LOG(debug) << "DONE : " << mV0sTmp[0].size() << " " << mCascadesTmp[0].size() << " " << m3bodyTmp[0].size(); + + svDebug->Close(); } //__________________________________________________________________ void SVertexer::init() { + svDebug = std::make_unique("svtxDebug.root", "recreate"); } //__________________________________________________________________ @@ -577,8 +598,10 @@ void SVertexer::buildT2V(const o2::globaltracking::RecoContainer& recoData) // a } //__________________________________________________________________ -bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, int iN, int ithread) +// bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, int iN, int ithread) +bool SVertexer::checkV0(const o2::globaltracking::RecoContainer& recoData, const TrackCand& seedP, const TrackCand& seedN, int iP, int iN, int ithread) { + 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); @@ -619,6 +642,13 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, fitterV0.setCollinear(true); } + auto mclabelP = recoData.getTrackMCLabel(seedP.gid); + auto mclabelN = recoData.getTrackMCLabel(seedN.gid); + (*svDebug) << "V0Tree" + //<< "PTrackGID.=" << seedP.gid << "NTrackGID.=" << seedN.gid + << "PTrackMCGID=" << mclabelP << "NTrackMCGID=" << mclabelN + << "\n"; + // feed DCAFitter int nCand = fitterV0.process(seedP, seedN); if (mSVParams->mTPCTrackPhotonTune && isTPConly) { @@ -638,14 +668,15 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, // check closeness to the beam-line float dxv0 = v0XYZ[0] - mMeanVertex.getX(), dyv0 = v0XYZ[1] - mMeanVertex.getY(), r2v0 = dxv0 * dxv0 + dyv0 * dyv0; if (r2v0 < mMinR2ToMeanVertex) { - LOG(debug) << "RejMinR2ToMeanVertex"; - return false; + if (!debugFor3BodyDecays) + 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; - return false; + if (!debugFor3BodyDecays) + return false; } const int cand = 0; if (!fitterV0.isPropagateTracksToVertexDone(cand) && !fitterV0.propagateTracksToVertex(cand)) { @@ -665,11 +696,13 @@ 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; - return false; + if (!debugFor3BodyDecays) + return false; } if (pV0[2] * pV0[2] / pt2V0 > mMaxTgl2V0) { // tgLambda cut LOG(debug) << "RejTgL " << pV0[2] * pV0[2] / pt2V0; - return false; + if (!debugFor3BodyDecays) + return false; } float p2V0 = pt2V0 + pV0[2] * pV0[2], ptV0 = std::sqrt(pt2V0); // apply mass selections @@ -711,6 +744,8 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, (!mSVParams->checkV0Hypothesis || good3bodyV0Hyp) && (pt2V0 > 0.5) && (!mSVParams->mSkipTPCOnly3Body || !isTPConly); + if (debugFor3BodyDecays) + checkFor3BodyDecays = true; // for debug bool rejectAfter3BodyCheck = false; // To reject v0s which can be 3-body decay candidates but not cascade or v0 bool checkForCascade = mEnableCascades && (!mSVParams->mSkipTPCOnlyCascade || !usesTPCOnly) && @@ -744,7 +779,8 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, if (checkFor3BodyDecays) { // use looser cuts for 3-body decay candidates if (dca2 > mMaxDCAXY2ToMeanVertex3bodyV0 || cosPAXY < mSVParams->minCosPAXYMeanVertex3bodyV0) { LOG(debug) << "Rej for 3 body decays DCAXY2: " << dca2 << " << cosPAXY: " << cosPAXY; - checkFor3BodyDecays = false; + if (!debugFor3BodyDecays) + checkFor3BodyDecays = false; } } @@ -802,11 +838,28 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, rejectIfNotCascade = true; } int nV0Ini = mV0sIdxTmp[ithread].size(); + // check 3 body decays if (checkFor3BodyDecays) { + + // 3body debug + (*svDebug) << "V0TreeAfterCut" + //<< "PTrackGID=" << seedP.gid << "NTrackGID=" << seedN.gid + << "PTrackMCGID=" << mclabelP << "NTrackMCGID=" << mclabelN + << "V0R=" << rv0 << "V0DrP=" << drv0P << "V0DrN=" << drv0N << "V0Pt=" << ptV0 << "V0TgLambda=" << pV0[2] * pV0[2] / pt2V0 + << "V0Dca=" << std::sqrt(dca2) << "V0CosXY=" << cosPAXY << "V0CosPA=" << bestCosPA + << "V0LambdaMass=" << mV0Hyps[2].calcMass(p2Pos, p2Neg, p2V0) << "V0AntiLambdaMass=" << mV0Hyps[3].calcMass(p2Pos, p2Neg, p2V0) + << "V0PTrackCharge=" << trPProp.getAbsCharge() << "V0NTrackCharge=" << trNProp.getAbsCharge() + << "V0PTrackPx=" << pP[0] << "V0PTrackPy=" << pP[1] << "V0PTrackPz=" << pP[2] + << "V0NTrackPx=" << pN[0] << "V0NTrackPy=" << pN[1] << "V0NTrackPz=" << pN[2] + << "V0H3LMass=" << mV0Hyps[4].calcMass(p2Pos, p2Neg, p2V0) << "V0AntiH3LMass=" << mV0Hyps[5].calcMass(p2Pos, p2Neg, p2V0) + << "\n"; + int n3bodyDecays = 0; - n3bodyDecays += check3bodyDecays(v0Idxnew, v0new, rv0, pV0, p2V0, iN, NEG, vlist, ithread); - n3bodyDecays += check3bodyDecays(v0Idxnew, v0new, rv0, pV0, p2V0, iP, POS, vlist, ithread); + // n3bodyDecays += check3bodyDecays(v0Idxnew, v0new, rv0, pV0, p2V0, iN, NEG, vlist, ithread); + // n3bodyDecays += check3bodyDecays(v0Idxnew, v0new, rv0, pV0, p2V0, iP, POS, vlist, ithread); + n3bodyDecays += check3bodyDecays(recoData, v0Idxnew, v0new, rv0, pV0, p2V0, iN, NEG, vlist, ithread); + n3bodyDecays += check3bodyDecays(recoData, v0Idxnew, v0new, rv0, pV0, p2V0, iP, POS, vlist, ithread); } if (rejectAfter3BodyCheck) { return false; @@ -1025,7 +1078,8 @@ int SVertexer::checkCascades(const V0Index& v0Idx, const V0& v0, float rv0, std: } //__________________________________________________________________ -int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, std::array pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread) +// int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, std::array pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread) +int SVertexer::check3bodyDecays(const o2::globaltracking::RecoContainer& recoData, const V0Index& v0Idx, const V0& v0, float rv0, std::array pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread) { // check last added V0 for belonging to cascade auto& fitter3body = mFitter3body[ithread]; @@ -1064,8 +1118,19 @@ int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, s decay3bodyVlist.setMax(v0Idx.getVertexID()); } + // for 3body debug + auto mclabel0 = recoData.getTrackMCLabel(v0Idx.getProngID(0)); + auto mclabel1 = recoData.getTrackMCLabel(v0Idx.getProngID(1)); + auto mclabel2 = recoData.getTrackMCLabel(bach.gid); + + (*svDebug) << "VtxTree" + //<< "Track0GID.=" << v0Idx.getProngID(0) << "Track1GID.=" << v0Idx.getProngID(1) << "Track2GID.=" << bach.gid + << "Track0MCGID=" << mclabel0 << "Track1MCGID=" << mclabel1 << "Track2MCGID=" << mclabel2 + << "\n"; + if (bach.getPt() < 0.6) { - continue; + if (!debugFor3BodyDecays) + continue; } int n3bodyVtx = fitter3body.process(v0.getProng(0), v0.getProng(1), bach); @@ -1078,7 +1143,8 @@ int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, s // make sure the 3 body vertex radius is close to that of the mean vertex float dxc = vertexXYZ[0] - mMeanVertex.getX(), dyc = vertexXYZ[1] - mMeanVertex.getY(), dzc = vertexXYZ[2] - mMeanVertex.getZ(), r2vertex = dxc * dxc + dyc * dyc; if (std::abs(rv0 - std::sqrt(r2vertex)) > mSVParams->maxRDiffV03body || r2vertex < mMinR2ToMeanVertex) { - continue; + if (!debugFor3BodyDecays) + continue; } float drvtxBach = std::sqrt(r2vertex) - bach.minR; if (drvtxBach > mSVParams->causalityRTolerance || drvtxBach < -mSVParams->maxV0ToProngsRDiff) { @@ -1145,7 +1211,8 @@ int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, s } } if (!goodHyp) { - continue; + if (!debugFor3BodyDecays || decay3bodyVtxID == -1) + continue; } const auto& decay3bodyPv = mPVertices[decay3bodyVtxID]; @@ -1154,13 +1221,26 @@ int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, s o2::dataformats::DCA dca; if (!trc.propagateToDCA(decay3bodyPv, fitter3body.getBz(), &dca, 5.) || std::abs(dca.getY()) > mSVParams->maxDCAXY3Body || std::abs(dca.getZ()) > mSVParams->maxDCAZ3Body) { - continue; + if (!debugFor3BodyDecays) + continue; } if (mSVParams->createFull3Bodies) { candidate3B.setCosPA(vtxCosPA); candidate3B.setDCA(fitter3body.getChi2AtPCACandidate()); m3bodyTmp[ithread].push_back(candidate3B); } + // for debug + float sqP0 = p0[0] * p0[0] + p0[1] * p0[1] + p0[2] * p0[2], sqP1 = p1[0] * p1[0] + p1[1] * p1[1] + p1[2] * p1[2], sqPBach = pbach[0] * pbach[0] + pbach[1] * pbach[1] + pbach[2] * pbach[2]; + float pt2Candidate = p3B[0] * p3B[0] + p3B[1] * p3B[1], p2Candidate = pt2Candidate + p3B[2] * p3B[2]; + float ptCandidate = std::sqrt(pt2Candidate); + (*svDebug) << "VtxTreeAfterCut" + //<< "Track0GID.=" << v0Idx.getProngID(0) << "Track1GID.=" << v0Idx.getProngID(1) << "Track2GID.=" << bach.gid + << "Track0MCGID=" << mclabel0 << "Track1MCGID=" << mclabel1 << "Track2MCGID=" << mclabel2 + << "VtxBachR=" << bach.minR << "VtxDrBach=" << drvtxBach << "VtxPt=" << ptCandidate << "VtxDiffR=" << std::abs(rv0 - std::sqrt(r2vertex)) << "VtxR=" << std::sqrt(r2vertex) + << "VtxTgLambda=" << p3B[2] * p3B[2] / pt2Candidate << "VtxCosPA=" << vtxCosPA << "VtxDcaY=" << dca.getY() << "VtxDcaZ=" << dca.getZ() + << "VtxH3LMass=" << m3bodyHyps[0].calcMass(sqP0, sqP1, sqPBach, ptCandidate) << "VtxAntiH3LMass=" << m3bodyHyps[1].calcMass(sqP0, sqP1, sqPBach, ptCandidate) + << "VtxMassMargin=" << m3bodyHyps[0].getMargin(ptCandidate) << "VtxMassSigma=" << m3bodyHyps[0].getSigma(ptCandidate) + << "\n"; m3bodyIdxTmp[ithread].emplace_back(decay3bodyVtxID, v0Idx.getProngID(0), v0Idx.getProngID(1), bach.gid); Decay3BodyIndex decay3bodyIdx(decay3bodyVtxID, v0Idx.getProngID(0), v0Idx.getProngID(1), bach.gid);