3030#include < Framework/runDataProcessing.h>
3131
3232#include < TMCProcess.h>
33+ #include < TPDGCode.h>
3334
3435#include < cstdint>
36+ #include < string>
3537#include < vector>
3638
3739using namespace o2 ;
3840using namespace o2 ::framework;
3941using aod::track::TrackSelectionFlags;
4042
43+ int getMotherCategory (int pdg)
44+ {
45+ switch (std::abs (pdg)) {
46+ case 310 :
47+ return 1 ; // K0S
48+ case 130 :
49+ return 2 ; // K0L
50+ case 3122 :
51+ return 3 ; // Lambda
52+
53+ case 3222 : // Sigma+
54+ case 3112 :
55+ return 4 ; // Sigma-
56+
57+ case 3312 :
58+ return 5 ; // Xi-
59+
60+ case 3334 :
61+ return 6 ; // Omega
62+
63+ case 211 :
64+ return 7 ; // Pion
65+ case 321 :
66+ return 8 ; // Kaon
67+ case 2112 :
68+ return 9 ; // Neutron
69+ case 13 :
70+ return 10 ; // Muon
71+
72+ default :
73+ return 11 ; // Other
74+ }
75+ }
76+
77+ static std::string getMotherName (int category)
78+ {
79+ static const std::vector<std::string> names = {
80+ " K0S" ,
81+ " K0L" ,
82+ " Lambda" ,
83+ " SigmaCharged" ,
84+ " XiCharged" ,
85+ " Omega" ,
86+ " Pi" ,
87+ " K" ,
88+ " n" ,
89+ " Muon" ,
90+ " Other" };
91+
92+ return names.at (category - 1 );
93+ }
94+
4195struct PccQa {
4296 HistogramRegistry histos;
4397 Service<o2::framework::O2DatabasePDG> pdg;
4498
99+ std::vector<std::shared_ptr<TH2 >> motherDCA;
100+ std::vector<std::shared_ptr<TH2 >> motherDCAWeighted;
101+ std::vector<std::shared_ptr<TProfile>> motherWeights;
102+
45103 static constexpr float MaxVtxZ = 10 .f;
46104
47105 Configurable<uint32_t > maxMult{" maxMult" , 100 , " max multiplicity" };
@@ -73,44 +131,76 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
73131void PccQa::init (InitContext const &)
74132{
75133
76- histos. add ( " eventCounter " , " " , kTH1D , {{ 1 , 0.5 , 1.5 }}) ;
77- const AxisSpec dcaAxis = {1000 , -1 ., 1 ., " #it{DCA}_{xy}" , " dca" };
134+ // const AxisSpec dcaAxis = {1000, -1. , 1., "#it{DCA}_{xy}", "dca"} ;
135+ const AxisSpec dcaAxis = {99 , -1 ., 1 ., " #it{DCA}_{xy}" , " dca" };
78136 std::vector<double > ptBinEdges = {0.2 , 0.4 , 0.6 , 0.8 , 1.0 , 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 };
79137 const AxisSpec ptAxis{ptBinEdges, " #it{p}_{T} (GeV/#it{c})" , " pt" };
80138
81- const int nBinsMult = 30 ;
82- const AxisSpec nchAxis = {nBinsMult, 0 ., static_cast <double >(maxMult), " #it{N}_{ch}" , " nch" };
139+ const int nBinsMult = 25 ;
140+ const AxisSpec nchAxis = {nBinsMult, 1 ., static_cast <double >(maxMult), " #it{N}_{ch}" , " nch" };
141+
142+ histos.add (" multDist" , " " , kTH1D , {nchAxis});
83143
84- histos.add (" DCAxyVsPt" , " " , kTH2D , {ptAxis, dcaAxis});
85144 histos.add (" DCAxyVsPtVsMult" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
86145
87146 if (doprocessMC) {
88- histos.add (" DCAxyVsPt_weighted" , " " , kTH2D , {ptAxis, dcaAxis});
89147 histos.add (" DCAxyVsPtVsMult_weighted" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
90148
91- histos.add (" prim/DCAxyVsPt" , " " , kTH2D , {ptAxis, dcaAxis});
92149 histos.add (" prim/DCAxyVsPtVsMult" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
93-
94150 histos.add (" prim/DCAxyVsPtVsMult_weighted" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
95- histos.add (" prim/DCAxyVsPt_weighted" , " " , kTH2D , {ptAxis, dcaAxis});
96151
97- histos.add (" sec/DCAxyVsPt" , " " , kTH2D , {ptAxis, dcaAxis});
98152 histos.add (" sec/DCAxyVsPtVsMult" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
99-
100- histos.add (" sec/DCAxyVsPt_weighted" , " " , kTH2D , {ptAxis, dcaAxis});
101153 histos.add (" sec/DCAxyVsPtVsMult_weighted" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
102154
103- histos.add (" sec/dec/DCAxyVsPt" , " " , kTH2D , {ptAxis, dcaAxis});
104155 histos.add (" sec/dec/DCAxyVsPtVsMult" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
105-
106- histos.add (" sec/dec/DCAxyVsPt_weighted" , " " , kTH2D , {ptAxis, dcaAxis});
107156 histos.add (" sec/dec/DCAxyVsPtVsMult_weighted" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
108157
109- histos.add (" sec/mat/DCAxyVsPt" , " " , kTH2D , {ptAxis, dcaAxis});
110158 histos.add (" sec/mat/DCAxyVsPtVsMult" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
111-
112- histos.add (" sec/mat/DCAxyVsPt_weighted" , " " , kTH2D , {ptAxis, dcaAxis});
113159 histos.add (" sec/mat/DCAxyVsPtVsMult_weighted" , " " , kTH3D , {nchAxis, ptAxis, dcaAxis});
160+
161+ histos.add (" prim/weights" , " " , kTProfile , {ptAxis});
162+ histos.add (" sec/weights" , " " , kTProfile , {ptAxis});
163+ histos.add (" sec/dec/weights" , " " , kTProfile , {ptAxis});
164+ histos.add (" sec/mat/weights" , " " , kTProfile , {ptAxis});
165+
166+ std::vector<std::string> labels = {
167+ " K^{0}_{S}" , " K^{0}_{L}" , " #Lambda" ,
168+ " #Sigma^{#pm}" ,
169+ " #Xi^{#pm}" ,
170+ " #Omega^{#pm}" ,
171+ " #pi^{#pm}" , " K^{#pm}" , " n" ,
172+ " #mu^{#pm}" ,
173+ " Other" };
174+
175+ const int nMotherCategories = static_cast <int >(labels.size ());
176+
177+ histos.add (
178+ " sec/dec/motherPDG" ,
179+ " Mother particles of decay secondaries" ,
180+ kTH1D ,
181+ {{nMotherCategories, 0.5 , nMotherCategories + 0.5 }});
182+
183+ auto hMother = histos.get <TH1 >(HIST (" sec/dec/motherPDG" ));
184+
185+ for (int i = 0 ; i < nMotherCategories; ++i) {
186+ hMother->GetXaxis ()->SetBinLabel (i + 1 , labels[i].c_str ());
187+ }
188+
189+ for (int category = 1 ; category <= nMotherCategories; ++category) {
190+ const auto mother = getMotherName (category);
191+
192+ motherDCA.push_back (histos.add <TH2 >(
193+ fmt::format (" sec/dec/{}/DCAxyVsPt" , mother),
194+ " " , kTH2D , {ptAxis, dcaAxis}));
195+
196+ motherDCAWeighted.push_back (histos.add <TH2 >(
197+ fmt::format (" sec/dec/{}/DCAxyVsPt_weighted" , mother),
198+ " " , kTH2D , {ptAxis, dcaAxis}));
199+
200+ motherWeights.push_back (histos.add <TProfile>(
201+ fmt::format (" sec/dec/{}/weights" , mother),
202+ " " , kTProfile , {ptAxis}));
203+ }
114204 }
115205}
116206
@@ -135,44 +225,109 @@ void PccQa::processMeas(const C& collision, const T& tracks)
135225 if ((std::abs (collision.posZ ()) > MaxVtxZ) || !collision.sel8 ()) {
136226 return ;
137227 }
138- histos.fill (HIST (" eventCounter" ), 1 );
228+
229+ auto MultCounterMeas = 0 ;
139230
140231 for (const auto & track : tracks) {
232+
141233 if (!TrackSelectionFlags::checkFlag (track.trackCutFlag (), TrackSelectionFlags::kGlobalTrackWoDCA )) {
142234 continue ;
143235 }
144- histos.fill (HIST (" DCAxyVsPt" ), track.pt (), track.dcaXY ());
145- histos.fill (HIST (" DCAxyVsPtVsMult" ), tracks.size (), track.pt (), track.dcaXY ());
236+
237+ MultCounterMeas++;
238+ }
239+
240+ /* LOGP(info, "MCB: --------------------------");
241+ LOGP(info, "MCB: ++ New event with Nch = {} ",MultCounterMeas);
242+ LOGP(info, "MCB: --------------------------"); */
243+
244+ histos.fill (HIST (" multDist" ), MultCounterMeas);
245+
246+ for (const auto & track : tracks) {
247+ if (!TrackSelectionFlags::checkFlag (track.trackCutFlag (), TrackSelectionFlags::kGlobalTrackWoDCA )) {
248+ continue ;
249+ }
250+
251+ histos.fill (HIST (" DCAxyVsPtVsMult" ), MultCounterMeas, track.pt (), track.dcaXY ());
146252
147253 if constexpr (IS_MC ) {
148254 if (!track.has_mcParticle ()) {
149255 continue ;
150256 }
257+
151258 const auto & particle = track.template mcParticle_as <ParticleTableMC>();
152259
153- histos.fill (HIST (" DCAxyVsPt_weighted" ), track.pt (), track.dcaXY (), particle.pccWeight ());
154- histos.fill (HIST (" DCAxyVsPtVsMult_weighted" ), tracks.size (), track.pt (), track.dcaXY (), particle.pccWeight ());
260+ histos.fill (HIST (" DCAxyVsPtVsMult_weighted" ), MultCounterMeas, track.pt (), track.dcaXY (), particle.pccWeight ());
155261
156262 if (particle.isPhysicalPrimary ()) {
157- histos.fill (HIST (" prim/DCAxyVsPt" ), track.pt (), track.dcaXY ());
158- histos.fill (HIST (" prim/DCAxyVsPt_weighted" ), track.pt (), track.dcaXY (), particle.pccWeight ());
159- histos.fill (HIST (" prim/DCAxyVsPtVsMult" ), tracks.size (), track.pt (), track.dcaXY ());
160- histos.fill (HIST (" prim/DCAxyVsPtVsMult_weighted" ), tracks.size (), track.pt (), track.dcaXY (), particle.pccWeight ());
263+
264+ histos.fill (HIST (" prim/weights" ), particle.pt (), particle.pccWeight ());
265+ histos.fill (HIST (" prim/DCAxyVsPtVsMult" ), MultCounterMeas, track.pt (), track.dcaXY ());
266+ histos.fill (HIST (" prim/DCAxyVsPtVsMult_weighted" ), MultCounterMeas, track.pt (), track.dcaXY (), particle.pccWeight ());
267+
161268 } else {
162- histos.fill (HIST (" sec/DCAxyVsPt" ), track.pt (), track.dcaXY ());
163- histos.fill (HIST (" sec/DCAxyVsPt_weighted" ), track.pt (), track.dcaXY (), particle.pccWeight ());
164- histos.fill (HIST (" sec/DCAxyVsPtVsMult" ), tracks.size (), track.pt (), track.dcaXY ());
165- histos.fill (HIST (" sec/DCAxyVsPtVsMult_weighted" ), tracks.size (), track.pt (), track.dcaXY (), particle.pccWeight ());
166- if (particle.getProcess () == TMCProcess::kPDecay ) {
167- histos.fill (HIST (" sec/dec/DCAxyVsPt" ), track.pt (), track.dcaXY ());
168- histos.fill (HIST (" sec/dec/DCAxyVsPt_weighted" ), track.pt (), track.dcaXY (), particle.pccWeight ());
169- histos.fill (HIST (" sec/dec/DCAxyVsPtVsMult" ), tracks.size (), track.pt (), track.dcaXY ());
170- histos.fill (HIST (" sec/dec/DCAxyVsPtVsMult_weighted" ), tracks.size (), track.pt (), track.dcaXY (), particle.pccWeight ());
171- } else if (particle.getProcess () == TMCProcess::kPHInhelastic || particle.getProcess () == TMCProcess::kPHadronic || particle.getProcess () == TMCProcess::kPHElastic ) {
172- histos.fill (HIST (" sec/mat/DCAxyVsPt" ), track.pt (), track.dcaXY ());
173- histos.fill (HIST (" sec/mat/DCAxyVsPt_weighted" ), track.pt (), track.dcaXY (), particle.pccWeight ());
174- histos.fill (HIST (" sec/mat/DCAxyVsPtVsMult" ), tracks.size (), track.pt (), track.dcaXY ());
175- histos.fill (HIST (" sec/mat/DCAxyVsPtVsMult_weighted" ), tracks.size (), track.pt (), track.dcaXY (), particle.pccWeight ());
269+
270+ if (!particle.has_mothers ()) {
271+ continue ;
272+ }
273+
274+ histos.fill (HIST (" sec/weights" ), particle.pt (), particle.pccWeight ());
275+
276+ histos.fill (HIST (" sec/DCAxyVsPtVsMult" ), MultCounterMeas, track.pt (), track.dcaXY ());
277+ histos.fill (HIST (" sec/DCAxyVsPtVsMult_weighted" ), MultCounterMeas, track.pt (), track.dcaXY (), particle.pccWeight ());
278+
279+ if (particle.getProcess () == TMCProcess::kPDecay || particle.getProcess () == TMCProcess::kPPrimary ) {
280+
281+ histos.fill (HIST (" sec/dec/DCAxyVsPtVsMult" ),
282+ MultCounterMeas, track.pt (), track.dcaXY ());
283+
284+ histos.fill (HIST (" sec/dec/DCAxyVsPtVsMult_weighted" ),
285+ MultCounterMeas, track.pt (), track.dcaXY (),
286+ particle.pccWeight ());
287+
288+ histos.fill (HIST (" sec/dec/weights" ),
289+ particle.pt (), particle.pccWeight ());
290+
291+ /* LOGP(info,
292+ " -> Particle={}, PDG={}, process={}, "
293+ "nMoms={}, 1. mom ID={}",
294+ particle.globalIndex(),
295+ particle.pdgCode(),
296+ particle.getProcess(),
297+ particle.mothersIds().size(),
298+ particle.mothersIds()[0]); */
299+
300+ // Identify the immediate mother
301+ // auto mother = particle.template mothers_first_as<aod::McParticles>();
302+ auto mother = particle.template mothers_first_as <ParticleTableMC>();
303+
304+ const int motherCategory = getMotherCategory (mother.pdgCode ());
305+ const std::string motherName = getMotherName (motherCategory);
306+
307+ histos.fill (
308+ HIST (" sec/dec/motherPDG" ),
309+ motherCategory);
310+
311+ const int i = motherCategory - 1 ;
312+
313+ /* LOGP(info,
314+ " --> Mother PDG={}",
315+ mother.pdgCode()
316+ ); */
317+
318+ motherDCA[i]->Fill (track.pt (), track.dcaXY ());
319+
320+ motherDCAWeighted[i]->Fill (
321+ track.pt (), track.dcaXY (), particle.pccWeight ());
322+
323+ motherWeights[i]->Fill (
324+ particle.pt (), particle.pccWeight ());
325+
326+ } else {
327+
328+ histos.fill (HIST (" sec/mat/weights" ), particle.pt (), particle.pccWeight ());
329+ histos.fill (HIST (" sec/mat/DCAxyVsPtVsMult" ), MultCounterMeas, track.pt (), track.dcaXY ());
330+ histos.fill (HIST (" sec/mat/DCAxyVsPtVsMult_weighted" ), MultCounterMeas, track.pt (), track.dcaXY (), particle.pccWeight ());
176331 }
177332 }
178333 }
0 commit comments