diff --git a/PWGLF/TableProducer/Nuspex/deuteronInTriggeredEvents.cxx b/PWGLF/TableProducer/Nuspex/deuteronInTriggeredEvents.cxx index 9d855c71fda..22cca32068a 100644 --- a/PWGLF/TableProducer/Nuspex/deuteronInTriggeredEvents.cxx +++ b/PWGLF/TableProducer/Nuspex/deuteronInTriggeredEvents.cxx @@ -266,6 +266,7 @@ struct DeuteronInTriggeredEvents { Produces nucleiTable; // For data Produces nucleiTableMCExtension; // For MC analysis Produces genEventMCSel; // For MC reco events + Produces nucleiTableMC; // For MC analysis as data Service ccdb; Service pdgDB; // For INELgt0 gen MC selection Zorro zorro; // Definition of Zorro: helpful for skimmed data @@ -1111,6 +1112,77 @@ struct DeuteronInTriggeredEvents { } } PROCESS_SWITCH(DeuteronInTriggeredEvents, processMC, "MC analysis", false); + + void processMCasData(soa::Join const& collisions, aod::McCollisions const& mcCollisions, soa::Join const& tracks, aod::McParticles const& particlesMC, aod::BCsWithTimestamps const&) + { + nuclei::candidates.clear(); + std::vector goodCollisions(mcCollisions.size(), false); + + for (const auto& collision : collisions) { + if (!eventSelectionWithHisto(collision)) { + continue; + } + + // Avoid unwanted memory leaks + if (!collision.has_mcCollision()) + continue; + + int mcId = collision.mcCollisionId(); + if (mcId < 0 || mcId >= static_cast(mcCollisions.size())) + continue; + + goodCollisions[collision.mcCollisionId()] = true; + const auto& slicedTracks = tracks.sliceBy(tracksPerCollisions, collision.globalIndex()); + fillDataInfo(collision, slicedTracks); + } + + std::vector isReconstructed(particlesMC.size(), false); + for (size_t i{0}; i < nuclei::candidates.size(); ++i) { + auto& c = nuclei::candidates[i]; + if (c.fillTree) { + auto label = tracks.iteratorAt(c.globalIndex); + + if (label.mcParticleId() < -1 || label.mcParticleId() >= particlesMC.size()) { + continue; + } + + auto particle = particlesMC.iteratorAt(label.mcParticleId()); + + int motherPdgCode = 0; + float motherDecRadius = -1; + isReconstructed[particle.globalIndex()] = true; + + if (particle.isPhysicalPrimary()) { + c.flags |= kIsPhysicalPrimary; + if (particle.has_mothers()) { + for (const auto& motherparticle : particle.mothers_as()) { + if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) { + c.flags |= kIsSecondaryFromWeakDecay; + motherPdgCode = motherparticle.pdgCode(); + motherDecRadius = std::hypot(particle.vx() - motherparticle.vx(), particle.vy() - motherparticle.vy()); + break; + } + } + } + } else if (particle.getProcess() == TMCProcess::kPDecay) { + c.flags |= kIsSecondaryFromWeakDecay; + for (const auto& motherparticle : particle.mothers_as()) { + motherPdgCode = motherparticle.pdgCode(); + motherDecRadius = std::hypot(particle.vx() - motherparticle.vx(), particle.vy() - motherparticle.vy()); + } + } else { + c.flags |= kIsSecondaryFromMaterial; + } + + isReconstructed[particle.globalIndex()] = true; + float absoDecL = computeAbsoDecL(particle); + + nucleiTableMC(c.pt, c.eta, c.phi, c.tpcInnerParam, c.beta, c.zVertex, c.nContrib, c.dcaXY, c.dcaZ, c.tpcSignal, c.itsChi2, c.tpcChi2, c.tofChi2, c.flags, c.tpcFindableCls, c.tpcCrossedRows, c.itsClsMap, c.tpcNCls, c.tpcNClsShared, c.clusterSizesITS, goodCollisions[particle.mcCollisionId()], particle.pt(), particle.eta(), particle.phi(), particle.pdgCode(), motherPdgCode, motherDecRadius, absoDecL); + } + } + } + + PROCESS_SWITCH(DeuteronInTriggeredEvents, processMCasData, "MC as data analysis", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)