Skip to content

Commit c094921

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

File tree

1 file changed

+198
-0
lines changed

1 file changed

+198
-0
lines changed

‎PWGCF/FemtoDream/FemtoDreamDetaDphiStar.h‎

Lines changed: 198 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,198 @@
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] = std::get<std::shared_ptr<TH2>>(mHistogramRegistry->add((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] = std::get<std::shared_ptr<TH2>>(mHistogramRegistry->add((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+
std::cout << 0 << " " << i << std::endl;
58+
histdetadpiRadii[0][i] = std::get<std::shared_ptr<TH2>>(mHistogramRegistryQA->add((static_cast<std::string>(histNamesRadii[0][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}}));
59+
}
60+
}
61+
}
62+
if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kV0) {
63+
for (int i = 0; i < 2; i++) {
64+
histdetadpi[0][i] = std::get<std::shared_ptr<TH2>>(mHistogramRegistry->add((static_cast<std::string>(histNames[0][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}}));
65+
histdetadpi[1][i] = std::get<std::shared_ptr<TH2>>(mHistogramRegistry->add((static_cast<std::string>(histNames[1][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}}));
66+
if (plotForEveryRadii) {
67+
for (int j = 0; j < 9; j++) {
68+
histdetadpiRadii[i][j] = std::get<std::shared_ptr<TH2>>(mHistogramRegistryQA->add((static_cast<std::string>(histNamesRadii[i][j])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -2.5, 2.5}, {100, -2.5, 2.5}}));
69+
}
70+
}
71+
}
72+
}
73+
}
74+
/// Check if pair is close or not
75+
template <typename Part, typename Parts>
76+
bool isClosePair(Part const& part1, Part const& part2, Parts const& particles)
77+
{
78+
if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kTrack) {
79+
/// Track-Track combination
80+
// check if provided particles are in agreement with the class instantiation
81+
if (part1.partType() != o2::aod::femtodreamparticle::ParticleType::kTrack || part2.partType() != o2::aod::femtodreamparticle::ParticleType::kTrack) {
82+
LOG(fatal) << "FemtoDreamDetaDphiStar: passed arguments don't agree with FemtoDreamDetaDphiStar instantiation! Please provide kTrack,kTrack candidates.";
83+
return false;
84+
}
85+
auto deta = part1.eta() - part2.eta();
86+
auto dphiAvg = AveragePhiStar(part1, part2, 0);
87+
histdetadpi[0][0]->Fill(deta, dphiAvg);
88+
if (pow(dphiAvg, 2) / pow(deltaPhiMax, 2) + pow(deta, 2) / pow(deltaEtaMax, 2) < 1.) {
89+
return false;
90+
} else {
91+
histdetadpi[0][1]->Fill(deta, dphiAvg);
92+
return true;
93+
}
94+
95+
} else if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kV0) {
96+
/// Track-V0 combination
97+
// check if provided particles are in agreement with the class instantiation
98+
if (part1.partType() != o2::aod::femtodreamparticle::ParticleType::kTrack || part2.partType() != o2::aod::femtodreamparticle::ParticleType::kV0) {
99+
LOG(fatal) << "FemtoDreamDetaDphiStar: passed arguments don't agree with FemtoDreamDetaDphiStar instantiation! Please provide kTrack,kV0 candidates.";
100+
return false;
101+
}
102+
103+
bool pass = true;
104+
for (int i = 0; i < 2; i++) {
105+
auto daughter = particles.begin() + part2.indices()[i];
106+
auto deta = part1.eta() - part2.eta();
107+
auto dphiAvg = AveragePhiStar(part1, daughter, i);
108+
histdetadpi[i][0]->Fill(deta, dphiAvg);
109+
if (pow(dphiAvg, 2) / pow(deltaPhiMax, 2) + pow(deta, 2) / pow(deltaEtaMax, 2) < 1.) {
110+
return false;
111+
} else {
112+
histdetadpi[i][1]->Fill(deta, dphiAvg);
113+
return true;
114+
}
115+
}
116+
return pass;
117+
} else {
118+
LOG(fatal) << "FemtoDreamPairCleaner: Combination of objects not defined - quitting!";
119+
return false;
120+
}
121+
}
122+
123+
private:
124+
HistogramRegistry* mHistogramRegistry = nullptr; ///< For main output
125+
HistogramRegistry* mHistogramRegistryQA = nullptr; ///< For QA output
126+
static constexpr std::string_view histNames[2][2] = {{"detadphidetadphi0Before_0", "detadphidetadphi0Before_1"},
127+
{"detadphidetadphi0After_0", "detadphidetadphi0After_1"}};
128+
static constexpr std::string_view histNamesRadii[2][9] = {{"detadphidetadphi0Before_0_0", "detadphidetadphi0Before_0_1", "detadphidetadphi0Before_0_2",
129+
"detadphidetadphi0Before_0_3", "detadphidetadphi0Before_0_4", "detadphidetadphi0Before_0_5",
130+
"detadphidetadphi0Before_0_6", "detadphidetadphi0Before_0_7", "detadphidetadphi0Before_0_8"},
131+
{"detadphidetadphi0Before_1_0", "detadphidetadphi0Before_1_1", "detadphidetadphi0Before_1_2",
132+
"detadphidetadphi0Before_1_3", "detadphidetadphi0Before_1_4", "detadphidetadphi0Before_1_5",
133+
"detadphidetadphi0Before_1_6", "detadphidetadphi0Before_1_7", "detadphidetadphi0Before_1_8"}};
134+
135+
static constexpr o2::aod::femtodreamparticle::ParticleType mPartOneType = partOne; ///< Type of particle 1
136+
static constexpr o2::aod::femtodreamparticle::ParticleType mPartTwoType = partTwo; ///< Type of particle 2
137+
138+
static constexpr float tmpRadiiTPC[9] = {85., 105., 125., 145., 165., 185., 205., 225., 245.};
139+
140+
static constexpr uint32_t kSignMinusMask = 1;
141+
static constexpr uint32_t kSignPlusMask = 1 << 1;
142+
static constexpr uint32_t kValue0 = 0;
143+
144+
float deltaPhiMax;
145+
float deltaEtaMax;
146+
float magfield;
147+
bool plotForEveryRadii = false;
148+
149+
/// Calculate phi at all required radii stored in tmpRadiiTPC
150+
template <typename T>
151+
void PhiAtRadiiTPC(const T& part, std::vector<float>& tmpVec)
152+
{
153+
154+
float phi0 = part.phi();
155+
// Start: Get the charge from cutcontainer using masks
156+
float charge;
157+
if ((part.cut() & kSignMinusMask) == kValue0 && (part.cut() & kSignPlusMask) == kValue0) {
158+
charge = 0;
159+
} else if ((part.cut() & kSignPlusMask) == kSignPlusMask) {
160+
charge = 1;
161+
} else if ((part.cut() & kSignMinusMask) == kSignMinusMask) {
162+
charge = -1;
163+
} else {
164+
LOG(fatal) << "FemtoDreamDetaDphiStar: Charge bits are set wrong!";
165+
}
166+
// End: Get the charge from cutcontainer using masks
167+
float pt = part.pt();
168+
for (size_t i = 0; i < 9; i++) {
169+
tmpVec.push_back(phi0 - std::asin(0.3 * charge * 0.1 * magfield * tmpRadiiTPC[i] * 0.01 / (2. * pt)));
170+
}
171+
}
172+
173+
/// Calculate average phi
174+
template <typename T>
175+
float AveragePhiStar(const T& part1, const T& part2, int iHist)
176+
{
177+
std::vector<float> tmpVec1;
178+
std::vector<float> tmpVec2;
179+
PhiAtRadiiTPC(part1, tmpVec1);
180+
PhiAtRadiiTPC(part2, tmpVec2);
181+
const int num = tmpVec1.size();
182+
float dPhiAvg = 0;
183+
for (size_t i = 0; i < num; i++) {
184+
float dphi = tmpVec1.at(i) - tmpVec2.at(i);
185+
dphi = TVector2::Phi_mpi_pi(dphi);
186+
dPhiAvg += dphi;
187+
if (plotForEveryRadii) {
188+
histdetadpiRadii[iHist][i]->Fill(part1.eta() - part2.eta(), dphi);
189+
}
190+
}
191+
return (dPhiAvg / (float)num);
192+
}
193+
};
194+
195+
} /* namespace femtoDream */
196+
} /* namespace o2::analysis */
197+
198+
#endif /* ANALYSIS_TASKS_PWGCF_FEMTODREAM_FEMTODREAMDETADPHISTAR_H_ */

0 commit comments

Comments
 (0)