Skip to content

Commit 106e00e

Browse files
GiorgioAlbertoLuciaGiorgioAlbertoLucia
authored
[PWGLF] now the task fills both matter and antimatter (#14172)
1 parent 18e3069 commit 106e00e

File tree

1 file changed

+171
-89
lines changed

1 file changed

+171
-89
lines changed

‎PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx‎

Lines changed: 171 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,10 @@
1818
#include "Common/DataModel/EventSelection.h"
1919
#include "Common/DataModel/TrackSelectionTables.h"
2020

21+
#include "CCDB/BasicCCDBManager.h"
22+
#include "DataFormatsParameters/GRPMagField.h"
23+
#include "DataFormatsParameters/GRPObject.h"
24+
#include "DetectorsBase/Propagator.h"
2125
#include "Framework/ASoA.h"
2226
#include "Framework/ASoAHelpers.h"
2327
#include "Framework/AnalysisDataModel.h"
@@ -37,6 +41,7 @@
3741
#include <TVector2.h>
3842
#include <TVector3.h>
3943

44+
#include <string>
4045
#include <vector>
4146

4247
using namespace std;
@@ -46,9 +51,20 @@ using namespace o2::framework::expressions;
4651
using namespace o2::constants::physics;
4752
using std::array;
4853

49-
using MCTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA, aod::McTrackLabels>;
54+
namespace
55+
{
56+
constexpr int DeuteronPDG = o2::constants::physics::kDeuteron;
57+
constexpr int HePDG = o2::constants::physics::kHelium3;
58+
constexpr int HypertritonPDG = o2::constants::physics::kHyperTriton;
59+
constexpr int HyperHelium4PDG = o2::constants::physics::kHyperHelium4;
60+
} // namespace
5061

5162
struct nucleiFromHypertritonMap {
63+
64+
using MCCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>;
65+
using MCTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::McTrackLabels>;
66+
Preslice<aod::McParticles> mMcParticlesPerCollision = o2::aod::mcparticle::mcCollisionId;
67+
5268
HistogramRegistry registryMC{
5369
"registryMC",
5470
{},
@@ -70,119 +86,185 @@ struct nucleiFromHypertritonMap {
7086
Configurable<float> max_dcaz{"max_dcaz", 0.05f, "Maximum DCAz"};
7187
Configurable<float> min_nsigmaTPC{"min_nsigmaTPC", -2.0f, "Minimum nsigma TPC"};
7288
Configurable<float> max_nsigmaTPC{"max_nsigmaTPC", +2.0f, "Maximum nsigma TPC"};
73-
Configurable<float> min_pt{"min_pt", 0.0f, "minimum pt of the tracks"};
74-
Configurable<float> max_pt{"max_pt", 10.0f, "maximum pt of the tracks"};
75-
Configurable<int> nbin_pt{"nbin_pt", 50, "number of pt bins"};
76-
Configurable<int> nbin_dca = {"nbin_dca", 50, "number of DCA bins"};
89+
90+
Configurable<float> settingCutVertex{"settingCutVertex", 10.0f, "Accepted z-vertex range"};
91+
ConfigurableAxis pt_axis{"pt_axis", {50, -10.f, 10.f}, ";signed #it{p}_{T} (GeV/#it{c});Counts"};
92+
7793
Configurable<bool> saveHelium{"saveHelium", false, "Save helium candidates"};
7894

79-
int AntideuteronPDG = -1000010020;
80-
int AntihePDG = -1000020030;
81-
int AntiHypertritonPDG = -1010010030;
82-
int AntiHyperHelium4PDG = -1010020040;
95+
Configurable<int> cfgMaterialCorrection{"cfgMaterialCorrection", static_cast<int>(o2::base::Propagator::MatCorrType::USEMatCorrLUT), "Type of material correction"};
96+
Configurable<std::string> cfgCCDBurl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
97+
Service<o2::ccdb::BasicCCDBManager> mCcdb;
98+
int mRunNumber = 0;
99+
float mBz = 0.f;
100+
101+
int mSelectedPDG = 0;
102+
o2::dataformats::DCA mDcaInfoCov;
103+
o2::track::TrackParametrizationWithError<float> mTrackParCov;
104+
o2::base::MatLayerCylSet* mLut = nullptr;
83105

84106
void init(InitContext const&)
85107
{
86-
registryMC.add("hypertritonPtGen", "hypertritonPtGen", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
108+
mCcdb->setURL(cfgCCDBurl);
109+
mCcdb->setCaching(true);
110+
mCcdb->setLocalObjectValidityChecking();
111+
mCcdb->setFatalWhenNull(false);
112+
mLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(mCcdb->get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
113+
114+
registryMC.add("hypertritonPtGen", "hypertritonPtGen", HistType::kTH1F, {pt_axis});
115+
87116
if (saveHelium) {
88-
registryMC.add("he3SecPtRec_from_hypertriton", "he3SecPtRec_from_hypertriton", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
89-
registryMC.add("he3SecPtGen_from_hypertriton", "he3SecPtGen_from_hypertriton", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
90-
registryMC.add("hyperHe4PtGen", "hyperHe4PtGen", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
91-
registryMC.add("he3SecPtRec_from_hyperHe4", "he3SecPtRec_from_hyperHe4", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
92-
registryMC.add("he3SecPtGen_from_hyperHe4", "he3SecPtGen_from_hyperHe4", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
93-
registryMC.add("he3PtRec", "he3PtRec", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
94-
registryMC.add("he3PtGen", "he3PtGen", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
117+
registryMC.add("he3SecPtRec_from_hypertriton", "he3SecPtRec_from_hypertriton", HistType::kTH1F, {pt_axis});
118+
registryMC.add("he3SecPtGen_from_hypertriton", "he3SecPtGen_from_hypertriton", HistType::kTH1F, {pt_axis});
119+
registryMC.add("hyperHe4PtGen", "hyperHe4PtGen", HistType::kTH1F, {pt_axis});
120+
registryMC.add("he3SecPtRec_from_hyperHe4", "he3SecPtRec_from_hyperHe4", HistType::kTH1F, {pt_axis});
121+
registryMC.add("he3SecPtGen_from_hyperHe4", "he3SecPtGen_from_hyperHe4", HistType::kTH1F, {pt_axis});
122+
registryMC.add("he3PtRec", "he3PtRec", HistType::kTH1F, {pt_axis});
123+
registryMC.add("he3PtGen", "he3PtGen", HistType::kTH1F, {pt_axis});
95124
} else {
96-
registryMC.add("deutSecPtRec_from_hypertriton", "deutSecPtRec_from_hypertriton", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
97-
registryMC.add("deutSecPtGen_from_hypertriton", "deutSecPtGen_from_hypertriton", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
98-
registryMC.add("deutPtRec", "deutPtRec", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
99-
registryMC.add("deutPtGen", "deutPtGen", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
125+
registryMC.add("deutSecPtRec_from_hypertriton", "deutSecPtRec_from_hypertriton", HistType::kTH1F, {pt_axis});
126+
registryMC.add("deutSecPtGen_from_hypertriton", "deutSecPtGen_from_hypertriton", HistType::kTH1F, {pt_axis});
127+
registryMC.add("deutPtRec", "deutPtRec", HistType::kTH1F, {pt_axis});
128+
registryMC.add("deutPtGen", "deutPtGen", HistType::kTH1F, {pt_axis});
100129
}
101-
}
102130

103-
void processMC(const aod::McParticles& mcParticles, const MCTracks& tracks)
104-
{
105-
int selectedPDG = 0;
106131
if (saveHelium) {
107-
selectedPDG = AntihePDG;
132+
mSelectedPDG = HePDG;
108133
} else {
109-
selectedPDG = AntideuteronPDG;
134+
mSelectedPDG = DeuteronPDG;
110135
}
136+
}
137+
138+
void initCCDB(const aod::BCsWithTimestamps::iterator& bc)
139+
{
140+
if (mRunNumber == bc.runNumber())
141+
return;
142+
143+
auto timestamp = bc.timestamp();
144+
mRunNumber = bc.runNumber();
145+
146+
o2::parameters::GRPMagField* grpmag = mCcdb->getForTimeStamp<o2::parameters::GRPMagField>("GLO/Config/GRPMagField", timestamp);
147+
o2::base::Propagator::initFieldFromGRP(grpmag);
148+
o2::base::Propagator::Instance()->setMatLUT(mLut);
149+
mBz = static_cast<float>(grpmag->getNominalL3Field());
150+
LOGF(info, "Retrieved GRP for timestamp %ull (%i) with magnetic field of %1.2f kZG", timestamp, mRunNumber, mBz);
151+
}
152+
153+
void processMC(const MCCollisions& collisions, const aod::BCsWithTimestamps&, const aod::McParticles& mcParticles, const MCTracks& tracks)
154+
{
155+
for (const auto& collision : collisions) {
156+
157+
if (!collision.sel8() || std::abs(collision.posZ()) > settingCutVertex) {
158+
return;
159+
}
160+
161+
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
162+
initCCDB(bc);
163+
164+
const int mcCollisionId = collision.mcCollisionId();
165+
auto mcParticlesThisCollision = mcParticles.sliceBy(mMcParticlesPerCollision, mcCollisionId);
166+
mcParticlesThisCollision.bindExternalIndices(&mcParticles);
167+
168+
for (const auto& mcparticle : mcParticlesThisCollision) {
169+
170+
if (((std::abs(mcparticle.pdgCode()) != HypertritonPDG && std::abs(mcparticle.pdgCode()) != HyperHelium4PDG) || !mcparticle.has_daughters()) && std::abs(mcparticle.pdgCode()) != mSelectedPDG) {
171+
continue;
172+
}
173+
174+
const float particleSign = mcparticle.pdgCode() < 0 ? -1. : 1.;
175+
176+
if (std::abs(mcparticle.pdgCode()) == HypertritonPDG) {
111177

112-
for (const auto& mcparticle : mcParticles) {
113-
if (((mcparticle.pdgCode() == AntiHypertritonPDG || mcparticle.pdgCode() == AntiHyperHelium4PDG) && mcparticle.has_daughters()) || mcparticle.pdgCode() == selectedPDG) {
114-
if (mcparticle.pdgCode() == AntiHypertritonPDG) {
115-
for (auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
116-
if (daughter.pdgCode() == selectedPDG) {
117-
registryMC.fill(HIST("hypertritonPtGen"), mcparticle.pt());
118-
if (saveHelium) {
119-
registryMC.fill(HIST("he3SecPtGen_from_hypertriton"), daughter.pt());
120-
} else {
121-
registryMC.fill(HIST("deutSecPtGen_from_hypertriton"), daughter.pt());
122-
}
178+
for (const auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
179+
if (std::abs(daughter.pdgCode()) != mSelectedPDG) {
180+
continue;
181+
}
182+
const float daughterSign = daughter.pdgCode() < 0 ? -1. : 1.;
183+
184+
registryMC.fill(HIST("hypertritonPtGen"), particleSign * mcparticle.pt());
185+
if (saveHelium) {
186+
registryMC.fill(HIST("he3SecPtGen_from_hypertriton"), daughterSign * daughter.pt());
187+
} else {
188+
registryMC.fill(HIST("deutSecPtGen_from_hypertriton"), daughterSign * daughter.pt());
123189
}
124190
}
125-
}
126-
if (mcparticle.pdgCode() == AntiHyperHelium4PDG) {
127-
for (auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
128-
if (daughter.pdgCode() == selectedPDG) {
129-
registryMC.fill(HIST("hyperHe4PtGen"), mcparticle.pt());
130-
if (saveHelium) {
131-
registryMC.fill(HIST("he3SecPtGen_from_hyperHe4"), daughter.pt());
132-
}
191+
} else if (std::abs(mcparticle.pdgCode()) == HyperHelium4PDG) {
192+
for (const auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
193+
if (std::abs(daughter.pdgCode()) != mSelectedPDG) {
194+
continue;
195+
}
196+
const float daughterSign = daughter.pdgCode() < 0 ? -1. : 1.;
197+
198+
registryMC.fill(HIST("hyperHe4PtGen"), particleSign * mcparticle.pt());
199+
if (saveHelium) {
200+
registryMC.fill(HIST("he3SecPtGen_from_hyperHe4"), daughterSign * daughter.pt());
133201
}
134202
}
203+
} else if (std::abs(mcparticle.pdgCode()) == HePDG && mcparticle.isPhysicalPrimary()) {
204+
registryMC.fill(HIST("he3PtGen"), particleSign * mcparticle.pt());
205+
} else if (std::abs(mcparticle.pdgCode()) == DeuteronPDG && mcparticle.isPhysicalPrimary()) {
206+
registryMC.fill(HIST("deutPtGen"), particleSign * mcparticle.pt());
135207
}
136-
if (mcparticle.pdgCode() == AntihePDG && mcparticle.isPhysicalPrimary()) {
137-
registryMC.fill(HIST("he3PtGen"), mcparticle.pt());
208+
}
209+
210+
const o2::math_utils::Point3D<float> collisionVertex{collision.posX(), collision.posY(), collision.posZ()};
211+
212+
for (const auto& track : tracks) {
213+
214+
if (!track.has_mcParticle()) {
215+
continue;
138216
}
139-
if (mcparticle.pdgCode() == AntideuteronPDG && mcparticle.isPhysicalPrimary()) {
140-
registryMC.fill(HIST("deutPtGen"), mcparticle.pt());
217+
218+
const auto& mcparticle = track.mcParticle();
219+
if (mcparticle.pdgCode() != mSelectedPDG) {
220+
continue;
141221
}
142-
}
143-
}
144222

145-
for (const auto& track : tracks) {
146-
if (!track.has_mcParticle()) {
147-
continue;
148-
}
149-
auto mcparticle = track.mcParticle();
150-
if (mcparticle.pdgCode() != selectedPDG) {
151-
continue;
152-
}
223+
mDcaInfoCov.set(999, 999, 999, 999, 999);
224+
setTrackParCov(track, mTrackParCov);
225+
mTrackParCov.setPID(track.pidForTracking());
226+
std::array<float, 2> dcaInfo;
227+
o2::base::Propagator::Instance()->propagateToDCA(collisionVertex, mTrackParCov, mBz, 2.f, static_cast<o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value), &dcaInfo);
153228

154-
if (track.itsNCls() < min_ITS_nClusters ||
155-
track.itsNClsInnerBarrel() < min_ITS_InnerBarrel_nClusters ||
156-
track.itsNClsInnerBarrel() > max_ITS_InnerBarrel_nClusters ||
157-
track.tpcNClsFound() < min_TPC_nClusters ||
158-
track.tpcNClsCrossedRows() < min_TPC_nCrossedRows ||
159-
track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() ||
160-
track.tpcChi2NCl() > 4.f ||
161-
track.tpcChi2NCl() < min_chi2_TPC ||
162-
track.eta() < min_eta || track.eta() > max_eta ||
163-
track.dcaXY() > max_dcaxy || track.dcaXY() < -max_dcaxy ||
164-
track.dcaZ() > max_dcaz || track.dcaZ() < -max_dcaz ||
165-
track.itsChi2NCl() > 36.f) {
166-
continue;
167-
}
168-
if (mcparticle.pdgCode() == AntideuteronPDG && mcparticle.isPhysicalPrimary()) {
169-
registryMC.fill(HIST("deutPtRec"), track.pt());
170-
}
171-
if (mcparticle.pdgCode() == AntihePDG && mcparticle.isPhysicalPrimary()) {
172-
registryMC.fill(HIST("he3PtRec"), 2 * track.pt());
173-
}
229+
const float dcaXY = dcaInfo[0];
230+
const float dcaZ = dcaInfo[1];
231+
232+
if (track.itsNCls() < min_ITS_nClusters ||
233+
track.itsNClsInnerBarrel() < min_ITS_InnerBarrel_nClusters ||
234+
track.itsNClsInnerBarrel() > max_ITS_InnerBarrel_nClusters ||
235+
track.tpcNClsFound() < min_TPC_nClusters ||
236+
track.tpcNClsCrossedRows() < min_TPC_nCrossedRows ||
237+
track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() ||
238+
track.tpcChi2NCl() > max_chi2_TPC ||
239+
track.tpcChi2NCl() < min_chi2_TPC ||
240+
track.eta() < min_eta || track.eta() > max_eta ||
241+
dcaXY > max_dcaxy || dcaXY < -max_dcaxy ||
242+
dcaZ > max_dcaz || dcaZ < -max_dcaz ||
243+
track.itsChi2NCl() > 36.f) {
244+
continue;
245+
}
246+
247+
const float particleSign = mcparticle.pdgCode() < 0 ? -1. : 1.;
174248

175-
for (auto& motherparticle : mcparticle.mothers_as<aod::McParticles>()) {
176-
if (motherparticle.pdgCode() == AntiHypertritonPDG || motherparticle.pdgCode() == AntiHyperHelium4PDG) {
177-
if (motherparticle.pdgCode() == AntiHypertritonPDG) {
178-
if (mcparticle.pdgCode() == AntihePDG) {
179-
registryMC.fill(HIST("he3SecPtRec_from_hypertriton"), 2 * track.pt());
249+
if (std::abs(mcparticle.pdgCode()) == DeuteronPDG && mcparticle.isPhysicalPrimary()) {
250+
registryMC.fill(HIST("deutPtRec"), particleSign * track.pt());
251+
}
252+
if (std::abs(mcparticle.pdgCode()) == HePDG && mcparticle.isPhysicalPrimary()) {
253+
registryMC.fill(HIST("he3PtRec"), particleSign * 2 * track.pt());
254+
}
255+
256+
for (const auto& motherparticle : mcparticle.mothers_as<aod::McParticles>()) {
257+
258+
if (std::abs(motherparticle.pdgCode()) == HypertritonPDG) {
259+
if (std::abs(mcparticle.pdgCode()) == HePDG) {
260+
registryMC.fill(HIST("he3SecPtRec_from_hypertriton"), 2 * particleSign * track.pt());
180261
} else {
181-
registryMC.fill(HIST("deutSecPtRec_from_hypertriton"), track.pt());
262+
registryMC.fill(HIST("deutSecPtRec_from_hypertriton"), particleSign * track.pt());
263+
}
264+
} else if (std::abs(motherparticle.pdgCode()) == HyperHelium4PDG) {
265+
if (std::abs(mcparticle.pdgCode()) == HePDG) {
266+
registryMC.fill(HIST("he3SecPtRec_from_hyperHe4"), 2 * particleSign * track.pt());
182267
}
183-
}
184-
if (motherparticle.pdgCode() == AntiHyperHelium4PDG) {
185-
registryMC.fill(HIST("he3SecPtRec_from_hyperHe4"), 2 * track.pt());
186268
}
187269
}
188270
}

0 commit comments

Comments
 (0)