Skip to content

Commit 26077ec

Browse files
victor-gonzalezVictor
victor-gonzalez
and
Victor
authored
[PWGCF] DptDptCorrelations, support for on the fly MC (#8348)
Co-authored-by: Victor <victor@cern.ch>
1 parent 2fcd170 commit 26077ec

File tree

1 file changed

+206
-4
lines changed

1 file changed

+206
-4
lines changed

‎PWGCF/TableProducer/dptdptfilter.cxx‎

Lines changed: 206 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -11,10 +11,13 @@
1111

1212
#include <cmath>
1313
#include <algorithm>
14+
#include <string>
15+
#include <vector>
1416

1517
#include "Framework/AnalysisTask.h"
1618
#include "Framework/AnalysisDataModel.h"
1719
#include "Framework/ASoAHelpers.h"
20+
#include "CommonConstants/PhysicsConstants.h"
1821
#include "Common/DataModel/EventSelection.h"
1922
#include "Common/DataModel/Centrality.h"
2023
#include "Common/Core/TableHelper.h"
@@ -29,6 +32,7 @@
2932
#include "Framework/RunningWorkflowInfo.h"
3033
#include <TROOT.h>
3134
#include <TDatabasePDG.h>
35+
#include <TPDGCode.h>
3236
#include <TParameter.h>
3337
#include <TList.h>
3438
#include <TDirectory.h>
@@ -163,6 +167,169 @@ std::vector<int> partMultNeg; // multiplicity of negative particles
163167

164168
using namespace dptdptfilter;
165169

170+
//////////////////////////////////////////////////////////////////////////////
171+
// Multiplicity in principle for on the fly generated events
172+
//////////////////////////////////////////////////////////////////////////////
173+
174+
struct Multiplicity {
175+
enum multest {
176+
kV0M,
177+
kCL1,
178+
kCL1GAP
179+
};
180+
181+
float getMultiplicityClass() { return multiplicityclass; }
182+
float getMultiplicity() { return multiplicity; }
183+
184+
multest classestimator = kV0M;
185+
186+
float multiplicityclass = -1.0;
187+
float multiplicity = 0.0;
188+
bool inelgth0 = false;
189+
int V0AM = 0;
190+
int V0CM = 0;
191+
int CL1M = 0;
192+
int CL1EtaGapM = 0;
193+
int dNchdEta = 0;
194+
int nPart = 0;
195+
TH1F* fhNPartTot = nullptr; ///< total number of particles analyzed
196+
TH1F* fhMultiplicity; ///< the multiplicity distribution
197+
TH2F* fhV0Multiplicity; ///< the V0M multiplicity histogram
198+
TH2F* fhCL1Multiplicity; ///< the CL1 multiplicity histogram
199+
TH2F* fhCL1EtaGapMultiplicity; ///< the CL1 with an eta gap multiplicity histogram
200+
const TH1* fhV0MMultPercentile; ///< the V0M Centrality / Multiplicity percentile estimation histogram
201+
const TH1* fhCL1MultPercentile; ///< the CL1 Centrality / Multiplicity percentile estimation histogram
202+
const TH1* fhCL1EtaGapMultPercentile; ///< the CL1 with an eta gap Centrality / Multiplicity percentile estimation histogram
203+
204+
void init(TList* hlist)
205+
{
206+
fhNPartTot = new TH1F("CollisionNpart", "Collision analyzed particles;number of particles;counts", 8000, -0.5, 8000 - 0.5);
207+
fhMultiplicity = new TH1F("CollisionMultiplicity", "Event multiplicity;multiplicity (%);counts", 101, -0.5, 101 - 0.5);
208+
fhV0Multiplicity = new TH2F("V0Multiplicity", "V0M;V0M;d#it{N}/d#eta;counts", 3000, -9.5, 3000 - 9.5, 2500, -9.5, 2500 - 9.5);
209+
fhCL1Multiplicity = new TH2F("CL1Multiplicity", "CL1M;CL1M;d#it{N}/d#eta;counts", 3000, -9.5, 3000 - 9.5, 2500, -9.5, 2500 - 9.5);
210+
fhCL1EtaGapMultiplicity = new TH2F("CL1EtaGapMultiplicity", "CL1M (excl |#eta|<0.8);CL1M;d#it{N}/d#eta;counts", 3000, -9.5, 3000 - 9.5, 2500, -9.5, 2500 - 9.5);
211+
212+
hlist->Add(fhNPartTot);
213+
hlist->Add(fhMultiplicity);
214+
hlist->Add(fhV0Multiplicity);
215+
hlist->Add(fhCL1Multiplicity);
216+
hlist->Add(fhCL1EtaGapMultiplicity);
217+
}
218+
219+
void setMultiplicityPercentiles(TList* list)
220+
{
221+
LOGF(info, "setMultiplicityPercentiles()", "From list %s", list->GetName());
222+
fhV0MMultPercentile = reinterpret_cast<TH1*>(list->FindObject("V0MCentMult"));
223+
fhCL1MultPercentile = reinterpret_cast<TH1*>(list->FindObject("CL1MCentMult"));
224+
fhCL1EtaGapMultPercentile = reinterpret_cast<TH1*>(list->FindObject("CL1EtaGapMCentMult"));
225+
226+
if (fhV0MMultPercentile == nullptr || fhCL1MultPercentile == nullptr || fhCL1EtaGapMultPercentile == nullptr) {
227+
LOGF(fatal, "setMultiplicityPercentiles()", "Percentiles histograms not correctly loaded. ABORTING!!!");
228+
return;
229+
}
230+
}
231+
232+
template <typename Particle>
233+
bool addParticleToMultiplicity(const Particle& p)
234+
{
235+
/* on the fly MC production */
236+
/* get event multiplicity according to the passed eta range */
237+
/* event multiplicity as number of primary charged particles */
238+
/* based on AliAnalysisTaskPhiCorrelations implementation */
239+
int pdgcode = std::abs(p.pdgCode());
240+
auto addTo = [](const Particle& p, int& est, float etamin, float etamax) {
241+
if (p.eta() < etamax && etamin < p.eta()) {
242+
est = est + 1;
243+
}
244+
};
245+
246+
/* pdg checks */
247+
switch (pdgcode) {
248+
case kPiPlus:
249+
case kKPlus:
250+
case kProton:
251+
/* not clear if we should use IsPhysicalPrimary here */
252+
/* TODO: adapt to FT0M Run 3 and other estimators */
253+
if (0.001 < p.pt() && p.pt() < 50.0) {
254+
if (p.eta() < 1.0 && -1.0 < p.eta()) {
255+
inelgth0 = true;
256+
}
257+
addTo(p, V0AM, 2.8, 5.1);
258+
addTo(p, V0CM, -3.7, -1.7);
259+
addTo(p, CL1M, -1.4, 1.4);
260+
addTo(p, CL1EtaGapM, -1.4, -0.8);
261+
addTo(p, CL1EtaGapM, 0.8, 1.4);
262+
addTo(p, dNchdEta, -0.5, 0.5);
263+
nPart++;
264+
}
265+
break;
266+
default:
267+
break;
268+
}
269+
return true;
270+
}
271+
272+
template <typename CollisionParticles>
273+
void extractMultiplicity(const CollisionParticles& particles)
274+
{
275+
multiplicityclass = 105;
276+
multiplicity = 0;
277+
inelgth0 = false;
278+
nPart = 0;
279+
V0AM = 0;
280+
V0CM = 0;
281+
CL1M = 0;
282+
CL1EtaGapM = 0;
283+
dNchdEta = 0;
284+
285+
for (auto particle : particles) {
286+
addParticleToMultiplicity(particle);
287+
}
288+
289+
if (inelgth0) {
290+
if (fhNPartTot != nullptr) {
291+
fhNPartTot->Fill(nPart);
292+
}
293+
if (fhV0Multiplicity != nullptr) {
294+
fhV0Multiplicity->Fill(V0AM + V0CM, dNchdEta);
295+
}
296+
if (fhCL1Multiplicity != nullptr) {
297+
fhCL1Multiplicity->Fill(CL1M, dNchdEta);
298+
}
299+
if (fhCL1EtaGapMultiplicity != nullptr) {
300+
fhCL1EtaGapMultiplicity->Fill(CL1EtaGapM, dNchdEta);
301+
}
302+
switch (classestimator) {
303+
case kV0M:
304+
if (fhV0MMultPercentile != nullptr) {
305+
multiplicityclass = fhV0MMultPercentile->GetBinContent(fhV0MMultPercentile->FindFixBin(V0AM + V0CM));
306+
multiplicity = V0AM + V0CM;
307+
}
308+
break;
309+
case kCL1:
310+
if (fhCL1MultPercentile != nullptr) {
311+
multiplicityclass = fhCL1MultPercentile->GetBinContent(fhCL1MultPercentile->FindFixBin(CL1M));
312+
multiplicity = CL1M;
313+
}
314+
break;
315+
case kCL1GAP:
316+
if (fhCL1EtaGapMultPercentile != nullptr) {
317+
multiplicityclass = fhCL1EtaGapMultPercentile->GetBinContent(fhCL1EtaGapMultPercentile->FindFixBin(CL1EtaGapM));
318+
multiplicity = CL1EtaGapM;
319+
}
320+
break;
321+
default:
322+
break;
323+
}
324+
fhMultiplicity->Fill(multiplicityclass);
325+
}
326+
}
327+
};
328+
329+
//////////////////////////////////////////////////////////////////////////////
330+
// The filter class
331+
//////////////////////////////////////////////////////////////////////////////
332+
166333
struct DptDptFilter {
167334
struct : ConfigurableGroup {
168335
Configurable<std::string> cfgCCDBUrl{"input_ccdburl", "http://ccdb-test.cern.ch:8080", "The CCDB url for the input file"};
@@ -189,6 +356,8 @@ struct DptDptFilter {
189356
Produces<aod::DptDptCFAcceptedTrueCollisions> acceptedtrueevents;
190357
Produces<aod::DptDptCFGenCollisionsInfo> gencollisionsinfo;
191358

359+
Multiplicity multiplicity;
360+
192361
void init(InitContext const&)
193362
{
194363
using namespace dptdptfilter;
@@ -212,6 +381,8 @@ struct DptDptFilter {
212381
/* the centrality/multiplicity estimation */
213382
if (doprocessWithoutCent || doprocessWithoutCentDetectorLevel || doprocessWithoutCentGeneratorLevel) {
214383
fCentMultEstimator = kNOCM;
384+
} else if (doprocessOnTheFlyGeneratorLevel) {
385+
fCentMultEstimator = kFV0A;
215386
} else {
216387
fCentMultEstimator = getCentMultEstimator(cfgCentMultEstimator);
217388
}
@@ -271,14 +442,20 @@ struct DptDptFilter {
271442

272443
fhTrueVertexZB = new TH1F("TrueVertexZB", "Vertex Z before (truth); z_{vtx}", 60, -15, 15);
273444
fhTrueVertexZA = new TH1F("TrueVertexZA", "Vertex Z (truth); z_{vtx}", zvtxbins, zvtxlow, zvtxup);
274-
fhTrueVertexZAA = new TH1F("TrueVertexZAA", "Vertex Z (truth rec associated); z_{vtx}", zvtxbins, zvtxlow, zvtxup);
445+
if (!doprocessOnTheFlyGeneratorLevel) {
446+
fhTrueVertexZAA = new TH1F("TrueVertexZAA", "Vertex Z (truth rec associated); z_{vtx}", zvtxbins, zvtxlow, zvtxup);
447+
}
275448

276449
/* add the hstograms to the output list */
277450
fOutputList->Add(fhTrueCentMultB);
278451
fOutputList->Add(fhTrueCentMultA);
279452
fOutputList->Add(fhTrueVertexZB);
280453
fOutputList->Add(fhTrueVertexZA);
281-
fOutputList->Add(fhTrueVertexZAA);
454+
if (doprocessOnTheFlyGeneratorLevel) {
455+
multiplicity.init(fOutputList);
456+
} else {
457+
fOutputList->Add(fhTrueVertexZAA);
458+
}
282459
}
283460
}
284461

@@ -304,7 +481,7 @@ struct DptDptFilter {
304481
PROCESS_SWITCH(DptDptFilter, processWithoutCentDetectorLevel, "Process MC detector level without centrality", false);
305482

306483
template <typename CollisionObject, typename ParticlesList>
307-
void processGenerated(CollisionObject const& mccollision, ParticlesList const& mcparticles, float centormult);
484+
bool processGenerated(CollisionObject const& mccollision, ParticlesList const& mcparticles, float centormult);
308485

309486
template <typename CollisionsGroup, typename AllCollisions>
310487
void processGeneratorLevel(aod::McCollision const& mccollision,
@@ -331,6 +508,10 @@ struct DptDptFilter {
331508
aod::CollisionsEvSel const& allcollisions);
332509
PROCESS_SWITCH(DptDptFilter, processWithoutCentGeneratorLevel, "Process generated without centrality", false);
333510

511+
void processOnTheFlyGeneratorLevel(aod::McCollision const& mccollision,
512+
aod::McParticles const& mcparticles);
513+
PROCESS_SWITCH(DptDptFilter, processOnTheFlyGeneratorLevel, "Process on the fly generated events", false);
514+
334515
void processVertexGenerated(aod::McCollisions const&);
335516
PROCESS_SWITCH(DptDptFilter, processVertexGenerated, "Process vertex generator level", false);
336517
};
@@ -398,7 +579,7 @@ void DptDptFilter::processWithoutCentDetectorLevel(aod::CollisionEvSel const& co
398579
}
399580

400581
template <typename CollisionObject, typename ParticlesList>
401-
void DptDptFilter::processGenerated(CollisionObject const& mccollision, ParticlesList const&, float centormult)
582+
bool DptDptFilter::processGenerated(CollisionObject const& mccollision, ParticlesList const&, float centormult)
402583
{
403584
using namespace dptdptfilter;
404585

@@ -411,6 +592,7 @@ void DptDptFilter::processGenerated(CollisionObject const& mccollision, Particle
411592
} else {
412593
gencollisionsinfo(acceptedevent, centormult);
413594
}
595+
return static_cast<bool>(acceptedevent);
414596
}
415597

416598
template <typename CollisionsGroup, typename AllCollisions>
@@ -471,6 +653,26 @@ void DptDptFilter::processWithoutCentGeneratorLevel(aod::McCollision const& mcco
471653
processGeneratorLevel(mccollision, collisions, mcparticles, allcollisions, 50.0);
472654
}
473655

656+
void DptDptFilter::processOnTheFlyGeneratorLevel(aod::McCollision const& mccollision,
657+
aod::McParticles const& mcparticles)
658+
{
659+
uint8_t acceptedEvent = uint8_t(false);
660+
fhTrueVertexZB->Fill(mccollision.posZ());
661+
/* we assign a default value for the time being */
662+
float centormult = 50.0f;
663+
if (IsEvtSelected(mccollision, centormult)) {
664+
acceptedEvent = true;
665+
multiplicity.extractMultiplicity(mcparticles);
666+
fhTrueVertexZA->Fill((mccollision.posZ()));
667+
centormult = multiplicity.getMultiplicityClass();
668+
}
669+
if (fullDerivedData) {
670+
acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), acceptedEvent, centormult);
671+
} else {
672+
gencollisionsinfo(acceptedEvent, centormult);
673+
}
674+
}
675+
474676
void DptDptFilter::processVertexGenerated(aod::McCollisions const& mccollisions)
475677
{
476678
for (aod::McCollision const& mccollision : mccollisions) {

0 commit comments

Comments
 (0)