Skip to content

Commit 867a9d4

Browse files
kgwizdzikgwizdzi
authored
[PWGCF] FU: fixing histos for prompt fraction; process to store only reco events in MC truth (#17064)
1 parent 67499f8 commit 867a9d4

2 files changed

Lines changed: 244 additions & 308 deletions

File tree

‎PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerMCTruthTask.cxx‎

Lines changed: 113 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,8 @@
1818
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h"
1919
#include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h"
2020

21+
#include "Common/CCDB/EventSelectionParams.h"
22+
#include "Common/CCDB/RCTSelectionFlags.h"
2123
#include "Common/CCDB/TriggerAliases.h"
2224
#include "Common/Core/RecoDecay.h"
2325
#include "Common/DataModel/EventSelection.h"
@@ -40,13 +42,15 @@
4042
#include <chrono>
4143
#include <cstddef>
4244
#include <cstdint>
45+
#include <set>
4346
#include <vector>
4447

4548
using namespace o2;
4649
using namespace o2::analysis::femto_universe;
4750
using namespace o2::framework;
4851
using namespace o2::framework::expressions;
4952
using namespace o2::constants::physics;
53+
using namespace o2::aod::rctsel;
5054

5155
namespace o2::aod
5256
{
@@ -76,22 +80,26 @@ int getRowDaughters(int daughID, T const& vecID)
7680
}
7781

7882
struct FemtoUniverseProducerMCTruthTask {
79-
int mRunNumber;
80-
float mMagField;
81-
Service<o2::ccdb::BasicCCDBManager> ccdb; /// Accessing the CCDB
83+
84+
float mMagField = 0.;
85+
Service<o2::ccdb::BasicCCDBManager> ccdb{}; /// Accessing the CCDB
8286

8387
// Tables being produced
8488
Produces<aod::FdCollisions> outputCollision;
8589
Produces<aod::FDParticles> outputParts;
8690
// Produces<aod::FDMCLabels> outputPartsMCLabels;
8791
// Produces<aod::FdMCParticles> outputPartsMC;
8892

93+
RCTFlagsChecker rctChecker;
94+
95+
static constexpr int kHadElasticScatt = 20; // kPHElastic - hadronic elastic scattering
96+
static constexpr int kHadInelasticScatt = 23; // kPHInhelastic - hadronic inelastic scattering
97+
8998
// Analysis configs
9099
Configurable<bool> confIsRun3{"confIsRun3", false, "Running on Run3 or pilot"};
91-
Configurable<bool> confIsMC{"confIsMC", false, "Running on MC; implemented only for Run3"};
92-
Configurable<bool> confIsForceGRP{"confIsForceGRP", false, "Set true if the magnetic field configuration is not available in the usual CCDB directory (e.g. for Run 2 converted data or unanchorad Monte Carlo)"};
93100
Configurable<std::vector<int>> confPDGCodes{"confPDGCodes", std::vector<int>{211, -211, 2212, -2212, 333}, "PDG of particles to be stored"};
94101
Configurable<bool> confAnalysisWithPID{"confAnalysisWithPID", true, "1: take only particles with specified PDG, 0: all particles"};
102+
Configurable<bool> confStoreMotherPDG{"confStoreMotherPDG", false, "Store mother's PDG in tempFitVar."};
95103

96104
/// Event cuts
97105
Configurable<float> confEvtZvtx{"confEvtZvtx", 10.f, "Evt sel: Max. z-Vertex (cm)"};
@@ -101,6 +109,7 @@ struct FemtoUniverseProducerMCTruthTask {
101109
Configurable<float> confCentFT0Min{"confCentFT0Min", 0.f, "Min CentFT0 value for centrality selection"};
102110
Configurable<float> confCentFT0Max{"confCentFT0Max", 200.f, "Max CentFT0 value for centrality selection"};
103111
Configurable<bool> confDoSpher{"confDoSpher", false, "Calculate sphericity. If false sphericity will take value of 2."};
112+
Configurable<bool> confIsCheckRCTFlags{"confIsCheckRCTFlags", true, "Use RCTFlags"};
104113

105114
// Track cuts
106115
struct : o2::framework::ConfigurableGroup {
@@ -113,21 +122,26 @@ struct FemtoUniverseProducerMCTruthTask {
113122
Configurable<float> yD0CandGenMax{"yD0CandGenMax", 0.5, "Rapidity cut for the D0/D0bar mesons"};
114123

115124
FemtoUniverseCollisionSelection colCuts;
125+
FemtoUniverseCollisionSelection recoCollCuts;
116126
FemtoUniverseTrackSelection trackCuts;
117127
HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::QAObject};
128+
HistogramRegistry qaRecoRegistry{"QARecoHistos", {}, OutputObjHandlingPolicy::QAObject};
118129

119130
void init(InitContext&)
120131
{
121-
if ((doprocessTrackMC) == false) {
122-
LOGF(fatal, "Neither processFullData nor processFullMC enabled. Please choose one.");
132+
if (!doprocessTrackMC || !doprocessTrackMcOnlyRecoColl) {
133+
LOGF(fatal, "Neither processTrackMC nor processTrackMcOnlyRecoColl enabled. Please choose one.");
123134
}
124-
135+
rctChecker.init("CBT_hadronPID", false, true);
136+
// MC Truth collisions
125137
colCuts.setCuts(confEvtZvtx, confEvtTriggerCheck, confEvtTriggerSel, confEvtOfflineCheck, confIsRun3, confCentFT0Min, confCentFT0Max);
126-
127138
colCuts.init(&qaRegistry);
139+
// MC Reco collisions
140+
recoCollCuts.setCuts(confEvtZvtx, confEvtTriggerCheck, confEvtTriggerSel, confEvtOfflineCheck, confIsRun3, confCentFT0Min, confCentFT0Max);
141+
recoCollCuts.init(&qaRecoRegistry);
142+
128143
trackCuts.init<aod::femtouniverseparticle::ParticleType::kTrack, aod::femtouniverseparticle::TrackType::kNoChild, aod::femtouniverseparticle::CutContainerType>(&qaRegistry);
129144

130-
mRunNumber = 0;
131145
mMagField = 0.0;
132146

133147
/// Initializing CCDB
@@ -159,6 +173,43 @@ struct FemtoUniverseProducerMCTruthTask {
159173
}
160174
}
161175

176+
template <typename CollisionType>
177+
bool checkMcRecoCollisions(CollisionType const& col)
178+
{
179+
if (confIsCheckRCTFlags && !rctChecker(col)) {
180+
return false;
181+
}
182+
183+
// check whether the basic event selection criteria are fulfilled
184+
// if the basic selection is NOT fulfilled MC Truth collision is rejected
185+
if (!recoCollCuts.isSelected(col)) {
186+
return false;
187+
}
188+
// additional checks on the reconstructed collisions
189+
if (col.selection_bit(aod::evsel::kNoSameBunchPileup) && col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) && col.selection_bit(aod::evsel::kIsVertexITSTPC)) {
190+
recoCollCuts.fillQA(col);
191+
return true;
192+
}
193+
return false;
194+
}
195+
196+
template <typename ParticleType>
197+
int getMotherPDG(ParticleType const& particle)
198+
{
199+
if (particle.isPhysicalPrimary()) {
200+
return 0;
201+
}
202+
if (particle.has_mothers()) {
203+
if (particle.getProcess() == kHadElasticScatt || particle.getProcess() == kHadInelasticScatt) { // treat particles from hadronic scattering (20, 23) as primary
204+
return 0;
205+
}
206+
auto motherparticlesMC = particle.template mothers_as<aod::McParticles>();
207+
const auto& motherparticleMC = motherparticlesMC.front();
208+
return motherparticleMC.pdgCode();
209+
}
210+
return 999;
211+
}
212+
162213
template <typename TrackType>
163214
void fillParticles(TrackType const& tracks)
164215
{
@@ -168,28 +219,27 @@ struct FemtoUniverseProducerMCTruthTask {
168219
/// if the most open selection criteria are not fulfilled there is no
169220
/// point looking further at the track
170221

171-
if (particle.pt() < ConfFilteringTracks.confPtLowFilterCut || particle.pt() > ConfFilteringTracks.confPtHighFilterCut)
222+
if (particle.pt() < ConfFilteringTracks.confPtLowFilterCut || particle.pt() > ConfFilteringTracks.confPtHighFilterCut) {
172223
continue;
224+
}
173225

174226
int pdgCode = particle.pdgCode();
175227

176228
if (confAnalysisWithPID) {
177229
bool pass = false;
178230
std::vector<int> tmpPDGCodes = confPDGCodes; // necessary due to some features of the Configurable
179231
for (auto const& pdg : tmpPDGCodes) {
180-
if (pdgCode == Pdg::kPhi) { // phi meson
181-
pass = true;
182-
} else if (std::abs(pdgCode) == Pdg::kD0) { // D0(bar) meson
183-
pass = true;
184-
} else if (pdgCode == Pdg::kDPlus) { // D+ meson
232+
if (pdgCode == Pdg::kPhi || // phi meson
233+
std::abs(pdgCode) == Pdg::kD0 || // D0(bar) meson
234+
pdgCode == Pdg::kDPlus || // D+ meson
235+
(pdg == pdgCode && particle.isPhysicalPrimary())) {
185236
pass = true;
186-
} else if (static_cast<int>(pdg) == static_cast<int>(pdgCode)) {
187-
if (particle.isPhysicalPrimary())
188-
pass = true;
237+
break; // Exit early once a match is found
189238
}
190239
}
191-
if (!pass)
240+
if (!pass) {
192241
continue;
242+
}
193243
}
194244

195245
// check if D0/D0bar mesons pass the rapidity cut
@@ -199,15 +249,10 @@ struct FemtoUniverseProducerMCTruthTask {
199249
if (std::abs(particle.pdgCode()) == Pdg::kD0) {
200250
if (std::abs(particle.y()) > yD0CandGenMax) {
201251
continue;
202-
} else {
203-
origin = RecoDecay::getCharmHadronOrigin(tracks, particle);
204-
}
205-
} else {
206-
if (std::abs(particle.eta()) > ConfFilteringTracks.confEtaFilterCut) {
207-
continue;
208-
} else {
209-
origin = -99;
210252
}
253+
origin = RecoDecay::getCharmHadronOrigin(tracks, particle);
254+
} else if (std::abs(particle.eta()) > ConfFilteringTracks.confEtaFilterCut) {
255+
continue;
211256
}
212257

213258
/// check if we end-up with the correct final state using MC info
@@ -216,7 +261,8 @@ struct FemtoUniverseProducerMCTruthTask {
216261
/// check if we have D0(bar) → π± K∓
217262
continue;
218263
}
219-
264+
// Getting the PDG code of the mother particle
265+
int motherPDG = confStoreMotherPDG ? getMotherPDG(particle) : particle.pdgCode();
220266
// we cannot use isSelectedMinimal since it takes Ncls
221267
// if (!trackCuts.isSelectedMinimal(track)) {
222268
// continue;
@@ -239,22 +285,54 @@ struct FemtoUniverseProducerMCTruthTask {
239285
pdgCode,
240286
childIDs,
241287
origin,
242-
-999.);
288+
motherPDG);
243289
}
244290
}
245291

246-
void
247-
processTrackMC(aod::McCollision const&,
248-
soa::SmallGroups<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels>> const& collisions,
249-
aod::McParticles const& mcParticles,
250-
aod::BCsWithTimestamps const&)
292+
void processTrackMC(aod::McCollision const&,
293+
soa::SmallGroups<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels>> const& collisions,
294+
aod::McParticles const& mcParticles,
295+
aod::BCsWithTimestamps const&)
251296
{
252297
// magnetic field for run not needed for mc truth
253298
// fill the tables
254299
fillCollisions(collisions, mcParticles);
255300
fillParticles(mcParticles);
256301
}
257302
PROCESS_SWITCH(FemtoUniverseProducerMCTruthTask, processTrackMC, "Provide MC data for track analysis", true);
303+
304+
Preslice<aod::McParticles> perMCCollision = aod::mcparticle::mcCollisionId;
305+
PresliceUnsorted<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels>> recoCollsPerMCColl = aod::mcparticle::mcCollisionId;
306+
307+
void processTrackMcOnlyRecoColl(aod::McCollisions const& mccols,
308+
soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels> const& collisions,
309+
aod::McParticles const& mcParticles,
310+
aod::BCsWithTimestamps const&)
311+
{
312+
// check on the reco collisions
313+
static std::set<int> mcColIds;
314+
mcColIds.clear();
315+
316+
for (const auto& col : collisions) {
317+
const auto colcheck = checkMcRecoCollisions(col);
318+
if (colcheck) {
319+
mcColIds.insert(col.mcCollisionId());
320+
}
321+
}
322+
// filling truth collisions and particles
323+
for (const auto& mccol : mccols) {
324+
if (!mcColIds.contains(mccol.globalIndex())) {
325+
continue;
326+
}
327+
auto groupedMCParticles = mcParticles.sliceBy(perMCCollision, mccol.globalIndex());
328+
auto groupedCollisions = collisions.sliceBy(recoCollsPerMCColl, mccol.globalIndex());
329+
// magnetic field for run not needed for mc truth
330+
// fill the tables
331+
fillCollisions(groupedCollisions, groupedMCParticles);
332+
fillParticles(groupedMCParticles);
333+
}
334+
}
335+
PROCESS_SWITCH(FemtoUniverseProducerMCTruthTask, processTrackMcOnlyRecoColl, "Provide MC data for track analysis from truth coll which were recontructed ", false);
258336
};
259337

260338
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)