Skip to content

Commit e8f9a73

Browse files
aalkinaalkin
authored
PWGMM: Update dndeta task (#438)
* Vertex cut is applied in post-processing * using partitions for mult. estimation * cleanup
1 parent ec50aa6 commit e8f9a73

1 file changed

Lines changed: 98 additions & 51 deletions

File tree

‎PWGMM/Tasks/dndeta.h‎

Lines changed: 98 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,8 @@
2424
#include "Common/DataModel/TrackSelectionTables.h"
2525
#include "Common/Core/MC.h"
2626

27+
#include "CommonConstants/MathConstants.h"
28+
2729
#include "TDatabasePDG.h"
2830

2931
namespace o2::pwgmm::multiplicity
@@ -61,11 +63,9 @@ template <uint8_t TRACKTYPE>
6163
struct PseudorapidityDensity {
6264
Service<TDatabasePDG> pdg;
6365

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

68+
Configurable<bool> useEvSel{"useEvSel", true, "use event selection"};
6969
Configurable<bool> useDCA{"useDCA", false, "use DCA cuts"};
7070
Configurable<float> maxDCAXY{"maxDCAXY", 2.4, "max allowed transverse DCA"};
7171
Configurable<float> maxDCAZ{"maxDCAZ", 3.2, "max allowed longitudal DCA"};
@@ -74,13 +74,12 @@ struct PseudorapidityDensity {
7474
{VARIABLE_WIDTH, 0., 0.01, 0.1, 0.5, 1, 5, 10, 15, 20, 30, 40, 50, 70, 100},
7575
"Centrality/multiplicity percentile binning"};
7676

77-
Configurable<bool> useEvSel{"useEvSel", true, "use event selection"};
78-
7977
HistogramRegistry registry{
8078
"registry",
8179
{
8280
{"EventsNtrkZvtx", "; N_{trk}; Z_{vtx}; events", {HistType::kTH2F, {{301, -0.5, 300.5}, {201, -20.1, 20.1}}}}, //
8381
{"TracksEtaZvtx", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}}, //
82+
{"TracksEtaZvtx_gt0", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}}, //
8483
{"TracksPhiEta", "; #varphi; #eta; tracks", {HistType::kTH2F, {{600, 0, 2 * M_PI}, {21, -2.1, 2.1}}}}, //
8584
{"EventSelection", ";status;events", {HistType::kTH1F, {{7, 0.5, 7.5}}}} //
8685
} //
@@ -92,23 +91,31 @@ struct PseudorapidityDensity {
9291
auto* x = hstat->GetXaxis();
9392
x->SetBinLabel(1, "All");
9493
x->SetBinLabel(2, "Selected");
95-
x->SetBinLabel(3, "Rejected");
96-
x->SetBinLabel(4, "Good BCs");
97-
x->SetBinLabel(5, "BCs with collisions");
98-
x->SetBinLabel(6, "BCs with selected collisions");
99-
x->SetBinLabel(7, "BCs with pile-up");
94+
x->SetBinLabel(3, "Selected INEL>0");
95+
x->SetBinLabel(4, "Rejected");
96+
x->SetBinLabel(5, "Good BCs");
97+
x->SetBinLabel(6, "BCs with collisions");
98+
x->SetBinLabel(7, "BCs with pile-up/splitting");
10099

101100
if (doprocessGen) {
102101
registry.add({"EventsNtrkZvtxGen", "; N_{trk}; Z_{vtx}; events", {HistType::kTH2F, {{301, -0.5, 300.5}, {201, -20.1, 20.1}}}});
102+
registry.add({"EventsNtrkZvtxGen_t", "; N_{part}; Z_{vtx}; events", {HistType::kTH2F, {{301, -0.5, 300.5}, {201, -20.1, 20.1}}}});
103103
registry.add({"TracksEtaZvtxGen", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}});
104+
registry.add({"TracksEtaZvtxGen_t", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}});
105+
registry.add({"TracksEtaZvtxGen_gt0", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}});
106+
registry.add({"TracksEtaZvtxGen_gt0t", "; #eta; Z_{vtx}; tracks", {HistType::kTH2F, {{21, -2.1, 2.1}, {201, -20.1, 20.1}}}});
107+
104108
registry.add({"TracksPhiEtaGen", "; #varphi; #eta; tracks", {HistType::kTH2F, {{600, 0, 2 * M_PI}, {21, -2.1, 2.1}}}});
105-
registry.add({"EventEfficiency", "; status; events", {HistType::kTH1F, {{3, 0.5, 3.5}}}});
109+
registry.add({"EventEfficiency", "; status; events", {HistType::kTH1F, {{5, 0.5, 5.5}}}});
110+
registry.add({"NotFoundEventZvtx", " ; Z_{vtx}", {HistType::kTH1F, {{201, -20.1, 20.1}}}});
106111

107112
auto heff = registry.get<TH1>(HIST("EventEfficiency"));
108113
x = heff->GetXaxis();
109114
x->SetBinLabel(1, "Generated");
110-
x->SetBinLabel(2, "Reconstructed");
111-
x->SetBinLabel(3, "Selected");
115+
x->SetBinLabel(2, "Generated INEL>0");
116+
x->SetBinLabel(3, "Reconstructed");
117+
x->SetBinLabel(4, "Selected");
118+
x->SetBinLabel(5, "Selected INEL>0");
112119
}
113120

114121
if (doprocessBinned) {
@@ -140,65 +147,73 @@ struct PseudorapidityDensity {
140147
return BCSelector<B, TRACKTYPE>(bc);
141148
}
142149

143-
expressions::Filter etaFilter = (aod::track::eta < etaMax) && (aod::track::eta > etaMin);
144-
expressions::Filter trackTypeFilter = (aod::track::trackType == TRACKTYPE);
145-
expressions::Filter DCAFilter = ifnode(useDCA.node(), nabs(aod::track::dcaXY) <= maxDCAXY && nabs(aod::track::dcaZ) <= maxDCAZ, framework::expressions::LiteralNode{true});
146-
expressions::Filter posZFilter = (aod::collision::posZ < vtxZMax) && (aod::collision::posZ > vtxZMin);
147-
expressions::Filter posZFilterMC = (aod::mccollision::posZ < vtxZMax) && (aod::mccollision::posZ > vtxZMin);
148-
149150
using FullBCs = soa::Join<aod::BCsWithTimestamps, aod::BcSels>;
150-
void processTagging(FullBCs const& bcs, soa::Filtered<soa::Join<aod::Collisions, aod::EvSels>> const& collisions, aod::FT0s const& ft0s)
151+
void processTagging(FullBCs const& bcs, soa::Join<aod::Collisions, aod::EvSels> const& collisions)
151152
{
152153
std::vector<typename std::decay_t<decltype(collisions)>::iterator> cols;
153154
for (auto& bc : bcs) {
154155
if (selectBC(bc)) {
155-
registry.fill(HIST("EventSelection"), 4);
156+
registry.fill(HIST("EventSelection"), 5.);
156157
cols.clear();
157158
for (auto& collision : collisions) {
158-
if ((collision.has_foundFT0()) && (collision.foundFT0().bcId() == bc.globalIndex())) {
159+
if (collision.has_foundBC()) {
160+
if (collision.foundBCId() == bc.globalIndex()) {
161+
cols.emplace_back(collision);
162+
}
163+
} else if (collision.bcId() == bc.globalIndex()) {
159164
cols.emplace_back(collision);
160165
}
161166
}
162-
LOGP(info, "BC {} has {} collisions", bc.globalBC(), cols.size());
167+
LOGP(debug, "BC {} has {} collisions", bc.globalBC(), cols.size());
163168
if (!cols.empty()) {
164-
registry.fill(HIST("EventSelection"), 5);
165-
if (std::any_of(cols.begin(), cols.end(), [&](auto const& x) { return select(x); })) {
166-
registry.fill(HIST("EventSelection"), 6);
167-
}
169+
registry.fill(HIST("EventSelection"), 6.);
168170
if (cols.size() > 1) {
169-
registry.fill(HIST("EventSelection"), 7);
171+
registry.fill(HIST("EventSelection"), 7.);
170172
}
171173
}
172174
}
173175
}
174176
}
175177

176178
PROCESS_SWITCH(PseudorapidityDensity<TRACKTYPE>, processTagging, "Collect event sample stats", false);
179+
expressions::Filter trackTypeFilter = (aod::track::trackType == TRACKTYPE);
180+
expressions::Filter DCAFilter = ifnode(useDCA.node(), nabs(aod::track::dcaXY) <= maxDCAXY && nabs(aod::track::dcaZ) <= maxDCAZ, framework::expressions::LiteralNode{true});
177181

178-
void process(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels>>::iterator const& collision, soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtended>> const& tracks)
182+
using Trks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtended>>;
183+
Partition<Trks> sample = nabs(aod::track::eta) < estimatorEta;
184+
185+
void process(soa::Join<aod::Collisions, aod::EvSels>::iterator const& collision, Trks const& tracks)
179186
{
180187
registry.fill(HIST("EventSelection"), 1.);
181188
if (select(collision)) {
182189
registry.fill(HIST("EventSelection"), 2.);
183190
auto z = collision.posZ();
184-
registry.fill(HIST("EventsNtrkZvtx"), tracks.size(), z);
191+
auto perCollisionSample = sample->sliceByCached(aod::track::collisionId, collision.globalIndex());
192+
if (perCollisionSample.size() > 0) {
193+
registry.fill(HIST("EventSelection"), 3.);
194+
}
195+
registry.fill(HIST("EventsNtrkZvtx"), perCollisionSample.size(), z);
185196
for (auto& track : tracks) {
186197
registry.fill(HIST("TracksEtaZvtx"), track.eta(), z);
187198
registry.fill(HIST("TracksPhiEta"), track.phi(), track.eta());
199+
if (perCollisionSample.size() > 0) {
200+
registry.fill(HIST("TracksEtaZvtx_gt0"), track.eta(), z);
201+
}
188202
}
189203
} else {
190-
registry.fill(HIST("EventSelection"), 3.);
204+
registry.fill(HIST("EventSelection"), 4.);
191205
}
192206
}
193207

194-
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)
208+
void processBinned(soa::Join<aod::Collisions, aod::EvSels, aod::CentV0Ms>::iterator const& collision, soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtended>> const& tracks)
195209
{
196210
auto p = collision.centV0M();
197211
registry.fill(HIST("EventSelectionBin"), 1., p);
198212
if (select(collision)) {
199213
registry.fill(HIST("EventSelectionBin"), 2., p);
200214
auto z = collision.posZ();
201-
registry.fill(HIST("EventsNtrkZvtxBin"), tracks.size(), z, p);
215+
auto localSample = sample->sliceByCached(aod::track::collisionId, collision.globalIndex());
216+
registry.fill(HIST("EventsNtrkZvtxBin"), localSample.size(), z, p);
202217
for (auto& track : tracks) {
203218
registry.fill(HIST("TracksEtaZvtxBin"), track.eta(), z, p);
204219
registry.fill(HIST("TracksPhiEtaBin"), track.phi(), track.eta(), p);
@@ -211,33 +226,65 @@ struct PseudorapidityDensity {
211226
PROCESS_SWITCH(PseudorapidityDensity<TRACKTYPE>, processBinned, "Process centrality/mult. percentile binned", false);
212227

213228
using Particles = aod::McParticles;
229+
expressions::Filter primaries = (aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary;
230+
Partition<soa::Filtered<Particles>> mcSample = nabs(aod::mcparticle::eta) < estimatorEta;
214231

215-
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)
232+
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*/)
216233
{
234+
auto perCollisionMCSample = mcSample->sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex());
235+
auto nCharged = 0;
236+
for (auto& particle : perCollisionMCSample) {
237+
auto charge = 0;
238+
auto p = pdg->GetParticle(particle.pdgCode());
239+
if (p != nullptr) {
240+
charge = (int)p->Charge();
241+
}
242+
if (charge != 0) {
243+
nCharged++;
244+
}
245+
}
246+
registry.fill(HIST("EventsNtrkZvtxGen_t"), nCharged, mcCollision.posZ());
217247
registry.fill(HIST("EventEfficiency"), 1.);
218-
for (auto& collision : collisions) {
248+
249+
if (nCharged > 0) {
219250
registry.fill(HIST("EventEfficiency"), 2.);
251+
}
252+
bool atLeastOne = false;
253+
bool atLeastOne_gt0 = false;
254+
LOGP(debug, "MC col {} has {} reco cols", mcCollision.globalIndex(), collisions.size());
255+
for (auto& collision : collisions) {
256+
registry.fill(HIST("EventEfficiency"), 3.);
220257
if (select(collision)) {
221-
registry.fill(HIST("EventEfficiency"), 3.);
222-
auto stracks = tracks.sliceBy(aod::track::collisionId, collision.globalIndex());
223-
registry.fill(HIST("EventsNtrkZvtxGen"), stracks.size(), mcCollision.posZ());
258+
atLeastOne = true;
259+
auto perCollisionSample = sample->sliceByCached(aod::track::collisionId, collision.globalIndex());
260+
registry.fill(HIST("EventEfficiency"), 4.);
261+
if (perCollisionSample.size() > 0) {
262+
atLeastOne_gt0 = true;
263+
registry.fill(HIST("EventEfficiency"), 5.);
264+
}
265+
registry.fill(HIST("EventsNtrkZvtxGen"), perCollisionSample.size(), collision.posZ());
224266
}
225267
}
268+
if (collisions.size() == 0) {
269+
registry.fill(HIST("NotFoundEventZvtx"), mcCollision.posZ());
270+
}
226271
for (auto& particle : particles) {
227272
auto p = pdg->GetParticle(particle.pdgCode());
228-
int charge = 0;
229-
if (p == nullptr) {
230-
// unknown particles will be skipped
231-
if (particle.pdgCode() > 1000000000) {
232-
LOGP(debug, "[{}] Nucleus with PDG code {}", particle.globalIndex(), particle.pdgCode());
233-
} else {
234-
LOGP(debug, "[{}] Unknown particle with PDG code {}", particle.globalIndex(), particle.pdgCode());
235-
}
236-
} else {
237-
charge = p->Charge();
273+
auto charge = 0;
274+
if (p != nullptr) {
275+
charge = (int)p->Charge();
238276
}
239-
if (charge != 0 && particle.isPhysicalPrimary() && (particle.eta() < etaMax) && (particle.eta() > etaMin)) {
240-
registry.fill(HIST("TracksEtaZvtxGen"), particle.eta(), mcCollision.posZ());
277+
if (charge != 0) {
278+
registry.fill(HIST("TracksEtaZvtxGen_t"), particle.eta(), mcCollision.posZ());
279+
if (perCollisionMCSample.size() > 0) {
280+
registry.fill(HIST("TracksEtaZvtxGen_gt0t"), particle.eta(), mcCollision.posZ());
281+
}
282+
if (atLeastOne) {
283+
registry.fill(HIST("TracksEtaZvtxGen"), particle.eta(), mcCollision.posZ());
284+
if (atLeastOne_gt0) {
285+
registry.fill(HIST("TracksEtaZvtxGen_gt0"), particle.eta(), mcCollision.posZ());
286+
}
287+
}
241288
registry.fill(HIST("TracksPhiEtaGen"), particle.phi(), particle.eta());
242289
}
243290
}

0 commit comments

Comments
 (0)