Skip to content

Commit 427181c

Browse files
lauraserLaura Serksnyte
lauraser
and
Laura Serksnyte
authored
Study of CPR - possible bug found (#5603)
* Study of CPR - possible bug found * Update femtoDreamDetaDphiStar.h Cpplint fix --------- Co-authored-by: Laura Serksnyte <laura.serksnyte@cern.ch>
1 parent 42ecac3 commit 427181c

File tree

3 files changed

+64
-28
lines changed

3 files changed

+64
-28
lines changed

‎PWGCF/FemtoDream/Core/femtoDreamDetaDphiStar.h‎

Lines changed: 34 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -41,32 +41,33 @@ class FemtoDreamDetaDphiStar
4141
/// Destructor
4242
virtual ~FemtoDreamDetaDphiStar() = default;
4343
/// Initialization of the histograms and setting required values
44-
void init(HistogramRegistry* registry, HistogramRegistry* registryQA, float ldeltaPhiMax, float ldeltaEtaMax, bool lplotForEveryRadii)
44+
void init(HistogramRegistry* registry, HistogramRegistry* registryQA, float ldeltaPhiMax, float ldeltaEtaMax, bool lplotForEveryRadii, int meORse = 0, bool oldversion = true)
4545
{
4646
deltaPhiMax = ldeltaPhiMax;
4747
deltaEtaMax = ldeltaEtaMax;
4848
plotForEveryRadii = lplotForEveryRadii;
49+
runOldVersion = oldversion;
4950
mHistogramRegistry = registry;
5051
mHistogramRegistryQA = registryQA;
5152

5253
if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kTrack) {
5354
std::string dirName = static_cast<std::string>(dirNames[0]);
54-
histdetadpi[0][0] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[0][0])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
55-
histdetadpi[0][1] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[1][0])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
55+
histdetadpi[0][0] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[0][0]) + static_cast<std::string>(histNameSEorME[meORse])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
56+
histdetadpi[0][1] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[1][0]) + static_cast<std::string>(histNameSEorME[meORse])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
5657
if (plotForEveryRadii) {
5758
for (int i = 0; i < 9; i++) {
58-
histdetadpiRadii[0][i] = mHistogramRegistryQA->add<TH2>((dirName + static_cast<std::string>(histNamesRadii[0][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
59+
histdetadpiRadii[0][i] = mHistogramRegistryQA->add<TH2>((dirName + static_cast<std::string>(histNamesRadii[0][i]) + static_cast<std::string>(histNameSEorME[meORse])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
5960
}
6061
}
6162
}
6263
if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kV0) {
6364
for (int i = 0; i < 2; i++) {
6465
std::string dirName = static_cast<std::string>(dirNames[1]);
65-
histdetadpi[i][0] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[0][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
66-
histdetadpi[i][1] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[1][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
66+
histdetadpi[i][0] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[0][i]) + static_cast<std::string>(histNameSEorME[meORse])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
67+
histdetadpi[i][1] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[1][i]) + static_cast<std::string>(histNameSEorME[meORse])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
6768
if (plotForEveryRadii) {
6869
for (int j = 0; j < 9; j++) {
69-
histdetadpiRadii[i][j] = mHistogramRegistryQA->add<TH2>((dirName + static_cast<std::string>(histNamesRadii[i][j])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
70+
histdetadpiRadii[i][j] = mHistogramRegistryQA->add<TH2>((dirName + static_cast<std::string>(histNamesRadii[i][j]) + static_cast<std::string>(histNameSEorME[meORse])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
7071
}
7172
}
7273
}
@@ -128,6 +129,8 @@ class FemtoDreamDetaDphiStar
128129
HistogramRegistry* mHistogramRegistryQA = nullptr; ///< For QA output
129130
static constexpr std::string_view dirNames[2] = {"kTrack_kTrack/", "kTrack_kV0/"};
130131

132+
static constexpr std::string_view histNameSEorME[3] = {"_SEandME", "_SE", "_ME"};
133+
131134
static constexpr std::string_view histNames[2][2] = {{"detadphidetadphi0Before_0", "detadphidetadphi0Before_1"},
132135
{"detadphidetadphi0After_0", "detadphidetadphi0After_1"}};
133136
static constexpr std::string_view histNamesRadii[2][9] = {{"detadphidetadphi0Before_0_0", "detadphidetadphi0Before_0_1", "detadphidetadphi0Before_0_2",
@@ -150,6 +153,9 @@ class FemtoDreamDetaDphiStar
150153
float deltaEtaMax;
151154
float magfield;
152155
bool plotForEveryRadii = false;
156+
// a possible bug was found, but this must be tested on hyperloop with larger statistics
157+
// possiboility to run old code is turned on so a proper comparison of both code versions can be done
158+
bool runOldVersion = true;
153159

154160
std::array<std::array<std::shared_ptr<TH2>, 2>, 2> histdetadpi{};
155161
std::array<std::array<std::shared_ptr<TH2>, 9>, 2> histdetadpiRadii{};
@@ -175,7 +181,18 @@ class FemtoDreamDetaDphiStar
175181
// End: Get the charge from cutcontainer using masks
176182
float pt = part.pt();
177183
for (size_t i = 0; i < 9; i++) {
178-
tmpVec.push_back(phi0 - std::asin(0.3 * charge * 0.1 * magfield * tmpRadiiTPC[i] * 0.01 / (2. * pt)));
184+
if (runOldVersion) {
185+
tmpVec.push_back(phi0 - std::asin(0.3 * charge * 0.1 * magfield * tmpRadiiTPC[i] * 0.01 / (2. * pt)));
186+
}
187+
if (!runOldVersion) {
188+
auto arg = 0.3 * charge * magfield * tmpRadiiTPC[i] * 0.01 / (2. * pt);
189+
// for very low pT particles, this value goes outside of range -1 to 1 at at large tpc radius; asin fails
190+
if (abs(arg) < 1) {
191+
tmpVec.push_back(phi0 - std::asin(0.3 * charge * magfield * tmpRadiiTPC[i] * 0.01 / (2. * pt)));
192+
} else {
193+
tmpVec.push_back(999);
194+
}
195+
}
179196
}
180197
}
181198

@@ -188,16 +205,23 @@ class FemtoDreamDetaDphiStar
188205
PhiAtRadiiTPC(part1, tmpVec1);
189206
PhiAtRadiiTPC(part2, tmpVec2);
190207
int num = tmpVec1.size();
208+
int meaningfulEntries = num;
191209
float dPhiAvg = 0;
210+
float dphi;
192211
for (int i = 0; i < num; i++) {
193-
float dphi = tmpVec1.at(i) - tmpVec2.at(i);
212+
if (tmpVec1.at(i) != 999 && tmpVec2.at(i) != 999) {
213+
dphi = tmpVec1.at(i) - tmpVec2.at(i);
214+
} else {
215+
dphi = 0;
216+
meaningfulEntries = meaningfulEntries - 1;
217+
}
194218
dphi = TVector2::Phi_mpi_pi(dphi);
195219
dPhiAvg += dphi;
196220
if (plotForEveryRadii) {
197221
histdetadpiRadii[iHist][i]->Fill(part1.eta() - part2.eta(), dphi);
198222
}
199223
}
200-
return dPhiAvg / num;
224+
return dPhiAvg / static_cast<float>(meaningfulEntries);
201225
}
202226
};
203227

‎PWGCF/FemtoDream/Tasks/femtoDreamTripletTaskTrackTrackTrack.cxx‎

Lines changed: 13 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,9 @@ struct femtoDreamTripletTaskTrackTrackTrack {
5151

5252
/// Particle selection part
5353

54+
// which CPR to use, old is with a possible bug and new is fixed
55+
Configurable<bool> ConfUseOLD_possiblyWrong_CPR{"ConfUseOLD_possiblyWrong_CPR", true, "Use for old CPR, which possibly has a bug. This is implemented only for debugging reasons to compare old and new code on hyperloop datasets."};
56+
5457
/// Table for both particles
5558
Configurable<float> ConfTracksInMixedEvent{"ConfTracksInMixedEvent", 1, "Number of tracks of interest, contained in the mixed event sample: 1 - only events with at least one track of interest are used in mixing; ...; 3 - only events with at least three track of interest are used in mixing. Max value is 3"};
5659
Configurable<float> ConfMaxpT{"ConfMaxpT", 4.05f, "Maximum transverse momentum of the particles"};
@@ -107,7 +110,8 @@ struct femtoDreamTripletTaskTrackTrackTrack {
107110
FemtoDreamContainerThreeBody<femtoDreamContainerThreeBody::EventType::same, femtoDreamContainerThreeBody::Observable::Q3> sameEventCont;
108111
FemtoDreamContainerThreeBody<femtoDreamContainerThreeBody::EventType::mixed, femtoDreamContainerThreeBody::Observable::Q3> mixedEventCont;
109112
FemtoDreamPairCleaner<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kTrack> pairCleaner;
110-
FemtoDreamDetaDphiStar<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kTrack> pairCloseRejection;
113+
FemtoDreamDetaDphiStar<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kTrack> pairCloseRejectionSE;
114+
FemtoDreamDetaDphiStar<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kTrack> pairCloseRejectionME;
111115
/// Histogram output
112116
HistogramRegistry qaRegistry{"TrackQA", {}, OutputObjHandlingPolicy::AnalysisObject};
113117
HistogramRegistry resultRegistry{"Correlations", {}, OutputObjHandlingPolicy::AnalysisObject};
@@ -134,7 +138,8 @@ struct femtoDreamTripletTaskTrackTrackTrack {
134138
mixedEventCont.setPDGCodes(ConfPDGCodePart, ConfPDGCodePart, ConfPDGCodePart);
135139
pairCleaner.init(&qaRegistry); // SERKSNYTE : later check if init should be updated to have 3 separate histos
136140
if (ConfIsCPR.value) {
137-
pairCloseRejection.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value);
141+
pairCloseRejectionSE.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value, 1, ConfUseOLD_possiblyWrong_CPR);
142+
pairCloseRejectionME.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value, 2, ConfUseOLD_possiblyWrong_CPR);
138143
}
139144

140145
// get masses
@@ -199,13 +204,13 @@ struct femtoDreamTripletTaskTrackTrackTrack {
199204
/// Now build the combinations
200205
for (auto& [p1, p2, p3] : combinations(CombinationsStrictlyUpperIndexPolicy(groupSelectedParts, groupSelectedParts, groupSelectedParts))) {
201206
if (ConfIsCPR.value) {
202-
if (pairCloseRejection.isClosePair(p1, p2, parts, magFieldTesla)) {
207+
if (pairCloseRejectionSE.isClosePair(p1, p2, parts, magFieldTesla)) {
203208
continue;
204209
}
205-
if (pairCloseRejection.isClosePair(p2, p3, parts, magFieldTesla)) {
210+
if (pairCloseRejectionSE.isClosePair(p2, p3, parts, magFieldTesla)) {
206211
continue;
207212
}
208-
if (pairCloseRejection.isClosePair(p1, p3, parts, magFieldTesla)) {
213+
if (pairCloseRejectionSE.isClosePair(p1, p3, parts, magFieldTesla)) {
209214
continue;
210215
}
211216
}
@@ -318,14 +323,14 @@ struct femtoDreamTripletTaskTrackTrackTrack {
318323
{
319324
for (auto& [p1, p2, p3] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo, groupPartsThree))) {
320325
if (ConfIsCPR.value) {
321-
if (pairCloseRejection.isClosePair(p1, p2, parts, magFieldTesla)) {
326+
if (pairCloseRejectionME.isClosePair(p1, p2, parts, magFieldTesla)) {
322327
continue;
323328
}
324-
if (pairCloseRejection.isClosePair(p2, p3, parts, magFieldTesla)) {
329+
if (pairCloseRejectionME.isClosePair(p2, p3, parts, magFieldTesla)) {
325330
continue;
326331
}
327332

328-
if (pairCloseRejection.isClosePair(p1, p3, parts, magFieldTesla)) {
333+
if (pairCloseRejectionME.isClosePair(p1, p3, parts, magFieldTesla)) {
329334
continue;
330335
}
331336
}

‎PWGCF/FemtoDream/Tasks/femtoDreamTripletTaskTrackTrackV0.cxx‎

Lines changed: 17 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,9 @@ struct femtoDreamTripletTaskTrackTrackV0 {
5151
Configurable<bool> ConfMixIfTVOPairPresent{"ConfMixIfTVOPairPresent", false, "Use for mixing only events which have a TV0 pair (at least one track and one V0)"};
5252
Configurable<bool> ConfMixIfTOrVOPartsPresent{"ConfMixIfTOrVOPartsPresent", false, "Use for mixing only events which have at least one particle of interest"};
5353

54+
// which CPR to use, old is with a possible bug and new is fixed
55+
Configurable<bool> ConfUseOLD_possiblyWrong_CPR{"ConfUseOLD_possiblyWrong_CPR", true, "Use for old CPR, which possibly has a bug. This is implemented only for debugging reasons to compare old and new code on hyperloop datasets."};
56+
5457
/// Track selection
5558
Configurable<int> ConfPDGCodePart{"ConfPDGCodePart", 2212, "Particle PDG code"};
5659
Configurable<o2::aod::femtodreamparticle::cutContainerType> ConfCutPart{"ConfCutPart", 5542474, "Track - Selection bit from cutCulator"};
@@ -151,8 +154,10 @@ struct femtoDreamTripletTaskTrackTrackV0 {
151154
FemtoDreamContainerThreeBody<femtoDreamContainerThreeBody::EventType::mixed, femtoDreamContainerThreeBody::Observable::Q3> mixedEventCont;
152155
FemtoDreamPairCleaner<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kTrack> pairCleanerTrackTrack;
153156
FemtoDreamPairCleaner<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kV0> pairCleanerTrackV0;
154-
FemtoDreamDetaDphiStar<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kTrack> pairCloseRejectionTrackTrack;
155-
FemtoDreamDetaDphiStar<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kV0> pairCloseRejectionTrackV0;
157+
FemtoDreamDetaDphiStar<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kTrack> pairCloseRejectionTrackTrackSE;
158+
FemtoDreamDetaDphiStar<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kV0> pairCloseRejectionTrackV0SE;
159+
FemtoDreamDetaDphiStar<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kTrack> pairCloseRejectionTrackTrackME;
160+
FemtoDreamDetaDphiStar<aod::femtodreamparticle::ParticleType::kTrack, aod::femtodreamparticle::ParticleType::kV0> pairCloseRejectionTrackV0ME;
156161
/// Histogram output
157162
HistogramRegistry qaRegistry{"TrackQA", {}, OutputObjHandlingPolicy::AnalysisObject};
158163
HistogramRegistry resultRegistry{"Correlations", {}, OutputObjHandlingPolicy::AnalysisObject};
@@ -189,8 +194,10 @@ struct femtoDreamTripletTaskTrackTrackV0 {
189194
pairCleanerTrackTrack.init(&qaRegistry);
190195
pairCleanerTrackV0.init(&qaRegistry);
191196
if (ConfIsCPR.value) {
192-
pairCloseRejectionTrackTrack.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value);
193-
pairCloseRejectionTrackV0.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value);
197+
pairCloseRejectionTrackTrackSE.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value, 1, ConfUseOLD_possiblyWrong_CPR);
198+
pairCloseRejectionTrackV0SE.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value, 1, ConfUseOLD_possiblyWrong_CPR);
199+
pairCloseRejectionTrackTrackME.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value, 2, ConfUseOLD_possiblyWrong_CPR);
200+
pairCloseRejectionTrackV0ME.init(&resultRegistry, &qaRegistry, ConfCPRdeltaPhiMax.value, ConfCPRdeltaEtaMax.value, ConfCPRPlotPerRadii.value, 2, ConfUseOLD_possiblyWrong_CPR);
194201
}
195202

196203
// get masses
@@ -298,13 +305,13 @@ struct femtoDreamTripletTaskTrackTrackV0 {
298305

299306
// Close pair rejection
300307
if (ConfIsCPR.value) {
301-
if (pairCloseRejectionTrackTrack.isClosePair(T1, T2, parts, magFieldTesla)) {
308+
if (pairCloseRejectionTrackTrackSE.isClosePair(T1, T2, parts, magFieldTesla)) {
302309
continue;
303310
}
304-
if (pairCloseRejectionTrackV0.isClosePair(T1, V0, parts, magFieldTesla)) {
311+
if (pairCloseRejectionTrackV0SE.isClosePair(T1, V0, parts, magFieldTesla)) {
305312
continue;
306313
}
307-
if (pairCloseRejectionTrackV0.isClosePair(T2, V0, parts, magFieldTesla)) {
314+
if (pairCloseRejectionTrackV0SE.isClosePair(T2, V0, parts, magFieldTesla)) {
308315
continue;
309316
}
310317
}
@@ -485,13 +492,13 @@ struct femtoDreamTripletTaskTrackTrackV0 {
485492

486493
// Close pair rejection
487494
if (ConfIsCPR.value) {
488-
if (pairCloseRejectionTrackTrack.isClosePair(T1, T2, parts, magFieldTesla)) {
495+
if (pairCloseRejectionTrackTrackME.isClosePair(T1, T2, parts, magFieldTesla)) {
489496
continue;
490497
}
491-
if (pairCloseRejectionTrackV0.isClosePair(T1, V0, parts, magFieldTesla)) {
498+
if (pairCloseRejectionTrackV0ME.isClosePair(T1, V0, parts, magFieldTesla)) {
492499
continue;
493500
}
494-
if (pairCloseRejectionTrackV0.isClosePair(T2, V0, parts, magFieldTesla)) {
501+
if (pairCloseRejectionTrackV0ME.isClosePair(T2, V0, parts, magFieldTesla)) {
495502
continue;
496503
}
497504
}

0 commit comments

Comments
 (0)