Skip to content

[PWGCF] FU: fixing histos for prompt fraction; process to store only reco events in MC truth #17064

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
kgwizdzi wants to merge 1 commit into
base: master
Choose a base branch
from
+235 −305
Draft

[PWGCF] FU: fixing histos for prompt fraction; process to store only reco events in MC truth #17064

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
145 changes: 110 additions & 35 deletions PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerMCTruthTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@
#include "PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h"
#include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h"

#include "Common/CCDB/EventSelectionParams.h"
#include "Common/CCDB/RCTSelectionFlags.h"
#include "Common/CCDB/TriggerAliases.h"
#include "Common/Core/RecoDecay.h"
#include "Common/DataModel/EventSelection.h"
Expand Down Expand Up @@ -47,6 +49,7 @@
using namespace o2::framework;
using namespace o2::framework::expressions;
using namespace o2::constants::physics;
using namespace o2::aod::rctsel;

namespace o2::aod
{
Expand Down Expand Up @@ -76,22 +79,23 @@
}

struct FemtoUniverseProducerMCTruthTask {
int mRunNumber;
float mMagField;
Service<o2::ccdb::BasicCCDBManager> ccdb; /// Accessing the CCDB

float mMagField = 0.;
Service<o2::ccdb::BasicCCDBManager> ccdb{}; /// Accessing the CCDB

// Tables being produced
Produces<aod::FdCollisions> outputCollision;
Produces<aod::FDParticles> outputParts;
// Produces<aod::FDMCLabels> outputPartsMCLabels;
// Produces<aod::FdMCParticles> outputPartsMC;

RCTFlagsChecker rctChecker;

// Analysis configs
Configurable<bool> confIsRun3{"confIsRun3", false, "Running on Run3 or pilot"};
Configurable<bool> confIsMC{"confIsMC", false, "Running on MC; implemented only for Run3"};
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)"};
Configurable<std::vector<int>> confPDGCodes{"confPDGCodes", std::vector<int>{211, -211, 2212, -2212, 333}, "PDG of particles to be stored"};
Configurable<bool> confAnalysisWithPID{"confAnalysisWithPID", true, "1: take only particles with specified PDG, 0: all particles"};
Configurable<bool> confStoreMotherPDG{"confStoreMotherPDG", false, "Store mother's PDG in tempFitVar."};

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

// Track cuts
struct : o2::framework::ConfigurableGroup {
Expand All @@ -113,21 +118,26 @@
Configurable<float> yD0CandGenMax{"yD0CandGenMax", 0.5, "Rapidity cut for the D0/D0bar mesons"};

FemtoUniverseCollisionSelection colCuts;
FemtoUniverseCollisionSelection recoCollCuts;
FemtoUniverseTrackSelection trackCuts;
HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::QAObject};
HistogramRegistry qaRecoRegistry{"QARecoHistos", {}, OutputObjHandlingPolicy::QAObject};

void init(InitContext&)
{
if ((doprocessTrackMC) == false) {
LOGF(fatal, "Neither processFullData nor processFullMC enabled. Please choose one.");
if (!doprocessTrackMC || !doprocessTrackMcOnlyRecoColl) {
LOGF(fatal, "Neither processTrackMC nor processTrackMcOnlyRecoColl enabled. Please choose one.");
}

rctChecker.init("CBT_hadronPID", false, true);
// MC Truth collisions
colCuts.setCuts(confEvtZvtx, confEvtTriggerCheck, confEvtTriggerSel, confEvtOfflineCheck, confIsRun3, confCentFT0Min, confCentFT0Max);

colCuts.init(&qaRegistry);
// MC Reco collisions
recoCollCuts.setCuts(confEvtZvtx, confEvtTriggerCheck, confEvtTriggerSel, confEvtOfflineCheck, confIsRun3, confCentFT0Min, confCentFT0Max);
recoCollCuts.init(&qaRecoRegistry);

trackCuts.init<aod::femtouniverseparticle::ParticleType::kTrack, aod::femtouniverseparticle::TrackType::kNoChild, aod::femtouniverseparticle::CutContainerType>(&qaRegistry);

mRunNumber = 0;
mMagField = 0.0;

/// Initializing CCDB
Expand Down Expand Up @@ -159,6 +169,44 @@
}
}

template <typename CollisionType>
bool checkMcRecoCollisions(CollisionType const& col)
{
if (confIsCheckRCTFlags && !rctChecker(col)) {
return false;
}

// check whether the basic event selection criteria are fulfilled
// if the basic selection is NOT fulfilled MC Truth collision is rejected
if (!recoCollCuts.isSelected(col)) {
return false;
}
// additional checks on the reconstructed collisions
if (col.selection_bit(aod::evsel::kNoSameBunchPileup) && col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) && col.selection_bit(aod::evsel::kIsVertexITSTPC)) {
recoCollCuts.fillQA(col);
return true;
} else {
return false;
}
}

template <typename ParticleType>
int getMotherPDG(ParticleType particle)
{
if (particle.isPhysicalPrimary()) {
return 0;
} else if (particle.has_mothers()) {
if (particle.getProcess() == 20 || particle.getProcess() == 23) { // treat particles from hadronic scattering (20, 23) as primary

Check failure on line 199 in PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerMCTruthTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
return 0;
}
auto motherparticlesMC = particle.template mothers_as<aod::McParticles>();
auto motherparticleMC = motherparticlesMC.front();
return motherparticleMC.pdgCode();
} else {
return 999;
}
}

template <typename TrackType>
void fillParticles(TrackType const& tracks)
{
Expand All @@ -168,28 +216,27 @@
/// if the most open selection criteria are not fulfilled there is no
/// point looking further at the track

if (particle.pt() < ConfFilteringTracks.confPtLowFilterCut || particle.pt() > ConfFilteringTracks.confPtHighFilterCut)
if (particle.pt() < ConfFilteringTracks.confPtLowFilterCut || particle.pt() > ConfFilteringTracks.confPtHighFilterCut) {
continue;
}

int pdgCode = particle.pdgCode();

if (confAnalysisWithPID) {
bool pass = false;
std::vector<int> tmpPDGCodes = confPDGCodes; // necessary due to some features of the Configurable
for (auto const& pdg : tmpPDGCodes) {
if (pdgCode == Pdg::kPhi) { // phi meson
if (pdgCode == Pdg::kPhi || // phi meson
std::abs(pdgCode) == Pdg::kD0 || // D0(bar) meson
pdgCode == Pdg::kDPlus || // D+ meson
(pdg == pdgCode && particle.isPhysicalPrimary())) {
pass = true;
} else if (std::abs(pdgCode) == Pdg::kD0) { // D0(bar) meson
pass = true;
} else if (pdgCode == Pdg::kDPlus) { // D+ meson
pass = true;
} else if (static_cast<int>(pdg) == static_cast<int>(pdgCode)) {
if (particle.isPhysicalPrimary())
pass = true;
break; // Exit early once a match is found
}
}
if (!pass)
if (!pass) {
continue;
}
}

// check if D0/D0bar mesons pass the rapidity cut
Expand All @@ -199,15 +246,10 @@
if (std::abs(particle.pdgCode()) == Pdg::kD0) {
if (std::abs(particle.y()) > yD0CandGenMax) {
continue;
} else {
origin = RecoDecay::getCharmHadronOrigin(tracks, particle);
}
} else {
if (std::abs(particle.eta()) > ConfFilteringTracks.confEtaFilterCut) {
continue;
} else {
origin = -99;
}
origin = RecoDecay::getCharmHadronOrigin(tracks, particle);
} else if (std::abs(particle.eta()) > ConfFilteringTracks.confEtaFilterCut) {
continue;
}

/// check if we end-up with the correct final state using MC info
Expand All @@ -216,7 +258,8 @@
/// check if we have D0(bar) → π± K∓
continue;
}

// Getting the PDG code of the mother particle
int motherPDG = confStoreMotherPDG ? getMotherPDG(particle) : particle.pdgCode();
// we cannot use isSelectedMinimal since it takes Ncls
// if (!trackCuts.isSelectedMinimal(track)) {
// continue;
Expand All @@ -239,22 +282,54 @@
pdgCode,
childIDs,
origin,
-999.);
motherPDG);
}
}

void
processTrackMC(aod::McCollision const&,
soa::SmallGroups<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels>> const& collisions,
aod::McParticles const& mcParticles,
aod::BCsWithTimestamps const&)
void processTrackMC(aod::McCollision const&,
soa::SmallGroups<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels>> const& collisions,
aod::McParticles const& mcParticles,
aod::BCsWithTimestamps const&)
{
// magnetic field for run not needed for mc truth
// fill the tables
fillCollisions(collisions, mcParticles);
fillParticles(mcParticles);
}
PROCESS_SWITCH(FemtoUniverseProducerMCTruthTask, processTrackMC, "Provide MC data for track analysis", true);

Preslice<aod::McParticles> perMCCollision = aod::mcparticle::mcCollisionId;
PresliceUnsorted<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels>> recoCollsPerMCColl = aod::mcparticle::mcCollisionId;

void processTrackMcOnlyRecoColl(aod::McCollisions const& mccols,
soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels> const& collisions,
aod::McParticles const& mcParticles,
aod::BCsWithTimestamps const&)
{
// check on the reco collisions
static std::set<int> mcColIds;
mcColIds.clear();

for (const auto& col : collisions) {
const auto colcheck = checkMcRecoCollisions(col);
if (colcheck) {
mcColIds.insert(col.mcCollisionId());
}
}
// filling truth collisions and particles
for (const auto& mccol : mccols) {
if (!mcColIds.contains(mccol.globalIndex())) {
continue;
}
auto groupedMCParticles = mcParticles.sliceBy(perMCCollision, mccol.globalIndex());
auto groupedCollisions = collisions.sliceBy(recoCollsPerMCColl, mccol.globalIndex());
// magnetic field for run not needed for mc truth
// fill the tables
fillCollisions(groupedCollisions, groupedMCParticles);
fillParticles(groupedMCParticles);
}
}
PROCESS_SWITCH(FemtoUniverseProducerMCTruthTask, processTrackMcOnlyRecoColl, "Provide MC data for track analysis from truth coll which were recontructed ", false);
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
Expand Down
Loading
Loading