Skip to content

Commit c691a15

Browse files
ariedel-cernariedel-cern
authored
[PWGCF] [PWGC] Fixes in femto framework (#14686)
1 parent 5c7dfda commit c691a15

File tree

3 files changed

+123
-121
lines changed

3 files changed

+123
-121
lines changed

‎PWGCF/Femto/Core/pairHistManager.h‎

Lines changed: 11 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -304,8 +304,6 @@ class PairHistManager
304304
{
305305
mPdgMass1 = o2::analysis::femto::utils::getMass(PdgParticle1);
306306
mPdgMass2 = o2::analysis::femto::utils::getMass(PdgParticle2);
307-
mAverageMass = (mPdgMass1 + mPdgMass2) / 2.f;
308-
mReducedMass = 2.f * (mPdgMass1 * mPdgMass2) / (mPdgMass1 + mPdgMass2);
309307
}
310308
void setCharge(int chargeAbsParticle1, int chargeAbsParticle2)
311309
{
@@ -596,23 +594,28 @@ class PairHistManager
596594

597595
float getKt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2)
598596
{
599-
double kt = 0.5 * (part1 + part2).Pt();
597+
auto sum = (part1 + part2);
598+
double kt = 0.5 * sum.Pt();
600599
return static_cast<float>(kt);
601600
}
602601

603602
float getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2)
604603
{
605604
auto sum = part1 + part2;
606605
double mt = 0;
606+
double averageMass = 0;
607+
double reducedMass = 0;
607608
switch (mMtType) {
608609
case modes::TransverseMassType::kAveragePdgMass:
609-
mt = std::hypot(0.5 * sum.Pt(), mAverageMass);
610+
averageMass = 0.5 * (part1.M() + part2.M());
611+
mt = std::hypot(0.5 * sum.Pt(), averageMass);
610612
break;
611613
case modes::TransverseMassType::kReducedPdgMass:
612-
mt = std::hypot(0.5 * sum.Pt(), mReducedMass);
614+
reducedMass = 2. * (part1.M() * part2.M()) / (part1.M() + part2.M());
615+
mt = std::hypot(0.5 * sum.Pt(), reducedMass);
613616
break;
614617
case modes::TransverseMassType::kMt4Vector:
615-
mt = sum.Mt() / 2.f;
618+
mt = 0.5 * sum.Mt();
616619
break;
617620
default:
618621
LOG(fatal) << "Invalid transverse mass type, breaking...";
@@ -634,12 +637,10 @@ class PairHistManager
634637
}
635638

636639
o2::framework::HistogramRegistry* mHistogramRegistry = nullptr;
637-
float mPdgMass1 = 0.f;
638-
float mPdgMass2 = 0.f;
640+
double mPdgMass1 = 0.;
641+
double mPdgMass2 = 0.;
639642

640643
modes::TransverseMassType mMtType = modes::TransverseMassType::kAveragePdgMass;
641-
double mAverageMass = 0.f;
642-
double mReducedMass = 0.f;
643644

644645
int mAbsCharge1 = 1;
645646
int mAbsCharge2 = 1;

‎PWGCF/Femto/Core/particleCleaner.h‎

Lines changed: 90 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,9 @@ struct ConfParticleCleaner : o2::framework::ConfigurableGroup {
3232
o2::framework::Configurable<bool> activate{"activate", false, "Activate particle cleaner"};
3333
o2::framework::Configurable<std::vector<int>> requiredPdgCodes{"requiredPdgCodes", {}, "Only consider particles with this exact pdg code (including the sign!)"};
3434
o2::framework::Configurable<std::vector<int>> rejectedPdgCodes{"rejectedPdgCodes", {}, "Reject particles with this exact pdg code (including the sign!)"};
35-
o2::framework::Configurable<bool> rejectedParticleWithoutMcInformation{"rejectedParticleWithoutMcInformation", true, "If true, all particles which have no associated MC information, are rejected by default"};
35+
o2::framework::Configurable<bool> rejectParticleWithoutMcParticle{"rejectParticleWithoutMcParticle", false, "If true, particles which have no associated MC information, are rejected"};
36+
o2::framework::Configurable<bool> rejectParticleWithoutMcMother{"rejectParticleWithoutMcMother", false, "If true, particles which have no associated mother are rejected"};
37+
o2::framework::Configurable<bool> rejectParticleWithoutMcPartonicMother{"rejectParticleWithoutMcPartonicMother", false, "If true, all particles which have no associated partonic mother"};
3638
o2::framework::Configurable<std::vector<int>> requiredMotherPdgCodes{"requiredMotherPdgCodes", {}, "Only consider particles whose mothers have one of the supplied pdg codes (inclduing the sign!)"};
3739
o2::framework::Configurable<std::vector<int>> rejectMotherPdgCodes{"rejectMotherPdgCodes", {}, "Only consider particles whose mothers do not have one of the supplied pdg codes (inclduing the sign!)"};
3840
o2::framework::Configurable<std::vector<int>> requiredPartonicMotherPdgCodes{"requiredPartonicMotherPdgCodes", {}, "Only consider particles whose partonic mothers have one of the supplied pdg codes (inclduing the sign!)"};
@@ -87,122 +89,141 @@ class ParticleCleaner
8789
{
8890
mActivate = confMpc.activate.value;
8991

90-
mRejectParticleWithoutMcInformation = confMpc.rejectedParticleWithoutMcInformation.value;
92+
mRejectParticleWithoutMcParticle = confMpc.rejectParticleWithoutMcParticle.value;
93+
mRejectParticleWithoutMcMother = confMpc.rejectParticleWithoutMcMother.value;
94+
mRejectParticleWithoutMcPartonicMother = confMpc.rejectParticleWithoutMcPartonicMother.value;
9195

9296
mRequiredPdgCodes = confMpc.requiredPdgCodes.value;
9397
mRejectedPdgCodes = confMpc.rejectedPdgCodes.value;
9498

9599
mRequiredMotherPdgCodes = confMpc.requiredMotherPdgCodes.value;
96-
mRejectMotherPdgCodes = confMpc.rejectMotherPdgCodes.value;
100+
mRejectedMotherPdgCodes = confMpc.rejectMotherPdgCodes.value;
97101

98102
mRequiredPartonicMotherPdgCodes = confMpc.requiredPartonicMotherPdgCodes.value;
99-
mRejectPartonicMotherPdgCodes = confMpc.rejectPartonicMotherPdgCodes.value;
103+
mRejectedPartonicMotherPdgCodes = confMpc.rejectPartonicMotherPdgCodes.value;
100104
}
101105

102106
template <typename T1, typename T2, typename T3, typename T4>
103-
bool isClean(T1 const& particle, T2 const& /*mcParticles*/, T3 const& /*mcMothers*/, T4 const& /*mcPartonicMothers*/)
107+
bool isClean(T1 const& particle,
108+
T2 const& /*mcParticles*/,
109+
T3 const& /*mcMothers*/,
110+
T4 const& /*mcPartonicMothers*/)
104111
{
105-
// check whether we apply cuts at all
106112
if (!mActivate) {
107113
return true;
108114
}
109-
// check if we even have an associated mc particle
115+
116+
bool hasRequiredPdgCode = true;
117+
bool hasRejectedPdgCode = false;
118+
119+
bool hasMotherWithRequiredPdgCode = true;
120+
bool hasMotherWithRejectedPdgCode = false;
121+
122+
bool hasPartonicMotherWithRequiredPdgCode = true;
123+
bool hasPartonicMotherWithRejectedPdgCode = false;
124+
125+
// MC particle
110126
if (!particle.has_fMcParticle()) {
111-
if (mRejectParticleWithoutMcInformation) {
127+
if (mRejectParticleWithoutMcParticle || !mRequiredPdgCodes.empty()) {
112128
return false;
113-
} else {
114-
return true;
115129
}
116-
}
117-
// perfrom cuts based on mc information of the particle itself
118-
auto mcParticle = particle.template fMcParticle_as<T2>();
119-
120-
// if list is empty, set it to true and skip the looop
121-
bool hasRequiredPdgCode = true;
122-
if (!mRequiredPdgCodes.empty()) {
123-
hasRequiredPdgCode = false;
124-
for (int const& pdgCode : mRequiredPdgCodes) {
125-
if (pdgCode == mcParticle.pdgCode()) {
126-
hasRequiredPdgCode = true;
127-
break;
130+
} else {
131+
auto mcParticle = particle.template fMcParticle_as<T2>();
132+
133+
if (!mRequiredPdgCodes.empty()) {
134+
hasRequiredPdgCode = false;
135+
for (int const& pdgCode : mRequiredPdgCodes) {
136+
if (pdgCode == mcParticle.pdgCode()) {
137+
hasRequiredPdgCode = true;
138+
break;
139+
}
128140
}
129141
}
130-
}
131142

132-
bool hasRejectedPdgCode = false;
133-
if (!mRejectedPdgCodes.empty()) {
134-
for (int const& pdgCode : mRejectedPdgCodes) {
135-
if (pdgCode == mcParticle.pdgCode()) {
136-
hasRejectedPdgCode = true;
137-
break;
143+
if (!mRejectedPdgCodes.empty()) {
144+
for (int const& pdgCode : mRejectedPdgCodes) {
145+
if (pdgCode == mcParticle.pdgCode()) {
146+
hasRejectedPdgCode = true;
147+
break;
148+
}
138149
}
139150
}
140151
}
141152

142-
// perfrom cuts based on mc information of the mothers
143-
auto mother = particle.template fMcMother_as<T3>();
144-
145-
// if list is empty, set it to true and skip the looop
146-
bool hasMotherWithRequiredPdgCode = true;
147-
if (!mRequiredMotherPdgCodes.empty()) {
148-
hasMotherWithRequiredPdgCode = false;
149-
for (int const& pdgCode : mRequiredMotherPdgCodes) {
150-
if (pdgCode == mother.pdgCode()) {
151-
hasMotherWithRequiredPdgCode = true;
152-
break;
153+
// MC mother
154+
if (!particle.has_fMcMother()) {
155+
if (mRejectParticleWithoutMcMother || !mRequiredMotherPdgCodes.empty()) {
156+
return false;
157+
}
158+
} else {
159+
auto mother = particle.template fMcMother_as<T3>();
160+
161+
if (!mRequiredMotherPdgCodes.empty()) {
162+
hasMotherWithRequiredPdgCode = false;
163+
for (int const& pdgCode : mRequiredMotherPdgCodes) {
164+
if (pdgCode == mother.pdgCode()) {
165+
hasMotherWithRequiredPdgCode = true;
166+
break;
167+
}
153168
}
154169
}
155-
}
156170

157-
bool hasMotherWithRejectedPdgCode = false;
158-
if (!mRejectMotherPdgCodes.empty()) {
159-
for (int const& pdgCode : mRejectMotherPdgCodes) {
160-
if (pdgCode == mother.pdgCode()) {
161-
hasMotherWithRejectedPdgCode = true;
162-
break;
171+
if (!mRejectedMotherPdgCodes.empty()) {
172+
for (int const& pdgCode : mRejectedMotherPdgCodes) {
173+
if (pdgCode == mother.pdgCode()) {
174+
hasMotherWithRejectedPdgCode = true;
175+
break;
176+
}
163177
}
164178
}
165179
}
166180

167-
// perfrom cuts based on mc information of the partonic mothers
168-
auto partonicMother = particle.template fMcPartMoth_as<T4>();
169-
170-
// if list is empty, set it to true and skip the looop
171-
bool hasPartonicMotherWithRequiredPdgCode = true;
172-
if (!mRequiredPartonicMotherPdgCodes.empty()) {
173-
hasPartonicMotherWithRequiredPdgCode = false;
174-
for (int const& pdgCode : mRequiredPartonicMotherPdgCodes) {
175-
if (pdgCode == partonicMother.pdgCode()) {
176-
hasPartonicMotherWithRequiredPdgCode = true;
177-
break;
181+
// MC partonic mother
182+
if (!particle.has_fMcPartMoth()) {
183+
if (mRejectParticleWithoutMcPartonicMother ||
184+
!mRequiredPartonicMotherPdgCodes.empty()) {
185+
return false;
186+
}
187+
} else {
188+
auto partonicMother = particle.template fMcPartMoth_as<T4>();
189+
190+
if (!mRequiredPartonicMotherPdgCodes.empty()) {
191+
hasPartonicMotherWithRequiredPdgCode = false;
192+
for (int const& pdgCode : mRequiredPartonicMotherPdgCodes) {
193+
if (pdgCode == partonicMother.pdgCode()) {
194+
hasPartonicMotherWithRequiredPdgCode = true;
195+
break;
196+
}
178197
}
179198
}
180-
}
181199

182-
bool hasPartonicMotherWithRejectedPdgCode = false;
183-
if (!mRejectPartonicMotherPdgCodes.empty()) {
184-
for (int const& pdgCode : mRejectPartonicMotherPdgCodes) {
185-
if (pdgCode == partonicMother.pdgCode()) {
186-
hasPartonicMotherWithRejectedPdgCode = true;
187-
break;
200+
if (!mRejectedPartonicMotherPdgCodes.empty()) {
201+
for (int const& pdgCode : mRejectedPartonicMotherPdgCodes) {
202+
if (pdgCode == partonicMother.pdgCode()) {
203+
hasPartonicMotherWithRejectedPdgCode = true;
204+
break;
205+
}
188206
}
189207
}
190208
}
191209

192210
return hasRequiredPdgCode && !hasRejectedPdgCode &&
193211
hasMotherWithRequiredPdgCode && !hasMotherWithRejectedPdgCode &&
194-
hasPartonicMotherWithRequiredPdgCode && !hasPartonicMotherWithRejectedPdgCode;
212+
hasPartonicMotherWithRequiredPdgCode &&
213+
!hasPartonicMotherWithRejectedPdgCode;
195214
}
196215

197216
private:
198217
bool mActivate = false;
199-
bool mRejectParticleWithoutMcInformation = true;
218+
bool mRejectParticleWithoutMcParticle = true;
219+
bool mRejectParticleWithoutMcMother = true;
220+
bool mRejectParticleWithoutMcPartonicMother = true;
200221
std::vector<int> mRequiredPdgCodes = {};
201222
std::vector<int> mRejectedPdgCodes = {};
202223
std::vector<int> mRequiredMotherPdgCodes{};
203-
std::vector<int> mRejectMotherPdgCodes{};
224+
std::vector<int> mRejectedMotherPdgCodes{};
204225
std::vector<int> mRequiredPartonicMotherPdgCodes{};
205-
std::vector<int> mRejectPartonicMotherPdgCodes{};
226+
std::vector<int> mRejectedPartonicMotherPdgCodes{};
206227
};
207228

208229
} // namespace particlecleaner

‎PWGCF/Femto/Core/tripletHistManager.h‎

Lines changed: 22 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
#include <Math/GenVector/Boost.h>
2828
#include <Math/Vector4D.h>
2929

30+
#include <algorithm>
3031
#include <array>
3132
#include <cmath>
3233
#include <map>
@@ -273,9 +274,9 @@ class TripletHistManager
273274
mParticle3 = ROOT::Math::PtEtaPhiMVector(mAbsCharge3 * particle3.pt(), particle3.eta(), particle3.phi(), mPdgMass3);
274275

275276
// set mT
276-
mMt = getMt(mParticle1, mPdgMass1, mParticle2, mPdgMass2, mParticle3, mPdgMass3);
277+
mMt = getMt(mParticle1, mParticle2, mParticle3);
277278
// set Q3
278-
mQ3 = getQ3(mParticle1, mPdgMass1, mParticle2, mPdgMass2, mParticle3, mPdgMass3);
279+
mQ3 = getQ3(mParticle1, mParticle2, mParticle3);
279280

280281
// if one of the particles has a mass getter, we cache the value for the filling later
281282
if constexpr (modes::hasMass(particleType1)) {
@@ -323,9 +324,9 @@ class TripletHistManager
323324
mTrueParticle3 = ROOT::Math::PtEtaPhiMVector(mAbsCharge3 * mcParticle3.pt(), mcParticle3.eta(), mcParticle3.phi(), mPdgMass3);
324325

325326
// set true mT
326-
mTrueMt = getMt(mTrueParticle1, mPdgMass1, mTrueParticle2, mPdgMass2, mTrueParticle3, mPdgMass3);
327+
mTrueMt = getMt(mTrueParticle1, mTrueParticle2, mTrueParticle3);
327328
// set true Q3
328-
mTrueQ3 = getQ3(mTrueParticle1, mPdgMass1, mTrueParticle2, mPdgMass2, mTrueParticle3, mPdgMass3);
329+
mTrueQ3 = getQ3(mTrueParticle1, mTrueParticle2, mTrueParticle3);
329330
}
330331

331332
template <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6>
@@ -382,26 +383,13 @@ class TripletHistManager
382383
float getQ3() const { return mQ3; }
383384

384385
private:
385-
template <typename T>
386-
T getqij(const T& vecparti, const T& vecpartj)
386+
float getQ3(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, ROOT::Math::PtEtaPhiMVector const& part3)
387387
{
388-
T trackSum = vecparti + vecpartj;
389-
T trackDifference = vecparti - vecpartj;
390-
float scaling = trackDifference.Dot(trackSum) / trackSum.Dot(trackSum);
391-
return trackDifference - scaling * trackSum;
392-
}
393-
394-
template <typename T>
395-
float getQ3(const T& part1, const float mass1, const T& part2, const float mass2, const T& part3, const float mass3)
396-
{
397-
ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1);
398-
ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2);
399-
ROOT::Math::PtEtaPhiMVector vecpart3(part3.pt(), part3.eta(), part3.phi(), mass3);
400-
auto q12 = getqij(vecpart1, vecpart2);
401-
auto q23 = getqij(vecpart2, vecpart3);
402-
auto q31 = getqij(vecpart3, vecpart1);
403-
float q = q12.M2() + q23.M2() + q31.M2();
404-
return std::sqrt(-q);
388+
double q12 = -(part1 - part2).M2();
389+
double q23 = -(part2 - part3).M2();
390+
double q31 = -(part3 - part1).M2();
391+
double q = q12 + q23 + q31;
392+
return static_cast<float>(std::sqrt(std::max(0., q))); // protection if rounding error give a value below 0
405393
}
406394

407395
void initAnalysis(std::map<TripletHist, std::vector<o2::framework::AxisSpec>> const& Specs)
@@ -497,22 +485,19 @@ class TripletHistManager
497485
}
498486
}
499487

500-
float getMt(ROOT::Math::PtEtaPhiMVector const& part1,
501-
float pdgMass1,
502-
ROOT::Math::PtEtaPhiMVector const& part2,
503-
float pdgMass2)
488+
double getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2)
504489
{
505490
auto sum = part1 + part2;
506-
float mt = 0;
507-
float averageMass = 0;
508-
float reducedMass = 0;
491+
double mt = 0.;
492+
double averageMass = 0.;
493+
double reducedMass = 0.;
509494
switch (mMtType) {
510495
case modes::TransverseMassType::kAveragePdgMass:
511-
averageMass = 0.5 * (pdgMass1 + pdgMass2);
496+
averageMass = 0.5 * (part1.M() + part2.M());
512497
mt = std::hypot(0.5 * sum.Pt(), averageMass);
513498
break;
514499
case modes::TransverseMassType::kReducedPdgMass:
515-
reducedMass = 2.f * (pdgMass1 * pdgMass2 / (pdgMass1 + pdgMass2));
500+
reducedMass = 2. * (part1.M() * part2.M() / (part1.M() + part2.M()));
516501
mt = std::hypot(0.5 * sum.Pt(), reducedMass);
517502
break;
518503
case modes::TransverseMassType::kMt4Vector:
@@ -524,20 +509,15 @@ class TripletHistManager
524509
return mt;
525510
}
526511

527-
float getMt(ROOT::Math::PtEtaPhiMVector const& part1,
528-
float mass1,
529-
ROOT::Math::PtEtaPhiMVector const& part2,
530-
float mass2,
531-
ROOT::Math::PtEtaPhiMVector const& part3,
532-
float mass3)
512+
float getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, ROOT::Math::PtEtaPhiMVector const& part3)
533513
{
534-
return (getMt(part1, mass1, part2, mass2) + getMt(part2, mass2, part3, mass3) + getMt(part1, mass1, part3, mass3)) / 3.;
514+
return static_cast<float>((getMt(part1, part2) + getMt(part2, part3) + getMt(part1, part3)) / 3.);
535515
}
536516

537517
o2::framework::HistogramRegistry* mHistogramRegistry = nullptr;
538-
float mPdgMass1 = 0.f;
539-
float mPdgMass2 = 0.f;
540-
float mPdgMass3 = 0.f;
518+
double mPdgMass1 = 0.;
519+
double mPdgMass2 = 0.;
520+
double mPdgMass3 = 0.;
541521

542522
modes::TransverseMassType mMtType = modes::TransverseMassType::kAveragePdgMass;
543523

0 commit comments

Comments
 (0)