Skip to content

Commit 90cc905

Browse files
ariedel-cernariedel-cern
committed
Feat: add track-track and track-v0 pairing task
1 parent 4a6f75a commit 90cc905

22 files changed

+632
-514
lines changed

‎PWGCF/Femto/Core/baseSelection.h‎

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@ class BaseSelection
6363
if (static_cast<size_t>(observableIndex) >= NumObservables) {
6464
LOG(fatal) << "Observable is not valid. Observable (index) has to be smaller than " << NumObservables;
6565
}
66-
if (skipMostPermissiveBit) {
66+
if (!selectionValues.empty() && skipMostPermissiveBit) {
6767
mNSelections += selectionValues.size() - 1;
6868
} else {
6969
mNSelections += selectionValues.size();
@@ -115,9 +115,6 @@ class BaseSelection
115115
/// \param observableIndex Index of the observable.
116116
void addSelection(int mode, int observableIndex)
117117
{
118-
if (static_cast<size_t>(observableIndex) >= NumObservables) {
119-
LOG(fatal) << "Observable is not valid. Observable (index) has to be smaller than " << NumObservables;
120-
}
121118
switch (mode) {
122119
case -1: // cut is optional and we store bit for the cut
123120
mNSelections += 1;

‎PWGCF/Femto/Core/closePairRejection.h‎

Lines changed: 83 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -16,21 +16,15 @@
1616
#ifndef PWGCF_FEMTO_CORE_CLOSEPAIRREJECTION_H_
1717
#define PWGCF_FEMTO_CORE_CLOSEPAIRREJECTION_H_
1818

19-
#include "PWGCF/Femto/Core/dataTypes.h"
2019
#include "PWGCF/Femto/Core/femtoUtils.h"
2120
#include "PWGCF/Femto/Core/histManager.h"
22-
#include "PWGCF/Femto/Core/modes.h"
23-
#include "PWGCF/Femto/DataModel/FemtoTables.h"
2421

2522
#include "Framework/HistogramRegistry.h"
2623

2724
#include <array>
2825
#include <map>
29-
#include <memory>
3026
#include <numeric>
3127
#include <string>
32-
#include <type_traits>
33-
#include <unordered_map>
3428
#include <vector>
3529

3630
namespace o2::analysis::femto
@@ -58,8 +52,8 @@ struct ConfCpr : o2::framework::ConfigurableGroup {
5852
o2::framework::Configurable<bool> on{"on", true, "Turn on CPR"};
5953
o2::framework::Configurable<float> detaMax{"detaMax", 0.01f, "Maximium deta"};
6054
o2::framework::Configurable<float> dphistarMax{"dphistarMax", 0.01f, "Maximum dphistar"};
61-
o2::framework::ConfigurableAxis binningDeta{"binningDeta", {{200, -0.2, 0.2}}, "deta"};
62-
o2::framework::ConfigurableAxis binningDphistar{"binningDphistar", {{200, -0.2, 0.2}}, "dphi"};
55+
o2::framework::ConfigurableAxis binningDeta{"binningDeta", {{500, -0.5, 0.5}}, "deta"};
56+
o2::framework::ConfigurableAxis binningDphistar{"binningDphistar", {{500, -0.5, 0.5}}, "dphi"};
6357
};
6458

6559
// tpc radii for computing phistar
@@ -69,10 +63,8 @@ constexpr std::array<float, kNradii> kTpcRadius = {85., 105., 125., 145., 165.,
6963
// directory names
7064
constexpr char PrefixTrackTrackSe[] = "CPR_TrackTrack/SE/";
7165
constexpr char PrefixTrackTrackMe[] = "CPR_TrackTrack/ME/";
72-
constexpr char PrefixTrackPosDauSe[] = "CPR_TrackPosDau/SE/";
73-
constexpr char PrefixTrackNegDauSe[] = "CPR_TrackNegDau/SE/";
74-
constexpr char PrefixTrackPosDauMe[] = "CPR_TrackPosDau/ME/";
75-
constexpr char PrefixTrackNegDauMe[] = "CPR_TrackNegDau/ME/";
66+
constexpr char PrefixTrackV0Se[] = "CPR_TrackV0/SE/";
67+
constexpr char PrefixTrackV0Me[] = "CPR_TrackV0/ME/";
7668

7769
// must be in sync with enum TrackVariables
7870
// the enum gives the correct index in the array
@@ -91,58 +83,73 @@ constexpr std::array<histmanager::HistInfo<CprHist>, kCprHistogramLast> HistTabl
9183
template <typename T>
9284
auto makeCprHistSpecMap(const T& confCpr)
9385
{
94-
std::map<CprHist, std::vector<framework::AxisSpec>> specs;
95-
for (int i = 0; i < kCprHistogramLast; ++i) {
96-
specs[static_cast<CprHist>(i)] = {confCpr.binningDeta, confCpr.binningDphistar};
97-
}
98-
return specs;
86+
return std::map<CprHist, std::vector<framework::AxisSpec>>{
87+
{kAverage, {confCpr.binningDeta, confCpr.binningDphistar}},
88+
{kRadius0, {confCpr.binningDeta, confCpr.binningDphistar}},
89+
{kRadius1, {confCpr.binningDeta, confCpr.binningDphistar}},
90+
{kRadius2, {confCpr.binningDeta, confCpr.binningDphistar}},
91+
{kRadius3, {confCpr.binningDeta, confCpr.binningDphistar}},
92+
{kRadius4, {confCpr.binningDeta, confCpr.binningDphistar}},
93+
{kRadius5, {confCpr.binningDeta, confCpr.binningDphistar}},
94+
{kRadius6, {confCpr.binningDeta, confCpr.binningDphistar}},
95+
{kRadius7, {confCpr.binningDeta, confCpr.binningDphistar}},
96+
{kRadius8, {confCpr.binningDeta, confCpr.binningDphistar}}};
9997
};
10098

10199
template <const char* prefix>
102100
class CloseTrackRejection
103101
{
104102
public:
105-
CloseTrackRejection() {}
103+
CloseTrackRejection() = default;
104+
virtual ~CloseTrackRejection() = default;
106105

107-
void init(o2::framework::HistogramRegistry* registry, std::map<CprHist, std::vector<o2::framework::AxisSpec>>& specs, float detaMax, float dphistarMax)
106+
void init(o2::framework::HistogramRegistry* registry, std::map<CprHist, std::vector<o2::framework::AxisSpec>>& specs, float detaMax, float dphistarMax, int chargeTrack1, int chargeTrack2)
108107
{
109-
mHistogramRegistry = registry;
110108
mDetaMax = detaMax;
111109
mDphistarMax = dphistarMax;
112110

113-
for (int i = 0; i < kCprHistogramLast; ++i) {
114-
auto hist = static_cast<CprHist>(i);
115-
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(hist, HistTable), GetHistDesc(hist, HistTable), GetHistType(hist, HistTable), {specs.at(hist)});
111+
if (mDetaMax < o2::constants::math::Epsilon || mDphistarMax < o2::constants::math::Epsilon) {
112+
LOG(fatal) << "Either DetaMax or DphistarMax are 0 or negative. Either turn off CPR or specify reasonable values. Breaking ...";
113+
}
114+
mChargeTrack1 = chargeTrack1;
115+
mChargeTrack2 = chargeTrack2;
116+
117+
if (utils::sign(mChargeTrack1) != utils::sign(mChargeTrack2)) {
118+
LOG(warn) << "CPR is truned on for tracks with opposite charge. Is this intended?";
116119
}
120+
121+
mHistogramRegistry = registry;
122+
123+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kAverage, HistTable), GetHistDesc(kAverage, HistTable), GetHistType(kAverage, HistTable), {specs.at(kAverage)});
124+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kRadius0, HistTable), GetHistDesc(kRadius0, HistTable), GetHistType(kRadius0, HistTable), {specs.at(kRadius0)});
125+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kRadius1, HistTable), GetHistDesc(kRadius1, HistTable), GetHistType(kRadius1, HistTable), {specs.at(kRadius1)});
126+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kRadius2, HistTable), GetHistDesc(kRadius2, HistTable), GetHistType(kRadius2, HistTable), {specs.at(kRadius2)});
127+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kRadius3, HistTable), GetHistDesc(kRadius3, HistTable), GetHistType(kRadius3, HistTable), {specs.at(kRadius3)});
128+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kRadius4, HistTable), GetHistDesc(kRadius4, HistTable), GetHistType(kRadius4, HistTable), {specs.at(kRadius4)});
129+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kRadius5, HistTable), GetHistDesc(kRadius5, HistTable), GetHistType(kRadius5, HistTable), {specs.at(kRadius5)});
130+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kRadius6, HistTable), GetHistDesc(kRadius6, HistTable), GetHistType(kRadius6, HistTable), {specs.at(kRadius6)});
131+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kRadius7, HistTable), GetHistDesc(kRadius7, HistTable), GetHistType(kRadius7, HistTable), {specs.at(kRadius7)});
132+
mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(kRadius8, HistTable), GetHistDesc(kRadius8, HistTable), GetHistType(kRadius8, HistTable), {specs.at(kRadius8)});
117133
}
118134

119135
void setMagField(float magField) { mMagField = magField; }
120136

121-
void reset()
137+
template <typename T1, typename T2>
138+
void compute(const T1& track1, const T2& track2)
122139
{
123-
mSameCharge = false;
140+
// reset values
124141
mAverageDphistar = 0.f;
125142
mDeta = 0.f;
126143
mDphistar.fill(0.f);
127-
}
128144

129-
template <typename T1, typename T2>
130-
void compute(const T1& track1, const T2& track2)
131-
{
132-
reset();
133-
if (track1.sign() != track2.sign()) {
134-
return;
135-
}
136-
mSameCharge = true;
137145
mDeta = track1.eta() - track2.eta();
138-
for (size_t i = 0; i < kTpcRadius.size(); ++i) {
139-
auto phistar1 = utils::dphistar(mMagField, kTpcRadius[i], track1.sign(), track1.pt(), track1.phi());
140-
auto phistar2 = utils::dphistar(mMagField, kTpcRadius[i], track2.sign(), track2.pt(), track2.phi());
141-
146+
for (size_t i = 0; i < kTpcRadius.size(); i++) {
147+
auto phistar1 = utils::dphistar(mMagField, kTpcRadius[i], mChargeTrack1, track1.pt(), track1.phi());
148+
auto phistar2 = utils::dphistar(mMagField, kTpcRadius[i], mChargeTrack2, track2.pt(), track2.phi());
142149
if (phistar1 && phistar2) {
143150
// if the calculation for one phistar fails, keep the default value, which is 0
144151
// this makes it more likelier for the pair to be rejected sind the averave will be biased towards lower values
145-
mDphistar[i] = phistar1.value() - phistar2.value();
152+
mDphistar.at(i) = phistar1.value() - phistar2.value();
146153
}
147154
}
148155
mAverageDphistar = std::accumulate(mDphistar.begin(), mDphistar.end(), 0.f) / mDphistar.size();
@@ -154,28 +161,31 @@ class CloseTrackRejection
154161
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kAverage, HistTable)), mDeta, mAverageDphistar);
155162

156163
// fill radii hists
157-
for (int i = 0; i < kNradii; ++i) {
158-
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(static_cast<CprHist>(kRadius0 + 1), HistTable)), mDeta, mDphistar.at(i));
159-
}
164+
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kRadius0, HistTable)), mDeta, mDphistar.at(0));
165+
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kRadius1, HistTable)), mDeta, mDphistar.at(1));
166+
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kRadius2, HistTable)), mDeta, mDphistar.at(2));
167+
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kRadius3, HistTable)), mDeta, mDphistar.at(3));
168+
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kRadius4, HistTable)), mDeta, mDphistar.at(4));
169+
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kRadius5, HistTable)), mDeta, mDphistar.at(5));
170+
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kRadius6, HistTable)), mDeta, mDphistar.at(6));
171+
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kRadius7, HistTable)), mDeta, mDphistar.at(7));
172+
mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kRadius8, HistTable)), mDeta, mDphistar.at(8));
160173
}
161174

162175
bool isClosePair() const
163176
{
164-
if (!mSameCharge) {
165-
return false;
166-
} else {
167-
return ((mAverageDphistar * mAverageDphistar) / (mDphistarMax * mDphistarMax) + (mDeta * mDeta) / (mDetaMax * mDetaMax)) < 1.f;
168-
}
177+
return std::hypot(mAverageDphistar / mDphistarMax, mDeta / mDetaMax) < 1.f;
169178
}
170179

171180
private:
181+
int mChargeTrack1 = 0;
182+
int mChargeTrack2 = 0;
172183
float mMagField = 0.f;
173184
float mAverageDphistar = 0.f;
174-
bool mSameCharge = false;
175185
float mDeta = 0.f;
176186
float mDetaMax = 0.f;
177187
float mDphistarMax = 0.f;
178-
std::array<float, kNradii> mDphistar = {};
188+
std::array<float, kNradii> mDphistar = {0.f};
179189

180190
o2::framework::HistogramRegistry* mHistogramRegistry = nullptr;
181191
};
@@ -184,10 +194,10 @@ template <const char* prefix>
184194
class ClosePairRejectionTrackTrack
185195
{
186196
public:
187-
void init(o2::framework::HistogramRegistry* registry, std::map<CprHist, std::vector<o2::framework::AxisSpec>>& specs, float detaMax, float dphistarMax, bool isActivated)
197+
void init(o2::framework::HistogramRegistry* registry, std::map<CprHist, std::vector<o2::framework::AxisSpec>>& specs, float detaMax, float dphistarMax, int signTrack1, int absChargeTrack1, int signTrack2, int AbsChargeTrack2, bool isActivated)
188198
{
189199
mIsActivated = isActivated;
190-
mCtr.init(registry, specs, detaMax, dphistarMax);
200+
mCtr.init(registry, specs, detaMax, dphistarMax, signTrack1 * absChargeTrack1, signTrack2 * AbsChargeTrack2);
191201
}
192202

193203
void setMagField(float magField) { mCtr.setMagField(magField); }
@@ -205,41 +215,48 @@ class ClosePairRejectionTrackTrack
205215
bool mIsActivated = true;
206216
};
207217

208-
template <const char* prefixPosDau, const char* prefixNegDau>
218+
template <const char* prefix>
209219
class ClosePairRejectionTrackV0 // can also be used for any particle type that has pos/neg daughters, like resonances
210220
{
211221
public:
212-
void init(o2::framework::HistogramRegistry* registry, std::map<CprHist, std::vector<o2::framework::AxisSpec>>& specs, float detaMax, float dphistarMax, bool isActivated)
222+
void init(o2::framework::HistogramRegistry* registry, std::map<CprHist, std::vector<o2::framework::AxisSpec>>& specs, float detaMax, float dphistarMax, int signTrack, int absChargeTrack, bool isActivated)
213223
{
214224
mIsActivated = isActivated;
215-
mCtrPosDau.init(registry, specs, detaMax, dphistarMax);
216-
mCtrNegDau.init(registry, specs, detaMax, dphistarMax);
225+
mSignTrack = signTrack;
226+
227+
// initialize CPR with charge of the track and the same charge for the daughter particle
228+
// absolute charge of the daughter track will be 1, so we can pass the sign directly
229+
mCtr.init(registry, specs, detaMax, dphistarMax, signTrack * absChargeTrack, signTrack);
217230
}
218231

219232
void setMagField(float magField)
220233
{
221-
mCtrPosDau.setMagField(magField);
222-
mCtrNegDau.setMagField(magField);
234+
mCtr.setMagField(magField);
223235
}
224236
template <typename T1, typename T2, typename T3>
225237
void setPair(const T1& track, const T2& v0, const T3 /*trackTable*/)
226238
{
227-
auto posDaughter = v0.template posDau_as<T3>();
228-
auto negDaughter = v0.template negDau_as<T3>();
229-
mCtrPosDau.compute(track, posDaughter);
230-
mCtrNegDau.compute(track, negDaughter);
239+
if (mSignTrack == 1) {
240+
auto daughter = v0.template posDau_as<T3>();
241+
mCtr.compute(track, daughter);
242+
} else if (mSignTrack == -1) {
243+
auto daughter = v0.template negDau_as<T3>();
244+
mCtr.compute(track, daughter);
245+
} else {
246+
LOG(fatal) << "CPR Track-V0: Wrong track sign";
247+
}
231248
}
232-
bool isClosePair() const { return mCtrPosDau.isClosePair() || mCtrNegDau.isClosePair(); }
249+
250+
bool isClosePair() const { return mCtr.isClosePair(); }
233251
void fill()
234252
{
235-
mCtrPosDau.fill();
236-
mCtrNegDau.fill();
253+
mCtr.fill();
237254
}
238255
bool isActivated() const { return mIsActivated; }
239256

240257
private:
241-
CloseTrackRejection<prefixPosDau> mCtrPosDau;
242-
CloseTrackRejection<prefixNegDau> mCtrNegDau;
258+
CloseTrackRejection<prefix> mCtr;
259+
int mSignTrack = 0;
243260
bool mIsActivated = true;
244261
};
245262

0 commit comments

Comments
 (0)