diff --git a/PWGUD/Tasks/upcRhoPrimeAnalysis.cxx b/PWGUD/Tasks/upcRhoPrimeAnalysis.cxx index 703f0989283..9878d96f54f 100644 --- a/PWGUD/Tasks/upcRhoPrimeAnalysis.cxx +++ b/PWGUD/Tasks/upcRhoPrimeAnalysis.cxx @@ -9,10 +9,10 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. /// -/// \brief Task for analysis of rho' in UPCs using UD tables (from SG producer). -/// \author Cesar Ramirez, cesar.ramirez@cern.ch +/// \file upcRhoPrimeAnalysis.cxx +/// \brief Task for analysis of rho prime in UPCs using UD tables (from SG producer). +/// \author Cesar Omar Ramirez Alvarez (cesar.ramirez@cern.ch), Autonomous University of Puebla -#include "PWGUD/Core/UPCTauCentralBarrelHelperRL.h" #include "PWGUD/DataModel/UDTables.h" #include @@ -24,64 +24,69 @@ #include #include #include +#include #include -#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include // IWYU pragma: keep #include #include +#include #include #include #include -#include -#include #include using namespace o2; using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace ROOT::Math; // Define UD tables using UDtracks = soa::Join; -using UDCollisions = soa::Join; +using UDCollisions = soa::Join; +using UDMcCollisions = aod::UDMcCollisions; +using UDMcParticles = aod::UDMcParticles; +using UDtracksMC = soa::Join; +using UDCollisionsMC = soa::Join; namespace o2::aod { namespace fourpi { -// Declare columns -DECLARE_SOA_COLUMN(RunNumber, runNumber, int32_t); // Run number for event identification -DECLARE_SOA_COLUMN(M, m, double); // Invariant mass of the system -DECLARE_SOA_COLUMN(Pt, pt, double); // Transverse momentum of the system -DECLARE_SOA_COLUMN(Eta, eta, double); // Pseudorapidity of the system -DECLARE_SOA_COLUMN(Phi, phi, double); // Azimuthal angle of the system +// Columns RD / MC reco +DECLARE_SOA_COLUMN(RunNumber, runNumber, int32_t); // Run number +DECLARE_SOA_COLUMN(M, m, double); // System invariant mass +DECLARE_SOA_COLUMN(Pt, pt, double); // System pT +DECLARE_SOA_COLUMN(Eta, eta, double); // System pseudorapidity +DECLARE_SOA_COLUMN(Phi, phi, double); // System azimuthal angle DECLARE_SOA_COLUMN(PosX, posX, double); // Vertex X position DECLARE_SOA_COLUMN(PosY, posY, double); // Vertex Y position DECLARE_SOA_COLUMN(PosZ, posZ, double); // Vertex Z position -DECLARE_SOA_COLUMN(TotalCharge, totalCharge, int); // Total charge of selected tracks +DECLARE_SOA_COLUMN(TotalCharge, totalCharge, int); // Real total charge of the 4 selected tracks DECLARE_SOA_COLUMN(TotalFT0AmplitudeA, totalFT0AmplitudeA, float); // FT0A amplitude DECLARE_SOA_COLUMN(TotalFT0AmplitudeC, totalFT0AmplitudeC, float); // FT0C amplitude DECLARE_SOA_COLUMN(TotalFV0AmplitudeA, totalFV0AmplitudeA, float); // FV0A amplitude -DECLARE_SOA_COLUMN(NumContrib, numContrib, int32_t); // Number of primary vertex contributors +DECLARE_SOA_COLUMN(NumContrib, numContrib, int32_t); // Number of PV contributors DECLARE_SOA_COLUMN(Sign, sign, std::vector); // Track charges -DECLARE_SOA_COLUMN(TrackPt, trackPt, std::vector); // Track pT values -DECLARE_SOA_COLUMN(TrackEta, trackEta, std::vector); // Track eta values -DECLARE_SOA_COLUMN(TrackPhi, trackPhi, std::vector); // Track phi values -DECLARE_SOA_COLUMN(TPCNSigmaEl, tpcNSigmaEl, std::vector); // TPC nσ for electrons -DECLARE_SOA_COLUMN(TPCNSigmaPi, tpcNSigmaPi, std::vector); // TPC nσ for pions -DECLARE_SOA_COLUMN(TPCNSigmaKa, tpcNSigmaKa, std::vector); // TPC nσ for kaons -DECLARE_SOA_COLUMN(TPCNSigmaPr, tpcNSigmaPr, std::vector); // TPC nσ for protons -DECLARE_SOA_COLUMN(TrackID, trackID, std::vector); // Track identifiers -DECLARE_SOA_COLUMN(IsReconstructedWithUPC, isReconstructedWithUPC, bool); // UPC mode reconstruction flag -DECLARE_SOA_COLUMN(TimeZNA, timeZNA, float); // ZNA timing -DECLARE_SOA_COLUMN(TimeZNC, timeZNC, float); // ZNC timing +DECLARE_SOA_COLUMN(TrackPt, trackPt, std::vector); // Track pT +DECLARE_SOA_COLUMN(TrackEta, trackEta, std::vector); // Track eta +DECLARE_SOA_COLUMN(TrackPhi, trackPhi, std::vector); // Track phi +DECLARE_SOA_COLUMN(TPCNSigmaEl, tpcNSigmaEl, std::vector); // TPC nSigma electron +DECLARE_SOA_COLUMN(TPCNSigmaPi, tpcNSigmaPi, std::vector); // TPC nSigma pion +DECLARE_SOA_COLUMN(TPCNSigmaKa, tpcNSigmaKa, std::vector); // TPC nSigma kaon +DECLARE_SOA_COLUMN(TPCNSigmaPr, tpcNSigmaPr, std::vector); // TPC nSigma proton +DECLARE_SOA_COLUMN(TrackID, trackID, std::vector); // Track index within the system +DECLARE_SOA_COLUMN(IsReconstructedWithUPC, isReconstructedWithUPC, bool); // UPC reconstruction mode flag +DECLARE_SOA_COLUMN(TimeZNA, timeZNA, float); // ZNA time +DECLARE_SOA_COLUMN(TimeZNC, timeZNC, float); // ZNC time DECLARE_SOA_COLUMN(EnergyCommonZNA, energyCommonZNA, float); // ZNA energy DECLARE_SOA_COLUMN(EnergyCommonZNC, energyCommonZNC, float); // ZNC energy -DECLARE_SOA_COLUMN(IsChargeZero, isChargeZero, bool); // Neutral system flag -DECLARE_SOA_COLUMN(OccupancyInTime, occupancyInTime, int); // Occupancy in time +DECLARE_SOA_COLUMN(IsChargeZero, isChargeZero, bool); // TotalCharge == 0 +DECLARE_SOA_COLUMN(OccupancyInTime, occupancyInTime, int); // Occupancy DECLARE_SOA_COLUMN(HadronicRate, hadronicRate, double); // Hadronic interaction rate } // namespace fourpi -// Define the output DECLARE_SOA_TABLE(SYSTEMTREE, "AOD", "SystemTree", fourpi::RunNumber, fourpi::M, fourpi::Pt, fourpi::Eta, fourpi::Phi, fourpi::PosX, fourpi::PosY, fourpi::PosZ, fourpi::TotalCharge, @@ -92,10 +97,63 @@ DECLARE_SOA_TABLE(SYSTEMTREE, "AOD", "SystemTree", fourpi::TrackID, fourpi::IsReconstructedWithUPC, fourpi::TimeZNA, fourpi::TimeZNC, fourpi::EnergyCommonZNA, fourpi::EnergyCommonZNC, fourpi::IsChargeZero, fourpi::OccupancyInTime, fourpi::HadronicRate); + +namespace mcgen4pi +{ +// Columns MC gen +DECLARE_SOA_COLUMN(McMotherPdg, mcMotherPdg, int); // Common mother PDG +DECLARE_SOA_COLUMN(McMotherPt, mcMotherPt, float); // Mother pT +DECLARE_SOA_COLUMN(McMotherPhi, mcMotherPhi, float); // Mother phi +DECLARE_SOA_COLUMN(McMotherMass, mcMotherMass, float); // Mother invariant mass +DECLARE_SOA_COLUMN(McMotherRapidity, mcMotherRapidity, float); // Mother rapidity +DECLARE_SOA_COLUMN(McTotalCharge, mcTotalCharge, int); // Real total charge of the 4 generated particles + +// Per-particle info +DECLARE_SOA_COLUMN(McTrackPdg, mcTrackPdg, int[4]); +DECLARE_SOA_COLUMN(McTrackPt, mcTrackPt, float[4]); +DECLARE_SOA_COLUMN(McTrackEta, mcTrackEta, float[4]); +DECLARE_SOA_COLUMN(McTrackPhi, mcTrackPhi, float[4]); +DECLARE_SOA_COLUMN(McTrackSign, mcTrackSign, int[4]); +DECLARE_SOA_COLUMN(McTrackIsPrimary, mcTrackIsPrimary, int[4]); + +// Generated vertex +DECLARE_SOA_COLUMN(McPosX, mcPosX, float); +DECLARE_SOA_COLUMN(McPosY, mcPosY, float); +DECLARE_SOA_COLUMN(McPosZ, mcPosZ, float); + +// Generated collision index +DECLARE_SOA_COLUMN(McCollisionIndex, mcCollisionIndex, int); + +// Real run number for MC gen +DECLARE_SOA_COLUMN(McRunNumber, mcRunNumber, int); + +DECLARE_SOA_COLUMN(RecoIndex, recoIndex, int); +} // namespace mcgen4pi + +// MC reco<->gen match table +DECLARE_SOA_TABLE(FourPiMcMatchTree, "AOD", "FOURPIMCMATCH", + fourpi::IsReconstructedWithUPC, + mcgen4pi::McMotherPdg, mcgen4pi::McMotherPt, mcgen4pi::McMotherPhi, + mcgen4pi::McMotherMass, mcgen4pi::McMotherRapidity, mcgen4pi::McTotalCharge, + mcgen4pi::McTrackPdg, mcgen4pi::McTrackPt, mcgen4pi::McTrackEta, mcgen4pi::McTrackPhi, + mcgen4pi::McTrackSign, mcgen4pi::McTrackIsPrimary, + mcgen4pi::McPosX, mcgen4pi::McPosY, mcgen4pi::McPosZ, mcgen4pi::McCollisionIndex, + mcgen4pi::McRunNumber, mcgen4pi::RecoIndex); + +// MC All generated collisions +DECLARE_SOA_TABLE(FourPiMcGenAllTree, "AOD", "FOURPIMCGALL", + mcgen4pi::McMotherPdg, mcgen4pi::McMotherPt, mcgen4pi::McMotherPhi, + mcgen4pi::McMotherMass, mcgen4pi::McMotherRapidity, mcgen4pi::McTotalCharge, + mcgen4pi::McTrackPdg, mcgen4pi::McTrackPt, mcgen4pi::McTrackEta, mcgen4pi::McTrackPhi, + mcgen4pi::McTrackSign, mcgen4pi::McTrackIsPrimary, + mcgen4pi::McPosX, mcgen4pi::McPosY, mcgen4pi::McPosZ, mcgen4pi::McCollisionIndex, + mcgen4pi::McRunNumber); } // namespace o2::aod struct upcRhoPrimeAnalysis { Produces systemTree; + Produces fourPiMcMatchTree; + Produces fourPiMcGenAllTree; // System selection configuration Configurable systemYCut{"systemYCut", 0.5, "Max Rapidity of rho prime"}; @@ -127,11 +185,16 @@ struct upcRhoPrimeAnalysis { Configurable dcaZcut{"dcaZcut", 2, "dcaZ cut"}; Configurable minTPCFindableClusters{"minTPCFindableClusters", 70, "Minimum number of findable TPC clusters"}; - // Define histogram registry + // Optional generatorId filter + Configurable genId{"genId", -1, "generator ID; -1 = no filter"}; + + // Define histogram registry for RD / MC reco HistogramRegistry registry{ "registry", {// Event flow histograms - {"Events/Flow", "Event flow;Cut;Counts", {HistType::kTH1F, {{9, 0, 9}}}}, + {"Events/Flow", "Event flow;Cut;Counts", {HistType::kTH1D, {{11, 0, 11}}}}, + {"Events/FlowDetailed", "Detailed event flow;Cut;Counts", {HistType::kTH1D, {{13, 0, 13}}}}, + {"Events/hRecoMode", "Reconstruction mode;;Counts", {HistType::kTH1D, {{2, 0, 2}}}}, {"Events/VertexZ", "Vertex Z;z (cm);Counts", {HistType::kTH1F, {{200, -20, 20}}}}, {"Events/NumContrib", "Number of contributors;N_{contrib};Counts", {HistType::kTH1F, {{100, 0, 100}}}}, {"Events/FV0Amplitude", "FV0 amplitude;Amplitude;Counts", {HistType::kTH1F, {{200, 0, 200}}}}, @@ -145,23 +208,148 @@ struct upcRhoPrimeAnalysis { {"Tracks/TPCNSigmaPi", "TPC n#sigma for #pi;n#sigma;Counts", {HistType::kTH1F, {{200, -10, 10}}}}, {"Tracks/TPCChi2NCl", "TPC #chi^{2}/N_{cls};#chi^{2}/N_{cls};Counts", {HistType::kTH1F, {{200, 0, 20}}}}, {"Tracks/ITSChi2NCl", "ITS #chi^{2}/N_{cls};#chi^{2}/N_{cls};Counts", {HistType::kTH1F, {{200, 0, 50}}}}, - {"Tracks/RejectionReasons", "Track rejection reasons;Reason;Counts", {HistType::kTH1F, {{12, 0, 12}}}}, + {"Tracks/RejectionReasons", "Track rejection reasons;Reason;Counts", {HistType::kTH1F, {{15, 0, 15}}}}, {"Tracks/DCASpectrum", "Track DCA spectrum;DCA (cm);Counts", {HistType::kTH1F, {{100, 0, 5}}}}, {"Tracks/ChargeDistribution", "Track charge distribution;Charge;Counts", {HistType::kTH1F, {{3, -1.5, 1.5}}}}, {"Tracks/TPCClusters", "TPC clusters findable;N_{clusters};Counts", {HistType::kTH1F, {{100, 0, 200}}}}, + {"Tracks/NGoodTracksPerEvent", "Good tracks per event (before the ==4 cut);N;Counts", {HistType::kTH1F, {{11, -0.5, 10.5}}}}, // System kinematics histograms {"System/hM", ";m (GeV/#it{c}^{2});counts", {HistType::kTH1F, {{1000, 0.0, 10.0}}}}, - {"System/hPt", ";p_{T} (GeV/#it{c});counts", {HistType::kTH1F, {{1000, 0.0, 10.0}}}}, + {"System/hPt", ";p_{T} (GeV/#it{c});counts", {HistType::kTH1F, {{1000, 0.0, 1.1}}}}, {"System/hEta", ";#eta;counts", {HistType::kTH1F, {{180, -0.9, 0.9}}}}, {"System/hPhi", ";#phi;counts", {HistType::kTH1F, {{180, 0.0, 6.28}}}}, {"System/hY", ";y;counts", {HistType::kTH1F, {{180, -0.9, 0.9}}}}, + {"System/hTotalChargeBefore", "Total charge before M/Pt/Y cuts;Q;counts", {HistType::kTH1F, {{9, -4.5, 4.5}}}}, + // -4=----, -2=---+ (3-,1+), 0=+-+- (2+,2-), +2=+++- (3+,1-), +4=++++ + {"System/hTotalCharge", "System total charge (4 tracks, after M/Pt/Y cuts);Q;counts", {HistType::kTH1F, {{9, -4.5, 4.5}}}}, + {"System/hMVsTotalChargeBefore", "Invariant mass vs charge combination (before M/Pt/Y cuts);m (GeV/#it{c}^{2});Q", {HistType::kTH2F, {{1000, 0.0, 10.0}, {9, -4.5, 4.5}}}}, + {"System/hMVsTotalCharge", "Invariant mass vs charge combination (after M/Pt/Y cuts);m (GeV/#it{c}^{2});Q", {HistType::kTH2F, {{1000, 0.0, 10.0}, {9, -4.5, 4.5}}}}, // Comparison histograms {"Cuts/MBefore", "Mass before cuts;m (GeV/c^{2});Counts", {HistType::kTH1F, {{1000, 0, 10}}}}, {"Cuts/MAfter", "Mass after cuts;m (GeV/c^{2});Counts", {HistType::kTH1F, {{1000, 0, 10}}}}, - {"Cuts/PtBefore", "p_{T} before cuts;p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0, 1}}}}, - {"Cuts/PtAfter", "p_{T} after cuts;p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0, 10}}}}}}; + {"Cuts/PtBefore", "p_{T} before cuts;p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0, 1.1}}}}, + {"Cuts/PtAfter", "p_{T} after cuts;p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0, 1.1}}}}}}; + + // Define histogram registry for MC analysis + HistogramRegistry mcRegistry{ + "mcRegistry", + {// Event-level MC histograms + {"MC/Events/hAllEvents", "All MC Events", {HistType::kTH1F, {{1, 0, 1}}}}, + {"MC/Events/hAccepted", "Accepted MC Events", {HistType::kTH1F, {{1, 0, 1}}}}, + {"MC/Events/hVertexZ", "MC Vertex Z;z (cm);Counts", {HistType::kTH1F, {{400, -20.0, 20.0}}}}, + {"MC/Events/hNPrimaries", "Number of primary particles;N;Counts", {HistType::kTH1F, {{10, -0.5, 9.5}}}}, + + // Track-level MC histograms + {"MC/Tracks/hPt", "MC Track p_{T};p_{T} (GeV/c);Counts", {HistType::kTH1F, {{200, 0.0, 2.0}}}}, + {"MC/Tracks/hEta", "MC Track #eta;#eta;Counts", {HistType::kTH1F, {{200, -2.0, 2.0}}}}, + {"MC/Tracks/hPhi", "MC Track #phi;#phi;Counts", {HistType::kTH1F, {{200, 0.0, 6.28}}}}, + {"MC/Tracks/hPdgCode", "PDG codes;PDG code;Counts", {HistType::kTH1F, {{2000, -1000, 1000}}}}, + + // System-level (mother) MC histograms + {"MC/System/hM", "MC Invariant Mass;m (GeV/c^{2});Counts", {HistType::kTH1F, {{1000, 0.0, 10.0}}}}, + {"MC/System/hPt", "MC p_{T};p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 1.1}}}}, + {"MC/System/hY", "MC Rapidity;y;Counts", {HistType::kTH1F, {{180, -0.9, 0.9}}}}, + {"MC/System/hMvsPt", "MC Mass vs p_{T};m (GeV/c^{2});p_{T} (GeV/c)", {HistType::kTH2F, {{1000, 0.0, 10.0}, {1000, 0.0, 1.1}}}}, + {"MC/System/hMvsY", "MC Mass vs Rapidity;m (GeV/c^{2});y", {HistType::kTH2F, {{1000, 0.0, 10.0}, {200, -2.0, 2.0}}}}, + {"MC/System/hTotalCharge", "Generated total charge;Q;Counts", {HistType::kTH1F, {{9, -4.5, 4.5}}}}, + + // Basic summary + {"MC/Summary/hEventCounter", "Generated vs reconstructed summary;;Counts", {HistType::kTH1D, {{2, 0, 2}}}}, + {"MC/Summary/hMatchStatus", "Reco<->Gen match status;;Counts", {HistType::kTH1D, {{3, 0, 3}}}}, + {"MC/Summary/hRecoMode", "Reconstruction mode;;Counts", {HistType::kTH1D, {{2, 0, 2}}}}, + {"MC/Summary/hRecoTotalCharge", "Reco total charge (MC events, after M/Pt/Y cuts);Q;Counts", {HistType::kTH1D, {{9, -4.5, 4.5}}}}, + + // Rho prime specific histograms + {"MC/RhoPrime/hFound", "Rho Prime Found;Found;Counts", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, + {"MC/RhoPrime/hMass", "Rho Prime Mass;m (GeV/c^{2});Counts", {HistType::kTH1F, {{1000, 0.0, 10.0}}}}, + {"MC/RhoPrime/hMassUPC", "Rho Prime Mass UPC;m (GeV/c^{2});Counts", {HistType::kTH1F, {{1000, 0.0, 10.0}}}}, + {"MC/RhoPrime/hMassSTD", "Rho Prime Mass STD;m (GeV/c^{2});Counts", {HistType::kTH1F, {{1000, 0.0, 10.0}}}}, + {"MC/RhoPrime/hPt", "Rho Prime p_{T};p_{T} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 1.1}}}}, + + // Control histograms + {"MC/Control/hNDaughtersOfMother", "Number of daughters of the common mother;N;Counts", {HistType::kTH1F, {{10, 0, 10}}}}, + {"MC/Control/hMotherPdg", "PDG of the common mother found;PDG;Counts", {HistType::kTH1F, {{200, 0, 40000}}}}, + {"MC/Control/hTracksWithMcParticle", "Reco tracks with an associated mcParticle (of 4);N;Counts", {HistType::kTH1F, {{5, -0.5, 4.5}}}}, + + // MC reco<->gen match histograms + {"MC/Match/hRecoEvents", "Reco events entering the match check;;Counts", {HistType::kTH1F, {{1, 0, 1}}}}, + {"MC/Match/hMatchedGenM", "Matched candidates - generated M;m_{gen} (GeV/c^{2});Counts", {HistType::kTH1F, {{1000, 0.0, 10.0}}}}, + {"MC/Match/hMatchedGenPt", "Matched candidates - generated p_{T};p_{T,gen} (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 1.1}}}}, + {"MC/Match/hMatchedGenY", "Matched candidates - generated y;y_{gen};Counts", {HistType::kTH1F, {{180, -0.9, 0.9}}}}, + {"MC/Match/hRecoVsGenM", "Reco vs Gen M;m_{gen} (GeV/c^{2});m_{reco} (GeV/c^{2})", {HistType::kTH2F, {{500, 0.0, 5.0}, {500, 0.0, 5.0}}}}}}; + + // Helper functions for kinematic calculations + static float pt(float px, float py) { return std::sqrt(px * px + py * py); } + + static float eta(float px, float py, float pz) + { + if (std::abs(pz) > 1e-10) { + float p = std::sqrt(px * px + py * py + pz * pz); + return 0.5f * std::log((p + pz) / (p - pz)); + } + return 0.0f; + } + + static float phi(float px, float py) + { + if (std::abs(px) > 1e-10 || std::abs(py) > 1e-10) { + return std::atan2(py, px); + } + return 0.0f; + } + + // Generic "common mother" + struct MotherInfo { + bool found = false; + int pdg = 0; + float px = 0, py = 0; + }; + + // 4 particles share the same direct mother + template + MotherInfo findCommonMotherGeneric(const std::vector& daughters) + { + MotherInfo info; + if (daughters.empty() || !daughters[0].has_mothers()) { + return info; + } + auto firstMothers = daughters[0].template mothers_as(); + if (firstMothers.begin() == firstMothers.end()) { + return info; + } + auto motherIt = firstMothers.begin(); + int64_t motherGlobalIndex = motherIt->globalIndex(); + for (size_t i = 1; i < daughters.size(); i++) { + if (!daughters[i].has_mothers()) { + return info; + } + auto iMothers = daughters[i].template mothers_as(); + if (iMothers.begin() == iMothers.end() || iMothers.begin()->globalIndex() != motherGlobalIndex) { + return info; + } + } + const auto& mother = *motherIt; + info.found = true; + info.pdg = mother.pdgCode(); + info.px = mother.px(); + info.py = mother.py(); + return info; + } + + // Charge sign from PDG code + static int signFromPdg(int pdgCode) + { + if (pdgCode == 211) { + return 1; + } + if (pdgCode == -211) { + return -1; + } + return (pdgCode > 0) ? 1 : (pdgCode < 0) ? -1 + : 0; + } void init(InitContext&) { @@ -173,32 +361,135 @@ struct upcRhoPrimeAnalysis { hFlow->GetXaxis()->SetBinLabel(4, "ITS ROFb cut"); hFlow->GetXaxis()->SetBinLabel(5, "TFB cut"); hFlow->GetXaxis()->SetBinLabel(6, "Gap Side cut"); - hFlow->GetXaxis()->SetBinLabel(7, "PV contrib cut"); - hFlow->GetXaxis()->SetBinLabel(8, "Z vtx cut"); - hFlow->GetXaxis()->SetBinLabel(9, "4 tracks cut"); + hFlow->GetXaxis()->SetBinLabel(7, "ZDC energy cut"); + hFlow->GetXaxis()->SetBinLabel(8, "PV contrib cut"); + hFlow->GetXaxis()->SetBinLabel(9, "Z vtx cut"); + hFlow->GetXaxis()->SetBinLabel(10, "4 tracks cut"); + hFlow->GetXaxis()->SetBinLabel(11, "System cuts (M,Pt,Y)"); + + auto hFlowDetailed = registry.get(HIST("Events/FlowDetailed")); + hFlowDetailed->GetXaxis()->SetBinLabel(1, "All events"); + hFlowDetailed->GetXaxis()->SetBinLabel(2, "vtxITSTPC"); + hFlowDetailed->GetXaxis()->SetBinLabel(3, "sbp"); + hFlowDetailed->GetXaxis()->SetBinLabel(4, "itsROFb"); + hFlowDetailed->GetXaxis()->SetBinLabel(5, "tfb"); + hFlowDetailed->GetXaxis()->SetBinLabel(6, "gapSide"); + hFlowDetailed->GetXaxis()->SetBinLabel(7, "FV0 < cut"); + hFlowDetailed->GetXaxis()->SetBinLabel(8, "FT0A < cut"); + hFlowDetailed->GetXaxis()->SetBinLabel(9, "FT0C < cut"); + hFlowDetailed->GetXaxis()->SetBinLabel(10, "ZDC energy < cut"); + hFlowDetailed->GetXaxis()->SetBinLabel(11, "numContrib == 4"); + hFlowDetailed->GetXaxis()->SetBinLabel(12, "posZ < cut"); + hFlowDetailed->GetXaxis()->SetBinLabel(13, "4 tracks"); + + // Readable labels for Events/hRecoMode + auto hRecoMode = registry.get(HIST("Events/hRecoMode")); + hRecoMode->GetXaxis()->SetBinLabel(1, "STD"); + hRecoMode->GetXaxis()->SetBinLabel(2, "UPC"); + + auto hChargeBefore = registry.get(HIST("System/hTotalChargeBefore")); + hChargeBefore->GetXaxis()->SetBinLabel(1, "----"); + hChargeBefore->GetXaxis()->SetBinLabel(3, "---+"); + hChargeBefore->GetXaxis()->SetBinLabel(5, "+-+-"); + hChargeBefore->GetXaxis()->SetBinLabel(7, "+++-"); + hChargeBefore->GetXaxis()->SetBinLabel(9, "++++"); + + auto hCharge = registry.get(HIST("System/hTotalCharge")); + hCharge->GetXaxis()->SetBinLabel(1, "----"); + hCharge->GetXaxis()->SetBinLabel(3, "---+"); + hCharge->GetXaxis()->SetBinLabel(5, "+-+-"); + hCharge->GetXaxis()->SetBinLabel(7, "+++-"); + hCharge->GetXaxis()->SetBinLabel(9, "++++"); + + auto hMChargeBefore = registry.get(HIST("System/hMVsTotalChargeBefore")); + hMChargeBefore->GetYaxis()->SetBinLabel(1, "----"); + hMChargeBefore->GetYaxis()->SetBinLabel(3, "---+"); + hMChargeBefore->GetYaxis()->SetBinLabel(5, "+-+-"); + hMChargeBefore->GetYaxis()->SetBinLabel(7, "+++-"); + hMChargeBefore->GetYaxis()->SetBinLabel(9, "++++"); + + auto hMCharge = registry.get(HIST("System/hMVsTotalCharge")); + hMCharge->GetYaxis()->SetBinLabel(1, "----"); + hMCharge->GetYaxis()->SetBinLabel(3, "---+"); + hMCharge->GetYaxis()->SetBinLabel(5, "+-+-"); + hMCharge->GetYaxis()->SetBinLabel(7, "+++-"); + hMCharge->GetYaxis()->SetBinLabel(9, "++++"); - // Configure track rejection reasons histogram labels auto hReject = registry.get(HIST("Tracks/RejectionReasons")); - hReject->GetXaxis()->SetBinLabel(1, "All Tracks"); - hReject->GetXaxis()->SetBinLabel(2, "PV Contributor"); - hReject->GetXaxis()->SetBinLabel(3, "Has ITS+TPC"); - hReject->GetXaxis()->SetBinLabel(4, "pT > 0.1 GeV/c"); - hReject->GetXaxis()->SetBinLabel(5, "TPC chi2/cluster"); - hReject->GetXaxis()->SetBinLabel(6, "ITS chi2/cluster"); - hReject->GetXaxis()->SetBinLabel(7, "TPC clusters findable"); - hReject->GetXaxis()->SetBinLabel(8, "TPC nSigmaPi"); - hReject->GetXaxis()->SetBinLabel(9, "Eta acceptance"); - hReject->GetXaxis()->SetBinLabel(10, "DCAz cut"); - hReject->GetXaxis()->SetBinLabel(11, "DCAxy cut"); - hReject->GetXaxis()->SetBinLabel(12, "Accepted Tracks"); + hReject->GetXaxis()->SetBinLabel(1, "All tracks"); + hReject->GetXaxis()->SetBinLabel(2, "isPVContributor"); + hReject->GetXaxis()->SetBinLabel(3, "hasITS && hasTPC"); + hReject->GetXaxis()->SetBinLabel(4, "pT > 0.1"); + hReject->GetXaxis()->SetBinLabel(5, "tpcChi2NCl < cut"); + hReject->GetXaxis()->SetBinLabel(6, "itsChi2NCl < cut"); + hReject->GetXaxis()->SetBinLabel(7, "tpcNClsFindable > cut"); + hReject->GetXaxis()->SetBinLabel(8, "|tpcNSigmaPi| < cut"); + hReject->GetXaxis()->SetBinLabel(9, "|eta| < cut"); + hReject->GetXaxis()->SetBinLabel(10, "|dcaZ| < cut"); + hReject->GetXaxis()->SetBinLabel(11, "|dcaXY| < cut"); + hReject->GetXaxis()->SetBinLabel(12, "Accepted tracks"); + + auto hPdg = mcRegistry.get(HIST("MC/Tracks/hPdgCode")); + if (hPdg) { + int bin211 = hPdg->GetXaxis()->FindBin(211); + int binM211 = hPdg->GetXaxis()->FindBin(-211); + int bin30113 = hPdg->GetXaxis()->FindBin(30113); + if (bin211 > 0 && bin211 <= hPdg->GetNbinsX()) + hPdg->GetXaxis()->SetBinLabel(bin211, "#pi^{+}"); + if (binM211 > 0 && binM211 <= hPdg->GetNbinsX()) + hPdg->GetXaxis()->SetBinLabel(binM211, "#pi^{-}"); + if (bin30113 > 0 && bin30113 <= hPdg->GetNbinsX()) + hPdg->GetXaxis()->SetBinLabel(bin30113, "rho'"); + } + + auto hRhoFound = mcRegistry.get(HIST("MC/RhoPrime/hFound")); + hRhoFound->GetXaxis()->SetBinLabel(1, "Not Found"); + hRhoFound->GetXaxis()->SetBinLabel(2, "Found"); + + auto hMcCharge = mcRegistry.get(HIST("MC/System/hTotalCharge")); + hMcCharge->GetXaxis()->SetBinLabel(1, "----"); + hMcCharge->GetXaxis()->SetBinLabel(3, "---+"); + hMcCharge->GetXaxis()->SetBinLabel(5, "+-+-"); + hMcCharge->GetXaxis()->SetBinLabel(7, "+++-"); + hMcCharge->GetXaxis()->SetBinLabel(9, "++++"); + + auto hMcSummary = mcRegistry.get(HIST("MC/Summary/hEventCounter")); + hMcSummary->GetXaxis()->SetBinLabel(1, "Generated (denominator, gen all)"); + hMcSummary->GetXaxis()->SetBinLabel(2, "Reconstructed + matched (numerator)"); + + auto hMatchStatus = mcRegistry.get(HIST("MC/Summary/hMatchStatus")); + hMatchStatus->GetXaxis()->SetBinLabel(1, "No mcParticle on some track"); + hMatchStatus->GetXaxis()->SetBinLabel(2, "Has mcParticle but no common mother"); + hMatchStatus->GetXaxis()->SetBinLabel(3, "Match found"); + + auto hMcRecoMode = mcRegistry.get(HIST("MC/Summary/hRecoMode")); + hMcRecoMode->GetXaxis()->SetBinLabel(1, "STD"); + hMcRecoMode->GetXaxis()->SetBinLabel(2, "UPC"); + + auto hMcRecoCharge = mcRegistry.get(HIST("MC/Summary/hRecoTotalCharge")); + hMcRecoCharge->GetXaxis()->SetBinLabel(1, "----"); + hMcRecoCharge->GetXaxis()->SetBinLabel(3, "---+"); + hMcRecoCharge->GetXaxis()->SetBinLabel(5, "+-+-"); + hMcRecoCharge->GetXaxis()->SetBinLabel(7, "+++-"); + hMcRecoCharge->GetXaxis()->SetBinLabel(9, "++++"); + + if (doprocessDataCols && (doprocessMcCols || doprocessMcGenAll)) { + LOGP(fatal, + "Invalid config: processDataCols must not run together with processMcCols/processMcGenAll " + "(it would duplicate SystemTree rows). Disable one of the two groups in your configuration.json."); + } } - void process(UDCollisions::iterator const& collision, UDtracks const& tracks) + // Collision + track selection + template + std::vector selectAndFillSystemTree(ColType const& collision, TracksType const& tracks, bool& outFilled) { - // Count all processed events + outFilled = false; + std::vector goodTracks; + registry.fill(HIST("Events/Flow"), 0); + registry.fill(HIST("Events/FlowDetailed"), 0); - // Fill basic event diagnostics registry.fill(HIST("Events/VertexZ"), collision.posZ()); registry.fill(HIST("Events/NumContrib"), collision.numContrib()); registry.fill(HIST("Events/FV0Amplitude"), collision.totalFV0AmplitudeA()); @@ -207,64 +498,84 @@ struct upcRhoPrimeAnalysis { registry.fill(HIST("Events/ZDCEnergy"), collision.energyCommonZNA()); registry.fill(HIST("Events/ZDCEnergy"), collision.energyCommonZNC()); - // Apply event selection cuts in sequence - if (collision.vtxITSTPC() != vtxITSTPCcut) - return; + if (collision.vtxITSTPC() != vtxITSTPCcut) { + return goodTracks; + } registry.fill(HIST("Events/Flow"), 1); + registry.fill(HIST("Events/FlowDetailed"), 1); - if (collision.sbp() != sbpCut) - return; + if (collision.sbp() != sbpCut) { + return goodTracks; + } registry.fill(HIST("Events/Flow"), 2); + registry.fill(HIST("Events/FlowDetailed"), 2); - if (collision.itsROFb() != itsROFbCut) - return; + if (collision.itsROFb() != itsROFbCut) { + return goodTracks; + } registry.fill(HIST("Events/Flow"), 3); + registry.fill(HIST("Events/FlowDetailed"), 3); - if (collision.tfb() != tfbCut) - return; + if (collision.tfb() != tfbCut) { + return goodTracks; + } registry.fill(HIST("Events/Flow"), 4); + registry.fill(HIST("Events/FlowDetailed"), 4); + + if (specifyGapSide && collision.gapSide() != gapSide) { + return goodTracks; + } - if (specifyGapSide && collision.gapSide() != gapSide) - return; - if (collision.totalFV0AmplitudeA() > fv0Cut) - return; - if (collision.totalFT0AmplitudeA() > ft0aCut) - return; - if (collision.totalFT0AmplitudeC() > ft0cCut) - return; - if (collision.energyCommonZNA() > zdcCut || collision.energyCommonZNC() > zdcCut) - return; registry.fill(HIST("Events/Flow"), 5); + registry.fill(HIST("Events/FlowDetailed"), 5); - if (collision.numContrib() != numPVContrib) - return; + if (collision.totalFV0AmplitudeA() > fv0Cut) { + return goodTracks; + } + registry.fill(HIST("Events/FlowDetailed"), 6); + + if (collision.totalFT0AmplitudeA() > ft0aCut) { + return goodTracks; + } + registry.fill(HIST("Events/FlowDetailed"), 7); + + if (collision.totalFT0AmplitudeC() > ft0cCut) { + return goodTracks; + } + registry.fill(HIST("Events/FlowDetailed"), 8); + + if (collision.energyCommonZNA() > zdcCut || collision.energyCommonZNC() > zdcCut) { + return goodTracks; + } registry.fill(HIST("Events/Flow"), 6); + registry.fill(HIST("Events/FlowDetailed"), 9); - if (std::abs(collision.posZ()) > vZCut) - return; + if (collision.numContrib() != numPVContrib) { + return goodTracks; + } registry.fill(HIST("Events/Flow"), 7); + registry.fill(HIST("Events/FlowDetailed"), 10); - std::vector posPions; - std::vector negPions; - posPions.reserve(2); - negPions.reserve(2); + if (std::abs(collision.posZ()) > vZCut) { + return goodTracks; + } + registry.fill(HIST("Events/Flow"), 8); + registry.fill(HIST("Events/FlowDetailed"), 11); - // Loop over all tracks in the event + // --- Track selection: up to 4 good tracks, charge combination --- + goodTracks.reserve(4); for (const auto& track : tracks) { - registry.fill(HIST("Tracks/RejectionReasons"), 0); // Count all tracks + registry.fill(HIST("Tracks/RejectionReasons"), 0); - // Track selection criteria applied in sequence: if (useOnlyPVtracks && !track.isPVContributor()) { registry.fill(HIST("Tracks/RejectionReasons"), 1); continue; } - if (!track.hasITS() || !track.hasTPC()) { registry.fill(HIST("Tracks/RejectionReasons"), 2); continue; } - // Fill track spectra registry.fill(HIST("Tracks/Pt"), track.pt()); registry.fill(HIST("Tracks/Eta"), eta(track.px(), track.py(), track.pz())); registry.fill(HIST("Tracks/TPCNSigmaPi"), track.tpcNSigmaPi()); @@ -278,7 +589,6 @@ struct upcRhoPrimeAnalysis { registry.fill(HIST("Tracks/RejectionReasons"), 3); continue; } - if (track.tpcChi2NCl() > tpcChi2NClsCut) { registry.fill(HIST("Tracks/RejectionReasons"), 4); continue; @@ -287,84 +597,71 @@ struct upcRhoPrimeAnalysis { registry.fill(HIST("Tracks/RejectionReasons"), 5); continue; } - if (track.tpcNClsFindable() < minTPCFindableClusters) { registry.fill(HIST("Tracks/RejectionReasons"), 6); continue; } - if (std::abs(track.tpcNSigmaPi()) > nSigmaTPCcut) { registry.fill(HIST("Tracks/RejectionReasons"), 7); continue; } - - float trackEta = eta(track.px(), track.py(), track.pz()); - if (std::abs(trackEta) > etaCut) { + if (std::abs(eta(track.px(), track.py(), track.pz())) > etaCut) { registry.fill(HIST("Tracks/RejectionReasons"), 8); continue; } - if (std::abs(track.dcaZ()) > dcaZcut) { registry.fill(HIST("Tracks/RejectionReasons"), 9); continue; } - float maxDCAxy = 0.0105 + 0.035 / std::pow(track.pt(), 1.1); - if (dcaXYcut == 0 && (std::fabs(track.dcaXY()) > maxDCAxy)) { - registry.fill(HIST("Tracks/RejectionReasons"), 10); - continue; - } else if (dcaXYcut != 0 && (std::fabs(track.dcaXY()) > dcaXYcut)) { + if (dcaXYcut == 0 && std::fabs(track.dcaXY()) > maxDCAxy) { registry.fill(HIST("Tracks/RejectionReasons"), 10); continue; } - // Track passed all selection criteria registry.fill(HIST("Tracks/RejectionReasons"), 11); - - if (track.sign() > 0 && posPions.size() < 2) { - posPions.push_back(track); - } else if (track.sign() < 0 && negPions.size() < 2) { - negPions.push_back(track); - } - - if (posPions.size() == 2 && negPions.size() == 2) + goodTracks.push_back(track); + if (goodTracks.size() == 4) { break; + } } - if (posPions.size() != 2 || negPions.size() != 2) { - return; - } - registry.fill(HIST("Events/Flow"), 8); + // Basic control + registry.fill(HIST("Tracks/NGoodTracksPerEvent"), goodTracks.size()); - std::vector selectedTracks; - selectedTracks.insert(selectedTracks.end(), posPions.begin(), posPions.end()); - selectedTracks.insert(selectedTracks.end(), negPions.begin(), negPions.end()); + if (goodTracks.size() != 4) { + return goodTracks; + } + registry.fill(HIST("Events/Flow"), 9); + registry.fill(HIST("Events/FlowDetailed"), 12); - // Reconstruct the 4-pion system - ROOT::Math::PxPyPzMVector fourPionSystem; - std::vector pionFourVectors; + // Real total charge + int totalCharge = 0; + for (const auto& track : goodTracks) { + totalCharge += track.sign(); + } + bool isChargeZero = (totalCharge == 0); + registry.fill(HIST("System/hTotalChargeBefore"), totalCharge); - for (const auto& track : selectedTracks) { - ROOT::Math::PxPyPzMVector pionVec( - track.px(), track.py(), track.pz(), - o2::constants::physics::MassPionCharged); - fourPionSystem += pionVec; - pionFourVectors.push_back(pionVec); + PxPyPzMVector fourPionSystem; + for (const auto& track : goodTracks) { + fourPionSystem += PxPyPzMVector(track.px(), track.py(), track.pz(), o2::constants::physics::MassPionCharged); } - // Fill pre-cut system histograms registry.fill(HIST("Cuts/MBefore"), fourPionSystem.M()); registry.fill(HIST("Cuts/PtBefore"), fourPionSystem.Pt()); + registry.fill(HIST("System/hMVsTotalChargeBefore"), fourPionSystem.M(), totalCharge); - // Apply system-level kinematic cuts - if (fourPionSystem.M() < systemMassMinCut || fourPionSystem.M() > systemMassMaxCut) - return; - if (fourPionSystem.Pt() > systemPtCut) - return; - if (std::abs(fourPionSystem.Rapidity()) > systemYCut) - return; + if (fourPionSystem.M() < systemMassMinCut || fourPionSystem.M() > systemMassMaxCut) { + return goodTracks; + } + if (fourPionSystem.Pt() > systemPtCut) { + return goodTracks; + } + if (std::abs(fourPionSystem.Rapidity()) > systemYCut) { + return goodTracks; + } - // Fill post-cut system histograms registry.fill(HIST("Cuts/MAfter"), fourPionSystem.M()); registry.fill(HIST("Cuts/PtAfter"), fourPionSystem.Pt()); registry.fill(HIST("System/hM"), fourPionSystem.M()); @@ -372,13 +669,14 @@ struct upcRhoPrimeAnalysis { registry.fill(HIST("System/hEta"), fourPionSystem.Eta()); registry.fill(HIST("System/hPhi"), fourPionSystem.Phi() + o2::constants::math::PI); registry.fill(HIST("System/hY"), fourPionSystem.Rapidity()); + registry.fill(HIST("System/hTotalCharge"), totalCharge); + registry.fill(HIST("System/hMVsTotalCharge"), fourPionSystem.M(), totalCharge); std::vector trackPts, trackEtas, trackPhis; std::vector trackSigns, trackIDs; std::vector tpcNSigmasEl, tpcNSigmasPi, tpcNSigmasKa, tpcNSigmasPr; - - for (size_t i = 0; i < selectedTracks.size(); i++) { - const auto& track = selectedTracks[i]; + for (size_t i = 0; i < goodTracks.size(); i++) { + const auto& track = goodTracks[i]; trackPts.push_back(track.pt()); trackEtas.push_back(eta(track.px(), track.py(), track.pz())); trackPhis.push_back(phi(track.px(), track.py())); @@ -387,48 +685,272 @@ struct upcRhoPrimeAnalysis { tpcNSigmasPi.push_back(track.tpcNSigmaPi()); tpcNSigmasKa.push_back(track.tpcNSigmaKa()); tpcNSigmasPr.push_back(track.tpcNSigmaPr()); - trackIDs.push_back(i); + trackIDs.push_back(static_cast(i)); } bool isReconstructedWithUPC = (collision.flags() == 1); - // Fill the output + registry.fill(HIST("Events/Flow"), 10); + registry.fill(HIST("Events/hRecoMode"), isReconstructedWithUPC ? 1 : 0); + systemTree( collision.runNumber(), - fourPionSystem.M(), - fourPionSystem.Pt(), - fourPionSystem.Rapidity(), - fourPionSystem.Phi(), - collision.posX(), - collision.posY(), - collision.posZ(), - 0, // Total charge = 0 - collision.totalFT0AmplitudeA(), - collision.totalFT0AmplitudeC(), - collision.totalFV0AmplitudeA(), + fourPionSystem.M(), fourPionSystem.Pt(), fourPionSystem.Rapidity(), fourPionSystem.Phi(), + collision.posX(), collision.posY(), collision.posZ(), + totalCharge, + collision.totalFT0AmplitudeA(), collision.totalFT0AmplitudeC(), collision.totalFV0AmplitudeA(), collision.numContrib(), - trackSigns, - trackPts, - trackEtas, - trackPhis, - tpcNSigmasEl, - tpcNSigmasPi, - tpcNSigmasKa, - tpcNSigmasPr, - trackIDs, + trackSigns, trackPts, trackEtas, trackPhis, + tpcNSigmasEl, tpcNSigmasPi, tpcNSigmasKa, tpcNSigmasPr, + trackIDs, isReconstructedWithUPC, + collision.timeZNA(), collision.timeZNC(), collision.energyCommonZNA(), collision.energyCommonZNC(), + isChargeZero, collision.occupancyInTime(), collision.hadronicRate()); + + outFilled = true; + return goodTracks; + } // end selectAndFillSystemTree + + int getMcRunNumber(aod::BCs const& bcs) + { + if (bcs.size() == 0) { + return -1; + } + auto bc = bcs.begin(); + return bc.runNumber(); + } + + // processDataCols: RD or MC reco + void processDataCols(UDCollisions::iterator const& collision, UDtracks const& tracks) + { + bool filled = false; + selectAndFillSystemTree(collision, tracks, filled); + } + PROCESS_SWITCH(upcRhoPrimeAnalysis, processDataCols, "Process real data or MC reco", true); + + // processMcCols: MC only + void processMcCols(UDCollisionsMC::iterator const& collision, UDtracksMC const& tracks, UDMcParticles const&, UDMcCollisions const&, aod::BCs const& bcs) + { + if (genId != -1) { + if (!collision.has_udMcCollision() || collision.template udMcCollision_as().generatorsID() != genId) { + return; + } + } + + bool filled = false; + auto selTrks = selectAndFillSystemTree(collision, tracks, filled); + if (!filled) { + return; + } + + bool isReconstructedWithUPC = (collision.flags() == 1); + + mcRegistry.fill(HIST("MC/Match/hRecoEvents"), 0); + mcRegistry.fill(HIST("MC/Summary/hRecoMode"), isReconstructedWithUPC ? 1 : 0); + + int recoTotalCharge = 0; + for (const auto& track : selTrks) { + recoTotalCharge += track.sign(); + } + mcRegistry.fill(HIST("MC/Summary/hRecoTotalCharge"), recoTotalCharge); + + // Defaults (pdg=0 / kinematics=-999 => "no match") + int motherPdg = 0; + float motherPt = -999, motherPhi = -999, motherMass = -999, motherRap = -999; + int mcTotalCharge = 0; + int trackPdgs[4] = {0, 0, 0, 0}; + float trackPts[4] = {-999, -999, -999, -999}; + float trackEtas[4] = {-999, -999, -999, -999}; + float trackPhis[4] = {-999, -999, -999, -999}; + int trackSigns[4] = {0, 0, 0, 0}; + int isPrimary[4] = {0, 0, 0, 0}; + float mcPosX = -999, mcPosY = -999, mcPosZ = -999; + int mcCollisionIndex = -1; + int mcRunNumber = getMcRunNumber(bcs); + + // Basic control + int nTracksWithMcParticle = 0; + for (const auto& track : selTrks) { + if (track.has_udMcParticle()) { + nTracksWithMcParticle++; + } + } + mcRegistry.fill(HIST("MC/Control/hTracksWithMcParticle"), nTracksWithMcParticle); + + // All 4 reco tracks need an associated mcParticle to look for the mother + std::vector mcParts; + bool allHaveMcParticle = true; + for (const auto& track : selTrks) { + if (!track.has_udMcParticle()) { + allHaveMcParticle = false; + break; + } + mcParts.push_back(track.template udMcParticle_as()); + } + + MotherInfo mi; + if (allHaveMcParticle) { + mi = findCommonMotherGeneric(mcParts); + } + + if (!allHaveMcParticle) { + mcRegistry.fill(HIST("MC/Summary/hMatchStatus"), 0); + } else if (!mi.found) { + mcRegistry.fill(HIST("MC/Summary/hMatchStatus"), 1); + } else { + mcRegistry.fill(HIST("MC/Summary/hMatchStatus"), 2); + mcRegistry.fill(HIST("MC/Summary/hEventCounter"), 1); + } + + if (allHaveMcParticle) { + for (size_t i = 0; i < mcParts.size() && i < 4; i++) { + trackPdgs[i] = mcParts[i].pdgCode(); + trackPts[i] = pt(mcParts[i].px(), mcParts[i].py()); + trackEtas[i] = eta(mcParts[i].px(), mcParts[i].py(), mcParts[i].pz()); + trackPhis[i] = phi(mcParts[i].px(), mcParts[i].py()); + trackSigns[i] = signFromPdg(trackPdgs[i]); + isPrimary[i] = mcParts[i].isPhysicalPrimary() ? 1 : 0; + mcTotalCharge += trackSigns[i]; + } + + auto mcCollision = mcParts[0].template udMcCollision_as(); + mcPosX = mcCollision.posX(); + mcPosY = mcCollision.posY(); + mcPosZ = mcCollision.posZ(); + mcCollisionIndex = static_cast(mcCollision.globalIndex()); + } + + if (mi.found) { + motherPdg = mi.pdg; + motherPt = pt(mi.px, mi.py); + motherPhi = phi(mi.px, mi.py); + PxPyPzMVector genSystem; + for (const auto& p : mcParts) { + genSystem += PxPyPzMVector(p.px(), p.py(), p.pz(), o2::constants::physics::MassPionCharged); + } + motherMass = genSystem.M(); + motherRap = genSystem.Rapidity(); + + mcRegistry.fill(HIST("MC/Control/hMotherPdg"), motherPdg); + mcRegistry.fill(HIST("MC/Match/hMatchedGenM"), motherMass); + mcRegistry.fill(HIST("MC/Match/hMatchedGenPt"), motherPt); + mcRegistry.fill(HIST("MC/Match/hMatchedGenY"), motherRap); + + PxPyPzMVector recoSystem; + for (const auto& track : selTrks) { + recoSystem += PxPyPzMVector(track.px(), track.py(), track.pz(), o2::constants::physics::MassPionCharged); + } + mcRegistry.fill(HIST("MC/Match/hRecoVsGenM"), motherMass, recoSystem.M()); + if (motherPdg == 30113) { // only when the mother found is actually the rho prime + if (isReconstructedWithUPC) { + mcRegistry.fill(HIST("MC/RhoPrime/hMassUPC"), recoSystem.M()); + } else { + mcRegistry.fill(HIST("MC/RhoPrime/hMassSTD"), recoSystem.M()); + } + } + } + + fourPiMcMatchTree( isReconstructedWithUPC, - collision.timeZNA(), - collision.timeZNC(), - collision.energyCommonZNA(), - collision.energyCommonZNC(), - true, // Always charge zero for our selection - collision.occupancyInTime(), - collision.hadronicRate()); + motherPdg, motherPt, motherPhi, motherMass, motherRap, mcTotalCharge, + trackPdgs, trackPts, trackEtas, trackPhis, trackSigns, isPrimary, + mcPosX, mcPosY, mcPosZ, mcCollisionIndex, mcRunNumber, + systemTree.lastIndex()); + } + PROCESS_SWITCH(upcRhoPrimeAnalysis, processMcCols, "Match MC reco<->gen (MC only)", false); + + // processMcGenAll: MC All generated collisions + void processMcGenAll(UDMcCollisions::iterator const& mcCollision, UDMcParticles const& mcParticles, aod::BCs const& bcs) + { + if (genId != -1 && mcCollision.generatorsID() != genId) { + return; + } + + int mcRunNumber = getMcRunNumber(bcs); + + mcRegistry.fill(HIST("MC/Events/hAllEvents"), 0); + mcRegistry.fill(HIST("MC/Events/hVertexZ"), mcCollision.posZ()); + + std::vector primaries; + for (const auto& part : mcParticles) { + if (part.isPhysicalPrimary()) { + primaries.push_back(part); + } + } + mcRegistry.fill(HIST("MC/Events/hNPrimaries"), primaries.size()); + + if (primaries.size() != 4) { + return; + } + mcRegistry.fill(HIST("MC/Summary/hEventCounter"), 0); + + int trackPdgs[4] = {0, 0, 0, 0}; + float trackPts[4] = {0, 0, 0, 0}; + float trackEtas[4] = {0, 0, 0, 0}; + float trackPhis[4] = {0, 0, 0, 0}; + int trackSigns[4] = {0, 0, 0, 0}; + int isPrimary[4] = {1, 1, 1, 1}; + int mcTotalCharge = 0; + + PxPyPzMVector genSystem; + for (size_t i = 0; i < 4; i++) { + const auto& p = primaries[i]; + trackPdgs[i] = p.pdgCode(); + trackPts[i] = pt(p.px(), p.py()); + trackEtas[i] = eta(p.px(), p.py(), p.pz()); + trackPhis[i] = phi(p.px(), p.py()); + trackSigns[i] = signFromPdg(trackPdgs[i]); + mcTotalCharge += trackSigns[i]; + + mcRegistry.fill(HIST("MC/Tracks/hPt"), trackPts[i]); + mcRegistry.fill(HIST("MC/Tracks/hEta"), trackEtas[i]); + mcRegistry.fill(HIST("MC/Tracks/hPhi"), trackPhis[i] + o2::constants::math::PI); // same 0-2pi convention as System/hPhi + mcRegistry.fill(HIST("MC/Tracks/hPdgCode"), trackPdgs[i]); + + genSystem += PxPyPzMVector(p.px(), p.py(), p.pz(), o2::constants::physics::MassPionCharged); + } + + // Generic common-mother search + int motherPdg = 0; + float motherPt = -999, motherPhi = -999, motherMass = -999, motherRap = -999; + MotherInfo mi = findCommonMotherGeneric(primaries); + mcRegistry.fill(HIST("MC/Control/hNDaughtersOfMother"), mi.found ? 4 : 0); + if (mi.found) { + motherPdg = mi.pdg; + motherPt = pt(mi.px, mi.py); + motherPhi = phi(mi.px, mi.py); + motherMass = genSystem.M(); + motherRap = genSystem.Rapidity(); + mcRegistry.fill(HIST("MC/Control/hMotherPdg"), motherPdg); + } + + mcRegistry.fill(HIST("MC/System/hM"), genSystem.M()); + mcRegistry.fill(HIST("MC/System/hPt"), genSystem.Pt()); + mcRegistry.fill(HIST("MC/System/hY"), genSystem.Rapidity()); + mcRegistry.fill(HIST("MC/System/hMvsPt"), genSystem.M(), genSystem.Pt()); + mcRegistry.fill(HIST("MC/System/hMvsY"), genSystem.M(), genSystem.Rapidity()); + mcRegistry.fill(HIST("MC/System/hTotalCharge"), mcTotalCharge); + + if (motherPdg == 30113) { + mcRegistry.fill(HIST("MC/RhoPrime/hFound"), 1); + mcRegistry.fill(HIST("MC/RhoPrime/hMass"), motherMass); + mcRegistry.fill(HIST("MC/RhoPrime/hPt"), motherPt); + } else { + mcRegistry.fill(HIST("MC/RhoPrime/hFound"), 0); + } + + fourPiMcGenAllTree( + motherPdg, motherPt, motherPhi, motherMass, motherRap, mcTotalCharge, + trackPdgs, trackPts, trackEtas, trackPhis, trackSigns, isPrimary, + mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), + static_cast(mcCollision.globalIndex()), mcRunNumber); + + mcRegistry.fill(HIST("MC/Events/hAccepted"), 0); } + PROCESS_SWITCH(upcRhoPrimeAnalysis, processMcGenAll, "All generated collisions (MC only)", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc, TaskName{"upc-rho-prime-analysis"})}; } diff --git a/PWGUD/Tasks/upcVmRof.cxx b/PWGUD/Tasks/upcVmRof.cxx index 5d112ad9ecb..00e3dc6ac7c 100644 --- a/PWGUD/Tasks/upcVmRof.cxx +++ b/PWGUD/Tasks/upcVmRof.cxx @@ -13,6 +13,7 @@ /// \brief analysis of UPC vector meson production ROF by ROF /// /// \author Guillermo Contreras (jesus.guillermo.contreras.nuno@cern.ch), Czech Technical University in Prague +/// \author Cesar Omar Ramirez Alvarez (cesar.ramirez@cern.ch), Autonomous University of Puebla #include "Common/CCDB/EventSelectionParams.h" #include "Common/DataModel/EventSelection.h" @@ -41,6 +42,7 @@ #include #include +#include #include #include #include @@ -48,6 +50,7 @@ #include #include #include +#include #include #include @@ -57,11 +60,14 @@ using namespace o2::framework::expressions; using BCsTSsSels = soa::Join; using ColSels = soa::Join; +using ColSelsMc = soa::Join; using TRKs = soa::Join; +using TRKsMc = soa::Join; using ColSel = ColSels::iterator; +using ColSelMc = ColSelsMc::iterator; namespace o2::aod { @@ -74,6 +80,7 @@ DECLARE_SOA_COLUMN(PosY, posY, float); DECLARE_SOA_COLUMN(PosZ, posZ, float); DECLARE_SOA_COLUMN(Chi2, chi2, float); DECLARE_SOA_COLUMN(LocalBC, localBC, int); +DECLARE_SOA_COLUMN(NearestBCB, nearestBCB, int); DECLARE_SOA_COLUMN(LocalTF, localTF, int); DECLARE_SOA_COLUMN(LocalROF, localROF, int); DECLARE_SOA_COLUMN(UpcFlag, upcFlag, int); @@ -140,9 +147,49 @@ DECLARE_SOA_COLUMN(HasTof3, hasTof3, int); DECLARE_SOA_COLUMN(HasTof4, hasTof4, int); } // namespace datarows +namespace mcgen +{ +DECLARE_SOA_COLUMN(McMotherPdg, mcMotherPdg, int); +DECLARE_SOA_COLUMN(McMotherPt, mcMotherPt, float); +DECLARE_SOA_COLUMN(McMotherPhi, mcMotherPhi, float); +DECLARE_SOA_COLUMN(McMotherMass, mcMotherMass, float); +DECLARE_SOA_COLUMN(McMotherRapidity, mcMotherRapidity, float); + +DECLARE_SOA_COLUMN(McPdg1, mcPdg1, int); +DECLARE_SOA_COLUMN(McPt1, mcPt1, float); +DECLARE_SOA_COLUMN(McEta1, mcEta1, float); +DECLARE_SOA_COLUMN(McPhi1, mcPhi1, float); +DECLARE_SOA_COLUMN(McQ1, mcQ1, int); +DECLARE_SOA_COLUMN(IsPrimary1, isPrimary1, int); +DECLARE_SOA_COLUMN(McPdg2, mcPdg2, int); +DECLARE_SOA_COLUMN(McPt2, mcPt2, float); +DECLARE_SOA_COLUMN(McEta2, mcEta2, float); +DECLARE_SOA_COLUMN(McPhi2, mcPhi2, float); +DECLARE_SOA_COLUMN(McQ2, mcQ2, int); +DECLARE_SOA_COLUMN(IsPrimary2, isPrimary2, int); +DECLARE_SOA_COLUMN(McPdg3, mcPdg3, int); +DECLARE_SOA_COLUMN(McPt3, mcPt3, float); +DECLARE_SOA_COLUMN(McEta3, mcEta3, float); +DECLARE_SOA_COLUMN(McPhi3, mcPhi3, float); +DECLARE_SOA_COLUMN(McQ3, mcQ3, int); +DECLARE_SOA_COLUMN(IsPrimary3, isPrimary3, int); +DECLARE_SOA_COLUMN(McPdg4, mcPdg4, int); +DECLARE_SOA_COLUMN(McPt4, mcPt4, float); +DECLARE_SOA_COLUMN(McEta4, mcEta4, float); +DECLARE_SOA_COLUMN(McPhi4, mcPhi4, float); +DECLARE_SOA_COLUMN(McQ4, mcQ4, int); +DECLARE_SOA_COLUMN(IsPrimary4, isPrimary4, int); + +DECLARE_SOA_COLUMN(McPosX, mcPosX, float); +DECLARE_SOA_COLUMN(McPosY, mcPosY, float); +DECLARE_SOA_COLUMN(McPosZ, mcPosZ, float); +DECLARE_SOA_COLUMN(McLocalBC, mcLocalBC, int); +DECLARE_SOA_COLUMN(RecoIndex, recoIndex, int); +} // namespace mcgen + DECLARE_SOA_TABLE(TwoTrkTable, "AOD", "TWOTRKTABLE", datarows::RunNumber, datarows::PosX, datarows::PosY, datarows::PosZ, datarows::Chi2, - datarows::LocalBC, datarows::LocalTF, datarows::LocalROF, datarows::UpcFlag, + datarows::LocalBC, datarows::NearestBCB, datarows::LocalTF, datarows::LocalROF, datarows::UpcFlag, datarows::AmplitudeFT0A, datarows::AmplitudeFT0C, datarows::AmplitudeFV0A, datarows::AmplitudeFDDA, datarows::AmplitudeFDDC, datarows::TimeFT0A, datarows::TimeFT0C, datarows::TimeFV0A, datarows::TimeFDDA, datarows::TimeFDDC, datarows::ChannelsFT0A, datarows::ChannelsFT0C, datarows::ChannelsFV0A, datarows::ChannelsFDDA, datarows::ChannelsFDDC, @@ -152,7 +199,7 @@ DECLARE_SOA_TABLE(TwoTrkTable, "AOD", "TWOTRKTABLE", datarows::HasTof1, datarows::HasTof2); DECLARE_SOA_TABLE(FourTrkTable, "AOD", "FOURTRKTABLE", datarows::RunNumber, datarows::PosX, datarows::PosY, datarows::PosZ, datarows::Chi2, - datarows::LocalBC, datarows::LocalTF, datarows::LocalROF, datarows::UpcFlag, + datarows::LocalBC, datarows::NearestBCB, datarows::LocalTF, datarows::LocalROF, datarows::UpcFlag, datarows::AmplitudeFT0A, datarows::AmplitudeFT0C, datarows::AmplitudeFV0A, datarows::AmplitudeFDDA, datarows::AmplitudeFDDC, datarows::TimeFT0A, datarows::TimeFT0C, datarows::TimeFV0A, datarows::TimeFDDA, datarows::TimeFDDC, datarows::ChannelsFT0A, datarows::ChannelsFT0C, datarows::ChannelsFV0A, datarows::ChannelsFDDA, datarows::ChannelsFDDC, @@ -162,6 +209,32 @@ DECLARE_SOA_TABLE(FourTrkTable, "AOD", "FOURTRKTABLE", datarows::Pt3, datarows::Eta3, datarows::Phi3, datarows::Q3, datarows::PidPion3, datarows::PidElectron3, datarows::PidKaon3, datarows::PidProton3, datarows::Pt4, datarows::Eta4, datarows::Phi4, datarows::Q4, datarows::PidPion4, datarows::PidElectron4, datarows::PidKaon4, datarows::PidProton4, datarows::HasTof1, datarows::HasTof2, datarows::HasTof3, datarows::HasTof4); + +DECLARE_SOA_TABLE(TwoTrkRecGenTable, "AOD", "TWOTRKRECGEN", + mcgen::McMotherPdg, mcgen::McMotherPt, mcgen::McMotherPhi, mcgen::McMotherMass, mcgen::McMotherRapidity, + mcgen::McPdg1, mcgen::McPt1, mcgen::McEta1, mcgen::McPhi1, mcgen::McQ1, mcgen::IsPrimary1, + mcgen::McPdg2, mcgen::McPt2, mcgen::McEta2, mcgen::McPhi2, mcgen::McQ2, mcgen::IsPrimary2, + mcgen::McPosX, mcgen::McPosY, mcgen::McPosZ, mcgen::McLocalBC, mcgen::RecoIndex); +DECLARE_SOA_TABLE(FourTrkRecGenTable, "AOD", "FOURTRKRECGEN", + mcgen::McMotherPdg, mcgen::McMotherPt, mcgen::McMotherPhi, mcgen::McMotherMass, mcgen::McMotherRapidity, + mcgen::McPdg1, mcgen::McPt1, mcgen::McEta1, mcgen::McPhi1, mcgen::McQ1, mcgen::IsPrimary1, + mcgen::McPdg2, mcgen::McPt2, mcgen::McEta2, mcgen::McPhi2, mcgen::McQ2, mcgen::IsPrimary2, + mcgen::McPdg3, mcgen::McPt3, mcgen::McEta3, mcgen::McPhi3, mcgen::McQ3, mcgen::IsPrimary3, + mcgen::McPdg4, mcgen::McPt4, mcgen::McEta4, mcgen::McPhi4, mcgen::McQ4, mcgen::IsPrimary4, + mcgen::McPosX, mcgen::McPosY, mcgen::McPosZ, mcgen::McLocalBC, mcgen::RecoIndex); + +DECLARE_SOA_TABLE(TwoTrkGenTable, "AOD", "TWOTRKGEN", + mcgen::McMotherPdg, mcgen::McMotherPt, mcgen::McMotherPhi, mcgen::McMotherMass, mcgen::McMotherRapidity, + mcgen::McPdg1, mcgen::McPt1, mcgen::McEta1, mcgen::McPhi1, mcgen::McQ1, mcgen::IsPrimary1, + mcgen::McPdg2, mcgen::McPt2, mcgen::McEta2, mcgen::McPhi2, mcgen::McQ2, mcgen::IsPrimary2, + mcgen::McPosX, mcgen::McPosY, mcgen::McPosZ, mcgen::McLocalBC); +DECLARE_SOA_TABLE(FourTrkGenTable, "AOD", "FOURTRKGEN", + mcgen::McMotherPdg, mcgen::McMotherPt, mcgen::McMotherPhi, mcgen::McMotherMass, mcgen::McMotherRapidity, + mcgen::McPdg1, mcgen::McPt1, mcgen::McEta1, mcgen::McPhi1, mcgen::McQ1, mcgen::IsPrimary1, + mcgen::McPdg2, mcgen::McPt2, mcgen::McEta2, mcgen::McPhi2, mcgen::McQ2, mcgen::IsPrimary2, + mcgen::McPdg3, mcgen::McPt3, mcgen::McEta3, mcgen::McPhi3, mcgen::McQ3, mcgen::IsPrimary3, + mcgen::McPdg4, mcgen::McPt4, mcgen::McEta4, mcgen::McPhi4, mcgen::McQ4, mcgen::IsPrimary4, + mcgen::McPosX, mcgen::McPosY, mcgen::McPosZ, mcgen::McLocalBC); } // namespace o2::aod struct UpcVmRof { @@ -169,6 +242,10 @@ struct UpcVmRof { // output Produces twoTrkTable; Produces fourTrkTable; + Produces twoTrkRecGenTable; + Produces fourTrkRecGenTable; + Produces twoTrkGenTable; + Produces fourTrkGenTable; // services Service ccdb{}; // access to database @@ -188,6 +265,9 @@ struct UpcVmRof { std::bitset bcPatternB; std::bitset bcPatternC; std::vector bcbIdx; + std::vector vecNearestBCB; + std::vector vecDistNearestBCB; + int nbcB = 0; // variables to store ITS ROF info @@ -230,6 +310,7 @@ struct UpcVmRof { Configurable minTrkTpcClusters{"minTrkTpcClusters", 70.0, "minimum number of TPC clusters associated to the track"}; Configurable maxTrkDcaZ{"maxTrkDcaZ", 2.0, "max DCA in z of track to vtx (cm)"}; Configurable tfPerBin{"tfPerBin", 10000, "timeframes per bin 1e4 means some 28 s"}; + Configurable genId{"genId", -1, "generator ID; -1 = no filter"}; //-------------------------------------------------------------------------------- // get ITS ROF info @@ -242,6 +323,63 @@ struct UpcVmRof { rofPerOrbit = static_cast(o2::constants::lhc::LHCMaxBunches / rofLength); } + //-------------------------------------------------------------------------------- + // find nearest bcb for this filling scheme + void setNearestBCB() + { + vecNearestBCB.clear(); + vecDistNearestBCB.clear(); + for (int i = 0; i < o2::constants::lhc::LHCMaxBunches; i++) { + // extend vector + vecNearestBCB.push_back(-1); + vecDistNearestBCB.push_back(-1); + // are we in a bc-b? + if (bcPatternB.test(i)) { + vecNearestBCB[i] = i; + vecDistNearestBCB[i] = 0; + continue; + } + // find nearest previous bc-b + int nearestLeft = -1; + for (int j = 1; j < o2::constants::lhc::LHCMaxBunches; j++) { + int k = i - j; + if (k < 0) + k = o2::constants::lhc::LHCMaxBunches + k; + if (bcPatternB.test(k)) { + nearestLeft = k; + break; + } + } // end nearestLeft loop + + // find nearest next bc-b + int nearesRight = -1; + for (int j = 1; j < o2::constants::lhc::LHCMaxBunches; j++) { + int k = i + j; + if (k > (o2::constants::lhc::LHCMaxBunches - 1)) + k = k - o2::constants::lhc::LHCMaxBunches; + if (bcPatternB.test(k)) { + nearesRight = k; + break; + } + } // end nearesRight loop + + // find nearest bc-b + int dLeft = i - nearestLeft; + // cppcheck-suppress knownConditionTrueFalse + if (dLeft < 0) + dLeft += o2::constants::lhc::LHCMaxBunches; + int dRight = nearesRight - i; + // cppcheck-suppress knownConditionTrueFalse + if (dRight < 0) + dRight += o2::constants::lhc::LHCMaxBunches; + int dMin = std::min(dLeft, dRight); + int nearest = ((dMin == dLeft) ? nearestLeft : nearesRight); + vecNearestBCB[i] = nearest; + int dist = ((dMin == dLeft) ? -dLeft : dRight); + vecDistNearestBCB[i] = dist; + } // end loop over bc + } + //-------------------------------------------------------------------------------- // get filling scheme void getFillingScheme() @@ -460,6 +598,109 @@ struct UpcVmRof { ccdb->setFatalWhenNull(false); } // end init() + //-------------------------------------------------------------------------------- + // helper function for MC: pt + static float mcPt(float px, float py) { return std::sqrt(px * px + py * py); } + + //-------------------------------------------------------------------------------- + // helper function for MC: pseudorapidity + static float mcEta(float px, float py, float pz) + { + float p = std::sqrt(px * px + py * py + pz * pz); + if (std::abs(p - std::abs(pz)) > 1e-10) { + return 0.5f * std::log((p + pz) / (p - pz)); + } + return 0.0f; + } + + //-------------------------------------------------------------------------------- + // helper function for MC: azimuth + static float mcPhi(float px, float py) + { + if (std::abs(px) > 1e-10 || std::abs(py) > 1e-10) { + return std::atan2(py, px); + } + return 0.0f; + } + + //-------------------------------------------------------------------------------- + // helper function for MC: rapidity + static float mcRapidity(float px, float py, float pz, float mass) + { + float energy = std::sqrt(px * px + py * py + pz * pz + mass * mass); + if (std::abs(energy - std::abs(pz)) > 1e-10) { + return 0.5f * std::log((energy + pz) / (energy - pz)); + } + return 0.0f; + } + + //-------------------------------------------------------------------------------- + // mother search + struct MotherInfo { + bool found = false; + int pdg = 0; + float px = 0, py = 0, pz = 0, mass = 0; + }; + + //-------------------------------------------------------------------------------- + // checks if all given daughters share the same direct mother: found=false if not + template + MotherInfo findCommonMother(aod::McParticles const& mcParticles, const std::vector& daughters) + { + MotherInfo info; + if (daughters.empty() || !daughters[0].has_mothers()) { + return info; + } + int64_t motherIdx = daughters[0].mothersIds()[0]; + for (size_t i = 1; i < daughters.size(); i++) { + if (!daughters[i].has_mothers() || daughters[i].mothersIds()[0] != motherIdx) { + return info; + } + } + auto mother = mcParticles.rawIteratorAt(motherIdx); + info.found = true; + info.pdg = mother.pdgCode(); + info.px = mother.px(); + info.py = mother.py(); + info.pz = mother.pz(); + info.mass = std::sqrt(mother.e() * mother.e() - mother.px() * mother.px() - mother.py() * mother.py() - mother.pz() * mother.pz()); + return info; + } + + //-------------------------------------------------------------------------------- + // for processMcGenAll: same mother search. + template + MotherInfo findCommonMotherFromIters(/* aod::McParticles const& mcParticles, */ const std::vector& daughters) + { + MotherInfo info; + if (daughters.empty() || !daughters[0].has_mothers()) { + return info; + } + auto firstMothers = daughters[0].template mothers_as(); + if (firstMothers.begin() == firstMothers.end()) { + return info; + } + auto motherIt = firstMothers.begin(); + int64_t motherGlobalIndex = motherIt->globalIndex(); + for (size_t i = 1; i < daughters.size(); i++) { + if (!daughters[i].has_mothers()) { + return info; + } + auto iMothers = daughters[i].template mothers_as(); + if (iMothers.begin() == iMothers.end() || iMothers.begin()->globalIndex() != motherGlobalIndex) { + return info; + } + } + const auto& mother = *motherIt; + info.found = true; + info.pdg = mother.pdgCode(); + info.px = mother.px(); + info.py = mother.py(); + info.pz = mother.pz(); + info.mass = std::sqrt(mother.e() * mother.e() - mother.px() * mother.px() - mother.py() * mother.py() - mother.pz() * mother.pz()); + return info; + } + //-------------------------------------------------------------------------------- // process BCs to get trigger information void processBCs(BCsTSsSels const& bcs) @@ -528,11 +769,17 @@ struct UpcVmRof { PROCESS_SWITCH(UpcVmRof, processBCs, "get BCs and trigger information", true); //-------------------------------------------------------------------------------- - // get collision information - void processCols(ColSel const& col, BCsTSsSels const&, TRKs const& tracks, - aod::FV0As const&, aod::FT0s const&, aod::FDDs const&, - aod::Zdcs const&) + // shared collision+track selection; used by both real data and MC + template + std::vector fillCollisionTables( + ColType const& col, BCsTSsSels const& bcs, TracksType const& tracks, + aod::FV0As const&, aod::FT0s const&, aod::FDDs const&, aod::Zdcs const&, + bool& outIsTwoBody, bool& outIsFourBody) { + outIsTwoBody = false; + outIsFourBody = false; + std::vector selTrks; + // get info for this bc auto bc = col.template foundBC_as(); if (runNumberCol != bc.runNumber()) { // new run @@ -540,6 +787,7 @@ struct UpcVmRof { getRunInfo(runNumberCol); getFillingScheme(); addColHistos(runNumberCol); + setNearestBCB(); } int64_t thisBC = getBcWithinOrbit(bc.globalBC()); int64_t thisTF = getTimeFrame(bc.globalBC()); @@ -547,18 +795,14 @@ struct UpcVmRof { // select collision if (!checkColFlags(col, runNumberCol)) { - return; + return selTrks; } - // accept only -B bcs - if (!bcPatternB.test(thisBC)) { - return; - } colTH1Pointers[Form("col/%d/colSel_H", runNumberCol)]->Fill(10); // select on zVtx if (std::abs(col.posZ()) > maxAbsPosZ) { - return; + return selTrks; } colTH1Pointers[Form("col/%d/colSel_H", runNumberCol)]->Fill(11); @@ -566,7 +810,7 @@ struct UpcVmRof { bool isTwoContributors = (col.numContrib() == NTrksTwoBody); bool isFourContributors = (col.numContrib() == NTrksFourBody); if (!isTwoContributors && !isFourContributors) { - return; + return selTrks; } if (isTwoContributors) { colTH1Pointers[Form("col/%d/colSel_H", runNumberCol)]->Fill(16); @@ -576,7 +820,6 @@ struct UpcVmRof { } // select tracks - std::vector selTrks; colTH1Pointers[Form("col/%d/trkSel_H", runNumberCol)]->Fill(0); for (const auto& track : tracks) { if (!track.isPVContributor()) { @@ -623,7 +866,7 @@ struct UpcVmRof { bool isTwoBody = isTwoContributors && (selTrks.size() == NTrksTwoBody); bool isFourBody = isFourContributors && (selTrks.size() == NTrksFourBody); if (!isTwoBody && !isFourBody) { - return; + return selTrks; } // selected events @@ -634,6 +877,40 @@ struct UpcVmRof { colTH1Pointers[Form("col/%d/colSel_H", runNumberCol)]->Fill(19); } + // find the bc row for the nearest bc if different from the current one + bool foundBCB = true; + auto nearbcb = bcs.iteratorAt(bc.globalIndex()); + if (vecDistNearestBCB[thisBC] > 0) { + foundBCB = false; + auto gidxBCB = bc.globalIndex() + vecDistNearestBCB[thisBC]; + if (gidxBCB <= bcs.iteratorAt(bcs.size() - 1).globalIndex()) { + while (nearbcb != bcs.end()) { + ++nearbcb; + if (nearbcb.globalIndex() > gidxBCB) { + break; + } else if (getBcWithinOrbit(nearbcb.globalBC()) == vecNearestBCB[thisBC]) { + foundBCB = true; + break; + } + } // end while + } + } // end search in future direction + if (vecDistNearestBCB[thisBC] < 0) { + foundBCB = false; + auto gidxBCB = bc.globalIndex() + vecDistNearestBCB[thisBC]; + if (gidxBCB >= bcs.iteratorAt(0).globalIndex()) { + while (nearbcb != bcs.iteratorAt(0)) { + --nearbcb; + if (nearbcb.globalIndex() < gidxBCB) { + break; + } else if (getBcWithinOrbit(nearbcb.globalBC()) == vecNearestBCB[thisBC]) { + foundBCB = true; + break; + } + } // end while + } + } // end search in past direction + // FT0 selection float aFT0A = 0; float aFT0C = 0; @@ -641,34 +918,34 @@ struct UpcVmRof { float tFT0C = 33; // default time to mark events without FT0 info int nFT0A = 0; int nFT0C = 0; - if (bc.has_foundFT0()) { + if (foundBCB && nearbcb.has_foundFT0()) { // a side - if (bc.foundFT0().isValidTimeA()) { // valid time - tFT0A = bc.foundFT0().timeA(); + if (nearbcb.foundFT0().isValidTimeA()) { // valid time + tFT0A = nearbcb.foundFT0().timeA(); if (std::abs(tFT0A) > maxAbsTimeFT0) { - return; + return selTrks; } colTH1Pointers[Form("col/%d/colSel_H", runNumberCol)]->Fill(12); - aFT0A = bc.foundFT0().sumAmpA(); + aFT0A = nearbcb.foundFT0().sumAmpA(); if (aFT0A > maxAmpFT0) { - return; + return selTrks; } colTH1Pointers[Form("col/%d/colSel_H", runNumberCol)]->Fill(13); - nFT0A = (bc.foundFT0().amplitudeA()).size(); + nFT0A = (nearbcb.foundFT0().amplitudeA()).size(); } // a side // c side - if (bc.foundFT0().isValidTimeC()) { // valid time - tFT0C = bc.foundFT0().timeC(); + if (nearbcb.foundFT0().isValidTimeC()) { // valid time + tFT0C = nearbcb.foundFT0().timeC(); if (std::abs(tFT0C) > maxAbsTimeFT0) { - return; + return selTrks; } colTH1Pointers[Form("col/%d/colSel_H", runNumberCol)]->Fill(14); - aFT0C = bc.foundFT0().sumAmpC(); + aFT0C = nearbcb.foundFT0().sumAmpC(); if (aFT0C > maxAmpFT0) { - return; + return selTrks; } colTH1Pointers[Form("col/%d/colSel_H", runNumberCol)]->Fill(15); - nFT0C = (bc.foundFT0().amplitudeC()).size(); + nFT0C = (nearbcb.foundFT0().amplitudeC()).size(); } // c side } // FT0 selection @@ -684,9 +961,9 @@ struct UpcVmRof { float aFV0A = 0; float tFV0A = 33; // default time to mark events without FV0 info int nFV0A = 0; - if (bc.has_foundFV0()) { - tFV0A = bc.foundFV0().time(); - auto v = bc.foundFV0().amplitude(); + if (foundBCB && nearbcb.has_foundFV0()) { + tFV0A = nearbcb.foundFV0().time(); + auto v = nearbcb.foundFV0().amplitude(); aFV0A = std::accumulate(v.begin(), v.end(), 0.f); nFV0A = v.size(); } // FV0A info @@ -698,9 +975,9 @@ struct UpcVmRof { float aFDDC = 0; float tFDDC = 33; // default time to mark events without FDD info int nFDDC = 0; - if (bc.has_foundFDD()) { - tFDDA = bc.foundFDD().timeA(); - auto vA = bc.foundFDD().chargeA(); + if (foundBCB && nearbcb.has_foundFDD()) { + tFDDA = nearbcb.foundFDD().timeA(); + auto vA = nearbcb.foundFDD().chargeA(); // channelPairs = {{0, 4}, {1, 5}, {2, 6}, {3, 7}}; if (vA[0] > 0 && vA[4] > 0) { aFDDA += 0.5 * (vA[0] + vA[4]); @@ -718,8 +995,8 @@ struct UpcVmRof { aFDDA += 0.5 * (vA[3] + vA[7]); nFDDA++; } - tFDDC = bc.foundFDD().timeC(); - auto vC = bc.foundFDD().chargeC(); + tFDDC = nearbcb.foundFDD().timeC(); + auto vC = nearbcb.foundFDD().chargeC(); // channelPairs = {{0, 4}, {1, 5}, {2, 6}, {3, 7}}; if (vC[0] > 0 && vC[4] > 0) { aFDDC += 0.5 * (vC[0] + vC[4]); @@ -744,11 +1021,11 @@ struct UpcVmRof { float tZNC = -999; // default time to mark events without ZN info float eZNA = -999; float eZNC = -999; - if (bc.has_zdc()) { - tZNA = (bc.zdc()).timeZNA(); - tZNC = (bc.zdc()).timeZNC(); - eZNA = (bc.zdc()).energyCommonZNA(); - eZNC = (bc.zdc()).energyCommonZNC(); + if (foundBCB && nearbcb.has_zdc()) { + tZNA = (nearbcb.zdc()).timeZNA(); + tZNC = (nearbcb.zdc()).timeZNC(); + eZNA = (nearbcb.zdc()).energyCommonZNA(); + eZNC = (nearbcb.zdc()).energyCommonZNC(); if (!std::isfinite(tZNA)) { tZNA = -999; } @@ -774,7 +1051,7 @@ struct UpcVmRof { tof[1] = 1; } colTH1Pointers[Form("col/%d/twoTrkTF_H", runNumberCol)]->Fill(thisTF); - twoTrkTable(runNumberCol, col.posX(), col.posY(), col.posZ(), col.chi2(), thisBC, thisTF, thisROF, recoFlag, + twoTrkTable(runNumberCol, col.posX(), col.posY(), col.posZ(), col.chi2(), thisBC, vecNearestBCB[thisBC], thisTF, thisROF, recoFlag, aFT0A, aFT0C, aFV0A, aFDDA, aFDDC, tFT0A, tFT0C, tFV0A, tFDDA, tFDDC, nFT0A, nFT0C, nFV0A, nFDDA, nFDDC, eZNA, eZNC, tZNA, tZNC, selTrks[0].pt(), selTrks[0].eta(), selTrks[0].phi(), selTrks[0].sign(), @@ -797,7 +1074,7 @@ struct UpcVmRof { tof[3] = 1; } colTH1Pointers[Form("col/%d/fourTrkTF_H", runNumberCol)]->Fill(thisTF); - fourTrkTable(runNumberCol, col.posX(), col.posY(), col.posZ(), col.chi2(), thisBC, thisTF, thisROF, recoFlag, + fourTrkTable(runNumberCol, col.posX(), col.posY(), col.posZ(), col.chi2(), thisBC, vecNearestBCB[thisBC], thisTF, thisROF, recoFlag, aFT0A, aFT0C, aFV0A, aFDDA, aFDDC, tFT0A, tFT0C, tFV0A, tFDDA, tFDDC, nFT0A, nFT0C, nFV0A, nFDDA, nFDDC, eZNA, eZNC, tZNA, tZNC, selTrks[0].pt(), selTrks[0].eta(), selTrks[0].phi(), selTrks[0].sign(), @@ -811,12 +1088,287 @@ struct UpcVmRof { tof[0], tof[1], tof[2], tof[3]); } - } // end processCol - PROCESS_SWITCH(UpcVmRof, processCols, "get collisions and track information", true); + outIsTwoBody = isTwoBody; + outIsFourBody = isFourBody; + return selTrks; + } // end fillCollisionTables + + //-------------------------------------------------------------------------------- + // process real data or MC at reconstruction level + void processDataCols(ColSel const& col, BCsTSsSels const& bcs, TRKs const& tracks, + aod::FV0As const& fv0s, aod::FT0s const& ft0s, aod::FDDs const& fdds, aod::Zdcs const& zdcs) + { + bool isTwoBody, isFourBody; + fillCollisionTables(col, bcs, tracks, fv0s, ft0s, fdds, zdcs, isTwoBody, isFourBody); + } // end processDataCols + PROCESS_SWITCH(UpcVmRof, processDataCols, "process real data or mc reco, no mc truth", true); + + //-------------------------------------------------------------------------------- + // MC only: fills the same reconstruction tables and the generated information of the reconstructed event + void processMcCols(ColSelMc const& col, BCsTSsSels const& bcs, TRKsMc const& tracks, + aod::FV0As const& fv0s, aod::FT0s const& ft0s, aod::FDDs const& fdds, aod::Zdcs const& zdcs, + aod::McCollisions const&, aod::McParticles const& mcParticles) + { + if (genId != -1) { + if (!col.has_mcCollision() || col.mcCollision().getGeneratorId() != genId) { + return; + } + } + + bool isTwoBody, isFourBody; + auto selTrks = fillCollisionTables(col, bcs, tracks, fv0s, ft0s, fdds, zdcs, isTwoBody, isFourBody); + if (!isTwoBody && !isFourBody) { + return; // event was not written to the reco tables either + } + + // MC gen + if (isTwoBody) { + auto mcPart1 = selTrks[0].has_mcParticle() ? std::optional(selTrks[0].mcParticle()) : std::nullopt; + auto mcPart2 = selTrks[1].has_mcParticle() ? std::optional(selTrks[1].mcParticle()) : std::nullopt; + + int pdg1 = 0, pdg2 = 0, sign1 = 0, sign2 = 0, isPrim1 = 0, isPrim2 = 0; + float pt1 = -999, eta1 = -999, phi1 = -999, pt2 = -999, eta2 = -999, phi2 = -999; + int motherPdg = 0; + float motherPt = -999, motherPhi = -999, motherMass = -999, motherRap = -999; + float mcPosX = -999, mcPosY = -999, mcPosZ = -999; + int mcLocalBC = -999; + + if (mcPart1 && mcPart2) { + pdg1 = mcPart1->pdgCode(); + pt1 = mcPt(mcPart1->px(), mcPart1->py()); + eta1 = mcEta(mcPart1->px(), mcPart1->py(), mcPart1->pz()); + phi1 = mcPhi(mcPart1->px(), mcPart1->py()); + sign1 = (pdg1 > 0) ? 1 : (pdg1 < 0) ? -1 + : 0; + isPrim1 = mcPart1->isPhysicalPrimary() ? 1 : 0; + + pdg2 = mcPart2->pdgCode(); + pt2 = mcPt(mcPart2->px(), mcPart2->py()); + eta2 = mcEta(mcPart2->px(), mcPart2->py(), mcPart2->pz()); + phi2 = mcPhi(mcPart2->px(), mcPart2->py()); + sign2 = (pdg2 > 0) ? 1 : (pdg2 < 0) ? -1 + : 0; + isPrim2 = mcPart2->isPhysicalPrimary() ? 1 : 0; + + auto mcCollision = mcPart1->mcCollision(); + mcPosX = mcCollision.posX(); + mcPosY = mcCollision.posY(); + mcPosZ = mcCollision.posZ(); + mcLocalBC = static_cast(getBcWithinOrbit(mcCollision.template bc_as().globalBC())); + + float sumPx = mcPart1->px() + mcPart2->px(); + float sumPy = mcPart1->py() + mcPart2->py(); + float sumPz = mcPart1->pz() + mcPart2->pz(); + float sumE = mcPart1->e() + mcPart2->e(); + motherPt = mcPt(sumPx, sumPy); + motherPhi = mcPhi(sumPx, sumPy); + motherMass = std::sqrt(sumE * sumE - sumPx * sumPx - sumPy * sumPy - sumPz * sumPz); + motherRap = mcRapidity(sumPx, sumPy, sumPz, motherMass); + + MotherInfo mi = findCommonMother(mcParticles, std::vector{*mcPart1, *mcPart2}); + if (mi.found) { + motherPdg = mi.pdg; + } + } + + twoTrkRecGenTable( + motherPdg, motherPt, motherPhi, motherMass, motherRap, + pdg1, pt1, eta1, phi1, sign1, isPrim1, + pdg2, pt2, eta2, phi2, sign2, isPrim2, + mcPosX, mcPosY, mcPosZ, mcLocalBC, + twoTrkTable.lastIndex() // row of TwoTrkTable this mc info belongs to + ); + } + if (isFourBody) { + auto mcPart1 = selTrks[0].has_mcParticle() ? std::optional(selTrks[0].mcParticle()) : std::nullopt; + auto mcPart2 = selTrks[1].has_mcParticle() ? std::optional(selTrks[1].mcParticle()) : std::nullopt; + auto mcPart3 = selTrks[2].has_mcParticle() ? std::optional(selTrks[2].mcParticle()) : std::nullopt; + auto mcPart4 = selTrks[3].has_mcParticle() ? std::optional(selTrks[3].mcParticle()) : std::nullopt; + + int pdg1 = 0, pdg2 = 0, pdg3 = 0, pdg4 = 0; + int sign1 = 0, sign2 = 0, sign3 = 0, sign4 = 0; + int isPrim1 = 0, isPrim2 = 0, isPrim3 = 0, isPrim4 = 0; + float pt1 = -999, eta1 = -999, phi1 = -999; + float pt2 = -999, eta2 = -999, phi2 = -999; + float pt3 = -999, eta3 = -999, phi3 = -999; + float pt4 = -999, eta4 = -999, phi4 = -999; + int motherPdg = 0; + float motherPt = -999, motherPhi = -999, motherMass = -999, motherRap = -999; + float mcPosX = -999, mcPosY = -999, mcPosZ = -999; + int mcLocalBC = -999; + + if (mcPart1 && mcPart2 && mcPart3 && mcPart4) { + pdg1 = mcPart1->pdgCode(); + pt1 = mcPt(mcPart1->px(), mcPart1->py()); + eta1 = mcEta(mcPart1->px(), mcPart1->py(), mcPart1->pz()); + phi1 = mcPhi(mcPart1->px(), mcPart1->py()); + sign1 = (pdg1 > 0) ? 1 : (pdg1 < 0) ? -1 + : 0; + isPrim1 = mcPart1->isPhysicalPrimary() ? 1 : 0; + + pdg2 = mcPart2->pdgCode(); + pt2 = mcPt(mcPart2->px(), mcPart2->py()); + eta2 = mcEta(mcPart2->px(), mcPart2->py(), mcPart2->pz()); + phi2 = mcPhi(mcPart2->px(), mcPart2->py()); + sign2 = (pdg2 > 0) ? 1 : (pdg2 < 0) ? -1 + : 0; + isPrim2 = mcPart2->isPhysicalPrimary() ? 1 : 0; + + pdg3 = mcPart3->pdgCode(); + pt3 = mcPt(mcPart3->px(), mcPart3->py()); + eta3 = mcEta(mcPart3->px(), mcPart3->py(), mcPart3->pz()); + phi3 = mcPhi(mcPart3->px(), mcPart3->py()); + sign3 = (pdg3 > 0) ? 1 : (pdg3 < 0) ? -1 + : 0; + isPrim3 = mcPart3->isPhysicalPrimary() ? 1 : 0; + + pdg4 = mcPart4->pdgCode(); + pt4 = mcPt(mcPart4->px(), mcPart4->py()); + eta4 = mcEta(mcPart4->px(), mcPart4->py(), mcPart4->pz()); + phi4 = mcPhi(mcPart4->px(), mcPart4->py()); + sign4 = (pdg4 > 0) ? 1 : (pdg4 < 0) ? -1 + : 0; + isPrim4 = mcPart4->isPhysicalPrimary() ? 1 : 0; + + auto mcCollision = mcPart1->mcCollision(); + mcPosX = mcCollision.posX(); + mcPosY = mcCollision.posY(); + mcPosZ = mcCollision.posZ(); + mcLocalBC = static_cast(getBcWithinOrbit(mcCollision.template bc_as().globalBC())); + + float sumPx = mcPart1->px() + mcPart2->px() + mcPart3->px() + mcPart4->px(); + float sumPy = mcPart1->py() + mcPart2->py() + mcPart3->py() + mcPart4->py(); + float sumPz = mcPart1->pz() + mcPart2->pz() + mcPart3->pz() + mcPart4->pz(); + float sumE = mcPart1->e() + mcPart2->e() + mcPart3->e() + mcPart4->e(); + motherPt = mcPt(sumPx, sumPy); + motherPhi = mcPhi(sumPx, sumPy); + motherMass = std::sqrt(sumE * sumE - sumPx * sumPx - sumPy * sumPy - sumPz * sumPz); + motherRap = mcRapidity(sumPx, sumPy, sumPz, motherMass); + + MotherInfo mi = findCommonMother(mcParticles, std::vector{*mcPart1, *mcPart2, *mcPart3, *mcPart4}); + if (mi.found) { + motherPdg = mi.pdg; + } + } + + fourTrkRecGenTable( + motherPdg, motherPt, motherPhi, motherMass, motherRap, + pdg1, pt1, eta1, phi1, sign1, isPrim1, + pdg2, pt2, eta2, phi2, sign2, isPrim2, + pdg3, pt3, eta3, phi3, sign3, isPrim3, + pdg4, pt4, eta4, phi4, sign4, isPrim4, + mcPosX, mcPosY, mcPosZ, mcLocalBC, + fourTrkTable.lastIndex() // row of FourTrkTable this mc info belongs to + ); + } + } // end processMcCols + PROCESS_SWITCH(UpcVmRof, processMcCols, "process mc truth matched to reco, mc only", false); + + //-------------------------------------------------------------------------------- + // Fill MC information for all generated events + void processMcGen(aod::McCollision const& mcCollision, BCsTSsSels const&, aod::McParticles const& mcParticles) + { + if (genId != -1 && mcCollision.getGeneratorId() != genId) { + return; + } + + // event counter + if (!colTH1Pointers["mcGen/genSel_H"]) { + colTH1Pointers["mcGen/genSel_H"] = colTH1Registry.add("mcGen/genSel_H", + "pure generator event counter; selID; Counter", + {HistType::kTH1D, {{3, -0.5, 2.5}}}); + } + colTH1Pointers["mcGen/genSel_H"]->Fill(0); + + std::vector primaries; + for (const auto& part : mcParticles) { + if (part.isPhysicalPrimary()) { + primaries.push_back(part); + } + } + + // local BC of the generated collision + int localBc = static_cast(getBcWithinOrbit(mcCollision.template bc_as().globalBC())); + + if (primaries.size() == NTrksTwoBody) { + auto p1 = primaries[0]; + auto p2 = primaries[1]; + + int pdg1 = p1.pdgCode(); + int pdg2 = p2.pdgCode(); + int sign1 = (pdg1 > 0) ? 1 : (pdg1 < 0) ? -1 + : 0; + int sign2 = (pdg2 > 0) ? 1 : (pdg2 < 0) ? -1 + : 0; + + // mother mass/momentum always from summed daughters + float sumPx = p1.px() + p2.px(); + float sumPy = p1.py() + p2.py(); + float sumPz = p1.pz() + p2.pz(); + float sumE = p1.e() + p2.e(); + float motherMass = std::sqrt(sumE * sumE - sumPx * sumPx - sumPy * sumPy - sumPz * sumPz); + + int motherPdg = 0; // 0 = no common mother confirmed in the gen tree + MotherInfo mi = findCommonMotherFromIters(/* mcParticles,*/ std::vector{p1, p2}); + if (mi.found) { + motherPdg = mi.pdg; + } + + colTH1Pointers["mcGen/genSel_H"]->Fill(1); + twoTrkGenTable( + motherPdg, mcPt(sumPx, sumPy), mcPhi(sumPx, sumPy), motherMass, mcRapidity(sumPx, sumPy, sumPz, motherMass), + pdg1, mcPt(p1.px(), p1.py()), mcEta(p1.px(), p1.py(), p1.pz()), mcPhi(p1.px(), p1.py()), sign1, 1, + pdg2, mcPt(p2.px(), p2.py()), mcEta(p2.px(), p2.py(), p2.pz()), mcPhi(p2.px(), p2.py()), sign2, 1, + mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), localBc); + } + + if (primaries.size() == NTrksFourBody) { + auto p1 = primaries[0]; + auto p2 = primaries[1]; + auto p3 = primaries[2]; + auto p4 = primaries[3]; + + int pdg1 = p1.pdgCode(); + int pdg2 = p2.pdgCode(); + int pdg3 = p3.pdgCode(); + int pdg4 = p4.pdgCode(); + int sign1 = (pdg1 > 0) ? 1 : (pdg1 < 0) ? -1 + : 0; + int sign2 = (pdg2 > 0) ? 1 : (pdg2 < 0) ? -1 + : 0; + int sign3 = (pdg3 > 0) ? 1 : (pdg3 < 0) ? -1 + : 0; + int sign4 = (pdg4 > 0) ? 1 : (pdg4 < 0) ? -1 + : 0; + + float sumPx = p1.px() + p2.px() + p3.px() + p4.px(); + float sumPy = p1.py() + p2.py() + p3.py() + p4.py(); + float sumPz = p1.pz() + p2.pz() + p3.pz() + p4.pz(); + float sumE = p1.e() + p2.e() + p3.e() + p4.e(); + float motherMass = std::sqrt(sumE * sumE - sumPx * sumPx - sumPy * sumPy - sumPz * sumPz); + + int motherPdg = 0; + MotherInfo mi = findCommonMotherFromIters(/* mcParticles, */ std::vector{p1, p2, p3, p4}); + if (mi.found) { + motherPdg = mi.pdg; + } + + colTH1Pointers["mcGen/genSel_H"]->Fill(2); + fourTrkGenTable( + motherPdg, mcPt(sumPx, sumPy), mcPhi(sumPx, sumPy), motherMass, mcRapidity(sumPx, sumPy, sumPz, motherMass), + pdg1, mcPt(p1.px(), p1.py()), mcEta(p1.px(), p1.py(), p1.pz()), mcPhi(p1.px(), p1.py()), sign1, 1, + pdg2, mcPt(p2.px(), p2.py()), mcEta(p2.px(), p2.py(), p2.pz()), mcPhi(p2.px(), p2.py()), sign2, 1, + pdg3, mcPt(p3.px(), p3.py()), mcEta(p3.px(), p3.py(), p3.pz()), mcPhi(p3.px(), p3.py()), sign3, 1, + pdg4, mcPt(p4.px(), p4.py()), mcEta(p4.px(), p4.py(), p4.pz()), mcPhi(p4.px(), p4.py()), sign4, 1, + mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), localBc); + } + } // end processMcGenAll + PROCESS_SWITCH(UpcVmRof, processMcGen, "process all generated collisions, mc only", false); }; // end of struct UpcVmRof WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) + { return WorkflowSpec{ adaptAnalysisTask(cfgc)};