diff --git a/PWGMM/Tasks/dndeta.h b/PWGMM/Tasks/dndeta.h index 08c1d2ee21e..be51f26bc86 100644 --- a/PWGMM/Tasks/dndeta.h +++ b/PWGMM/Tasks/dndeta.h @@ -24,6 +24,8 @@ #include "Common/DataModel/TrackSelectionTables.h" #include "Common/Core/MC.h" +#include "CommonConstants/MathConstants.h" + #include "TDatabasePDG.h" namespace o2::pwgmm::multiplicity @@ -61,11 +63,9 @@ template struct PseudorapidityDensity { Service pdg; - Configurable etaMax{"etaMax", 2.0, "max eta value"}; - Configurable etaMin{"etaMin", -2.0, "min eta value"}; - Configurable vtxZMax{"vtxZMax", 15, "max z vertex"}; - Configurable vtxZMin{"vtxZMin", -15, "min z vertex"}; + Configurable estimatorEta{"estimatorEta", 1.0, "eta range for INEL>0 sample definition"}; + Configurable useEvSel{"useEvSel", true, "use event selection"}; Configurable useDCA{"useDCA", false, "use DCA cuts"}; Configurable maxDCAXY{"maxDCAXY", 2.4, "max allowed transverse DCA"}; Configurable maxDCAZ{"maxDCAZ", 3.2, "max allowed longitudal DCA"}; @@ -74,13 +74,12 @@ struct PseudorapidityDensity { {VARIABLE_WIDTH, 0., 0.01, 0.1, 0.5, 1, 5, 10, 15, 20, 30, 40, 50, 70, 100}, "Centrality/multiplicity percentile binning"}; - Configurable useEvSel{"useEvSel", true, "use event selection"}; - HistogramRegistry registry{ "registry", { {"EventsNtrkZvtx", "; N_{trk}; Z_{vtx}; events", {HistType::kTH2F, {{301, -0.5, 300.5}, {201, -20.1, 20.1}}}}, // {"TracksEtaZvtx", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}}, // + {"TracksEtaZvtx_gt0", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}}, // {"TracksPhiEta", "; #varphi; #eta; tracks", {HistType::kTH2F, {{600, 0, 2 * M_PI}, {21, -2.1, 2.1}}}}, // {"EventSelection", ";status;events", {HistType::kTH1F, {{7, 0.5, 7.5}}}} // } // @@ -92,23 +91,31 @@ struct PseudorapidityDensity { auto* x = hstat->GetXaxis(); x->SetBinLabel(1, "All"); x->SetBinLabel(2, "Selected"); - x->SetBinLabel(3, "Rejected"); - x->SetBinLabel(4, "Good BCs"); - x->SetBinLabel(5, "BCs with collisions"); - x->SetBinLabel(6, "BCs with selected collisions"); - x->SetBinLabel(7, "BCs with pile-up"); + x->SetBinLabel(3, "Selected INEL>0"); + x->SetBinLabel(4, "Rejected"); + x->SetBinLabel(5, "Good BCs"); + x->SetBinLabel(6, "BCs with collisions"); + x->SetBinLabel(7, "BCs with pile-up/splitting"); if (doprocessGen) { registry.add({"EventsNtrkZvtxGen", "; N_{trk}; Z_{vtx}; events", {HistType::kTH2F, {{301, -0.5, 300.5}, {201, -20.1, 20.1}}}}); + registry.add({"EventsNtrkZvtxGen_t", "; N_{part}; Z_{vtx}; events", {HistType::kTH2F, {{301, -0.5, 300.5}, {201, -20.1, 20.1}}}}); registry.add({"TracksEtaZvtxGen", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}}); + registry.add({"TracksEtaZvtxGen_t", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}}); + registry.add({"TracksEtaZvtxGen_gt0", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}}); + registry.add({"TracksEtaZvtxGen_gt0t", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}}); + registry.add({"TracksPhiEtaGen", "; #varphi; #eta; tracks", {HistType::kTH2F, {{600, 0, 2 * M_PI}, {21, -2.1, 2.1}}}}); - registry.add({"EventEfficiency", "; status; events", {HistType::kTH1F, {{3, 0.5, 3.5}}}}); + registry.add({"EventEfficiency", "; status; events", {HistType::kTH1F, {{5, 0.5, 5.5}}}}); + registry.add({"NotFoundEventZvtx", " ; Z_{vtx}", {HistType::kTH1F, {{201, -20.1, 20.1}}}}); auto heff = registry.get(HIST("EventEfficiency")); x = heff->GetXaxis(); x->SetBinLabel(1, "Generated"); - x->SetBinLabel(2, "Reconstructed"); - x->SetBinLabel(3, "Selected"); + x->SetBinLabel(2, "Generated INEL>0"); + x->SetBinLabel(3, "Reconstructed"); + x->SetBinLabel(4, "Selected"); + x->SetBinLabel(5, "Selected INEL>0"); } if (doprocessBinned) { @@ -140,33 +147,28 @@ struct PseudorapidityDensity { return BCSelector(bc); } - expressions::Filter etaFilter = (aod::track::eta < etaMax) && (aod::track::eta > etaMin); - expressions::Filter trackTypeFilter = (aod::track::trackType == TRACKTYPE); - expressions::Filter DCAFilter = ifnode(useDCA.node(), nabs(aod::track::dcaXY) <= maxDCAXY && nabs(aod::track::dcaZ) <= maxDCAZ, framework::expressions::LiteralNode{true}); - expressions::Filter posZFilter = (aod::collision::posZ < vtxZMax) && (aod::collision::posZ > vtxZMin); - expressions::Filter posZFilterMC = (aod::mccollision::posZ < vtxZMax) && (aod::mccollision::posZ > vtxZMin); - using FullBCs = soa::Join; - void processTagging(FullBCs const& bcs, soa::Filtered> const& collisions, aod::FT0s const& ft0s) + void processTagging(FullBCs const& bcs, soa::Join const& collisions) { std::vector::iterator> cols; for (auto& bc : bcs) { if (selectBC(bc)) { - registry.fill(HIST("EventSelection"), 4); + registry.fill(HIST("EventSelection"), 5.); cols.clear(); for (auto& collision : collisions) { - if ((collision.has_foundFT0()) && (collision.foundFT0().bcId() == bc.globalIndex())) { + if (collision.has_foundBC()) { + if (collision.foundBCId() == bc.globalIndex()) { + cols.emplace_back(collision); + } + } else if (collision.bcId() == bc.globalIndex()) { cols.emplace_back(collision); } } - LOGP(info, "BC {} has {} collisions", bc.globalBC(), cols.size()); + LOGP(debug, "BC {} has {} collisions", bc.globalBC(), cols.size()); if (!cols.empty()) { - registry.fill(HIST("EventSelection"), 5); - if (std::any_of(cols.begin(), cols.end(), [&](auto const& x) { return select(x); })) { - registry.fill(HIST("EventSelection"), 6); - } + registry.fill(HIST("EventSelection"), 6.); if (cols.size() > 1) { - registry.fill(HIST("EventSelection"), 7); + registry.fill(HIST("EventSelection"), 7.); } } } @@ -174,31 +176,44 @@ struct PseudorapidityDensity { } PROCESS_SWITCH(PseudorapidityDensity, processTagging, "Collect event sample stats", false); + expressions::Filter trackTypeFilter = (aod::track::trackType == TRACKTYPE); + expressions::Filter DCAFilter = ifnode(useDCA.node(), nabs(aod::track::dcaXY) <= maxDCAXY && nabs(aod::track::dcaZ) <= maxDCAZ, framework::expressions::LiteralNode{true}); - void process(soa::Filtered>::iterator const& collision, soa::Filtered> const& tracks) + using Trks = soa::Filtered>; + Partition sample = nabs(aod::track::eta) < estimatorEta; + + void process(soa::Join::iterator const& collision, Trks const& tracks) { registry.fill(HIST("EventSelection"), 1.); if (select(collision)) { registry.fill(HIST("EventSelection"), 2.); auto z = collision.posZ(); - registry.fill(HIST("EventsNtrkZvtx"), tracks.size(), z); + auto perCollisionSample = sample->sliceByCached(aod::track::collisionId, collision.globalIndex()); + if (perCollisionSample.size() > 0) { + registry.fill(HIST("EventSelection"), 3.); + } + registry.fill(HIST("EventsNtrkZvtx"), perCollisionSample.size(), z); for (auto& track : tracks) { registry.fill(HIST("TracksEtaZvtx"), track.eta(), z); registry.fill(HIST("TracksPhiEta"), track.phi(), track.eta()); + if (perCollisionSample.size() > 0) { + registry.fill(HIST("TracksEtaZvtx_gt0"), track.eta(), z); + } } } else { - registry.fill(HIST("EventSelection"), 3.); + registry.fill(HIST("EventSelection"), 4.); } } - void processBinned(soa::Filtered>::iterator const& collision, soa::Filtered> const& tracks) + void processBinned(soa::Join::iterator const& collision, soa::Filtered> const& tracks) { auto p = collision.centV0M(); registry.fill(HIST("EventSelectionBin"), 1., p); if (select(collision)) { registry.fill(HIST("EventSelectionBin"), 2., p); auto z = collision.posZ(); - registry.fill(HIST("EventsNtrkZvtxBin"), tracks.size(), z, p); + auto localSample = sample->sliceByCached(aod::track::collisionId, collision.globalIndex()); + registry.fill(HIST("EventsNtrkZvtxBin"), localSample.size(), z, p); for (auto& track : tracks) { registry.fill(HIST("TracksEtaZvtxBin"), track.eta(), z, p); registry.fill(HIST("TracksPhiEtaBin"), track.phi(), track.eta(), p); @@ -211,33 +226,65 @@ struct PseudorapidityDensity { PROCESS_SWITCH(PseudorapidityDensity, processBinned, "Process centrality/mult. percentile binned", false); using Particles = aod::McParticles; + expressions::Filter primaries = (aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary; + Partition> mcSample = nabs(aod::mcparticle::eta) < estimatorEta; - void processGen(soa::Filtered::iterator const& mcCollision, o2::soa::SmallGroups> const& collisions, Particles const& particles, soa::Filtered> const& tracks) + void processGen(aod::McCollisions::iterator const& mcCollision, o2::soa::SmallGroups> const& collisions, soa::Filtered const& particles, soa::Filtered> const& /*tracks*/) { + auto perCollisionMCSample = mcSample->sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex()); + auto nCharged = 0; + for (auto& particle : perCollisionMCSample) { + auto charge = 0; + auto p = pdg->GetParticle(particle.pdgCode()); + if (p != nullptr) { + charge = (int)p->Charge(); + } + if (charge != 0) { + nCharged++; + } + } + registry.fill(HIST("EventsNtrkZvtxGen_t"), nCharged, mcCollision.posZ()); registry.fill(HIST("EventEfficiency"), 1.); - for (auto& collision : collisions) { + + if (nCharged > 0) { registry.fill(HIST("EventEfficiency"), 2.); + } + bool atLeastOne = false; + bool atLeastOne_gt0 = false; + LOGP(debug, "MC col {} has {} reco cols", mcCollision.globalIndex(), collisions.size()); + for (auto& collision : collisions) { + registry.fill(HIST("EventEfficiency"), 3.); if (select(collision)) { - registry.fill(HIST("EventEfficiency"), 3.); - auto stracks = tracks.sliceBy(aod::track::collisionId, collision.globalIndex()); - registry.fill(HIST("EventsNtrkZvtxGen"), stracks.size(), mcCollision.posZ()); + atLeastOne = true; + auto perCollisionSample = sample->sliceByCached(aod::track::collisionId, collision.globalIndex()); + registry.fill(HIST("EventEfficiency"), 4.); + if (perCollisionSample.size() > 0) { + atLeastOne_gt0 = true; + registry.fill(HIST("EventEfficiency"), 5.); + } + registry.fill(HIST("EventsNtrkZvtxGen"), perCollisionSample.size(), collision.posZ()); } } + if (collisions.size() == 0) { + registry.fill(HIST("NotFoundEventZvtx"), mcCollision.posZ()); + } for (auto& particle : particles) { auto p = pdg->GetParticle(particle.pdgCode()); - int charge = 0; - if (p == nullptr) { - // unknown particles will be skipped - if (particle.pdgCode() > 1000000000) { - LOGP(debug, "[{}] Nucleus with PDG code {}", particle.globalIndex(), particle.pdgCode()); - } else { - LOGP(debug, "[{}] Unknown particle with PDG code {}", particle.globalIndex(), particle.pdgCode()); - } - } else { - charge = p->Charge(); + auto charge = 0; + if (p != nullptr) { + charge = (int)p->Charge(); } - if (charge != 0 && particle.isPhysicalPrimary() && (particle.eta() < etaMax) && (particle.eta() > etaMin)) { - registry.fill(HIST("TracksEtaZvtxGen"), particle.eta(), mcCollision.posZ()); + if (charge != 0) { + registry.fill(HIST("TracksEtaZvtxGen_t"), particle.eta(), mcCollision.posZ()); + if (perCollisionMCSample.size() > 0) { + registry.fill(HIST("TracksEtaZvtxGen_gt0t"), particle.eta(), mcCollision.posZ()); + } + if (atLeastOne) { + registry.fill(HIST("TracksEtaZvtxGen"), particle.eta(), mcCollision.posZ()); + if (atLeastOne_gt0) { + registry.fill(HIST("TracksEtaZvtxGen_gt0"), particle.eta(), mcCollision.posZ()); + } + } registry.fill(HIST("TracksPhiEtaGen"), particle.phi(), particle.eta()); } }