Skip to content

Commit 953f08e

Browse files
zchochulzchochul
committed
FemtoWorld_v.21.07
1 parent 74c2f14 commit 953f08e

9 files changed

Lines changed: 1137 additions & 0 deletions

‎PWGCF/FemtoWorld/CMakeLists.txt‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
# Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
# All rights not expressly granted are reserved.
4+
#
5+
# This software is distributed under the terms of the GNU General Public
6+
# License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
#
8+
# In applying this license CERN does not waive the privileges and immunities
9+
# granted to it by virtue of its status as an Intergovernmental Organization
10+
# or submit itself to any jurisdiction.
11+
12+
o2physics_add_dpl_workflow(femto-world-pair-track-track
13+
SOURCES femtoWorldPairTaskTrackTrack.cxx
14+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
15+
COMPONENT_NAME Analysis)

‎PWGCF/FemtoWorld/FemtoWorldContainer.h‎

Lines changed: 140 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,140 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
/// \file FemtoWorldContainer.h
13+
/// \brief Definition of the FemtoWorldContainer
14+
/// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de
15+
/// \author Valentina Mantovani Sarti, valentina.mantovani-sarti@tum.de
16+
17+
#ifndef ANALYSIS_TASKS_PWGCF_FEMTOWORLD_FEMTOWORDLCONTAINER_H_
18+
#define ANALYSIS_TASKS_PWGCF_FEMTOWORLD_FEMTOWORDLCONTAINER_H_
19+
20+
#include "Framework/HistogramRegistry.h"
21+
#include "FemtoWorldMath.h"
22+
23+
#include "Math/Vector4D.h"
24+
#include "TMath.h"
25+
#include "TDatabasePDG.h"
26+
27+
using namespace o2::framework;
28+
29+
namespace o2::analysis::femtoWorld
30+
{
31+
32+
namespace femtoWorldContainer
33+
{
34+
/// Femtoscopic observable to be computed
35+
enum Observable { kstar ///< kstar
36+
};
37+
38+
/// Type of the event processind
39+
enum EventType { same, ///< Pair from same event
40+
mixed ///< Pair from mixed event
41+
};
42+
}; // namespace femtoWorldContainer
43+
44+
/// \class FemtoWorldContainer
45+
/// \brief Container for all histogramming related to the correlation function. The two
46+
/// particles of the pair are passed here, and the correlation function and QA histograms
47+
/// are filled according to the specified observable
48+
/// \tparam eventType Type of the event (same/mixed)
49+
/// \tparam obs Observable to be computed (k*/Q_inv/...)
50+
template <femtoWorldContainer::EventType eventType, femtoWorldContainer::Observable obs>
51+
class FemtoWorldContainer
52+
{
53+
public:
54+
/// Destructor
55+
virtual ~FemtoWorldContainer() = default;
56+
57+
/// Initializes histograms for the task
58+
/// \tparam T Type of the configurable for the axis configuration
59+
/// \param registry Histogram registry to be passed
60+
/// \param kstarBins k* binning for the histograms
61+
/// \param multBins multiplicity binning for the histograms
62+
/// \param kTBins kT binning for the histograms
63+
/// \param mTBins mT binning for the histograms
64+
template <typename T>
65+
void init(HistogramRegistry* registry, T& kstarBins, T& multBins, T& kTBins, T& mTBins)
66+
{
67+
mHistogramRegistry = registry;
68+
std::string femtoObs;
69+
if constexpr (mFemtoObs == femtoWorldContainer::Observable::kstar) {
70+
femtoObs = "#it{k*} (GeV/#it{c})";
71+
}
72+
std::vector<double> tmpVecMult = multBins;
73+
framework::AxisSpec multAxis = {tmpVecMult, "Multiplicity"};
74+
framework::AxisSpec femtoObsAxis = {kstarBins, femtoObs.c_str()};
75+
framework::AxisSpec kTAxis = {kTBins, "#it{k}_{T} (GeV/#it{c})"};
76+
framework::AxisSpec mTAxis = {mTBins, "#it{m}_{T} (GeV/#it{c}^{2})"};
77+
78+
std::string folderName = static_cast<std::string>(mFolderSuffix[mEventType]);
79+
mHistogramRegistry->add((folderName + "relPairDist").c_str(), ("; " + femtoObs + "; Entries").c_str(), kTH1F, {femtoObsAxis});
80+
mHistogramRegistry->add((folderName + "relPairkT").c_str(), "; #it{k}_{T} (GeV/#it{c}); Entries", kTH1F, {kTAxis});
81+
mHistogramRegistry->add((folderName + "relPairkstarkT").c_str(), ("; " + femtoObs + "; #it{k}_{T} (GeV/#it{c})").c_str(), kTH2F, {femtoObsAxis, kTAxis});
82+
mHistogramRegistry->add((folderName + "relPairkstarmT").c_str(), ("; " + femtoObs + "; #it{m}_{T} (GeV/#it{c}^{2})").c_str(), kTH2F, {femtoObsAxis, mTAxis});
83+
mHistogramRegistry->add((folderName + "relPairkstarMult").c_str(), ("; " + femtoObs + "; Multiplicity").c_str(), kTH2F, {femtoObsAxis, multAxis});
84+
mHistogramRegistry->add((folderName + "kstarPtPart1").c_str(), ("; " + femtoObs + "; #it{p} _{T} Particle 1 (GeV/#it{c})").c_str(), kTH2F, {femtoObsAxis, {375, 0., 7.5}});
85+
mHistogramRegistry->add((folderName + "kstarPtPart2").c_str(), ("; " + femtoObs + "; #it{p} _{T} Particle 2 (GeV/#it{c})").c_str(), kTH2F, {femtoObsAxis, {375, 0., 7.5}});
86+
mHistogramRegistry->add((folderName + "MultPtPart1").c_str(), "; #it{p} _{T} Particle 1 (GeV/#it{c}); Multiplicity", kTH2F, {{375, 0., 7.5}, multAxis});
87+
mHistogramRegistry->add((folderName + "MultPtPart2").c_str(), "; #it{p} _{T} Particle 2 (GeV/#it{c}); Multiplicity", kTH2F, {{375, 0., 7.5}, multAxis});
88+
mHistogramRegistry->add((folderName + "PtPart1PtPart2").c_str(), "; #it{p} _{T} Particle 1 (GeV/#it{c}); #it{p} _{T} Particle 2 (GeV/#it{c})", kTH2F, {{375, 0., 7.5}, {375, 0., 7.5}});
89+
}
90+
91+
/// Set the PDG codes of the two particles involved
92+
/// \param pdg1 PDG code of particle one
93+
/// \param pdg2 PDG code of particle two
94+
void setPDGCodes(const int pdg1, const int pdg2)
95+
{
96+
mMassOne = TDatabasePDG::Instance()->GetParticle(pdg1)->Mass();
97+
mMassTwo = TDatabasePDG::Instance()->GetParticle(pdg2)->Mass();
98+
}
99+
100+
/// Pass a pair to the container and compute all the relevant observables
101+
/// \tparam T type of the femtoworldparticle
102+
/// \param part1 Particle one
103+
/// \param part2 Particle two
104+
/// \param mult Multiplicity of the event
105+
template <typename T>
106+
void setPair(T const& part1, T const& part2, const int mult)
107+
{
108+
float femtoObs;
109+
if constexpr (mFemtoObs == femtoWorldContainer::Observable::kstar) {
110+
femtoObs = FemtoWorldMath::getkstar(part1, mMassOne, part2, mMassTwo);
111+
}
112+
const float kT = FemtoWorldMath::getkT(part1, mMassOne, part2, mMassTwo);
113+
const float mT = FemtoWorldMath::getmT(part1, mMassOne, part2, mMassTwo);
114+
115+
if (mHistogramRegistry) {
116+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairDist"), femtoObs);
117+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairkT"), kT);
118+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairkstarkT"), femtoObs, kT);
119+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairkstarmT"), femtoObs, mT);
120+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairkstarMult"), femtoObs, mult);
121+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("kstarPtPart1"), femtoObs, part1.pt());
122+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("kstarPtPart2"), femtoObs, part2.pt());
123+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("MultPtPart1"), part1.pt(), mult);
124+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("MultPtPart2"), part2.pt(), mult);
125+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("PtPart1PtPart2"), part1.pt(), part2.pt());
126+
}
127+
}
128+
129+
protected:
130+
HistogramRegistry* mHistogramRegistry = nullptr; ///< For QA output
131+
static constexpr std::string_view mFolderSuffix[2] = {"SameEvent/", "MixedEvent/"}; ///< Folder naming for the output according to mEventType
132+
static constexpr femtoWorldContainer::Observable mFemtoObs = obs; ///< Femtoscopic observable to be computed (according to femtoWorldContainer::Observable)
133+
static constexpr int mEventType = eventType; ///< Type of the event (same/mixed, according to femtoWorldContainer::EventType)
134+
float mMassOne = 0.f; ///< PDG mass of particle 1
135+
float mMassTwo = 0.f; ///< PDG mass of particle 2
136+
};
137+
138+
} // namespace o2::analysis::femtoWorld
139+
140+
#endif /* ANALYSIS_TASKS_PWGCF_FEMTOWORLD_FEMTOWORDLCONTAINER_H_ */

‎PWGCF/FemtoWorld/FemtoWorldDetaDphiStar.h‎

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

0 commit comments

Comments
 (0)