Skip to content

Commit 4c13d9e

Browse files
lauraserlauraser
committed
Close Pair Rejection for Femto
1 parent cc6e44a commit 4c13d9e

File tree

1 file changed

+197
-0
lines changed

1 file changed

+197
-0
lines changed

‎PWGCF/FemtoDream/FemtoDreamDetaDphiStar.h‎

Lines changed: 197 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
11+
/// \file FemtoDreamDetaDphiStar.h
12+
/// \brief FemtoDreamDetaDphiStar - Checks particles for the close pair rejection.
13+
/// \author Laura Serksnyte, TU München, laura.serksnyte@tum.de
14+
15+
#ifndef ANALYSIS_TASKS_PWGCF_FEMTODREAM_FEMTODREAMDETADPHISTAR_H_
16+
#define ANALYSIS_TASKS_PWGCF_FEMTODREAM_FEMTODREAMDETADPHISTAR_H_
17+
18+
#include "PWGCF/DataModel/FemtoDerived.h"
19+
#include "Framework/HistogramRegistry.h"
20+
#include <string>
21+
22+
namespace o2::analysis
23+
{
24+
namespace femtoDream
25+
{
26+
namespace
27+
{
28+
static std::array<std::array<std::shared_ptr<TH2>, 2>, 2> histdetadpi{};
29+
static std::array<std::array<std::shared_ptr<TH2>, 9>, 2> histdetadpiRadii{};
30+
} // namespace
31+
32+
/// \class FemtoDreamDetaDphiStar
33+
/// \brief Class to check particles for the close pair rejection.
34+
/// \tparam partOne Type of particle 1 (Track/V0/Cascade/...)
35+
/// \tparam partTwo Type of particle 2 (Track/V0/Cascade/...)
36+
template <o2::aod::femtodreamparticle::ParticleType partOne, o2::aod::femtodreamparticle::ParticleType partTwo>
37+
class FemtoDreamDetaDphiStar
38+
{
39+
public:
40+
/// Destructor
41+
virtual ~FemtoDreamDetaDphiStar() = default;
42+
/// Initalization of the histograms and setting required values
43+
void init(HistogramRegistry* registry, HistogramRegistry* registryQA, float ldeltaPhiMax, float ldeltaEtaMax, float lmagfield, bool lplotForEveryRadii)
44+
{
45+
deltaPhiMax = ldeltaPhiMax;
46+
deltaEtaMax = ldeltaEtaMax;
47+
magfield = lmagfield;
48+
plotForEveryRadii = lplotForEveryRadii;
49+
mHistogramRegistry = registry;
50+
mHistogramRegistryQA = registryQA;
51+
52+
if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kTrack) {
53+
histdetadpi[0][0] = mHistogramRegistry->add<TH2>((static_cast<std::string>(histNames[0][0])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}});
54+
histdetadpi[0][1] = mHistogramRegistry->add<TH2>((static_cast<std::string>(histNames[1][0])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}});
55+
if (plotForEveryRadii) {
56+
for (int i = 0; i < 9; i++) {
57+
histdetadpiRadii[0][i] = mHistogramRegistryQA->add<TH2>((static_cast<std::string>(histNamesRadii[0][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}});
58+
}
59+
}
60+
}
61+
if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kV0) {
62+
for (int i = 0; i < 2; i++) {
63+
histdetadpi[0][i] = mHistogramRegistry->add<TH2>((static_cast<std::string>(histNames[0][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}});
64+
histdetadpi[1][i] = mHistogramRegistry->add<TH2>((static_cast<std::string>(histNames[1][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}});
65+
if (plotForEveryRadii) {
66+
for (int j = 0; j < 9; j++) {
67+
histdetadpiRadii[i][j] = mHistogramRegistryQA->add<TH2>((static_cast<std::string>(histNamesRadii[i][j])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}});
68+
}
69+
}
70+
}
71+
}
72+
}
73+
/// Check if pair is close or not
74+
template <typename Part, typename Parts>
75+
bool isClosePair(Part const& part1, Part const& part2, Parts const& particles)
76+
{
77+
if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kTrack) {
78+
/// Track-Track combination
79+
// check if provided particles are in agreement with the class instantiation
80+
if (part1.partType() != o2::aod::femtodreamparticle::ParticleType::kTrack || part2.partType() != o2::aod::femtodreamparticle::ParticleType::kTrack) {
81+
LOG(fatal) << "FemtoDreamDetaDphiStar: passed arguments don't agree with FemtoDreamDetaDphiStar instantiation! Please provide kTrack,kTrack candidates.";
82+
return false;
83+
}
84+
auto deta = part1.eta() - part2.eta();
85+
auto dphiAvg = AveragePhiStar(part1, part2, 0);
86+
histdetadpi[0][0]->Fill(deta, dphiAvg);
87+
if (pow(dphiAvg, 2) / pow(deltaPhiMax, 2) + pow(deta, 2) / pow(deltaEtaMax, 2) < 1.) {
88+
return false;
89+
} else {
90+
histdetadpi[0][1]->Fill(deta, dphiAvg);
91+
return true;
92+
}
93+
94+
} else if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kV0) {
95+
/// Track-V0 combination
96+
// check if provided particles are in agreement with the class instantiation
97+
if (part1.partType() != o2::aod::femtodreamparticle::ParticleType::kTrack || part2.partType() != o2::aod::femtodreamparticle::ParticleType::kV0) {
98+
LOG(fatal) << "FemtoDreamDetaDphiStar: passed arguments don't agree with FemtoDreamDetaDphiStar instantiation! Please provide kTrack,kV0 candidates.";
99+
return false;
100+
}
101+
102+
bool pass = true;
103+
for (int i = 0; i < 2; i++) {
104+
auto daughter = particles.begin() + part2.indices()[i];
105+
auto deta = part1.eta() - part2.eta();
106+
auto dphiAvg = AveragePhiStar(part1, daughter, i);
107+
histdetadpi[i][0]->Fill(deta, dphiAvg);
108+
if (pow(dphiAvg, 2) / pow(deltaPhiMax, 2) + pow(deta, 2) / pow(deltaEtaMax, 2) < 1.) {
109+
return false;
110+
} else {
111+
histdetadpi[i][1]->Fill(deta, dphiAvg);
112+
return true;
113+
}
114+
}
115+
return pass;
116+
} else {
117+
LOG(fatal) << "FemtoDreamPairCleaner: Combination of objects not defined - quitting!";
118+
return false;
119+
}
120+
}
121+
122+
private:
123+
HistogramRegistry* mHistogramRegistry = nullptr; ///< For main output
124+
HistogramRegistry* mHistogramRegistryQA = nullptr; ///< For QA output
125+
static constexpr std::string_view histNames[2][2] = {{"detadphidetadphi0Before_0", "detadphidetadphi0Before_1"},
126+
{"detadphidetadphi0After_0", "detadphidetadphi0After_1"}};
127+
static constexpr std::string_view histNamesRadii[2][9] = {{"detadphidetadphi0Before_0_0", "detadphidetadphi0Before_0_1", "detadphidetadphi0Before_0_2",
128+
"detadphidetadphi0Before_0_3", "detadphidetadphi0Before_0_4", "detadphidetadphi0Before_0_5",
129+
"detadphidetadphi0Before_0_6", "detadphidetadphi0Before_0_7", "detadphidetadphi0Before_0_8"},
130+
{"detadphidetadphi0Before_1_0", "detadphidetadphi0Before_1_1", "detadphidetadphi0Before_1_2",
131+
"detadphidetadphi0Before_1_3", "detadphidetadphi0Before_1_4", "detadphidetadphi0Before_1_5",
132+
"detadphidetadphi0Before_1_6", "detadphidetadphi0Before_1_7", "detadphidetadphi0Before_1_8"}};
133+
134+
static constexpr o2::aod::femtodreamparticle::ParticleType mPartOneType = partOne; ///< Type of particle 1
135+
static constexpr o2::aod::femtodreamparticle::ParticleType mPartTwoType = partTwo; ///< Type of particle 2
136+
137+
static constexpr float tmpRadiiTPC[9] = {85., 105., 125., 145., 165., 185., 205., 225., 245.};
138+
139+
static constexpr uint32_t kSignMinusMask = 1;
140+
static constexpr uint32_t kSignPlusMask = 1 << 1;
141+
static constexpr uint32_t kValue0 = 0;
142+
143+
float deltaPhiMax;
144+
float deltaEtaMax;
145+
float magfield;
146+
bool plotForEveryRadii = false;
147+
148+
/// Calculate phi at all required radii stored in tmpRadiiTPC
149+
template <typename T>
150+
void PhiAtRadiiTPC(const T& part, std::vector<float>& tmpVec)
151+
{
152+
153+
float phi0 = part.phi();
154+
// Start: Get the charge from cutcontainer using masks
155+
float charge;
156+
if ((part.cut() & kSignMinusMask) == kValue0 && (part.cut() & kSignPlusMask) == kValue0) {
157+
charge = 0;
158+
} else if ((part.cut() & kSignPlusMask) == kSignPlusMask) {
159+
charge = 1;
160+
} else if ((part.cut() & kSignMinusMask) == kSignMinusMask) {
161+
charge = -1;
162+
} else {
163+
LOG(fatal) << "FemtoDreamDetaDphiStar: Charge bits are set wrong!";
164+
}
165+
// End: Get the charge from cutcontainer using masks
166+
float pt = part.pt();
167+
for (size_t i = 0; i < 9; i++) {
168+
tmpVec.push_back(phi0 - std::asin(0.3 * charge * 0.1 * magfield * tmpRadiiTPC[i] * 0.01 / (2. * pt)));
169+
}
170+
}
171+
172+
/// Calculate average phi
173+
template <typename T>
174+
float AveragePhiStar(const T& part1, const T& part2, int iHist)
175+
{
176+
std::vector<float> tmpVec1;
177+
std::vector<float> tmpVec2;
178+
PhiAtRadiiTPC(part1, tmpVec1);
179+
PhiAtRadiiTPC(part2, tmpVec2);
180+
const int num = tmpVec1.size();
181+
float dPhiAvg = 0;
182+
for (size_t i = 0; i < num; i++) {
183+
float dphi = tmpVec1.at(i) - tmpVec2.at(i);
184+
dphi = TVector2::Phi_mpi_pi(dphi);
185+
dPhiAvg += dphi;
186+
if (plotForEveryRadii) {
187+
histdetadpiRadii[iHist][i]->Fill(part1.eta() - part2.eta(), dphi);
188+
}
189+
}
190+
return (dPhiAvg / (float)num);
191+
}
192+
};
193+
194+
} /* namespace femtoDream */
195+
} /* namespace o2::analysis */
196+
197+
#endif /* ANALYSIS_TASKS_PWGCF_FEMTODREAM_FEMTODREAMDETADPHISTAR_H_ */

0 commit comments

Comments
 (0)