Skip to content

PWGMM: Update dndeta task #438

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
aalkin merged 3 commits into from
Jan 11, 2022
Merged

PWGMM: Update dndeta task #438

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5a9461d
Vertex cut is applied in post-processing
aalkin Dec 14, 2021
6e4b905
using partitions for mult. estimation
aalkin Dec 17, 2021
9f49d51
cleanup
aalkin Jan 11, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
149 changes: 98 additions & 51 deletions PWGMM/Tasks/dndeta.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -61,11 +63,9 @@ template <uint8_t TRACKTYPE>
struct PseudorapidityDensity {
Service<TDatabasePDG> pdg;

Configurable<float> etaMax{"etaMax", 2.0, "max eta value"};
Configurable<float> etaMin{"etaMin", -2.0, "min eta value"};
Configurable<float> vtxZMax{"vtxZMax", 15, "max z vertex"};
Configurable<float> vtxZMin{"vtxZMin", -15, "min z vertex"};
Configurable<float> estimatorEta{"estimatorEta", 1.0, "eta range for INEL>0 sample definition"};

Configurable<bool> useEvSel{"useEvSel", true, "use event selection"};
Configurable<bool> useDCA{"useDCA", false, "use DCA cuts"};
Configurable<float> maxDCAXY{"maxDCAXY", 2.4, "max allowed transverse DCA"};
Configurable<float> maxDCAZ{"maxDCAZ", 3.2, "max allowed longitudal DCA"};
Expand All @@ -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<bool> 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}}}} //
} //
Expand All @@ -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<TH1>(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) {
Expand Down Expand Up @@ -140,65 +147,73 @@ struct PseudorapidityDensity {
return BCSelector<B, TRACKTYPE>(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<aod::BCsWithTimestamps, aod::BcSels>;
void processTagging(FullBCs const& bcs, soa::Filtered<soa::Join<aod::Collisions, aod::EvSels>> const& collisions, aod::FT0s const& ft0s)
void processTagging(FullBCs const& bcs, soa::Join<aod::Collisions, aod::EvSels> const& collisions)
{
std::vector<typename std::decay_t<decltype(collisions)>::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.);
}
}
}
}
}

PROCESS_SWITCH(PseudorapidityDensity<TRACKTYPE>, 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<soa::Join<aod::Collisions, aod::EvSels>>::iterator const& collision, soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtended>> const& tracks)
using Trks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtended>>;
Partition<Trks> sample = nabs(aod::track::eta) < estimatorEta;

void process(soa::Join<aod::Collisions, aod::EvSels>::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<soa::Join<aod::Collisions, aod::EvSels, aod::CentV0Ms>>::iterator const& collision, soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtended>> const& tracks)
void processBinned(soa::Join<aod::Collisions, aod::EvSels, aod::CentV0Ms>::iterator const& collision, soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtended>> 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);
Expand All @@ -211,33 +226,65 @@ struct PseudorapidityDensity {
PROCESS_SWITCH(PseudorapidityDensity<TRACKTYPE>, 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<soa::Filtered<Particles>> mcSample = nabs(aod::mcparticle::eta) < estimatorEta;

void processGen(soa::Filtered<aod::McCollisions>::iterator const& mcCollision, o2::soa::SmallGroups<soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>> const& collisions, Particles const& particles, soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtended>> const& tracks)
void processGen(aod::McCollisions::iterator const& mcCollision, o2::soa::SmallGroups<soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>> const& collisions, soa::Filtered<Particles> const& particles, soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtended>> 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());
}
}
Expand Down