Skip to content

svertexer: adaptive mass windows #12136

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
shahor02 merged 9 commits into from
Nov 2, 2023
Merged

svertexer: adaptive mass windows #12136

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
f13b248
svertexer: adaptive mass windows
Oct 24, 2023
93573a8
One more step: Lucia's changes
Oct 25, 2023
9e3673d
Tighter windows when doing decay chains
Oct 25, 2023
71ee9e2
Please consider the following formatting changes
alibuild Oct 25, 2023
3096aaa
Merge pull request #15 from alibuild/alibot-cleanup-12136
ddobrigk Oct 25, 2023
5318239
Merge branch 'AliceO2Group:dev' into svertmass
ddobrigk Oct 30, 2023
c0f98e3
Latest bugfixes
Oct 31, 2023
71db1b9
Please consider the following formatting changes
alibuild Oct 31, 2023
a1ddf3e
Merge pull request #16 from alibuild/alibot-cleanup-12136
ddobrigk Nov 1, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 46 additions & 10 deletions Detectors/Vertexing/include/DetectorsVertexing/SVertexHypothesis.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,15 +30,20 @@ class SVertexHypothesis

public:
using PID = o2::track::PID;
enum PIDParams { SigmaM, // sigma of mass res at 0 pt
NSigmaM, // number of sigmas of mass res
MarginM, // additive safety margin in mass cut
CPt }; // pT dependence of mass resolution parameterized as mSigma*(1+mC1*pt);
enum PIDParams { SigmaM, // sigma of mass res at 0 pt
NSigmaM, // number of sigmas of mass res
MarginM, // additive safety margin in mass cut
NSigmaTightM, // number of sigmas of mass res when doing tight cut around mass (V0s used in cascades)
MarginTightM, // additive safety margin in mass cut when doing tight cut around mass (V0s used in cascades)
CPt, // pT dependence of mass resolution parameterized as mSigma*(1+mC1*pt);
CPt1,
CPt2,
CPt3 }; // pT dependence of mass resolution of Cascade parameterized as CPt+CPt1*pt +CPt2*TMath::Exp(-CPt3*pt);

static constexpr int NPIDParams = 4;
static constexpr int NPIDParams = 9;

void set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float cpt, float bz = 0.f);
void set(PID v0, PID ppos, PID pneg, const float pars[NPIDParams], float bz = 0.f);
void set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float nSigTight, float marginTight, float cpt, float cpt1, float cpt2, float cpt3, float bz = 0.f, float maxSigma = 0.01);
void set(PID v0, PID ppos, PID pneg, const float pars[NPIDParams], float bz = 0.f, float maxSigma = 0.01);

float getMassV0Hyp() const { return PID::getMass(mPIDV0); }
float getMassPosProng() const { return PID::getMass(mPIDPosProng); }
Expand All @@ -57,14 +62,44 @@ class SVertexHypothesis
{ // check if given mass and pt is matching to hypothesis
return check(calcMass(p2Pos, p2Neg, p2V0), ptV0);
}

bool check(float mass, float pt) const
{ // check if given mass and pt is matching to hypothesis
return std::abs(mass - getMassV0Hyp()) < getMargin(pt);
}

bool checkTight(float p2Pos, float p2Neg, float p2V0, float ptV0) const
{ // check if given mass and pt is matching to hypothesis
return checkTight(calcMass(p2Pos, p2Neg, p2V0), ptV0);
}
bool checkTight(float mass, float pt) const
{ // check if given mass and pt is matching to hypothesis
return std::abs(mass - getMassV0Hyp()) < getMarginTight(pt);
}

float getSigmaV0Cascade(float pt) const { return mPars[CPt] + mPars[CPt1] * pt + mPars[CPt2] * std::exp(-mPars[CPt3] * pt); }
float getSigma(float pt) const { return mPars[SigmaM] * (1.f + mPars[CPt] * pt); }
float getMargin(float pt) const { return mPars[NSigmaM] * getSigma(pt) + mPars[MarginM]; }
float getMargin(float pt, bool tight = false) const
{
int idxNsigma = NSigmaM;
int idxMargin = MarginM;
if (tight) { // move to indices for tight variables in case asked to do so (tighter peak cuts for decay chains)
idxNsigma = NSigmaTightM;
idxMargin = MarginTightM;
}
if (mPIDV0 == PID::XiMinus || mPIDV0 == PID::OmegaMinus) { // case for cascades, antiparticles included
float sigmaV0Cascade = getSigmaV0Cascade(pt);
if (sigmaV0Cascade > maxSigma) { // insuring that at low pt one gets reasonable width as the parametrisation function may explode to unphysical values
return mPars[idxNsigma] * maxSigma + mPars[idxMargin];
} else {
return mPars[idxNsigma] * sigmaV0Cascade + mPars[idxMargin];
}
} else if (mPIDV0 == PID::K0 || mPIDV0 == PID::Lambda) { // case for V0s, AntiLambda is included in PID::Lambda
return mPars[idxNsigma] * getSigmaV0Cascade(pt) + mPars[idxMargin];
} else {
return mPars[idxNsigma] * getSigma(pt) + mPars[idxMargin]; // case for HyperTriton and Hyperhydrog4
}
}
float getMarginTight(float pt) const { return getMargin(pt, true); }

private:
float getMass2PosProng() const { return PID::getMass2(mPIDPosProng); }
Expand All @@ -76,8 +111,9 @@ class SVertexHypothesis

public: // to be deleted
std::array<float, NPIDParams> mPars{};
float maxSigma;

ClassDefNV(SVertexHypothesis, 1);
ClassDefNV(SVertexHypothesis, 2);
};

class SVertex3Hypothesis
Expand Down
15 changes: 8 additions & 7 deletions Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,16 +96,17 @@ struct SVertexerParams : public o2::conf::ConfigurableParamHelper<SVertexerParam

// cuts on different V0 PID params
bool checkV0Hypothesis = true;
float pidCutsPhoton[SVertexHypothesis::NPIDParams] = {0.001, 20, 0.60, 0.0}; // Photon
float pidCutsK0[SVertexHypothesis::NPIDParams] = {0.003, 20, 0.07, 0.5}; // K0
float pidCutsLambda[SVertexHypothesis::NPIDParams] = {0.001, 20, 0.07, 0.5}; // Lambda
float pidCutsHTriton[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5}; // HyperTriton
float pidCutsHhydrog4[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5}; // Hyperhydrog4 - Need to update
float pidCutsPhoton[SVertexHypothesis::NPIDParams] = {0.001, 20, 0.60, 20, 0.60, 0.0, 0.0, 0.0, 0.0}; // Photon
Copy link
Collaborator

@shahor02 shahor02 Oct 25, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is it normal to have 600 MeV margin on the photon (old setting)?

Copy link
Contributor Author

@ddobrigk ddobrigk Oct 31, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a very good question but I did not touch it. I think the problem is that we currently apply logical OR to 7 types of mass hypotheses, so if one of them has very open cuts, the whole concept of windows will not be particularly helpful. I guess this is a question more to the PWG-EM conveners... Do we know where the 600 MeV/c2 comes from in the first place?

Copy link
Collaborator

@shahor02 shahor02 Oct 31, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No but @f3sch is now studying conversions, perhaps he can extract realistic reconstructed mass range for photons

float pidCutsK0[SVertexHypothesis::NPIDParams] = {0., 20, 0., 5.0, 0.0, 2.84798e-03, 9.84206e-04, 3.31951e-03, 2.39438}; // K0
float pidCutsLambda[SVertexHypothesis::NPIDParams] = {0., 20, 0., 5.0, 0.0, 1.09004e-03, 2.62291e-04, 8.93179e-03, 2.83121}; // Lambda
float pidCutsHTriton[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 14, 0.0, 0.5, 0.0, 0.0, 0.0}; // HyperTriton
float pidCutsHhydrog4[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 14, 0.0, 0.5, 0.0, 0.0, 0.0}; // Hyperhydrog4 - Need to update
//
// cuts on different Cascade PID params
bool checkCascadeHypothesis = true;
float pidCutsXiMinus[SVertexHypothesis::NPIDParams] = {0.001, 20, 0.07, 0.5}; // XiMinus
float pidCutsOmegaMinus[SVertexHypothesis::NPIDParams] = {0.001, 20, 0.07, 0.5}; // OmegaMinus
float pidCutsXiMinus[SVertexHypothesis::NPIDParams] = {0.0, 10, 0.0, 4.0, 0.0, 1.56315e-03, 2.23279e-04, 2.75136e-02, 3.309}; // XiMinus
float pidCutsOmegaMinus[SVertexHypothesis::NPIDParams] = {0.0, 10, 0.0, 4.0, 0.0, 1.43572e-03, 6.94416e-04, 2.13534e+05, 1.48889e+01}; // OmegaMinus
float maximalCascadeWidth = 0.006;
//
// cuts on different 3 body PID params
bool check3bodyHypothesis = true;
Expand Down
16 changes: 12 additions & 4 deletions Detectors/Vertexing/src/SVertexHypothesis.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,21 +17,29 @@

using namespace o2::vertexing;

void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float cpt, float bz)
void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float nSigTight, float marginTight, float cpt, float cpt1, float cpt2, float cpt3, float bz, float maxSigmaInput)
{
mPIDV0 = v0;
mPIDPosProng = ppos;
mPIDNegProng = pneg;
mPars[SigmaM] = sig;
mPars[NSigmaM] = nSig;
mPars[MarginM] = margin;
mPars[NSigmaTightM] = nSigTight;
mPars[MarginTightM] = marginTight;
mPars[CPt] = cpt;
mPars[CPt1] = cpt1;
mPars[CPt2] = cpt2;
mPars[CPt3] = cpt3;
maxSigma = maxSigmaInput;
float absBz{std::abs(bz)};
mPars[CPt] = absBz > 1e-3 ? cpt * 5.0066791 / absBz : 0.; // assume that pT dependent sigma is linear with B
if (cpt3 < 1)
mPars[CPt] = absBz > 1e-3 ? cpt * 5.0066791 / absBz : 0.; // assume that pT dependent sigma is linear with B; case for HyperTriton and Hyperhydrog4
}

void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, const float pars[NPIDParams], float bz)
void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, const float pars[NPIDParams], float bz, float maxSigmaInput)
{
set(v0, ppos, pneg, pars[SigmaM], pars[NSigmaM], pars[MarginM], pars[CPt], bz);
set(v0, ppos, pneg, pars[SigmaM], pars[NSigmaM], pars[MarginM], pars[NSigmaTightM], pars[MarginTightM], pars[CPt], pars[CPt1], pars[CPt2], pars[CPt3], bz, maxSigmaInput);
}

void SVertex3Hypothesis::set(PID v0, PID ppos, PID pneg, PID pbach, float sig, float nSig, float margin, float cpt, float bz)
Expand Down
16 changes: 11 additions & 5 deletions Detectors/Vertexing/src/SVertexer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -283,8 +283,8 @@ void SVertexer::updateTimeDependentParams()
mV0Hyps[HypV0::AntiHyperTriton].set(PID::HyperTriton, PID::Pion, PID::Helium3, mSVParams->pidCutsHTriton, bz);
mV0Hyps[HypV0::Hyperhydrog4].set(PID::Hyperhydrog4, PID::Alpha, PID::Pion, mSVParams->pidCutsHhydrog4, bz);
mV0Hyps[HypV0::AntiHyperhydrog4].set(PID::Hyperhydrog4, PID::Pion, PID::Alpha, mSVParams->pidCutsHhydrog4, bz);
mCascHyps[HypCascade::XiMinus].set(PID::XiMinus, PID::Lambda, PID::Pion, mSVParams->pidCutsXiMinus, bz);
mCascHyps[HypCascade::OmegaMinus].set(PID::OmegaMinus, PID::Lambda, PID::Kaon, mSVParams->pidCutsOmegaMinus, bz);
mCascHyps[HypCascade::XiMinus].set(PID::XiMinus, PID::Lambda, PID::Pion, mSVParams->pidCutsXiMinus, bz, mSVParams->maximalCascadeWidth);
mCascHyps[HypCascade::OmegaMinus].set(PID::OmegaMinus, PID::Lambda, PID::Kaon, mSVParams->pidCutsOmegaMinus, bz, mSVParams->maximalCascadeWidth);

m3bodyHyps[Hyp3body::H3L3body].set(PID::HyperTriton, PID::Proton, PID::Pion, PID::Deuteron, mSVParams->pidCutsH3L3body, bz);
m3bodyHyps[Hyp3body::AntiH3L3body].set(PID::HyperTriton, PID::Pion, PID::Proton, PID::Deuteron, mSVParams->pidCutsH3L3body, bz);
Expand Down Expand Up @@ -567,6 +567,12 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
goodHyp = hypCheckStatus[ipid] = true;
}
}
// check tight lambda mass only
bool goodLamForCascade = false, goodALamForCascade = false;
if (mV0Hyps[Lambda].checkTight(p2Pos, p2Neg, p2V0, ptV0))
goodLamForCascade = true;
if (mV0Hyps[AntiLambda].checkTight(p2Pos, p2Neg, p2V0, ptV0))
goodALamForCascade = true;

// apply mass selections for 3-body decay
bool good3bodyV0Hyp = false;
Expand All @@ -581,7 +587,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
// we want to reconstruct the 3 body decay of hypernuclei starting from the V0 of a proton and a pion (e.g. H3L->d + (p + pi-), or He4L->He3 + (p + pi-)))
bool checkFor3BodyDecays = mEnable3BodyDecays && (!mSVParams->checkV0Hypothesis || good3bodyV0Hyp) && (pt2V0 > 0.5);
bool rejectAfter3BodyCheck = false; // To reject v0s which can be 3-body decay candidates but not cascade or v0
bool checkForCascade = mEnableCascades && r2v0 < mMaxR2ToMeanVertexCascV0 && (!mSVParams->checkV0Hypothesis || (hypCheckStatus[HypV0::Lambda] || hypCheckStatus[HypV0::AntiLambda]));
bool checkForCascade = mEnableCascades && r2v0 < mMaxR2ToMeanVertexCascV0 && (!mSVParams->checkV0Hypothesis || (goodLamForCascade || goodALamForCascade));
bool rejectIfNotCascade = false;

if (!goodHyp && mSVParams->checkV0Hypothesis) {
Expand Down Expand Up @@ -673,10 +679,10 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
// check cascades
int nCascIni = mCascadesIdxTmp[ithread].size(), nV0Used = 0; // number of times this particular v0 (with assigned PV) was used (not counting using its clones with other PV)
if (checkForCascade) {
if (hypCheckStatus[HypV0::Lambda] || !mSVParams->checkCascadeHypothesis) {
if (goodLamForCascade || !mSVParams->checkCascadeHypothesis) {
nV0Used += checkCascades(v0Idxnew, v0new, rv0, pV0, p2V0, iN, NEG, vlist, ithread);
}
if (hypCheckStatus[HypV0::AntiLambda] || !mSVParams->checkCascadeHypothesis) {
if (goodALamForCascade || !mSVParams->checkCascadeHypothesis) {
nV0Used += checkCascades(v0Idxnew, v0new, rv0, pV0, p2V0, iP, POS, vlist, ithread);
}
}
Expand Down