Skip to content

Commit 5e782d0

Browse files
lauraserLaura Serksnytealibuild
authored
First code version for three-body task (#3609)
* First code version for three-body task * Please consider the following formatting changes --------- Co-authored-by: Laura Serksnyte <laura.serksnyte@cern.ch> Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent a12c4cd commit 5e782d0

4 files changed

Lines changed: 660 additions & 6 deletions

File tree

‎PWGCF/FemtoDream/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,11 @@ o2physics_add_dpl_workflow(femtodream-pair-track-track
2424
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
2525
COMPONENT_NAME Analysis)
2626

27+
o2physics_add_dpl_workflow(femtodream-pair-track-track-track
28+
SOURCES femtoDreamPairTaskTrackTrackTrack.cxx
29+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
30+
COMPONENT_NAME Analysis)
31+
2732
o2physics_add_dpl_workflow(femtodream-debug-track
2833
SOURCES femtoDreamDebugTrack.cxx
2934
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore

‎PWGCF/FemtoDream/FemtoDreamContainerThreeBody.h‎

Lines changed: 227 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,227 @@
1+
// Copyright 2019-2022 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 FemtoDreamContainer.h
13+
/// \brief Definition of the FemtoDreamContainer
14+
/// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de
15+
/// \author Valentina Mantovani Sarti, valentina.mantovani-sarti@tum.de
16+
/// \author Georgios Mantzaridis, TU München, georgios.mantzaridis@tum.de
17+
/// \author Anton Riedel, TU München, anton.riedel@tum.de
18+
/// \author Laura Serksnyte, TU München, laura.serksnyte@tum.de
19+
20+
#ifndef PWGCF_FEMTODREAM_FEMTODREAMCONTAINERTHREEBODY_H_
21+
#define PWGCF_FEMTODREAM_FEMTODREAMCONTAINERTHREEBODY_H_
22+
23+
#include <fairlogger/Logger.h>
24+
#include <vector>
25+
#include <string>
26+
27+
#include "Framework/HistogramRegistry.h"
28+
#include "FemtoDreamMath.h"
29+
30+
#include "Math/Vector4D.h"
31+
#include "TMath.h"
32+
#include "TDatabasePDG.h"
33+
34+
using namespace o2::framework;
35+
36+
namespace o2::analysis::femtoDream
37+
{
38+
39+
namespace femtoDreamContainerThreeBody
40+
{
41+
/// Femtoscopic observable to be computed
42+
enum Observable { Q3 ///< Q3
43+
};
44+
45+
/// Type of the event processind
46+
enum EventType { same, ///< Pair from same event
47+
mixed ///< Pair from mixed event
48+
};
49+
}; // namespace femtoDreamContainerThreeBody
50+
51+
/// \class FemtoDreamContainerThreeBody
52+
/// \brief Container for all histogramming related to the correlation function. The two
53+
/// particles of the pair are passed here, and the correlation function and QA histograms
54+
/// are filled according to the specified observable
55+
/// \tparam eventType Type of the event (same/mixed)
56+
/// \tparam obs Observable to be computed (Q_3/...)
57+
template <femtoDreamContainerThreeBody::EventType eventType, femtoDreamContainerThreeBody::Observable obs>
58+
class FemtoDreamContainerThreeBody
59+
{
60+
public:
61+
/// Destructor
62+
virtual ~FemtoDreamContainerThreeBody() = default;
63+
64+
/// Initializes histograms for the task in case of three-body femtoscopy
65+
/// Called by init both in case of reconstructed data/ Monte Carlo, and for Monte Carlo Truth
66+
/// \tparam T type of the axis Object
67+
/// \param folderName Name of the directory in the output file (no suffix for reconstructed data/ Monte Carlo; "_MC" for Monte Carlo Truth)
68+
/// \param femtoObs Title of the femto observable axis
69+
/// \param femtoObsAxis axis object for the femto observable axis
70+
/// \param multAxis axis object for the multiplicity axis
71+
template <typename T>
72+
void init_base(std::string folderName, std::string femtoObs, T femtoObsAxis, T multAxis)
73+
{
74+
75+
mHistogramRegistry->add((folderName + "/relTripletDist").c_str(), ("; " + femtoObs + "; Entries").c_str(), kTH1F, {femtoObsAxis});
76+
mHistogramRegistry->add((folderName + "/relTripletQ3Mult").c_str(), ("; " + femtoObs + "; Multiplicity").c_str(), kTH2F, {femtoObsAxis, multAxis});
77+
}
78+
79+
/// Initializes specialized Monte Carlo truth histograms for the task in case of three-body femtoscopy
80+
/// internal function called by init only in case of Monte Carlo truth
81+
/// \tparam T type of the xxis Object
82+
/// \param folderName Name of the directory in the output file (no suffix for reconstructed data/ Monte Carlo; "_MC" for Monte Carlo Truth)
83+
/// \param femtoObsAxis axis object for the femto observable axis
84+
template <typename T>
85+
void init_MC(std::string folderName, std::string femtoObs, T femtoObsAxis, T multAxis)
86+
{
87+
mHistogramRegistry->add((folderName + "/relPairDist_ReconNoFake").c_str(), ("; " + femtoObs + "; Entries").c_str(), kTH1F, {femtoObsAxis});
88+
mHistogramRegistry->add((folderName + "/relPairQ3Mult_ReconNoFake").c_str(), ("; " + femtoObs + "; Multiplicity").c_str(), kTH2F, {femtoObsAxis, multAxis});
89+
mHistogramRegistry->add((folderName + "/hNoMCtruthTripletCounter").c_str(), "; Counter; Entries", kTH1I, {{1, 0, 1}});
90+
mHistogramRegistry->add((folderName + "/hFakeTripletCounter").c_str(), "; Counter; Entries", kTH1I, {{1, 0, 1}});
91+
mHistogramRegistry->add((folderName + "/Q3_resolution").c_str(), "; #it{Q}_{3} reconstructed (GeV/#it{c}); #it{Q}_{3} truth (GeV/#it{c})", kTH2F, {femtoObsAxis, femtoObsAxis});
92+
}
93+
94+
/// Templated function to initialize the histograms for the task in case of three-body femtoscopy
95+
/// Always calls init_base to initialize the histograms for data/ Monte Carlo reconstructed
96+
/// In case of Monte Carlo, calls init_base again for Monte Carlo truth and the specialized function init_MC for additional histogramms
97+
/// \tparam T type of the configurable for the axis configuration
98+
/// \param registry Histogram registry to be passed
99+
/// \param Q3Bins Q3 binning for the histograms
100+
/// \param multBins multiplicity binning for the histograms
101+
/// \param isMC add Monte Carlo truth histograms to the output file
102+
template <typename T>
103+
void init(HistogramRegistry* registry, T& Q3Bins, T& multBins, bool isMC)
104+
{
105+
mHistogramRegistry = registry;
106+
std::string femtoObs;
107+
if constexpr (mFemtoObs == femtoDreamContainerThreeBody::Observable::Q3) {
108+
femtoObs = "#it{Q}_{3} (GeV/#it{c})";
109+
}
110+
std::vector<double> tmpVecMult = multBins;
111+
framework::AxisSpec multAxis = {tmpVecMult, "Multiplicity"};
112+
framework::AxisSpec femtoObsAxis = {Q3Bins, femtoObs.c_str()};
113+
114+
std::string folderName = static_cast<std::string>(mFolderSuffix[mEventType]) + static_cast<std::string>(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kRecon]);
115+
116+
init_base(folderName, femtoObs, femtoObsAxis, multAxis);
117+
if (isMC) {
118+
folderName = static_cast<std::string>(mFolderSuffix[mEventType]) + static_cast<std::string>(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]);
119+
init_base(folderName, femtoObs, femtoObsAxis, multAxis);
120+
init_MC(folderName, femtoObs, femtoObsAxis, multAxis);
121+
}
122+
}
123+
124+
/// Set the PDG codes of the three particles involved
125+
/// \param pdg1 PDG code of particle one
126+
/// \param pdg2 PDG code of particle two
127+
/// \param pdg3 PDG code of particle three
128+
void setPDGCodes(const int pdg1, const int pdg2, const int pdg3)
129+
{
130+
mMassOne = TDatabasePDG::Instance()->GetParticle(pdg1)->Mass();
131+
mMassTwo = TDatabasePDG::Instance()->GetParticle(pdg2)->Mass();
132+
mMassThree = TDatabasePDG::Instance()->GetParticle(pdg3)->Mass();
133+
mPDGOne = pdg1;
134+
mPDGTwo = pdg2;
135+
mPDGThree = pdg3;
136+
}
137+
138+
/// Pass a triplet to the container and compute all the relevant observables
139+
/// Called by setTriplet both in case of data/ and Monte Carlo reconstructed and for Monte Carlo truth
140+
/// \tparam T type of the femtodreamparticle
141+
/// \param part1 Particle one
142+
/// \param part2 Particle two
143+
/// \param part3 Particle three
144+
/// \param mult Multiplicity of the event
145+
template <o2::aod::femtodreamMCparticle::MCType mc, typename T>
146+
void setTriplet_base(const float femtoObs, T const& part1, T const& part2, T const& part3, const int mult)
147+
{
148+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/relTripletDist"), femtoObs);
149+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/relTripletQ3Mult"), femtoObs, mult);
150+
}
151+
152+
/// Called by setTriplet only in case of Monte Carlo truth
153+
/// Fills MC truth specific histogramms:
154+
/// - Q3 distribution plots with RECONSTRUCTED information but ONLY for non-fake candidates; needed for purity calculations of tracks
155+
/// - Q3 resolution matrix
156+
/// Note: Standard histograms with MC truth information are filled with the setPair_base function
157+
/// \param part1 Particle one
158+
/// \param part2 Particle two
159+
/// \param part3 Particle three
160+
/// \param mult Multiplicity of the event
161+
void setTriplet_MC(const float femtoObsMC, const float femtoObs, const int mult)
162+
{
163+
if (mHistogramRegistry) {
164+
// Fill the Q3 distributions with the reconstructed information but only for particles with the right PDG code
165+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/relPairDist_ReconNoFake"), femtoObs);
166+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/relPairkstarMult_ReconNoFake"), femtoObs, mult);
167+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/kstar_resolution"), femtoObsMC, femtoObs);
168+
}
169+
}
170+
171+
/// Templated function to handle data/ Monte Carlo reconstructed and Monte Carlo truth
172+
/// Always calls setTriplet_base to compute the observables with reconstructed data
173+
/// In case of Monte Carlo, calls setTriplet_base with MC info and specialized function setTriplet_MC for additional histogramms
174+
/// \tparam T type of the femtodreamparticle
175+
/// \param part1 Particle one
176+
/// \param part2 Particle two
177+
/// \param part3 Particle three
178+
/// \param mult Multiplicity of the event
179+
template <bool isMC, typename T>
180+
void setTriplet(T const& part1, T const& part2, T const& part3, const int mult)
181+
{
182+
float femtoObs, femtoObsMC;
183+
// Calculate femto observable and the mT with reconstructed information
184+
if constexpr (mFemtoObs == femtoDreamContainerThreeBody::Observable::Q3) {
185+
femtoObs = FemtoDreamMath::getQ3(part1, mMassOne, part2, mMassTwo, part3, mMassThree);
186+
}
187+
188+
if (mHistogramRegistry) {
189+
setTriplet_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(femtoObs, part1, part2, part3, mult);
190+
191+
if constexpr (isMC) {
192+
if (part1.has_fdMCParticle() && part2.has_fdMCParticle() && part3.has_fdMCParticle()) {
193+
// calculate the femto observable with MC truth information
194+
if constexpr (mFemtoObs == femtoDreamContainerThreeBody::Observable::Q3) {
195+
femtoObsMC = FemtoDreamMath::getQ3(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, part3.fdMCParticle(), mMassThree);
196+
}
197+
198+
if (abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo && abs(part3.fdMCParticle().pdgMCTruth()) == mPDGThree) { // Note: all triplet-histogramms are filled with MC truth information ONLY in case of non-fake candidates
199+
setTriplet_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(femtoObsMC, part1.fdMCParticle(), part2.fdMCParticle(), part3.fdMCParticle(), mult);
200+
setTriplet_MC(femtoObsMC, femtoObs, mult);
201+
} else {
202+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakeTripletCounter"), 0);
203+
}
204+
205+
} else {
206+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hNoMCtruthTripletCounter"), 0);
207+
}
208+
}
209+
}
210+
}
211+
212+
protected:
213+
HistogramRegistry* mHistogramRegistry = nullptr; ///< For QA output
214+
static constexpr std::string_view mFolderSuffix[2] = {"SameEvent", "MixedEvent"}; ///< Folder naming for the output according to mEventType
215+
static constexpr femtoDreamContainerThreeBody::Observable mFemtoObs = obs; ///< Femtoscopic observable to be computed (according to femtoDreamContainerThreeBody::Observable)
216+
static constexpr int mEventType = eventType; ///< Type of the event (same/mixed, according to femtoDreamContainerThreeBody::EventType)
217+
float mMassOne = 0.f; ///< PDG mass of particle 1
218+
float mMassTwo = 0.f; ///< PDG mass of particle 1
219+
float mMassThree = 0.f; ///< PDG mass of particle 3; if relevant
220+
int mPDGOne = 0; ///< PDG code of particle 1
221+
int mPDGTwo = 0; ///< PDG code of particle 2
222+
int mPDGThree = 0; ///< PDG code of particle 3; if relevant
223+
};
224+
225+
} // namespace o2::analysis::femtoDream
226+
227+
#endif // PWGCF_FEMTODREAM_FEMTODREAMCONTAINERTHREEBODY_H_

‎PWGCF/FemtoDream/FemtoDreamMath.h‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -92,13 +92,14 @@ class FemtoDreamMath
9292
template <typename T>
9393
static float getQ3(const T& part1, const float mass1, const T& part2, const float mass2, const T& part3, const float mass3)
9494
{
95-
float E1 = sqrt(pow(part1.px(), 2) + pow(part1.py(), 2) + pow(part1.pz(), 2) + pow(mass1, 2));
96-
float E2 = sqrt(pow(part2.px(), 2) + pow(part2.py(), 2) + pow(part2.pz(), 2) + pow(mass2, 2));
97-
float E3 = sqrt(pow(part3.px(), 2) + pow(part3.py(), 2) + pow(part3.pz(), 2) + pow(mass3, 2));
9895

99-
const ROOT::Math::PxPyPzEVector vecpart1(part1.px(), part1.py(), part1.pz(), E1);
100-
const ROOT::Math::PxPyPzEVector vecpart2(part2.px(), part2.py(), part2.pz(), E2);
101-
const ROOT::Math::PxPyPzEVector vecpart3(part3.px(), part3.py(), part3.pz(), E3);
96+
const ROOT::Math::PtEtaPhiMVector vecpart01(part1.pt(), part1.eta(), part1.phi(), mass1);
97+
const ROOT::Math::PtEtaPhiMVector vecpart02(part2.pt(), part2.eta(), part2.phi(), mass2);
98+
const ROOT::Math::PtEtaPhiMVector vecpart03(part3.pt(), part3.eta(), part3.phi(), mass3);
99+
100+
const ROOT::Math::PxPyPzEVector vecpart1(vecpart01);
101+
const ROOT::Math::PxPyPzEVector vecpart2(vecpart02);
102+
const ROOT::Math::PxPyPzEVector vecpart3(vecpart03);
102103

103104
ROOT::Math::PxPyPzEVector q12 = getqij(vecpart1, vecpart2);
104105
ROOT::Math::PxPyPzEVector q23 = getqij(vecpart2, vecpart3);

0 commit comments

Comments
 (0)