Skip to content

Commit

Permalink
Add debug output
Browse files Browse the repository at this point in the history
  • Loading branch information
wang-yuanzhe committed Nov 1, 2024
1 parent 5933a9d commit c0d5b96
Show file tree
Hide file tree
Showing 2 changed files with 111 additions and 17 deletions.
18 changes: 16 additions & 2 deletions Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 <TTree.h>
#include <TFile.h>

namespace o2
{
namespace tpc
Expand Down Expand Up @@ -153,9 +161,11 @@ class SVertexer
private:
template <class TVI, class TCI, class T3I, class TR>
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<float, 3> pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread);
int check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, std::array<float, 3> pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread);
//int check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, std::array<float, 3> 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<float, 3> 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();
Expand Down Expand Up @@ -224,6 +234,10 @@ class SVertexer
bool mEnableCascades = true;
bool mEnable3BodyDecays = false;
bool mUseMC = false;

//For Debug
std::unique_ptr<o2::utils::TreeStreamRedirector> svDebug;
bool debugFor3BodyDecays = true;
};

} // namespace vertexing
Expand Down
110 changes: 95 additions & 15 deletions Detectors/Vertexing/src/SVertexer.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand Down Expand Up @@ -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<o2::utils::TreeStreamRedirector>("svtxDebug.root", "recreate");
}

//__________________________________________________________________
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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) {
Expand All @@ -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)) {
Expand All @@ -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
Expand Down Expand Up @@ -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) &&
Expand Down Expand Up @@ -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;
}
}

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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<float, 3> pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread)
// int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, std::array<float, 3> 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<float, 3> pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread)
{
// check last added V0 for belonging to cascade
auto& fitter3body = mFitter3body[ithread];
Expand Down Expand Up @@ -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);
Expand All @@ -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) {
Expand Down Expand Up @@ -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];
Expand All @@ -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);
Expand Down

0 comments on commit c0d5b96

Please sign in to comment.