Skip to content

Commit cc8d9a9

Browse files
GiorgioAlbertoLuciaGiorgioAlbertoLucia
committed
now the task fills both matter and antimatter
1 parent cfc23d1 commit cc8d9a9

File tree

1 file changed

+170
-89
lines changed

1 file changed

+170
-89
lines changed

‎PWGLF/Tasks/Nuspex/nucleiFromHypertritonMap.cxx‎

Lines changed: 170 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"
@@ -46,9 +50,20 @@ using namespace o2::framework::expressions;
4650
using namespace o2::constants::physics;
4751
using std::array;
4852

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

5161
struct nucleiFromHypertritonMap {
62+
63+
using MCCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>;
64+
using MCTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::McTrackLabels>;
65+
Preslice<aod::McParticles> mMcParticlesPerCollision = o2::aod::mcparticle::mcCollisionId;
66+
5267
HistogramRegistry registryMC{
5368
"registryMC",
5469
{},
@@ -70,119 +85,185 @@ struct nucleiFromHypertritonMap {
7085
Configurable<float> max_dcaz{"max_dcaz", 0.05f, "Maximum DCAz"};
7186
Configurable<float> min_nsigmaTPC{"min_nsigmaTPC", -2.0f, "Minimum nsigma TPC"};
7287
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"};
88+
89+
Configurable<float> settingCutVertex{"settingCutVertex", 10.0f, "Accepted z-vertex range"};
90+
ConfigurableAxis pt_axis{"pt_axis", {50, -10.f, 10.f}, ";signed #it{p}_{T} (GeV/#it{c});Counts"};
91+
7792
Configurable<bool> saveHelium{"saveHelium", false, "Save helium candidates"};
7893

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

84105
void init(InitContext const&)
85106
{
86-
registryMC.add("hypertritonPtGen", "hypertritonPtGen", HistType::kTH1F, {{nbin_pt, min_pt, max_pt, "p_{T} (GeV/c)"}});
107+
mCcdb->setURL(cfgCCDBurl);
108+
mCcdb->setCaching(true);
109+
mCcdb->setLocalObjectValidityChecking();
110+
mCcdb->setFatalWhenNull(false);
111+
mLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(mCcdb->get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
112+
113+
registryMC.add("hypertritonPtGen", "hypertritonPtGen", HistType::kTH1F, {pt_axis});
114+
87115
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)"}});
116+
registryMC.add("he3SecPtRec_from_hypertriton", "he3SecPtRec_from_hypertriton", HistType::kTH1F, {pt_axis});
117+
registryMC.add("he3SecPtGen_from_hypertriton", "he3SecPtGen_from_hypertriton", HistType::kTH1F, {pt_axis});
118+
registryMC.add("hyperHe4PtGen", "hyperHe4PtGen", HistType::kTH1F, {pt_axis});
119+
registryMC.add("he3SecPtRec_from_hyperHe4", "he3SecPtRec_from_hyperHe4", HistType::kTH1F, {pt_axis});
120+
registryMC.add("he3SecPtGen_from_hyperHe4", "he3SecPtGen_from_hyperHe4", HistType::kTH1F, {pt_axis});
121+
registryMC.add("he3PtRec", "he3PtRec", HistType::kTH1F, {pt_axis});
122+
registryMC.add("he3PtGen", "he3PtGen", HistType::kTH1F, {pt_axis});
95123
} 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)"}});
124+
registryMC.add("deutSecPtRec_from_hypertriton", "deutSecPtRec_from_hypertriton", HistType::kTH1F, {pt_axis});
125+
registryMC.add("deutSecPtGen_from_hypertriton", "deutSecPtGen_from_hypertriton", HistType::kTH1F, {pt_axis});
126+
registryMC.add("deutPtRec", "deutPtRec", HistType::kTH1F, {pt_axis});
127+
registryMC.add("deutPtGen", "deutPtGen", HistType::kTH1F, {pt_axis});
100128
}
101-
}
102129

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

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-
}
177+
for (const auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
178+
if (std::abs(daughter.pdgCode()) != mSelectedPDG) {
179+
continue;
180+
}
181+
const float daughterSign = daughter.pdgCode() < 0 ? -1. : 1.;
182+
183+
registryMC.fill(HIST("hypertritonPtGen"), particleSign * mcparticle.pt());
184+
if (saveHelium) {
185+
registryMC.fill(HIST("he3SecPtGen_from_hypertriton"), daughterSign * daughter.pt());
186+
} else {
187+
registryMC.fill(HIST("deutSecPtGen_from_hypertriton"), daughterSign * daughter.pt());
123188
}
124189
}
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-
}
190+
} else if (std::abs(mcparticle.pdgCode()) == HyperHelium4PDG) {
191+
for (const auto& daughter : mcparticle.daughters_as<aod::McParticles>()) {
192+
if (std::abs(daughter.pdgCode()) != mSelectedPDG) {
193+
continue;
194+
}
195+
const float daughterSign = daughter.pdgCode() < 0 ? -1. : 1.;
196+
197+
registryMC.fill(HIST("hyperHe4PtGen"), particleSign * mcparticle.pt());
198+
if (saveHelium) {
199+
registryMC.fill(HIST("he3SecPtGen_from_hyperHe4"), daughterSign * daughter.pt());
133200
}
134201
}
202+
} else if (std::abs(mcparticle.pdgCode()) == HePDG && mcparticle.isPhysicalPrimary()) {
203+
registryMC.fill(HIST("he3PtGen"), particleSign * mcparticle.pt());
204+
} else if (std::abs(mcparticle.pdgCode()) == DeuteronPDG && mcparticle.isPhysicalPrimary()) {
205+
registryMC.fill(HIST("deutPtGen"), particleSign * mcparticle.pt());
135206
}
136-
if (mcparticle.pdgCode() == AntihePDG && mcparticle.isPhysicalPrimary()) {
137-
registryMC.fill(HIST("he3PtGen"), mcparticle.pt());
207+
}
208+
209+
const o2::math_utils::Point3D<float> collisionVertex{collision.posX(), collision.posY(), collision.posZ()};
210+
211+
for (const auto& track : tracks) {
212+
213+
if (!track.has_mcParticle()) {
214+
continue;
138215
}
139-
if (mcparticle.pdgCode() == AntideuteronPDG && mcparticle.isPhysicalPrimary()) {
140-
registryMC.fill(HIST("deutPtGen"), mcparticle.pt());
216+
217+
const auto& mcparticle = track.mcParticle();
218+
if (mcparticle.pdgCode() != mSelectedPDG) {
219+
continue;
141220
}
142-
}
143-
}
144221

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-
}
222+
mDcaInfoCov.set(999, 999, 999, 999, 999);
223+
setTrackParCov(track, mTrackParCov);
224+
mTrackParCov.setPID(track.pidForTracking());
225+
std::array<float, 2> dcaInfo;
226+
o2::base::Propagator::Instance()->propagateToDCA(collisionVertex, mTrackParCov, mBz, 2.f, static_cast<o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value), &dcaInfo);
153227

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-
}
228+
const float dcaXY = dcaInfo[0];
229+
const float dcaZ = dcaInfo[1];
230+
231+
if (track.itsNCls() < min_ITS_nClusters ||
232+
track.itsNClsInnerBarrel() < min_ITS_InnerBarrel_nClusters ||
233+
track.itsNClsInnerBarrel() > max_ITS_InnerBarrel_nClusters ||
234+
track.tpcNClsFound() < min_TPC_nClusters ||
235+
track.tpcNClsCrossedRows() < min_TPC_nCrossedRows ||
236+
track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() ||
237+
track.tpcChi2NCl() > max_chi2_TPC ||
238+
track.tpcChi2NCl() < min_chi2_TPC ||
239+
track.eta() < min_eta || track.eta() > max_eta ||
240+
dcaXY > max_dcaxy || dcaXY < -max_dcaxy ||
241+
dcaZ > max_dcaz || dcaZ < -max_dcaz ||
242+
track.itsChi2NCl() > 36.f) {
243+
continue;
244+
}
245+
246+
const float particleSign = mcparticle.pdgCode() < 0 ? -1. : 1.;
174247

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());
248+
if (std::abs(mcparticle.pdgCode()) == DeuteronPDG && mcparticle.isPhysicalPrimary()) {
249+
registryMC.fill(HIST("deutPtRec"), particleSign * track.pt());
250+
}
251+
if (std::abs(mcparticle.pdgCode()) == HePDG && mcparticle.isPhysicalPrimary()) {
252+
registryMC.fill(HIST("he3PtRec"), particleSign * 2 * track.pt());
253+
}
254+
255+
for (const auto& motherparticle : mcparticle.mothers_as<aod::McParticles>()) {
256+
257+
if (std::abs(motherparticle.pdgCode()) == HypertritonPDG) {
258+
if (std::abs(mcparticle.pdgCode()) == HePDG) {
259+
registryMC.fill(HIST("he3SecPtRec_from_hypertriton"), 2 * particleSign * track.pt());
180260
} else {
181-
registryMC.fill(HIST("deutSecPtRec_from_hypertriton"), track.pt());
261+
registryMC.fill(HIST("deutSecPtRec_from_hypertriton"), particleSign * track.pt());
262+
}
263+
} else if (std::abs(motherparticle.pdgCode()) == HyperHelium4PDG) {
264+
if (std::abs(mcparticle.pdgCode()) == HePDG) {
265+
registryMC.fill(HIST("he3SecPtRec_from_hyperHe4"), 2 * particleSign * track.pt());
182266
}
183-
}
184-
if (motherparticle.pdgCode() == AntiHyperHelium4PDG) {
185-
registryMC.fill(HIST("he3SecPtRec_from_hyperHe4"), 2 * track.pt());
186267
}
187268
}
188269
}

0 commit comments

Comments
 (0)