diff --git a/PWGLF/Tasks/Nuspex/multiplicityPt.cxx b/PWGLF/Tasks/Nuspex/multiplicityPt.cxx index 8a8b5e3c895..6edc236a27f 100644 --- a/PWGLF/Tasks/Nuspex/multiplicityPt.cxx +++ b/PWGLF/Tasks/Nuspex/multiplicityPt.cxx @@ -84,7 +84,7 @@ struct MultiplicityPt { Configurable cfgCutEtaMax{"cfgCutEtaMax", 0.8f, "Max eta range for tracks"}; Configurable cfgCutEtaMin{"cfgCutEtaMin", -0.8f, "Min eta range for tracks"}; Configurable cfgCutNsigma{"cfgCutNsigma", 3.0f, "nsigma cut range for tracks"}; - + Configurable enablePIDHistograms{"enablePIDHistograms", true, "Enable PID histograms"}; Configurable useCustomTrackCuts{"useCustomTrackCuts", true, "Flag to use custom track cuts"}; Configurable itsPattern{"itsPattern", 0, "0 = Run3ITSibAny, 1 = Run3ITSallAny, 2 = Run3ITSall7Layers, 3 = Run3ITSibTwo"}; @@ -134,9 +134,9 @@ struct MultiplicityPt { Preslice perMCCol = aod::mcparticle::mcCollisionId; enum ParticleSpecies : int { - PartPion = 0, - PartKaon = 1, - PartProton = 2, + PartPion = 0, + PartKaon = 1, + PartProton = 2, }; int getMagneticField(uint64_t timestamp) @@ -156,8 +156,10 @@ struct MultiplicityPt { float getTransformedPhi(const float phi, const int charge, const float magField) const { float transformedPhi = phi; - if (magField < 0) transformedPhi = o2::constants::math::TwoPI - transformedPhi; - if (charge < 0) transformedPhi = o2::constants::math::TwoPI - transformedPhi; + if (magField < 0) + transformedPhi = o2::constants::math::TwoPI - transformedPhi; + if (charge < 0) + transformedPhi = o2::constants::math::TwoPI - transformedPhi; transformedPhi += o2::constants::math::PI / 18.0f; transformedPhi = std::fmod(transformedPhi, o2::constants::math::PI / 9.0f); return transformedPhi; @@ -166,18 +168,23 @@ struct MultiplicityPt { template bool passedPhiCut(const TrackType& track, float magField) const { - if (!applyPhiCut.value) return true; - if (track.pt() < pTthresholdPhiCut.value) return true; + if (!applyPhiCut.value) + return true; + if (track.pt() < pTthresholdPhiCut.value) + return true; float phi = track.phi(); int charge = track.sign(); - if (magField < 0) phi = o2::constants::math::TwoPI - phi; - if (charge < 0) phi = o2::constants::math::TwoPI - phi; + if (magField < 0) + phi = o2::constants::math::TwoPI - phi; + if (charge < 0) + phi = o2::constants::math::TwoPI - phi; phi += o2::constants::math::PI / 18.0f; phi = std::fmod(phi, o2::constants::math::PI / 9.0f); - if (phi < fphiCutHigh->Eval(track.pt()) && phi > fphiCutLow->Eval(track.pt())) return false; + if (phi < fphiCutHigh->Eval(track.pt()) && phi > fphiCutLow->Eval(track.pt())) + return false; return true; } @@ -187,20 +194,30 @@ struct MultiplicityPt { int count = 0; for (const auto& particle : particles) { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (!pdgParticle || pdgParticle->Charge() == 0.) continue; - if (!particle.isPhysicalPrimary()) continue; - if (std::abs(particle.eta()) > etaMax) continue; - if (particle.pt() < ptMin) continue; + if (!pdgParticle || pdgParticle->Charge() == 0.) + continue; + if (!particle.isPhysicalPrimary()) + continue; + if (std::abs(particle.eta()) > etaMax) + continue; + if (particle.pt() < ptMin) + continue; count++; } return count; } template - bool passedNClTPCFoundCut(const T& trk) const { return !nClTPCFoundCut.value || trk.tpcNClsFound() >= minTPCNClsFound.value; } + bool passedNClTPCFoundCut(const T& trk) const + { + return !nClTPCFoundCut.value || trk.tpcNClsFound() >= minTPCNClsFound.value; + } template - bool passedNClTPCPIDCut(const T& trk) const { return !nClTPCPIDCut.value || trk.tpcNClsPID() >= minTPCNClsPID.value; } + bool passedNClTPCPIDCut(const T& trk) const + { + return !nClTPCPIDCut.value || trk.tpcNClsPID() >= minTPCNClsPID.value; + } template bool passesCutWoDCA(TrackType const& track) const @@ -208,8 +225,10 @@ struct MultiplicityPt { if (useCustomTrackCuts.value) { for (int i = 0; i < static_cast(TrackSelection::TrackCuts::kNCuts); i++) { if (i == static_cast(TrackSelection::TrackCuts::kDCAxy) || - i == static_cast(TrackSelection::TrackCuts::kDCAz)) continue; - if (!customTrackCuts.IsSelected(track, static_cast(i))) return false; + i == static_cast(TrackSelection::TrackCuts::kDCAz)) + continue; + if (!customTrackCuts.IsSelected(track, static_cast(i))) + return false; } return true; } @@ -220,7 +239,8 @@ struct MultiplicityPt { bool passesDCAxyCut(TrackType const& track) const { if (useCustomTrackCuts.value) { - if (!passesCutWoDCA(track)) return false; + if (!passesCutWoDCA(track)) + return false; constexpr float DcaXYConst = 0.0105f; constexpr float DcaXYPtScale = 0.0350f; constexpr float DcaXYPtPower = 1.1f; @@ -233,13 +253,20 @@ struct MultiplicityPt { template bool passesTrackSelection(TrackType const& track, float magField = 0) const { - if (track.eta() < cfgCutEtaMin.value || track.eta() > cfgCutEtaMax.value) return false; - if (track.tpcChi2NCl() < minChi2PerClusterTPC.value || track.tpcChi2NCl() > maxChi2PerClusterTPC.value) return false; - if (!passesCutWoDCA(track)) return false; - if (!passesDCAxyCut(track)) return false; - if (!passedNClTPCFoundCut(track)) return false; - if (!passedNClTPCPIDCut(track)) return false; - if (!passedPhiCut(track, magField)) return false; + if (track.eta() < cfgCutEtaMin.value || track.eta() > cfgCutEtaMax.value) + return false; + if (track.tpcChi2NCl() < minChi2PerClusterTPC.value || track.tpcChi2NCl() > maxChi2PerClusterTPC.value) + return false; + if (!passesCutWoDCA(track)) + return false; + if (!passesDCAxyCut(track)) + return false; + if (!passedNClTPCFoundCut(track)) + return false; + if (!passedNClTPCPIDCut(track)) + return false; + if (!passedPhiCut(track, magField)) + return false; return true; } @@ -248,7 +275,7 @@ struct MultiplicityPt { { float nsigmaTPC = 0.f; float cutValue = 0.f; - + if (species == PartPion) { nsigmaTPC = track.tpcNSigmaPi(); cutValue = cfgCutNsigmaPi.value; @@ -259,7 +286,7 @@ struct MultiplicityPt { nsigmaTPC = track.tpcNSigmaPr(); cutValue = cfgCutNsigmaPr.value; } - + return std::abs(nsigmaTPC) < cutValue; } @@ -267,10 +294,14 @@ struct MultiplicityPt { bool isGoodPrimary(ParticleType const& particle) const { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (!pdgParticle || pdgParticle->Charge() == 0.) return false; - if (!particle.isPhysicalPrimary()) return false; - if (std::abs(particle.eta()) >= cfgCutEtaMax.value) return false; - if (particle.pt() < cfgTrkLowPtCut.value) return false; + if (!pdgParticle || pdgParticle->Charge() == 0.) + return false; + if (!particle.isPhysicalPrimary()) + return false; + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) + return false; + if (particle.pt() < cfgTrkLowPtCut.value) + return false; return true; } @@ -322,20 +353,22 @@ void MultiplicityPt::init(InitContext const&) std::vector centBinningStd = {0., 1., 5., 10., 15., 20., 30., 40., 50., 60., 70., 80., 90., 100.}; AxisSpec centAxis = {centBinningStd, "FT0M Centrality (%)"}; - + std::vector centBinningFine; - for (int i = 0; i <= CentBinMax; i++) centBinningFine.push_back(static_cast(i)); + for (int i = 0; i <= CentBinMax; i++) + centBinningFine.push_back(static_cast(i)); AxisSpec centFineAxis = {centBinningFine, "FT0M Centrality (%)"}; std::vector multBins; - for (int i = 0; i <= MultBinMax; i++) multBins.push_back(static_cast(i)); + for (int i = 0; i <= MultBinMax; i++) + multBins.push_back(static_cast(i)); AxisSpec multAxis = {multBins, "N_{ch}^{gen} (|#eta|<0.8)"}; std::vector recoMultBins; - for (int i = 0; i <= RecMultBinMax; i++) recoMultBins.push_back(static_cast(i)); + for (int i = 0; i <= RecMultBinMax; i++) + recoMultBins.push_back(static_cast(i)); AxisSpec recoMultAxis = {recoMultBins, "N_{ch}^{reco}"}; - ue.add("Centrality/hCentRaw", "Raw FT0M Centrality;Centrality (%);Counts", HistType::kTH1D, {centFineAxis}); ue.add("Centrality/hCentAfterVtx", "Centrality after vertex cut;Centrality (%);Counts", HistType::kTH1D, {centFineAxis}); ue.add("Centrality/hCentAfterAll", "Centrality after all cuts;Centrality (%);Counts", HistType::kTH1D, {centFineAxis}); @@ -368,13 +401,11 @@ void MultiplicityPt::init(InitContext const&) ue.add("MC/EventLoss/GenMultVsCent", "Generated charged particles vs FT0M centrality;FT0M Centrality (%);N_{ch}^{gen}", HistType::kTH2D, {centAxis, multAxis}); ue.add("MC/EventLoss/GenMultVsCent_Selected", "Generated vs FT0M centrality (selected events);FT0M Centrality (%);N_{ch}^{gen}", HistType::kTH2D, {centAxis, multAxis}); - ue.add("Inclusive/hPtGen", "Generated primaries (physics selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Inclusive/hPtReco", "All reconstructed tracks (track selection only);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Inclusive/hPtPrimReco", "Reconstructed primaries (MC matched);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Inclusive/hPtSecReco", "Reconstructed secondaries (MC matched);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtPrimGenAll", "All generated #pi^{#pm} (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Pion/hPtGenINEL", "Generated #pi^{#pm} in INEL>0 events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Pion/hPtGenRecoEvent", "Generated #pi^{#pm} in reconstructed events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); @@ -382,16 +413,14 @@ void MultiplicityPt::init(InitContext const&) ue.add("Pion/hPtReco", "Reconstructed #pi^{#pm} (MC matched, any status);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Pion/hPtPrimReco", "Reconstructed primary #pi^{#pm} (MC matched);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Pion/hPtMeasured", "Measured #pi^{#pm} (PID selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtGenRecoEvent_Mult", "Generated #pi^{#pm} in reconstructed events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); ue.add("Pion/hPtGenINEL_Mult", "Generated #pi^{#pm} in INEL>0 events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); - + if (enablePIDHistograms) { ue.add("Pion/hNsigmaTPC", "TPC n#sigma #pi^{#pm};#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}", HistType::kTH2D, {ptAxis, {200, -10, 10}}); } - ue.add("Kaon/hPtPrimGenAll", "All generated K^{#pm} (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Kaon/hPtGenINEL", "Generated K^{#pm} in INEL>0 events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Kaon/hPtGenRecoEvent", "Generated K^{#pm} in reconstructed events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); @@ -402,32 +431,29 @@ void MultiplicityPt::init(InitContext const&) ue.add("Kaon/hPtGenRecoEvent_Mult", "Generated K^{#pm} in reconstructed events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); ue.add("Kaon/hPtGenINEL_Mult", "Generated K^{#pm} in INEL>0 events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); - + if (enablePIDHistograms) { ue.add("Kaon/hNsigmaTPC", "TPC n#sigma K^{#pm};#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}", HistType::kTH2D, {ptAxis, {200, -10, 10}}); } - ue.add("Proton/hPtPrimGenAll", "All generated p+#bar{p} (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Proton/hPtPrimGenAll", "All generated p+#bar{p} (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Proton/hPtGenINEL", "Generated p+#bar{p} in INEL>0 events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Proton/hPtGenRecoEvent", "Generated p+#bar{p} in reconstructed events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Proton/hPtGen", "Generated p+#bar{p} (physics selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Proton/hPtReco", "Reconstructed p+#bar{p} (MC matched, any status);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Proton/hPtPrimReco", "Reconstructed primary p+#bar{p} (MC matched);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); ue.add("Proton/hPtMeasured", "Measured p+#bar{p} (PID selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - - + ue.add("Proton/hPtGenRecoEvent_Mult", "Generated p+#bar{p} in reconstructed events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); ue.add("Proton/hPtGenINEL_Mult", "Generated p+#bar{p} in INEL>0 events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); - + if (enablePIDHistograms) { ue.add("Proton/hNsigmaTPC", "TPC n#sigma p+#bar{p};#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}", HistType::kTH2D, {ptAxis, {200, -10, 10}}); } - ue.add("hEventsReco_Cent", "Reconstructed events vs centrality;FT0M Centrality (%);Counts", HistType::kTH1D, {centAxis}); ue.add("hEventsINEL_Cent", "INEL>0 events vs centrality;FT0M Centrality (%);Counts", HistType::kTH1D, {centAxis}); - ue.add("hNclFoundTPC", "Number of TPC found clusters", HistType::kTH1D, {{200, 0, 200}}); ue.add("hNclPIDTPC", "Number of TPC PID clusters", HistType::kTH1D, {{200, 0, 200}}); ue.add("hEta", "Track eta;#eta;Counts", HistType::kTH1D, {{20, -0.8, 0.8}}); @@ -446,12 +472,11 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, BCsRun3 const& /*bcs*/) { - std::map mcCollisionToNch; std::set physicsSelectedMCCollisions; std::map mcCollisionToINELClass; std::set inel0MCCollisions; - std::map mcCollToCentFromReco; // MC collision ID -> centrality from reco + std::map mcCollToCentFromReco; // MC collision ID -> centrality from reco ue.fill(HIST("MC/GenRecoCollisions"), 1.f, mcCollisions.size()); @@ -463,14 +488,18 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, bool hasParticleInAcceptance = false; for (const auto& particle : particlesInCollision) { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (!pdgParticle || pdgParticle->Charge() == 0.) continue; - if (!particle.isPhysicalPrimary()) continue; - if (std::abs(particle.eta()) >= cfgCutEtaMax.value) continue; - if (particle.pt() < cfgTrkLowPtCut.value) continue; + if (!pdgParticle || pdgParticle->Charge() == 0.) + continue; + if (!particle.isPhysicalPrimary()) + continue; + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) + continue; + if (particle.pt() < cfgTrkLowPtCut.value) + continue; hasParticleInAcceptance = true; break; } - + if (hasParticleInAcceptance) { inel0MCCollisions.insert(mcCollId); } @@ -484,24 +513,30 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, mcCollisionToINELClass[mcCollId] = inelClass; bool physicsSelected = (std::abs(mcCollision.posZ()) <= cfgCutVertex.value); - if (cfgINELCut.value == INELgt0 && !inel0) physicsSelected = false; - if (cfgINELCut.value == INELgt1 && !inel1) physicsSelected = false; + if (cfgINELCut.value == INELgt0 && !inel0) + physicsSelected = false; + if (cfgINELCut.value == INELgt1 && !inel1) + physicsSelected = false; if (physicsSelected) { physicsSelectedMCCollisions.insert(mcCollId); - if (inel0) ue.fill(HIST("MC/GenRecoCollisions"), 2.f); - if (inel1) ue.fill(HIST("MC/GenRecoCollisions"), 3.f); + if (inel0) + ue.fill(HIST("MC/GenRecoCollisions"), 2.f); + if (inel1) + ue.fill(HIST("MC/GenRecoCollisions"), 3.f); } // Fill generated particle spectra for (const auto& particle : particlesInCollision) { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (!pdgParticle || pdgParticle->Charge() == 0.) continue; - if (!particle.isPhysicalPrimary()) continue; - + if (!pdgParticle || pdgParticle->Charge() == 0.) + continue; + if (!particle.isPhysicalPrimary()) + continue; + int pdgCode = std::abs(particle.pdgCode()); float pt = particle.pt(); - + // Fill hPtPrimGenAll for ALL generated particles (NO CUTS) if (pdgCode == PDG_t::kPiPlus) { ue.fill(HIST("Pion/hPtPrimGenAll"), pt); @@ -510,12 +545,14 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, } else if (pdgCode == PDG_t::kProton) { ue.fill(HIST("Proton/hPtPrimGenAll"), pt); } - + // Fill hPtGenINEL for particles in INEL>0 events (with acceptance cuts) if (hasParticleInAcceptance) { - if (std::abs(particle.eta()) >= cfgCutEtaMax.value) continue; - if (particle.pt() < cfgTrkLowPtCut.value) continue; - + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) + continue; + if (particle.pt() < cfgTrkLowPtCut.value) + continue; + if (pdgCode == PDG_t::kPiPlus) { ue.fill(HIST("Pion/hPtGenINEL"), pt); ue.fill(HIST("Pion/hPtGenINEL_Mult"), pt, nGenCharged); @@ -527,15 +564,17 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, ue.fill(HIST("Proton/hPtGenINEL_Mult"), pt, nGenCharged); } } - + // Apply acceptance cuts for physics-selected spectra - if (std::abs(particle.eta()) >= cfgCutEtaMax.value) continue; - if (particle.pt() < cfgTrkLowPtCut.value) continue; - + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) + continue; + if (particle.pt() < cfgTrkLowPtCut.value) + continue; + // Fill generated spectra (physics selected only) if (physicsSelected) { ue.fill(HIST("Inclusive/hPtGen"), pt); - + if (pdgCode == PDG_t::kPiPlus) { ue.fill(HIST("Pion/hPtGen"), pt); } else if (pdgCode == PDG_t::kKPlus) { @@ -543,14 +582,10 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, } else if (pdgCode == PDG_t::kProton) { ue.fill(HIST("Proton/hPtGen"), pt); } - - } } } - - std::map recoToMcMap; std::map recoToCentMap; @@ -570,62 +605,67 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, ue.fill(HIST("Centrality/hCentRaw"), centValue); } - - std::set reconstructedMCCollisions; std::set selectedMCCollisions; for (const auto& collision : collisions) { ue.fill(HIST("CutFlow/hCutStats"), 1); - + int64_t collId = collision.globalIndex(); // MC matching auto mcIt = recoToMcMap.find(collId); - if (mcIt == recoToMcMap.end()) continue; + if (mcIt == recoToMcMap.end()) + continue; ue.fill(HIST("CutFlow/hCutStats"), 2); - + int64_t mcCollId = mcIt->second; auto nchIt = mcCollisionToNch.find(mcCollId); - if (nchIt == mcCollisionToNch.end()) continue; + if (nchIt == mcCollisionToNch.end()) + continue; int nGenCharged = nchIt->second; - + auto inelIt = mcCollisionToINELClass.find(mcCollId); int inelClass = (inelIt != mcCollisionToINELClass.end()) ? inelIt->second : 0; // Centrality auto centIt = recoToCentMap.find(collId); - if (centIt == recoToCentMap.end()) continue; + if (centIt == recoToCentMap.end()) + continue; ue.fill(HIST("CutFlow/hCutStats"), 3); - + float cent = centIt->second; - if (cent < 0 || cent > CentBinMax) continue; - + if (cent < 0 || cent > CentBinMax) + continue; + // Store centrality for this MC collision (used later for MC truth plots) mcCollToCentFromReco[mcCollId] = cent; // Vertex cut bool passVertex = std::abs(collision.posZ()) <= cfgCutVertex.value; - if (!passVertex) continue; + if (!passVertex) + continue; ue.fill(HIST("CutFlow/hCutStats"), 4); ue.fill(HIST("Centrality/hCentAfterVtx"), cent); // INEL cut bool passINEL = true; - if (cfgINELCut.value == INELgt0 && inelClass < INELgt0) passINEL = false; - if (cfgINELCut.value == INELgt1 && inelClass < INELgt1) passINEL = false; - if (!passINEL) continue; + if (cfgINELCut.value == INELgt0 && inelClass < INELgt0) + passINEL = false; + if (cfgINELCut.value == INELgt1 && inelClass < INELgt1) + passINEL = false; + if (!passINEL) + continue; // Event passed all cuts ue.fill(HIST("CutFlow/hCutStats"), 5); ue.fill(HIST("Centrality/hCentAfterAll"), cent); - ue.fill(HIST("MC/EventLoss/GenMultVsCent_Selected"), cent, nGenCharged); ue.fill(HIST("Centrality/hCentVsMult"), cent, nGenCharged); ue.fill(HIST("Centrality/hMultVsCent"), nGenCharged, cent); ue.fill(HIST("hvtxZ"), collision.posZ()); - + selectedMCCollisions.insert(mcCollId); reconstructedMCCollisions.insert(mcCollId); @@ -636,27 +676,29 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, magField = getMagneticField(bc.timestamp()); } - int nTracksInEvent = 0; - + for (const auto& track : tracks) { - if (!track.has_collision()) continue; - if (track.collisionId() != collId) continue; - - if (!passesTrackSelection(track, magField)) continue; - + if (!track.has_collision()) + continue; + if (track.collisionId() != collId) + continue; + + if (!passesTrackSelection(track, magField)) + continue; + nTracksInEvent++; - + ue.fill(HIST("hNclFoundTPC"), track.tpcNClsFound()); ue.fill(HIST("hNclPIDTPC"), track.tpcNClsPID()); ue.fill(HIST("hEta"), track.eta()); ue.fill(HIST("hPhi"), track.phi()); ue.fill(HIST("Inclusive/hPtReco"), track.pt()); - + if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); int pdgCode = std::abs(particle.pdgCode()); - + if (pdgCode == PDG_t::kPiPlus) { ue.fill(HIST("Pion/hPtReco"), track.pt()); if (particle.isPhysicalPrimary()) { @@ -683,21 +725,21 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, } } } - + if (passesPIDSelection(track, PartPion)) { ue.fill(HIST("Pion/hPtMeasured"), track.pt()); if (enablePIDHistograms) { ue.fill(HIST("Pion/hNsigmaTPC"), track.pt(), track.tpcNSigmaPi()); } } - + if (passesPIDSelection(track, PartKaon)) { ue.fill(HIST("Kaon/hPtMeasured"), track.pt()); if (enablePIDHistograms) { ue.fill(HIST("Kaon/hNsigmaTPC"), track.pt(), track.tpcNSigmaKa()); } } - + if (passesPIDSelection(track, PartProton)) { ue.fill(HIST("Proton/hPtMeasured"), track.pt()); if (enablePIDHistograms) { @@ -705,36 +747,40 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, } } } - + ue.fill(HIST("Centrality/hRecoMultVsCent"), cent, nTracksInEvent); } - - for (const auto& mcCollision : mcCollisions) { int64_t mcCollId = mcCollision.globalIndex(); - - if (reconstructedMCCollisions.find(mcCollId) == reconstructedMCCollisions.end()) continue; - + + if (reconstructedMCCollisions.find(mcCollId) == reconstructedMCCollisions.end()) + continue; + auto centIt = mcCollToCentFromReco.find(mcCollId); - if (centIt == mcCollToCentFromReco.end()) continue; - + if (centIt == mcCollToCentFromReco.end()) + continue; + auto nchIt = mcCollisionToNch.find(mcCollId); int nGenCharged = (nchIt != mcCollisionToNch.end()) ? nchIt->second : 0; - + auto particlesInCollision = particles.sliceBy(perMCCol, mcCollId); - + for (const auto& particle : particlesInCollision) { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (!pdgParticle || pdgParticle->Charge() == 0.) continue; - if (!particle.isPhysicalPrimary()) continue; - - if (std::abs(particle.eta()) >= cfgCutEtaMax.value) continue; - if (particle.pt() < cfgTrkLowPtCut.value) continue; - + if (!pdgParticle || pdgParticle->Charge() == 0.) + continue; + if (!particle.isPhysicalPrimary()) + continue; + + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) + continue; + if (particle.pt() < cfgTrkLowPtCut.value) + continue; + int pdgCode = std::abs(particle.pdgCode()); float pt = particle.pt(); - + if (pdgCode == PDG_t::kPiPlus) { ue.fill(HIST("Pion/hPtGenRecoEvent"), pt); ue.fill(HIST("Pion/hPtGenRecoEvent_Mult"), pt, nGenCharged); @@ -748,7 +794,6 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, } } - for (const auto& mcCollId : inel0MCCollisions) { auto centIt = mcCollToCentFromReco.find(mcCollId); if (centIt != mcCollToCentFromReco.end()) { @@ -759,8 +804,6 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, } } - - for (const auto& mcCollId : reconstructedMCCollisions) { auto centIt = mcCollToCentFromReco.find(mcCollId); if (centIt != mcCollToCentFromReco.end()) { @@ -768,12 +811,10 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, } } - - ue.fill(HIST("hEventLossBreakdown"), 1.f, physicsSelectedMCCollisions.size()); ue.fill(HIST("hEventLossBreakdown"), 2.f, reconstructedMCCollisions.size()); ue.fill(HIST("hEventLossBreakdown"), 3.f, selectedMCCollisions.size()); - + ue.fill(HIST("MC/GenRecoCollisions"), 4.f, reconstructedMCCollisions.size()); ue.fill(HIST("MC/GenRecoCollisions"), 5.f, selectedMCCollisions.size()); @@ -782,7 +823,7 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, LOG(info) << "Physics selected MC collisions: " << physicsSelectedMCCollisions.size(); LOG(info) << "Reconstructed collisions: " << reconstructedMCCollisions.size(); LOG(info) << "Selected collisions: " << selectedMCCollisions.size(); - + if (physicsSelectedMCCollisions.size() > 0) { float efficiency = 100.f * selectedMCCollisions.size() / physicsSelectedMCCollisions.size(); LOG(info) << "Final efficiency: " << efficiency << "%"; @@ -793,5 +834,4 @@ void MultiplicityPt::processData(CollisionTableData::iterator const& /*collision TrackTableData const& /*tracks*/, BCsRun3 const& /*bcs*/) { - }