Skip to content

Commit c70f9cc

Browse files
fcatalan92fcatalan92
committed
Add Dplus analysis variables and update Dplus task histograms
1 parent 52d7f21 commit c70f9cc

File tree

3 files changed

+146
-37
lines changed

3 files changed

+146
-37
lines changed

‎Analysis/Core/include/Analysis/RecoDecay.h‎

Lines changed: 84 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -101,7 +101,7 @@ class RecoDecay
101101
}
102102

103103
/// Calculates scalar product of vectors.
104-
/// \note Promotes numbers to double before squaring to avoid precision loss in float multiplication.
104+
/// \note Promotes numbers to double to avoid precision loss in float multiplication.
105105
/// \param N dimension
106106
/// \param vec1,vec2 vectors
107107
/// \return scalar product
@@ -115,6 +115,19 @@ class RecoDecay
115115
return res;
116116
}
117117

118+
/// FIXME: probably cross and dot products should be in some utility class
119+
/// Calculates cross product of vectors in three dimensions.
120+
/// \note Promotes numbers to double to avoid precision loss in float multiplication.
121+
/// \param vec1,vec2 vectors
122+
/// \return cross-product vector
123+
template <typename T, typename U>
124+
static array<double, 3> crossProd(const array<T, 3>& vec1, const array<U, 3>& vec2)
125+
{
126+
return array<double, 3>{((double)vec1[1] * (double)vec2[2]) - ((double)vec1[2] * (double)vec2[1]),
127+
((double)vec1[2] * (double)vec2[0]) - ((double)vec1[0] * (double)vec2[2]),
128+
((double)vec1[0] * (double)vec2[1]) - ((double)vec1[1] * (double)vec2[0])};
129+
}
130+
118131
/// Calculates magnitude squared of a vector.
119132
/// \param N dimension
120133
/// \param vec vector
@@ -388,6 +401,76 @@ class RecoDecay
388401
return std::sqrt(M2(args...));
389402
}
390403

404+
// Calculation of topological quantities
405+
406+
/// Calculates impact parameter in the bending plane of the particle w.r.t. a point
407+
/// \param point {x, y, z} position of the point
408+
/// \param posSV {x, y, z} position of the secondary vertex
409+
/// \param mom {x, y, z} particle momentum array
410+
/// \return impact parameter in {x, y}
411+
template <typename T, typename U, typename V>
412+
static double ImpParXY(const T& point, const U& posSV, const array<V, 3>& mom)
413+
{
414+
// Ported from AliAODRecoDecay::ImpParXY
415+
auto flightLineXY = array{posSV[0] - point[0], posSV[1] - point[1]};
416+
auto k = dotProd(flightLineXY, array{mom[0], mom[1]}) / Pt2(mom);
417+
auto dx = flightLineXY[0] - k * (double)mom[0];
418+
auto dy = flightLineXY[1] - k * (double)mom[1];
419+
auto absImpPar = sqrtSumOfSquares(dx, dy);
420+
auto flightLine = array{posSV[0] - point[0], posSV[1] - point[1], posSV[2] - point[2]};
421+
auto cross = crossProd(mom, flightLine);
422+
return (cross[2] > 0. ? absImpPar : -1. * absImpPar);
423+
}
424+
425+
/// Calculates the difference between measured and expected track impact parameter
426+
/// normalized to its uncertainty
427+
/// \param decLenXY decay lenght in {x, y} plane
428+
/// \param errDecLenXY error on decay lenght in {x, y} plane
429+
/// \param momMother {x, y, z} or {x, y} candidate momentum array
430+
/// \param impParProng prong impact parameter
431+
/// \param errImpParProng error on prong impact parameter
432+
/// \param momProng {x, y, z} or {x, y} prong momentum array
433+
/// \return normalized difference between expected and observed impact parameter
434+
template <std::size_t N, std::size_t M, typename T, typename U, typename V, typename W, typename X, typename Y>
435+
static double normImpParMeasMinusExpProng(T decLenXY, U errDecLenXY, const array<V, N>& momMother, W impParProng,
436+
X errImpParProng, const array<Y, M>& momProng)
437+
{
438+
// Ported from AliAODRecoDecayHF::Getd0MeasMinusExpProng adding normalization directly in the function
439+
auto sinThetaP = ((double)momProng[0] * (double)momMother[1] - (double)momProng[1] * (double)momMother[0]) /
440+
(Pt(momProng) * Pt(momMother));
441+
auto diff = impParProng - (double)decLenXY * sinThetaP;
442+
auto errImpParExpProng = (double)errDecLenXY * sinThetaP;
443+
auto errDiff = sqrtSumOfSquares(errImpParProng, errImpParExpProng);
444+
return (errDiff > 0. ? diff / errDiff : 0.);
445+
}
446+
447+
/// Calculates maximum normalized difference between measured and expected impact parameter of candidate prongs
448+
/// \param posPV {x, y, z} or {x, y} position of primary vertex
449+
/// \param posSV {x, y, z} or {x, y} position of secondary vertex
450+
/// \param errDecLenXY error on decay lenght in {x, y} plane
451+
/// \param momMother {x, y, z} or {x, y} candidate momentum array
452+
/// \param arrImpPar array of prong impact parameters
453+
/// \param arrErrImpPar array of errors on prong impact parameter (same order as arrImpPar)
454+
/// \param momMom array of {x, y, z} or {x, y} prong momenta (same order as arrImpPar)
455+
/// \return maximum normalized difference between expected and observed impact parameters
456+
template <std::size_t N, std::size_t M, std::size_t K, typename T, typename U, typename V, typename W, typename X,
457+
typename Y, typename Z>
458+
static double maxNormalisedDeltaIP(const T& posPV, const U& posSV, V errDecLenXY, const array<W, M>& momMother,
459+
const array<X, N>& arrImpPar, const array<Y, N>& arrErrImpPar,
460+
const array<array<Z, K>, N>& arrMom)
461+
{
462+
auto decLenXY = distanceXY(posPV, posSV);
463+
double maxNormDeltaIP{0.};
464+
for (auto iProng = 0; iProng < N; ++iProng) {
465+
auto prongNormDeltaIP = normImpParMeasMinusExpProng(decLenXY, errDecLenXY, momMother, arrImpPar[iProng],
466+
arrErrImpPar[iProng], arrMom[iProng]);
467+
if (std::abs(prongNormDeltaIP) > std::abs(maxNormDeltaIP)) {
468+
maxNormDeltaIP = prongNormDeltaIP;
469+
}
470+
}
471+
return maxNormDeltaIP;
472+
}
473+
391474
/// Returns particle mass based on PDG code.
392475
/// \param pdg PDG code
393476
/// \return particle mass

‎Analysis/DataModel/include/Analysis/HFSecondaryVertex.h‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -119,6 +119,7 @@ DECLARE_SOA_COLUMN(ErrorDecayLengthXY, errorDecayLengthXY, float);
119119
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}); });
120120
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}); });
121121
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); });
122+
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}); });
122123
} // namespace hf_cand
123124

124125
// specific 2-prong decay properties
@@ -131,6 +132,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProduct, impactParameterProduct, [](fl
131132
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); });
132133
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); });
133134
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); });
135+
DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProngSqSum, impactParameterProngSqSum, [](float impParProng0, float impParProng1) { return RecoDecay::sumOfSquares(impParProng0, impParProng1); });
136+
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}}); });
134137
// MC matching result:
135138
// - bit 0: D0(bar) → π± K∓
136139
DECLARE_SOA_COLUMN(FlagMCMatchRec, flagMCMatchRec, uint8_t); // reconstruction level
@@ -218,6 +221,7 @@ DECLARE_SOA_TABLE(HfCandProng2Base, "AOD", "HFCANDP2BASE",
218221
hf_cand_prong2::M2<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1>,
219222
hf_cand_prong2::ImpactParameterProduct<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1>,
220223
hf_cand_prong2::CosThetaStar<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1>,
224+
hf_cand_prong2::ImpactParameterProngSqSum<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1>,
221225
/* dynamic columns that use candidate momentum components */
222226
hf_cand::Pt<hf_cand_prong2::Px, hf_cand_prong2::Py>,
223227
hf_cand::Pt2<hf_cand_prong2::Px, hf_cand_prong2::Py>,
@@ -227,6 +231,8 @@ DECLARE_SOA_TABLE(HfCandProng2Base, "AOD", "HFCANDP2BASE",
227231
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>,
228232
hf_cand::CPAXY<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand_prong2::Px, hf_cand_prong2::Py>,
229233
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>,
234+
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>,
235+
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>,
230236
hf_cand::Eta<hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
231237
hf_cand::Y<hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
232238
hf_cand::E<hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
@@ -254,6 +260,8 @@ DECLARE_SOA_EXPRESSION_COLUMN(Py, py, float, 1.f * aod::hf_cand::pyProng0 + 1.f
254260
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);
255261
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); });
256262
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); });
263+
DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProngSqSum, impactParameterProngSqSum, [](float impParProng0, float impParProng1, float impParProng2) { return RecoDecay::sumOfSquares(impParProng0, impParProng1, impParProng2); });
264+
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}}); });
257265
// MC matching result:
258266
// - bit 0: D± → π± K∓ π±
259267
// - bit 1: Lc± → p± K∓ π±
@@ -329,6 +337,7 @@ DECLARE_SOA_TABLE(HfCandProng3Base, "AOD", "HFCANDP3BASE",
329337
/* dynamic columns */
330338
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>,
331339
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>,
340+
hf_cand_prong3::ImpactParameterProngSqSum<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ImpactParameter2>,
332341
/* prong 2 */
333342
hf_cand::ImpactParameterNormalised2<hf_cand::ImpactParameter2, hf_cand::ErrorImpactParameter2>,
334343
hf_cand::PtProng2<hf_cand::PxProng2, hf_cand::PyProng2>,
@@ -343,6 +352,8 @@ DECLARE_SOA_TABLE(HfCandProng3Base, "AOD", "HFCANDP3BASE",
343352
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>,
344353
hf_cand::CPAXY<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand_prong3::Px, hf_cand_prong3::Py>,
345354
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>,
355+
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>,
356+
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>,
346357
hf_cand::Eta<hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
347358
hf_cand::Y<hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
348359
hf_cand::E<hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,

0 commit comments

Comments
 (0)