Skip to content

Commit

Permalink
PWGEM/PhotonMeson: update dalitz for MC (AliceO2Group#3541)
Browse files Browse the repository at this point in the history
  • Loading branch information
@dsekihat
dsekihat authored Oct 3, 2023
1 parent 8f2e52c commit 3bd0163
Show file tree
Hide file tree
Showing 9 changed files with 472 additions and 106 deletions.
11 changes: 10 additions & 1 deletion PWGEM/PhotonMeson/DataModel/gammaTables.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,9 +133,18 @@ DECLARE_SOA_COLUMN(McMask, mcMask, uint16_t);
// NOTE: MC labels. This table has one entry for each reconstructed track (joinable with the track tables)
DECLARE_SOA_TABLE(EMMCParticleLabels, "AOD", "EMMCPARLABEL", //!
emmcparticlelabel::EMMCParticleId, emmcparticlelabel::McMask);

using EMMCParticleLabel = EMMCParticleLabels::iterator;

namespace emprimarytrackmclabel
{
DECLARE_SOA_INDEX_COLUMN(EMMCParticle, emmcparticle); //!
DECLARE_SOA_COLUMN(McMask, mcMask, uint16_t);
} // namespace emprimarytrackmclabel

DECLARE_SOA_TABLE(EMPrimaryTrackMCLabels, "AOD", "EMPRMTRKMCLABEL", //!
emprimarytrackmclabel::EMMCParticleId, emprimarytrackmclabel::McMask);
using EMPrimaryTrackMCLabel = EMPrimaryTrackMCLabels::iterator;

namespace v0leg
{
DECLARE_SOA_INDEX_COLUMN(Collision, collision); //!
Expand Down
84 changes: 45 additions & 39 deletions PWGEM/PhotonMeson/TableProducer/createEMReducedMCEvent.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,9 @@ struct createEMReducedMCEvent {
Produces<o2::aod::EMReducedMCEventLabels> mceventlabels;
Produces<o2::aod::EMMCParticles> emmcparticles;
Produces<o2::aod::EMMCParticleLabels> emmcparticlelabels;
Produces<o2::aod::EMPrimaryTrackMCLabels> emprimarytrackmclabels;
Produces<o2::aod::V0KFEMReducedEventIds> v0kfeventid;
Produces<o2::aod::DalitzEEEMReducedEventIds> dalitzeventid;
Produces<o2::aod::PHOSEMReducedEventIds> phoseventid;
Produces<o2::aod::EMCEMReducedEventIds> emceventid;

Expand All @@ -59,13 +61,15 @@ struct createEMReducedMCEvent {

Preslice<aod::McParticles> perMcCollision = aod::mcparticle::mcCollisionId;
Preslice<aod::V0Photons> perCollision_pcm = aod::v0photon::collisionId;
Preslice<aod::DalitzEEs> perCollision_dalitz = aod::dalitzee::collisionId;
Preslice<aod::EMPrimaryTracks> perCollision_emprmtrk = aod::emprimarytrack::collisionId;
Preslice<aod::PHOSClusters> perCollision_phos = aod::skimmedcluster::collisionId;
Preslice<aod::SkimEMCClusters> perCollision_emc = aod::skimmedcluster::collisionId;

using MyCollisions = soa::Join<aod::Collisions, aod::Mults, aod::EvSels, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentNTPVs, aod::McCollisionLabels>;

template <uint8_t system, typename TTracks, typename TPCMs, typename TPCMLegs, typename TPHOSs, typename TEMCs>
void skimmingMC(MyCollisions const& collisions, aod::BCs const&, aod::McCollisions const&, aod::McParticles const& mcTracks, TTracks const&, TPCMs const& v0photons, TPCMLegs const& v0legs, TPHOSs const& phosclusters, TEMCs const& emcclusters)
template <uint8_t system, typename TTracks, typename TPCMs, typename TPCMLegs, typename TPHOSs, typename TEMCs, typename TDielectorns, typename TEMPrimaryTracks>
void skimmingMC(MyCollisions const& collisions, aod::BCs const&, aod::McCollisions const&, aod::McParticles const& mcTracks, TTracks const& o2tracks, TPCMs const& v0photons, TPCMLegs const& v0legs, TPHOSs const& phosclusters, TEMCs const& emcclusters, TDielectorns const& dileptons, TEMPrimaryTracks const& emprimarytracks)
{
// temporary variables used for the indexing of the skimmed MC stack
std::map<uint64_t, int> fNewLabels;
Expand All @@ -80,6 +84,7 @@ struct createEMReducedMCEvent {
int ng_pcm = 0;
int ng_phos = 0;
int ng_emc = 0;
int ng_dilepton = 0;

// TODO: investigate the collisions without corresponding mcCollision
if (!collision.has_mcCollision()) {
Expand All @@ -98,6 +103,12 @@ struct createEMReducedMCEvent {
for (int iv0 = 0; iv0 < v0photons_coll.size(); iv0++) {
v0kfeventid(events.lastIndex() + 1);
}

auto dileptons_coll = dileptons.sliceBy(perCollision_dalitz, collision.globalIndex());
ng_dilepton = v0photons_coll.size();
for (int iee = 0; iee < dileptons_coll.size(); iee++) {
dalitzeventid(events.lastIndex() + 1);
}
}
if constexpr (static_cast<bool>(system & kPHOS)) {
auto phos_coll = phosclusters.sliceBy(perCollision_phos, collision.globalIndex());
Expand Down Expand Up @@ -191,9 +202,6 @@ struct createEMReducedMCEvent {

for (auto& leg : {pos, ele}) { // be carefull of order {pos, ele}!
auto o2track = leg.template track_as<TracksMC>();
// if (!o2track.has_mcParticle()) {
// continue; // If no MC particle is found, skip the track
// }
auto mctrack = o2track.template mcParticle_as<aod::McParticles>();

// if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
Expand All @@ -207,6 +215,26 @@ struct createEMReducedMCEvent {
emmcparticlelabels(fNewLabels.find(mctrack.index())->second, o2track.mcMask());
} // end of leg loop
} // end of v0 loop

// for dalitz ee
auto emprimarytracks_coll = emprimarytracks.sliceBy(perCollision_emprmtrk, collision.globalIndex());
for (auto& emprimarytrack : emprimarytracks_coll) {
auto o2track = o2tracks.iteratorAt(emprimarytrack.trackId());
if (!o2track.has_mcParticle()) {
continue; // If no MC particle is found, skip the dilepton
}
auto mctrack = o2track.template mcParticle_as<aod::McParticles>();

// if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
if (!(fNewLabels.find(mctrack.globalIndex()) != fNewLabels.end())) {
fNewLabels[mctrack.globalIndex()] = fCounters[0];
fNewLabelsReversed[fCounters[0]] = mctrack.globalIndex();
// fMCFlags[mctrack.globalIndex()] = mcflags;
fEventIdx[mctrack.globalIndex()] = fEventLabels.find(mcCollision.globalIndex())->second;
fCounters[0]++;
}
emprimarytrackmclabels(fNewLabels.find(mctrack.index())->second, o2track.mcMask());
} // end of em primary track loop
}

} // end of collision loop
Expand All @@ -218,28 +246,6 @@ struct createEMReducedMCEvent {

std::vector<int> mothers;
if (mctrack.has_mothers()) {
////auto mp = mctrack.template mothers_first_as<TMCs>();
// int motherid = mctrack.mothersIds()[0];//first mother index
// while(motherid > -1) {
// if (motherid < mcTracks.size()) { // protect against bad mother indices. why is this needed?
// auto mp = mcTracks.iteratorAt(motherid);

// if (fNewLabels.find(mp.globalIndex()) != fNewLabels.end()) {
// mothers.push_back(fNewLabels.find(mp.globalIndex())->second);
// }

// if(mp.has_mothers()){
// motherid = mp.mothersIds()[0];//first mother index
// } else {
// motherid = -999;
// }
// }
// else{
// std::cout << "Mother label (" << motherid << ") exceeds the McParticles size (" << mcTracks.size() << ")" << std::endl;
// std::cout << " Check the MC generator" << std::endl;
// }
//}

for (auto& m : mctrack.mothersIds()) {
if (m < mcTracks.size()) { // protect against bad mother indices
if (fNewLabels.find(m) != fNewLabels.end()) {
Expand Down Expand Up @@ -289,40 +295,40 @@ struct createEMReducedMCEvent {
fCounters[1] = 0;
} // end of skimmingMC

void processMC_PCM(soa::SmallGroups<MyCollisions> const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs)
void processMC_PCM(soa::SmallGroups<MyCollisions> const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& o2tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::DalitzEEs const& dileptons, aod::EMPrimaryTracks const& emprimarytracks)
{
skimmingMC<kPCM>(collisions, bcs, mccollisions, mcTracks, tracks, v0photons, v0legs, nullptr, nullptr);
skimmingMC<kPCM>(collisions, bcs, mccollisions, mcTracks, o2tracks, v0photons, v0legs, nullptr, nullptr, dileptons, emprimarytracks);
}
void processMC_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, aod::PHOSClusters const& phosclusters)
{
skimmingMC<kPHOS>(collisions, bcs, mccollisions, mcTracks, nullptr, nullptr, nullptr, phosclusters, nullptr);
skimmingMC<kPHOS>(collisions, bcs, mccollisions, mcTracks, nullptr, nullptr, nullptr, phosclusters, nullptr, nullptr, nullptr);
}
void processMC_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, aod::SkimEMCClusters const& emcclusters)
{
skimmingMC<kEMC>(collisions, bcs, mccollisions, mcTracks, nullptr, nullptr, nullptr, nullptr, emcclusters);
skimmingMC<kEMC>(collisions, bcs, mccollisions, mcTracks, nullptr, nullptr, nullptr, nullptr, emcclusters, nullptr, nullptr);
}
void processMC_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::PHOSClusters const& phosclusters)
void processMC_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& o2tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::PHOSClusters const& phosclusters, aod::DalitzEEs const& dileptons, aod::EMPrimaryTracks const& emprimarytracks)
{
const uint8_t sysflag = kPCM | kPHOS;
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, tracks, v0photons, v0legs, phosclusters, nullptr);
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, o2tracks, v0photons, v0legs, phosclusters, nullptr, dileptons, emprimarytracks);
}
void processMC_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::SkimEMCClusters const& emcclusters)
void processMC_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& o2tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::SkimEMCClusters const& emcclusters, aod::DalitzEEs const& dileptons, aod::EMPrimaryTracks const& emprimarytracks)
{
const uint8_t sysflag = kPCM | kEMC;
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, tracks, v0photons, v0legs, nullptr, emcclusters);
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, o2tracks, v0photons, v0legs, nullptr, emcclusters, dileptons, emprimarytracks);
}
void processMC_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
{
const uint8_t sysflag = kPHOS | kEMC;
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, nullptr, nullptr, nullptr, phosclusters, emcclusters);
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, nullptr, nullptr, nullptr, phosclusters, emcclusters, nullptr, nullptr);
}
void processMC_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
void processMC_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& o2tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters, aod::DalitzEEs const& dileptons, aod::EMPrimaryTracks const& emprimarytracks)
{
const uint8_t sysflag = kPCM | kPHOS | kEMC;
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, tracks, v0photons, v0legs, phosclusters, emcclusters);
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, o2tracks, v0photons, v0legs, phosclusters, emcclusters, dileptons, emprimarytracks);
}

void processDummy(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs) {}
void processDummy(MyCollisions const& collisions) {}

PROCESS_SWITCH(createEMReducedMCEvent, processMC_PCM, "create em mc event table for PCM", false);
PROCESS_SWITCH(createEMReducedMCEvent, processMC_PHOS, "create em mc event table for PHOS", false);
Expand Down
103 changes: 64 additions & 39 deletions PWGEM/PhotonMeson/TableProducer/skimmerDalitzEE.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,9 +127,15 @@ struct skimmerDalitzEE {
mRunNumber = bc.runNumber();
}

template <typename TTrack>
template <bool isMC, typename TTrack>
bool checkTrack(TTrack const& track)
{
if constexpr (isMC) {
if (!track.has_mcParticle()) {
return false;
}
}

if (!track.hasITS() || !track.hasTPC()) {
return false;
}
Expand Down Expand Up @@ -159,12 +165,12 @@ struct skimmerDalitzEE {
return true;
}

template <EM_EEPairType pairtype, typename TCollision, typename TTracks1, typename TTracks2>
template <bool isMC, EM_EEPairType pairtype, typename TCollision, typename TTracks1, typename TTracks2>
void fillPairTable(TCollision const& collision, TTracks1 const& tracks1, TTracks2 const& tracks2)
{
if constexpr (pairtype == EM_EEPairType::kULS) { // ULS
for (auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(tracks1, tracks2))) {
if (!checkTrack(t1) || !checkTrack(t2)) {
if (!checkTrack<isMC>(t1) || !checkTrack<isMC>(t2)) {
continue;
}
fRegistry.fill(HIST("hNpairs"), static_cast<int>(pairtype));
Expand All @@ -184,7 +190,7 @@ struct skimmerDalitzEE {
} // end of pairing loop
} else { // LS
for (auto& [t1, t2] : combinations(CombinationsStrictlyUpperIndexPolicy(tracks1, tracks2))) {
if (!checkTrack(t1) || !checkTrack(t2)) {
if (!checkTrack<isMC>(t1) || !checkTrack<isMC>(t2)) {
continue;
}
fRegistry.fill(HIST("hNpairs"), static_cast<int>(pairtype));
Expand All @@ -204,11 +210,11 @@ struct skimmerDalitzEE {
}
}

template <typename TTracks>
template <bool isMC, typename TTracks>
void fillTrackTable(TTracks const& tracks)
{
for (auto& track : tracks) {
if (!checkTrack(track)) {
if (!checkTrack<isMC>(track)) {
continue;
}

Expand Down Expand Up @@ -249,13 +255,13 @@ struct skimmerDalitzEE {
// store track info which belongs to pairs. (i.e. only 1 track per event does not enter pair analysis.)
int npos = 0, nneg = 0;
for (auto& ptrack : posTracks_per_coll) {
if (!checkTrack(ptrack)) {
if (!checkTrack<false>(ptrack)) {
continue;
}
npos++;
}
for (auto& ntrack : negTracks_per_coll) {
if (!checkTrack(ntrack)) {
if (!checkTrack<false>(ntrack)) {
continue;
}
nneg++;
Expand All @@ -265,45 +271,64 @@ struct skimmerDalitzEE {
continue;
}

fillTrackTable(posTracks_per_coll);
fillTrackTable(negTracks_per_coll);
fillTrackTable<false>(posTracks_per_coll);
fillTrackTable<false>(negTracks_per_coll);

fillPairTable<EM_EEPairType::kULS>(collision, posTracks_per_coll, negTracks_per_coll); // ULS
fillPairTable<EM_EEPairType::kLSpp>(collision, posTracks_per_coll, posTracks_per_coll); // LS++
fillPairTable<EM_EEPairType::kLSnn>(collision, negTracks_per_coll, negTracks_per_coll); // LS--
fillPairTable<false, EM_EEPairType::kULS>(collision, posTracks_per_coll, negTracks_per_coll); // ULS
fillPairTable<false, EM_EEPairType::kLSpp>(collision, posTracks_per_coll, posTracks_per_coll); // LS++
fillPairTable<false, EM_EEPairType::kLSnn>(collision, negTracks_per_coll, negTracks_per_coll); // LS--

} // end of collision loop
fNewLabels.clear();
fCounter = 0;
}
PROCESS_SWITCH(skimmerDalitzEE, processRec, "process reconstructed info only", true);

// void processMC(soa::Join<aod::McCollisionLabels, aod::Collisions> const& collisions,
// aod::McCollisions const&,
// aod::BCsWithTimestamps const& bcs,
// aod::V0Datas const& theV0s,
// MyTracksMC const& theTracks)
// {
// for (auto& collision : collisions) {
//
// if (!collision.has_mcCollision()) {
// continue;
// }
// auto mcCollision = collision.mcCollision();
// auto bc = collision.bc_as<aod::BCsWithTimestamps>();
//
// auto lGroupedV0s = theV0s.sliceBy(perCollision, collision.globalIndex());
// for (auto& v0 : lGroupedV0s) {
// auto pos = v0.template posTrack_as<MyTracksMC>(); // positive daughter
// auto ele = v0.template negTrack_as<MyTracksMC>(); // negative daughter
//
// if (!ele.has_mcParticle() || !pos.has_mcParticle()) {
// continue; // If no MC particle is found, skip the v0
// }
// }
// }
// }
// PROCESS_SWITCH(skimmerDalitzEE, processMC, "process reconstructed and MC info ", false);
using MyFilteredTracksMC = soa::Filtered<MyTracksMC>;
Partition<MyFilteredTracksMC> posTracksMC = o2::aod::track::signed1Pt > 0.f;
Partition<MyFilteredTracksMC> negTracksMC = o2::aod::track::signed1Pt < 0.f;
void processMC(soa::Join<aod::McCollisionLabels, aod::Collisions> const& collisions, aod::McCollisions const&, aod::BCsWithTimestamps const&, MyFilteredTracksMC const& tracks)
{
for (auto& collision : collisions) {
if (!collision.has_mcCollision()) {
continue;
}
// auto mcCollision = collision.mcCollision();
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
initCCDB(bc);
auto posTracks_per_coll = posTracksMC->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache);
auto negTracks_per_coll = negTracksMC->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache);

// store track info which belongs to pairs. (i.e. only 1 track per event does not enter pair analysis.)
int npos = 0, nneg = 0;
for (auto& ptrack : posTracks_per_coll) {
if (!checkTrack<true>(ptrack)) {
continue;
}
npos++;
}
for (auto& ntrack : negTracks_per_coll) {
if (!checkTrack<true>(ntrack)) {
continue;
}
nneg++;
}

if (npos + nneg < 2) {
continue;
}

fillTrackTable<true>(posTracks_per_coll);
fillTrackTable<true>(negTracks_per_coll);

fillPairTable<true, EM_EEPairType::kULS>(collision, posTracks_per_coll, negTracks_per_coll); // ULS
fillPairTable<true, EM_EEPairType::kLSpp>(collision, posTracks_per_coll, posTracks_per_coll); // LS++
fillPairTable<true, EM_EEPairType::kLSnn>(collision, negTracks_per_coll, negTracks_per_coll); // LS--
} // end of collision loop
fNewLabels.clear();
fCounter = 0;
}
PROCESS_SWITCH(skimmerDalitzEE, processMC, "process reconstructed and MC info ", false);
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
Expand Down
4 changes: 4 additions & 0 deletions PWGEM/PhotonMeson/Tasks/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,10 @@ o2physics_add_dpl_workflow(dalitz-qc
PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore
COMPONENT_NAME Analysis)

o2physics_add_dpl_workflow(dalitz-qc-mc
SOURCES dalitzQCMC.cxx
PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore
COMPONENT_NAME Analysis)

o2physics_add_dpl_workflow(phos-qc
SOURCES phosQC.cxx
Expand Down
Loading

0 comments on commit 3bd0163

Please sign in to comment.