Skip to content

Commit 34b0c98

Browse files
njacazionjacazio
authored
[PWGLF] K0s ME add MC (#9732)
1 parent ab1c268 commit 34b0c98

File tree

1 file changed

+105
-24
lines changed

1 file changed

+105
-24
lines changed

‎PWGLF/Tasks/Strangeness/k0_mixed_events.cxx‎

Lines changed: 105 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,8 @@
3737

3838
#include "PWGCF/Femto3D/DataModel/singletrackselector.h"
3939
#include "PWGCF/Femto3D/Core/femto3dPairTask.h"
40+
#include "Common/DataModel/Centrality.h"
41+
#include "PWGLF/DataModel/mcCentrality.h"
4042

4143
using namespace o2;
4244
using namespace o2::soa;
@@ -64,6 +66,7 @@ class ResoPair : public MyFemtoPair
6466
bool isClosePair() const { return MyFemtoPair::IsClosePair(mDeltaEta, mDeltaPhi, mRadius); }
6567
void setEtaDiff(const float deta) { mDeltaEta = deta; }
6668
void setPhiStarDiff(const float dphi) { mDeltaPhi = dphi; }
69+
void setRadius(const float r) { mRadius = r; }
6770
void setPair(trkType const& first, trkType const& second)
6871
{
6972
MyFemtoPair::SetPair(first, second);
@@ -121,18 +124,20 @@ struct K0MixedEvents {
121124
Configurable<int> _particlePDGtoReject{"particlePDGtoRejectFromSecond", 0, "applied only if the particles are non-identical and only to the second particle in the pair!!!"};
122125
Configurable<std::vector<float>> _rejectWithinNsigmaTOF{"rejectWithinNsigmaTOF", std::vector<float>{-0.0f, 0.0f}, "TOF rejection Nsigma range for the particle specified with PDG to be rejected"};
123126

124-
Configurable<float> _deta{"deta", 0.01, "minimum allowed defference in eta between two tracks in a pair"};
125-
Configurable<float> _dphi{"dphi", 0.01, "minimum allowed defference in phi_star between two tracks in a pair"};
127+
Configurable<float> _deta{"deta", 1, "minimum allowed defference in eta between two tracks in a pair"};
128+
Configurable<float> _dphi{"dphi", 1, "minimum allowed defference in phi_star between two tracks in a pair"};
126129
Configurable<float> _radiusTPC{"radiusTPC", 1.2, "TPC radius to calculate phi_star for"};
127130

131+
Configurable<bool> useCentralityInvMass{"useCentralityInvMass", true, "Use the centrality vs inv. mass plots"};
128132
Configurable<bool> doMixedEvent{"doMixedEvent", false, "Do the mixed event"};
129133
Configurable<int> _multbinwidth{"multbinwidth", 50, "width of multiplicity bins within which the mixing is done"};
130134
Configurable<int> _vertexbinwidth{"vertexbinwidth", 2, "width of vertexZ bins within which the mixing is done"};
131135

132136
// Binnings
133-
ConfigurableAxis CFkStarBinning{"CFkStarBinning", {500, 0.4, 0.6}, "k* binning of the CF (Nbins, lowlimit, uplimit)"};
137+
ConfigurableAxis invMassBinning{"invMassBinning", {500, 0.4, 0.6}, "k* binning of the CF (Nbins, lowlimit, uplimit)"};
134138
ConfigurableAxis ptBinning{"ptBinning", {1000, 0.f, 10.f}, "pT binning (Nbins, lowlimit, uplimit)"};
135139
ConfigurableAxis dcaXyBinning{"dcaXyBinning", {100, -1.f, 1.f}, "dcaXY binning (Nbins, lowlimit, uplimit)"};
140+
ConfigurableAxis multPercentileBinning{"multPercentileBinning", {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 70.0, 100.0}, "Binning in multiplicity percentile"};
136141

137142
bool IsIdentical;
138143

@@ -180,43 +185,65 @@ struct K0MixedEvents {
180185
Pair->SetIdentical(IsIdentical);
181186
Pair->SetPDG1(_particlePDG_1);
182187
Pair->SetPDG2(_particlePDG_2);
183-
Pair->setEtaDiff(1);
188+
Pair->setEtaDiff(_deta);
189+
Pair->setPhiStarDiff(_dphi);
190+
Pair->setRadius(_radiusTPC);
184191

185192
TPCcuts_1 = std::make_pair(_particlePDG_1, _tpcNSigma_1);
186193
TOFcuts_1 = std::make_pair(_particlePDG_1, _tofNSigma_1);
187194
TPCcuts_2 = std::make_pair(_particlePDG_2, _tpcNSigma_2);
188195
TOFcuts_2 = std::make_pair(_particlePDG_2, _tofNSigma_2);
189196

190-
const AxisSpec invMassAxis{CFkStarBinning, "Inv. mass (GeV/c^{2})"};
197+
const AxisSpec invMassAxis{invMassBinning, "Inv. mass (GeV/c^{2})"};
191198
const AxisSpec ptAxis{ptBinning, "#it{p}_{T} (GeV/c)"};
192199
const AxisSpec dcaXyAxis{dcaXyBinning, "DCA_{xy} (cm)"};
200+
const AxisSpec multPercentileAxis{multPercentileBinning, "Mult. Perc."};
193201

194202
registry.add("Trks", "Trks", kTH1D, {{2, 0.5, 2.5, "Tracks"}});
195-
registry.add("VTXc", "VTXc", kTH1F, {{100, -20., 20., "vtx"}});
196-
registry.add("VTX", "VTX", kTH1F, {{100, -20., 20., "vtx"}});
197-
registry.add("SEcand", "SEcand", kTH1F, {{2, 0.5, 2.5}});
198-
registry.add("SE", "SE", kTH1F, {invMassAxis});
199-
registry.add("ME", "ME", kTH1F, {invMassAxis});
200-
registry.add("SEvsPt", "SEvsPt", kTH2D, {invMassAxis, ptAxis});
201-
registry.add("MEvsPt", "MEvsPt", kTH2D, {invMassAxis, ptAxis});
203+
registry.add("VTXc", "VTXc", kTH1D, {{100, -20., 20., "vtx"}});
204+
registry.add("VTX", "VTX", kTH1D, {{100, -20., 20., "vtx"}});
205+
registry.add("SEcand", "SEcand", kTH1D, {{2, 0.5, 2.5}});
206+
registry.add("SE", "SE", kTH1D, {invMassAxis});
207+
registry.add("ME", "ME", kTH1D, {invMassAxis});
208+
if (useCentralityInvMass) {
209+
registry.add("SEvsPt", "SEvsPt", kTH3F, {invMassAxis, ptAxis, multPercentileAxis});
210+
} else {
211+
registry.add("SEvsPt", "SEvsPt", kTH2D, {invMassAxis, ptAxis});
212+
}
213+
if (doMixedEvent) {
214+
if (useCentralityInvMass) {
215+
registry.add("MEvsPt", "MEvsPt", kTH3F, {invMassAxis, ptAxis, multPercentileAxis});
216+
} else {
217+
registry.add("MEvsPt", "MEvsPt", kTH2D, {invMassAxis, ptAxis});
218+
}
219+
}
202220
registry.add("eta", Form("eta_%i", _particlePDG_1.value), kTH2F, {ptAxis, {100, -10., 10., "#eta"}});
203-
registry.add("p_first", Form("p_%i", _particlePDG_1.value), kTH1F, {ptAxis});
221+
registry.add("p_first", Form("p_%i", _particlePDG_1.value), kTH1D, {ptAxis});
204222
registry.add("dcaXY_first", Form("dca_%i", _particlePDG_1.value), kTH2F, {ptAxis, dcaXyAxis});
205223
registry.add("nsigmaTOF_first", Form("nsigmaTOF_%i", _particlePDG_1.value), kTH2F, {ptAxis, {100, -10., 10., Form("N#sigma_{TOF}(%s))", pdgToSymbol(_particlePDG_1))}});
206224
registry.add("nsigmaTPC_first", Form("nsigmaTPC_%i", _particlePDG_1.value), kTH2F, {ptAxis, {100, -10., 10., Form("N#sigma_{TPC}(%s))", pdgToSymbol(_particlePDG_1))}});
207225
registry.add("rapidity_first", Form("rapidity_%i", _particlePDG_1.value), kTH2F, {ptAxis, {100, -10., 10., Form("y(%s)", pdgToSymbol(_particlePDG_1))}});
208226

209227
if (!IsIdentical) {
210-
registry.add("p_second", Form("p_%i", _particlePDG_2.value), kTH1F, {ptAxis});
228+
registry.add("p_second", Form("p_%i", _particlePDG_2.value), kTH1D, {ptAxis});
211229
registry.add("dcaXY_second", Form("dca_%i", _particlePDG_2.value), kTH2F, {ptAxis, dcaXyAxis});
212230
registry.add("nsigmaTOF_second", Form("nsigmaTOF_%i", _particlePDG_2.value), kTH2F, {ptAxis, {100, -10., 10., Form("N#sigma_{TOF}(%s))", pdgToSymbol(_particlePDG_2))}});
213231
registry.add("nsigmaTPC_second", Form("nsigmaTPC_%i", _particlePDG_2.value), kTH2F, {ptAxis, {100, -10., 10., Form("N#sigma_{TPC}(%s))", pdgToSymbol(_particlePDG_2))}});
214232
registry.add("rapidity_second", Form("rapidity_%i", _particlePDG_2.value), kTH2F, {ptAxis, {100, -10., 10., Form("y(%s)", pdgToSymbol(_particlePDG_2))}});
215233
}
234+
235+
if (!doprocessMCReco) {
236+
return;
237+
}
238+
if (useCentralityInvMass) {
239+
registry.add("MC/generatedInRecoEvs", "generatedInRecoEvs", kTH2D, {ptAxis, multPercentileAxis});
240+
} else {
241+
registry.add("MC/generatedInRecoEvs", "generatedInRecoEvs", kTH2D, {ptAxis});
242+
}
216243
}
217244

218245
template <typename Type>
219-
void mixTracks(Type const& tracks)
246+
void mixTracks(Type const& tracks, const float centrality)
220247
{ // template for identical particles from the same collision
221248

222249
LOG(debug) << "Mixing tracks of the same event";
@@ -233,14 +260,18 @@ struct K0MixedEvents {
233260
continue;
234261
}
235262
registry.fill(HIST("SEcand"), 2.f);
236-
registry.fill(HIST("SE"), Pair->getInvMass()); // close pair rejection and fillig the SE histo
237-
registry.fill(HIST("SEvsPt"), Pair->getInvMass(), Pair->getPt()); // close pair rejection and fillig the SE histo
263+
registry.fill(HIST("SE"), Pair->getInvMass()); // close pair rejection and fillig the SE histo
264+
if (useCentralityInvMass) {
265+
registry.fill(HIST("SEvsPt"), Pair->getInvMass(), Pair->getPt(), centrality); // close pair rejection and fillig the SE histo
266+
} else {
267+
registry.fill(HIST("SEvsPt"), Pair->getInvMass(), Pair->getPt()); // close pair rejection and fillig the SE histo
268+
}
238269
}
239270
}
240271
}
241272

242273
template <bool isSameEvent = false, typename Type>
243-
void mixTracks(Type const& tracks1, Type const& tracks2)
274+
void mixTracks(Type const& tracks1, Type const& tracks2, const float centrality)
244275
{
245276
LOG(debug) << "Mixing tracks of two different events";
246277
for (auto trk1 : tracks1) {
@@ -260,10 +291,18 @@ struct K0MixedEvents {
260291
if constexpr (isSameEvent) {
261292
registry.fill(HIST("SEcand"), 2.f);
262293
registry.fill(HIST("SE"), Pair->getInvMass());
263-
registry.fill(HIST("SEvsPt"), Pair->getInvMass(), Pair->getPt());
294+
if (useCentralityInvMass) {
295+
registry.fill(HIST("SEvsPt"), Pair->getInvMass(), Pair->getPt(), centrality);
296+
} else {
297+
registry.fill(HIST("SEvsPt"), Pair->getInvMass(), Pair->getPt());
298+
}
264299
} else {
265300
registry.fill(HIST("ME"), Pair->getInvMass());
266-
registry.fill(HIST("MEvsPt"), Pair->getInvMass(), Pair->getPt());
301+
if (useCentralityInvMass) {
302+
registry.fill(HIST("MEvsPt"), Pair->getInvMass(), Pair->getPt(), centrality);
303+
} else {
304+
registry.fill(HIST("MEvsPt"), Pair->getInvMass(), Pair->getPt());
305+
}
267306
}
268307
}
269308
}
@@ -406,7 +445,7 @@ struct K0MixedEvents {
406445
Pair->SetMagField1(col1->magField());
407446
Pair->SetMagField2(col1->magField());
408447

409-
mixTracks(selectedtracks_1[col1->index()]); // mixing SE identical
448+
mixTracks(selectedtracks_1[col1->index()], col1->multPerc()); // mixing SE identical
410449
if (!doMixedEvent) {
411450
continue;
412451
}
@@ -416,7 +455,7 @@ struct K0MixedEvents {
416455
auto col2 = (i->second)[indx2];
417456

418457
Pair->SetMagField2(col2->magField());
419-
mixTracks(selectedtracks_1[col1->index()], selectedtracks_1[col2->index()]); // mixing ME identical
458+
mixTracks(selectedtracks_1[col1->index()], selectedtracks_1[col2->index()], col1->multPerc()); // mixing ME identical
420459
}
421460
}
422461
}
@@ -434,7 +473,7 @@ struct K0MixedEvents {
434473
Pair->SetMagField1(col1->magField());
435474
Pair->SetMagField2(col1->magField());
436475

437-
mixTracks<true>(selectedtracks_1[col1->index()], selectedtracks_2[col1->index()]); // mixing SE non-identical
476+
mixTracks<true>(selectedtracks_1[col1->index()], selectedtracks_2[col1->index()], col1->multPerc()); // mixing SE non-identical
438477
if (!doMixedEvent) {
439478
continue;
440479
}
@@ -444,7 +483,7 @@ struct K0MixedEvents {
444483
auto col2 = (i->second)[indx2];
445484

446485
Pair->SetMagField2(col2->magField());
447-
mixTracks(selectedtracks_1[col1->index()], selectedtracks_2[col2->index()]); // mixing ME non-identical
486+
mixTracks(selectedtracks_1[col1->index()], selectedtracks_2[col2->index()], col1->multPerc()); // mixing ME non-identical
448487
}
449488
}
450489
}
@@ -466,6 +505,48 @@ struct K0MixedEvents {
466505
(i->second).clear();
467506
mixbins.clear();
468507
}
508+
509+
Filter eventFilter = (aod::evsel::sel8 == true && (nabs(o2::aod::collision::posZ) < _vertexZ));
510+
511+
using RecoMCCollisions = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::CentFT0Ms>;
512+
using GenMCCollisions = soa::Join<aod::McCollisions, aod::McCentFT0Ms>;
513+
514+
// Service<o2::framework::O2DatabasePDG> pdgDB;
515+
Preslice<aod::McParticles> perMCCol = aod::mcparticle::mcCollisionId;
516+
SliceCache cache;
517+
void processMCReco(soa::Filtered<RecoMCCollisions> const& collisions,
518+
GenMCCollisions const&,
519+
aod::McParticles const& mcParticles)
520+
{
521+
for (const auto& col : collisions) {
522+
if (!col.has_mcCollision()) {
523+
continue;
524+
}
525+
const auto& mcCollision = col.mcCollision_as<GenMCCollisions>();
526+
const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache);
527+
for (const auto& mcParticle : particlesInCollision) {
528+
switch (mcParticle.pdgCode()) {
529+
case 310:
530+
break;
531+
default:
532+
continue;
533+
}
534+
if (mcParticle.pdgCode() != 310) {
535+
LOG(fatal) << "Fatal in PDG";
536+
}
537+
if (std::abs(mcParticle.y()) > 0.5) {
538+
continue;
539+
}
540+
if (useCentralityInvMass) {
541+
registry.fill(HIST("MC/generatedInRecoEvs"), mcParticle.pt(), col.centFT0M());
542+
} else {
543+
registry.fill(HIST("MC/generatedInRecoEvs"), mcParticle.pt());
544+
}
545+
}
546+
}
547+
}
548+
549+
PROCESS_SWITCH(K0MixedEvents, processMCReco, "process mc", false);
469550
};
470551

471552
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)