Skip to content

PWGHF: Add Dplus variables and histograms #4797

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
jgrosseo merged 1 commit into from
Dec 3, 2020
Merged

PWGHF: Add Dplus variables and histograms #4797

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
85 changes: 84 additions & 1 deletion Analysis/Core/include/Analysis/RecoDecay.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ class RecoDecay
}

/// Calculates scalar product of vectors.
/// \note Promotes numbers to double before squaring to avoid precision loss in float multiplication.
/// \note Promotes numbers to double to avoid precision loss in float multiplication.
/// \param N dimension
/// \param vec1,vec2 vectors
/// \return scalar product
Expand All @@ -115,6 +115,19 @@ class RecoDecay
return res;
}

/// FIXME: probably cross and dot products should be in some utility class
/// Calculates cross product of vectors in three dimensions.
/// \note Promotes numbers to double to avoid precision loss in float multiplication.
/// \param vec1,vec2 vectors
/// \return cross-product vector
template <typename T, typename U>
static array<double, 3> crossProd(const array<T, 3>& vec1, const array<U, 3>& vec2)
{
return array<double, 3>{((double)vec1[1] * (double)vec2[2]) - ((double)vec1[2] * (double)vec2[1]),
((double)vec1[2] * (double)vec2[0]) - ((double)vec1[0] * (double)vec2[2]),
((double)vec1[0] * (double)vec2[1]) - ((double)vec1[1] * (double)vec2[0])};
}

/// Calculates magnitude squared of a vector.
/// \param N dimension
/// \param vec vector
Expand Down Expand Up @@ -388,6 +401,76 @@ class RecoDecay
return std::sqrt(M2(args...));
}

// Calculation of topological quantities

/// Calculates impact parameter in the bending plane of the particle w.r.t. a point
/// \param point {x, y, z} position of the point
/// \param posSV {x, y, z} position of the secondary vertex
/// \param mom {x, y, z} particle momentum array
/// \return impact parameter in {x, y}
template <typename T, typename U, typename V>
static double ImpParXY(const T& point, const U& posSV, const array<V, 3>& mom)
{
// Ported from AliAODRecoDecay::ImpParXY
auto flightLineXY = array{posSV[0] - point[0], posSV[1] - point[1]};
auto k = dotProd(flightLineXY, array{mom[0], mom[1]}) / Pt2(mom);
auto dx = flightLineXY[0] - k * (double)mom[0];
auto dy = flightLineXY[1] - k * (double)mom[1];
auto absImpPar = sqrtSumOfSquares(dx, dy);
auto flightLine = array{posSV[0] - point[0], posSV[1] - point[1], posSV[2] - point[2]};
auto cross = crossProd(mom, flightLine);
return (cross[2] > 0. ? absImpPar : -1. * absImpPar);
}

/// Calculates the difference between measured and expected track impact parameter
/// normalized to its uncertainty
/// \param decLenXY decay lenght in {x, y} plane
/// \param errDecLenXY error on decay lenght in {x, y} plane
/// \param momMother {x, y, z} or {x, y} candidate momentum array
/// \param impParProng prong impact parameter
/// \param errImpParProng error on prong impact parameter
/// \param momProng {x, y, z} or {x, y} prong momentum array
/// \return normalized difference between expected and observed impact parameter
template <std::size_t N, std::size_t M, typename T, typename U, typename V, typename W, typename X, typename Y>
static double normImpParMeasMinusExpProng(T decLenXY, U errDecLenXY, const array<V, N>& momMother, W impParProng,
X errImpParProng, const array<Y, M>& momProng)
{
// Ported from AliAODRecoDecayHF::Getd0MeasMinusExpProng adding normalization directly in the function
auto sinThetaP = ((double)momProng[0] * (double)momMother[1] - (double)momProng[1] * (double)momMother[0]) /
(Pt(momProng) * Pt(momMother));
auto diff = impParProng - (double)decLenXY * sinThetaP;
auto errImpParExpProng = (double)errDecLenXY * sinThetaP;
auto errDiff = sqrtSumOfSquares(errImpParProng, errImpParExpProng);
return (errDiff > 0. ? diff / errDiff : 0.);
}

/// Calculates maximum normalized difference between measured and expected impact parameter of candidate prongs
/// \param posPV {x, y, z} or {x, y} position of primary vertex
/// \param posSV {x, y, z} or {x, y} position of secondary vertex
/// \param errDecLenXY error on decay lenght in {x, y} plane
/// \param momMother {x, y, z} or {x, y} candidate momentum array
/// \param arrImpPar array of prong impact parameters
/// \param arrErrImpPar array of errors on prong impact parameter (same order as arrImpPar)
/// \param momMom array of {x, y, z} or {x, y} prong momenta (same order as arrImpPar)
/// \return maximum normalized difference between expected and observed impact parameters
template <std::size_t N, std::size_t M, std::size_t K, typename T, typename U, typename V, typename W, typename X,
typename Y, typename Z>
static double maxNormalisedDeltaIP(const T& posPV, const U& posSV, V errDecLenXY, const array<W, M>& momMother,
const array<X, N>& arrImpPar, const array<Y, N>& arrErrImpPar,
const array<array<Z, K>, N>& arrMom)
{
auto decLenXY = distanceXY(posPV, posSV);
double maxNormDeltaIP{0.};
for (auto iProng = 0; iProng < N; ++iProng) {
auto prongNormDeltaIP = normImpParMeasMinusExpProng(decLenXY, errDecLenXY, momMother, arrImpPar[iProng],
arrErrImpPar[iProng], arrMom[iProng]);
if (std::abs(prongNormDeltaIP) > std::abs(maxNormDeltaIP)) {
maxNormDeltaIP = prongNormDeltaIP;
}
}
return maxNormDeltaIP;
}

/// Returns particle mass based on PDG code.
/// \param pdg PDG code
/// \return particle mass
Expand Down
11 changes: 11 additions & 0 deletions Analysis/DataModel/include/Analysis/HFSecondaryVertex.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,7 @@ DECLARE_SOA_COLUMN(ErrorDecayLengthXY, errorDecayLengthXY, float);
DECLARE_SOA_DYNAMIC_COLUMN(CPA, cpa, [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz) { return RecoDecay::CPA(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}, array{px, py, pz}); });
DECLARE_SOA_DYNAMIC_COLUMN(CPAXY, cpaXY, [](float xVtxP, float yVtxP, float xVtxS, float yVtxS, float px, float py) { return RecoDecay::CPAXY(array{xVtxP, yVtxP}, array{xVtxS, yVtxS}, array{px, py}); });
DECLARE_SOA_DYNAMIC_COLUMN(Ct, ct, [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz, double m) { return RecoDecay::Ct(array{px, py, pz}, RecoDecay::distance(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}), m); });
DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterXY, impactParameterXY, [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz) { return RecoDecay::ImpParXY(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}, array{px, py, pz}); });
} // namespace hf_cand

// specific 2-prong decay properties
Expand All @@ -131,6 +132,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProduct, impactParameterProduct, [](fl
DECLARE_SOA_DYNAMIC_COLUMN(M, m, [](float px0, float py0, float pz0, float px1, float py1, float pz1, const array<double, 2>& m) { return RecoDecay::M(array{array{px0, py0, pz0}, array{px1, py1, pz1}}, m); });
DECLARE_SOA_DYNAMIC_COLUMN(M2, m2, [](float px0, float py0, float pz0, float px1, float py1, float pz1, const array<double, 2>& m) { return RecoDecay::M2(array{array{px0, py0, pz0}, array{px1, py1, pz1}}, m); });
DECLARE_SOA_DYNAMIC_COLUMN(CosThetaStar, cosThetaStar, [](float px0, float py0, float pz0, float px1, float py1, float pz1, const array<double, 2>& m, double mTot, int iProng) { return RecoDecay::CosThetaStar(array{array{px0, py0, pz0}, array{px1, py1, pz1}}, m, mTot, iProng); });
DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProngSqSum, impactParameterProngSqSum, [](float impParProng0, float impParProng1) { return RecoDecay::sumOfSquares(impParProng0, impParProng1); });
DECLARE_SOA_DYNAMIC_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, [](float xVtxP, float yVtxP, float xVtxS, float yVtxS, float errDlxy, float pxM, float pyM, float ip0, float errIp0, float ip1, float errIp1, float px0, float py0, float px1, float py1) { return RecoDecay::maxNormalisedDeltaIP(array{xVtxP, yVtxP}, array{xVtxS, yVtxS}, errDlxy, array{pxM, pyM}, array{ip0, ip1}, array{errIp0, errIp1}, array{array{px0, py0}, array{px1, py1}}); });
// MC matching result:
// - bit 0: D0(bar) → π± K∓
DECLARE_SOA_COLUMN(FlagMCMatchRec, flagMCMatchRec, uint8_t); // reconstruction level
Expand Down Expand Up @@ -218,6 +221,7 @@ DECLARE_SOA_TABLE(HfCandProng2Base, "AOD", "HFCANDP2BASE",
hf_cand_prong2::M2<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1>,
hf_cand_prong2::ImpactParameterProduct<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1>,
hf_cand_prong2::CosThetaStar<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1>,
hf_cand_prong2::ImpactParameterProngSqSum<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1>,
/* dynamic columns that use candidate momentum components */
hf_cand::Pt<hf_cand_prong2::Px, hf_cand_prong2::Py>,
hf_cand::Pt2<hf_cand_prong2::Px, hf_cand_prong2::Py>,
Expand All @@ -227,6 +231,8 @@ DECLARE_SOA_TABLE(HfCandProng2Base, "AOD", "HFCANDP2BASE",
hf_cand::CPA<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
hf_cand::CPAXY<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand_prong2::Px, hf_cand_prong2::Py>,
hf_cand::Ct<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
hf_cand::ImpactParameterXY<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
hf_cand_prong2::MaxNormalisedDeltaIP<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ErrorDecayLengthXY, hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand::ImpactParameter0, hf_cand::ErrorImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ErrorImpactParameter1, hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PxProng1, hf_cand::PyProng1>,
hf_cand::Eta<hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
hf_cand::Y<hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
hf_cand::E<hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
Expand Down Expand Up @@ -254,6 +260,8 @@ DECLARE_SOA_EXPRESSION_COLUMN(Py, py, float, 1.f * aod::hf_cand::pyProng0 + 1.f
DECLARE_SOA_EXPRESSION_COLUMN(Pz, pz, float, 1.f * aod::hf_cand::pzProng0 + 1.f * aod::hf_cand::pzProng1 + 1.f * aod::hf_cand::pzProng2);
DECLARE_SOA_DYNAMIC_COLUMN(M, m, [](float px0, float py0, float pz0, float px1, float py1, float pz1, float px2, float py2, float pz2, const array<double, 3>& m) { return RecoDecay::M(array{array{px0, py0, pz0}, array{px1, py1, pz1}, array{px2, py2, pz2}}, m); });
DECLARE_SOA_DYNAMIC_COLUMN(M2, m2, [](float px0, float py0, float pz0, float px1, float py1, float pz1, float px2, float py2, float pz2, const array<double, 3>& m) { return RecoDecay::M2(array{array{px0, py0, pz0}, array{px1, py1, pz1}, array{px2, py2, pz2}}, m); });
DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProngSqSum, impactParameterProngSqSum, [](float impParProng0, float impParProng1, float impParProng2) { return RecoDecay::sumOfSquares(impParProng0, impParProng1, impParProng2); });
DECLARE_SOA_DYNAMIC_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, [](float xVtxP, float yVtxP, float xVtxS, float yVtxS, float errDlxy, float pxM, float pyM, float ip0, float errIp0, float ip1, float errIp1, float ip2, float errIp2, float px0, float py0, float px1, float py1, float px2, float py2) { return RecoDecay::maxNormalisedDeltaIP(array{xVtxP, yVtxP}, array{xVtxS, yVtxS}, errDlxy, array{pxM, pyM}, array{ip0, ip1, ip2}, array{errIp0, errIp1, errIp2}, array{array{px0, py0}, array{px1, py1}, array{px2, py2}}); });
// MC matching result:
// - bit 0: D± → π± K∓ π±
// - bit 1: Lc± → p± K∓ π±
Expand Down Expand Up @@ -329,6 +337,7 @@ DECLARE_SOA_TABLE(HfCandProng3Base, "AOD", "HFCANDP3BASE",
/* dynamic columns */
hf_cand_prong3::M<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, hf_cand::PxProng2, hf_cand::PyProng2, hf_cand::PzProng2>,
hf_cand_prong3::M2<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, hf_cand::PxProng2, hf_cand::PyProng2, hf_cand::PzProng2>,
hf_cand_prong3::ImpactParameterProngSqSum<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ImpactParameter2>,
/* prong 2 */
hf_cand::ImpactParameterNormalised2<hf_cand::ImpactParameter2, hf_cand::ErrorImpactParameter2>,
hf_cand::PtProng2<hf_cand::PxProng2, hf_cand::PyProng2>,
Expand All @@ -343,6 +352,8 @@ DECLARE_SOA_TABLE(HfCandProng3Base, "AOD", "HFCANDP3BASE",
hf_cand::CPA<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
hf_cand::CPAXY<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand_prong3::Px, hf_cand_prong3::Py>,
hf_cand::Ct<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
hf_cand::ImpactParameterXY<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
hf_cand_prong3::MaxNormalisedDeltaIP<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ErrorDecayLengthXY, hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand::ImpactParameter0, hf_cand::ErrorImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ErrorImpactParameter1, hf_cand::ImpactParameter2, hf_cand::ErrorImpactParameter2, hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PxProng2, hf_cand::PyProng2>,
hf_cand::Eta<hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
hf_cand::Y<hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
hf_cand::E<hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
Expand Down
Loading