Skip to content

Commit a77a3a3

Browse files
limhyunji999alibuild
limhyunji999
and
alibuild
authored
[PWGLF] modify INEL>0 part (#14477)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent c25f15d commit a77a3a3

File tree

1 file changed

+116
-57
lines changed

1 file changed

+116
-57
lines changed

‎PWGLF/Tasks/Resonances/rho770analysis.cxx‎

Lines changed: 116 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -12,28 +12,41 @@
1212
/// \file rho770analysis.cxx
1313
/// \brief rho(770)0 analysis in pp 13 & 13.6 TeV
1414
/// \author Hyunji Lim (hyunji.lim@cern.ch)
15-
/// \since 03/12/2025
15+
/// \since 14/01/2026
1616

1717
#include "PWGLF/DataModel/LFResonanceTables.h"
18+
#include "PWGLF/DataModel/mcCentrality.h"
19+
#include "PWGLF/Utils/inelGt.h"
1820

21+
#include "Common/Core/RecoDecay.h"
1922
#include "Common/DataModel/Centrality.h"
2023
#include "Common/DataModel/EventSelection.h"
24+
#include "Common/DataModel/PIDResponse.h"
2125

2226
#include "CommonConstants/PhysicsConstants.h"
2327
#include "DataFormatsParameters/GRPObject.h"
2428
#include "Framework/ASoAHelpers.h"
2529
#include "Framework/AnalysisTask.h"
30+
#include "Framework/O2DatabasePDGPlugin.h"
2631
#include "Framework/runDataProcessing.h"
2732
#include <Framework/Configurable.h>
2833

29-
#include "TVector2.h"
30-
#include <TLorentzVector.h>
34+
#include "Math/Vector4D.h"
3135

3236
using namespace o2;
3337
using namespace o2::framework;
3438
using namespace o2::framework::expressions;
3539
using namespace o2::soa;
3640
using namespace o2::constants::physics;
41+
using LorentzVectorPxPyPzMVector = ROOT::Math::PxPyPzMVector;
42+
using LorentzVectorPxPyPzEVector = ROOT::Math::PxPyPzEVector;
43+
44+
enum class TrackPIDMode : int {
45+
TPCOrTOF = 0,
46+
OnlyTPC = 1,
47+
Combined = 2,
48+
TPCWithTOFVeto = 3
49+
};
3750

3851
struct rho770analysis {
3952
SliceCache cache;
@@ -62,17 +75,18 @@ struct rho770analysis {
6275
// kEtaRange)
6376
Configurable<bool> cfgPVContributor{"cfgPVContributor", true, "PV contributor track selection"}; // PV Contriuibutor
6477
Configurable<bool> cfgGlobalTrack{"cfgGlobalTrack", false, "Global track selection"}; // kGoldenChi2 | kDCAxy | kDCAz
65-
Configurable<int> cfgTPCcluster{"cfgTPCcluster", 0, "Number of TPC cluster"};
66-
Configurable<bool> cfgUseTPCRefit{"cfgUseTPCRefit", false, "Require TPC Refit"};
78+
Configurable<int> cfgTPCcluster{"cfgTPCcluster", 1, "Number of TPC cluster"};
79+
Configurable<bool> cfgUseTPCRefit{"cfgUseTPCRefit", false, "Require TPC Refit"}; // refit is included in global track selection
6780
Configurable<bool> cfgUseITSRefit{"cfgUseITSRefit", false, "Require ITS Refit"};
6881
Configurable<bool> cfgHasTOF{"cfgHasTOF", false, "Require TOF"};
82+
Configurable<int> cfgTPCRows{"cfgTPCRows", 80, "Minimum Number of TPC Crossed Rows "};
6983

7084
// PID
7185
Configurable<double> cMaxTOFnSigmaPion{"cMaxTOFnSigmaPion", 3.0, "TOF nSigma cut for Pion"}; // TOF
7286
Configurable<double> cMaxTPCnSigmaPion{"cMaxTPCnSigmaPion", 5.0, "TPC nSigma cut for Pion"}; // TPC
73-
Configurable<double> cMaxTPCnSigmaPionnoTOF{"cMaxTPCnSigmaPionnoTOF", 2.0, "TPC nSigma cut for Pion in no TOF case"}; // TPC
87+
Configurable<double> cMaxTPCnSigmaPionnoTOF{"cMaxTPCnSigmaPionnoTOF", 3.0, "TPC nSigma cut for Pion in no TOF case"}; // TPC
7488
Configurable<double> nsigmaCutCombinedPion{"nsigmaCutCombinedPion", 3.0, "Combined nSigma cut for Pion"};
75-
Configurable<int> selectType{"selectType", 0, "PID selection type"};
89+
Configurable<int> selectTypeInt{"selectTypeInt", static_cast<int>(TrackPIDMode::OnlyTPC), "Pion PID selection mode: 0=TPCOrTOF, 1=OnlyTPC, 2=Combined, 3=TPCWithTOFVeto"};
7690

7791
// Axis
7892
ConfigurableAxis massAxis{"massAxis", {400, 0.2, 2.2}, "Invariant mass axis"};
@@ -99,6 +113,8 @@ struct rho770analysis {
99113
histos.add("hInvMass_Kstar_LSpp", "Kstar ++ invariant mass", {HistType::kTHnSparseF, {massKstarAxis, ptAxis, centAxis}});
100114
histos.add("hInvMass_Kstar_LSmm", "Kstar -- invariant mass", {HistType::kTHnSparseF, {massKstarAxis, ptAxis, centAxis}});
101115

116+
histos.add("QA/hMultiplicityPercent", "Multiplicity percentile of collision", HistType::kTH1F, {centAxis});
117+
102118
histos.add("QA/Nsigma_TPC_BF", "", {HistType::kTH2F, {pTqaAxis, pidqaAxis}});
103119
histos.add("QA/Nsigma_TOF_BF", "", {HistType::kTH2F, {pTqaAxis, pidqaAxis}});
104120
histos.add("QA/TPC_TOF_BF", "", {HistType::kTH2F, {pidqaAxis, pidqaAxis}});
@@ -128,7 +144,7 @@ struct rho770analysis {
128144
{
129145
if (std::abs(track.pt()) < cfgMinPt)
130146
return false;
131-
if (std::fabs(track.eta()) > cfgMaxEta)
147+
if (std::abs(track.eta()) > cfgMaxEta)
132148
return false;
133149
if (std::abs(track.dcaXY()) > cfgMaxDCArToPVcut)
134150
return false;
@@ -150,26 +166,30 @@ struct rho770analysis {
150166
return false;
151167
if (cfgGlobalTrack && !track.isGlobalTrack())
152168
return false;
169+
if (track.tpcNClsCrossedRows() < cfgTPCRows)
170+
return false;
153171

154172
return true;
155173
}
156174

157175
template <typename TrackType>
158176
bool selPion(const TrackType track)
159177
{
160-
if (selectType == 0) {
178+
const auto mode = static_cast<TrackPIDMode>(selectTypeInt.value);
179+
180+
if (mode == TrackPIDMode::TPCOrTOF) { // TPC or TOF
161181
if (std::fabs(track.tpcNSigmaPi()) >= cMaxTPCnSigmaPion || std::fabs(track.tofNSigmaPi()) >= cMaxTOFnSigmaPion)
162182
return false;
163183
}
164-
if (selectType == 1) {
184+
if (mode == TrackPIDMode::OnlyTPC) { // only TPC
165185
if (std::fabs(track.tpcNSigmaPi()) >= cMaxTPCnSigmaPion)
166186
return false;
167187
}
168-
if (selectType == 2) {
188+
if (mode == TrackPIDMode::Combined) { // combining
169189
if (track.tpcNSigmaPi() * track.tpcNSigmaPi() + track.tofNSigmaPi() * track.tofNSigmaPi() >= nsigmaCutCombinedPion * nsigmaCutCombinedPion)
170190
return false;
171191
}
172-
if (selectType == 3) {
192+
if (mode == TrackPIDMode::TPCWithTOFVeto) { // TPC TOF veto
173193
if (track.hasTOF()) {
174194
if (std::fabs(track.tpcNSigmaPi()) >= cMaxTPCnSigmaPion || std::fabs(track.tofNSigmaPi()) >= cMaxTOFnSigmaPion)
175195
return false;
@@ -184,19 +204,21 @@ struct rho770analysis {
184204
template <typename TrackType>
185205
bool selKaon(const TrackType track)
186206
{
187-
if (selectType == 0) {
207+
const auto mode = static_cast<TrackPIDMode>(selectTypeInt.value);
208+
209+
if (mode == TrackPIDMode::TPCOrTOF) {
188210
if (std::fabs(track.tpcNSigmaKa()) >= cMaxTPCnSigmaPion || std::fabs(track.tofNSigmaKa()) >= cMaxTOFnSigmaPion)
189211
return false;
190212
}
191-
if (selectType == 1) {
213+
if (mode == TrackPIDMode::OnlyTPC) {
192214
if (std::fabs(track.tpcNSigmaKa()) >= cMaxTPCnSigmaPion)
193215
return false;
194216
}
195-
if (selectType == 2) {
217+
if (mode == TrackPIDMode::Combined) {
196218
if (track.tpcNSigmaKa() * track.tpcNSigmaKa() + track.tofNSigmaKa() * track.tofNSigmaKa() >= nsigmaCutCombinedPion * nsigmaCutCombinedPion)
197219
return false;
198220
}
199-
if (selectType == 3) {
221+
if (mode == TrackPIDMode::TPCWithTOFVeto) {
200222
if (track.hasTOF()) {
201223
if (std::fabs(track.tpcNSigmaKa()) >= cMaxTPCnSigmaPion || std::fabs(track.tofNSigmaKa()) >= cMaxTOFnSigmaPion)
202224
return false;
@@ -208,10 +230,12 @@ struct rho770analysis {
208230
return true;
209231
}
210232

211-
template <bool IsMC, typename CollisionType, typename TracksType>
212-
void fillHistograms(const CollisionType& collision, const TracksType& dTracks)
233+
template <bool IsMC, typename CollisionType, typename CenMult, typename TracksType>
234+
void fillHistograms(const CollisionType& /*collision*/, const CenMult& multiplicity, const TracksType& dTracks)
213235
{
214-
TLorentzVector part1, part2, reco;
236+
histos.fill(HIST("QA/hMultiplicityPercent"), multiplicity);
237+
238+
LorentzVectorPxPyPzMVector part1, part2, reco;
215239
for (const auto& [trk1, trk2] : combinations(CombinationsUpperIndexPolicy(dTracks, dTracks))) {
216240

217241
if (trk1.index() == trk2.index()) {
@@ -232,106 +256,141 @@ struct rho770analysis {
232256
continue;
233257

234258
if (selPion(trk1) && selPion(trk2)) {
235-
part1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massPi);
236-
part2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massPi);
259+
260+
part1 = LorentzVectorPxPyPzMVector(trk1.px(), trk1.py(), trk1.pz(), massPi);
261+
part2 = LorentzVectorPxPyPzMVector(trk2.px(), trk2.py(), trk2.pz(), massPi);
237262
reco = part1 + part2;
238263

239264
if (reco.Rapidity() > cfgMaxRap || reco.Rapidity() < cfgMinRap)
240265
continue;
241266

242267
if (trk1.sign() * trk2.sign() < 0) {
243-
histos.fill(HIST("hInvMass_rho770_US"), reco.M(), reco.Pt(), collision.cent());
244-
histos.fill(HIST("hInvMass_K0s_US"), reco.M(), reco.Pt(), collision.cent());
268+
histos.fill(HIST("hInvMass_rho770_US"), reco.M(), reco.Pt(), multiplicity);
269+
histos.fill(HIST("hInvMass_K0s_US"), reco.M(), reco.Pt(), multiplicity);
245270

246271
if constexpr (IsMC) {
247272
if (trk1.motherId() != trk2.motherId())
248273
continue;
249-
if (std::abs(trk1.pdgCode()) == 211 && std::abs(trk2.pdgCode()) == 211) {
250-
if (std::abs(trk1.motherPDG()) == 113) {
251-
histos.fill(HIST("MCL/hpT_rho770_REC"), reco.M(), reco.Pt(), collision.cent());
252-
} else if (std::abs(trk1.motherPDG()) == 223) {
253-
histos.fill(HIST("MCL/hpT_omega_REC"), reco.M(), reco.Pt(), collision.cent());
254-
} else if (std::abs(trk1.motherPDG()) == 310) {
255-
histos.fill(HIST("MCL/hpT_K0s_REC"), reco.M(), reco.Pt(), collision.cent());
256-
histos.fill(HIST("MCL/hpT_K0s_pipi_REC"), reco.M(), reco.Pt(), collision.cent());
274+
if (std::abs(trk1.pdgCode()) == kPiPlus && std::abs(trk2.pdgCode()) == kPiPlus) {
275+
if (std::abs(trk1.motherPDG()) == kRho770_0) {
276+
histos.fill(HIST("MCL/hpT_rho770_REC"), reco.M(), reco.Pt(), multiplicity);
277+
} else if (std::abs(trk1.motherPDG()) == kOmega) {
278+
histos.fill(HIST("MCL/hpT_omega_REC"), reco.M(), reco.Pt(), multiplicity);
279+
} else if (std::abs(trk1.motherPDG()) == kK0Short) {
280+
histos.fill(HIST("MCL/hpT_K0s_REC"), reco.M(), reco.Pt(), multiplicity);
281+
histos.fill(HIST("MCL/hpT_K0s_pipi_REC"), reco.M(), reco.Pt(), multiplicity);
257282
}
258-
} else if ((std::abs(trk1.pdgCode()) == 211 && std::abs(trk2.pdgCode()) == 321) || (std::abs(trk1.pdgCode()) == 321 && std::abs(trk2.pdgCode()) == 211)) {
259-
if (std::abs(trk1.motherPDG()) == 313) {
260-
histos.fill(HIST("MCL/hpT_Kstar_REC"), reco.M(), reco.Pt(), collision.cent());
283+
} else if ((std::abs(trk1.pdgCode()) == kPiPlus && std::abs(trk2.pdgCode()) == kKPlus) || (std::abs(trk1.pdgCode()) == kKPlus && std::abs(trk2.pdgCode()) == kPiPlus)) {
284+
if (std::abs(trk1.motherPDG()) == kK0Star892) {
285+
histos.fill(HIST("MCL/hpT_Kstar_REC"), reco.M(), reco.Pt(), multiplicity);
261286
}
262287
}
263288
}
264289
} else if (trk1.sign() > 0 && trk2.sign() > 0) {
265-
histos.fill(HIST("hInvMass_rho770_LSpp"), reco.M(), reco.Pt(), collision.cent());
266-
histos.fill(HIST("hInvMass_K0s_LSpp"), reco.M(), reco.Pt(), collision.cent());
290+
histos.fill(HIST("hInvMass_rho770_LSpp"), reco.M(), reco.Pt(), multiplicity);
291+
histos.fill(HIST("hInvMass_K0s_LSpp"), reco.M(), reco.Pt(), multiplicity);
267292
} else if (trk1.sign() < 0 && trk2.sign() < 0) {
268-
histos.fill(HIST("hInvMass_rho770_LSmm"), reco.M(), reco.Pt(), collision.cent());
269-
histos.fill(HIST("hInvMass_K0s_LSmm"), reco.M(), reco.Pt(), collision.cent());
293+
histos.fill(HIST("hInvMass_rho770_LSmm"), reco.M(), reco.Pt(), multiplicity);
294+
histos.fill(HIST("hInvMass_K0s_LSmm"), reco.M(), reco.Pt(), multiplicity);
270295
}
271296
}
272297

273298
if ((selPion(trk1) && selKaon(trk2)) || (selKaon(trk1) && selPion(trk2))) {
274299
if (selPion(trk1) && selKaon(trk2)) {
275-
part1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massPi);
276-
part2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massKa);
300+
part1 = LorentzVectorPxPyPzMVector(trk1.px(), trk1.py(), trk1.pz(), massPi);
301+
part2 = LorentzVectorPxPyPzMVector(trk2.px(), trk2.py(), trk2.pz(), massKa);
277302
} else if (selKaon(trk1) && selPion(trk2)) {
278-
part1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa);
279-
part2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massPi);
303+
part1 = LorentzVectorPxPyPzMVector(trk1.px(), trk1.py(), trk1.pz(), massKa);
304+
part2 = LorentzVectorPxPyPzMVector(trk2.px(), trk2.py(), trk2.pz(), massPi);
280305
}
281306
reco = part1 + part2;
282307

283308
if (reco.Rapidity() > cfgMaxRap || reco.Rapidity() < cfgMinRap)
284309
continue;
285310

286311
if (trk1.sign() * trk2.sign() < 0) {
287-
histos.fill(HIST("hInvMass_Kstar_US"), reco.M(), reco.Pt(), collision.cent());
312+
histos.fill(HIST("hInvMass_Kstar_US"), reco.M(), reco.Pt(), multiplicity);
288313

289314
if constexpr (IsMC) {
290315
if (trk1.motherId() != trk2.motherId())
291316
continue;
292-
if ((std::abs(trk1.pdgCode()) == 211 && std::abs(trk2.pdgCode()) == 321) || (std::abs(trk1.pdgCode()) == 321 && std::abs(trk2.pdgCode()) == 211)) {
293-
if (std::abs(trk1.motherPDG()) == 313) {
294-
histos.fill(HIST("MCL/hpT_Kstar_Kpi_REC"), reco.M(), reco.Pt(), collision.cent());
317+
if ((std::abs(trk1.pdgCode()) == kPiPlus && std::abs(trk2.pdgCode()) == kKPlus) || (std::abs(trk1.pdgCode()) == kKPlus && std::abs(trk2.pdgCode()) == kPiPlus)) {
318+
if (std::abs(trk1.motherPDG()) == kK0Star892) {
319+
histos.fill(HIST("MCL/hpT_Kstar_Kpi_REC"), reco.M(), reco.Pt(), multiplicity);
295320
}
296321
}
297322
}
298323
} else if (trk1.sign() > 0 && trk2.sign() > 0) {
299-
histos.fill(HIST("hInvMass_Kstar_LSpp"), reco.M(), reco.Pt(), collision.cent());
324+
histos.fill(HIST("hInvMass_Kstar_LSpp"), reco.M(), reco.Pt(), multiplicity);
300325
} else if (trk1.sign() < 0 && trk2.sign() < 0) {
301-
histos.fill(HIST("hInvMass_Kstar_LSmm"), reco.M(), reco.Pt(), collision.cent());
326+
histos.fill(HIST("hInvMass_Kstar_LSmm"), reco.M(), reco.Pt(), multiplicity);
302327
}
303328
}
304329
}
305330
}
306331

307332
void processData(aod::ResoCollision const& collision, aod::ResoTracks const& resotracks)
308333
{
309-
fillHistograms<false>(collision, resotracks);
334+
auto multiplicity = collision.cent();
335+
fillHistograms<false>(collision, multiplicity, resotracks);
310336
}
311337
PROCESS_SWITCH(rho770analysis, processData, "Process Event for data", true);
312338

313-
void processMCLight(aod::ResoCollision const& collision, soa::Join<aod::ResoTracks, aod::ResoMCTracks> const& resotracks)
339+
void processMCLight(ResoMCCols::iterator const& collision,
340+
aod::ResoCollisionColls const& collisionIndex,
341+
soa::Join<aod::ResoCollisionCandidatesMC, aod::PVMults> const& collisionsMC,
342+
soa::Join<aod::ResoTracks, aod::ResoMCTracks> const& resotracks,
343+
soa::Join<aod::McCollisions, aod::McCentFT0Ms> const&)
314344
{
315-
fillHistograms<true>(collision, resotracks);
345+
float multiplicity;
346+
auto linkRow = collisionIndex.iteratorAt(collision.globalIndex());
347+
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision
348+
349+
auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision
350+
351+
if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision
352+
return;
353+
354+
auto mcColl = coll.mcCollision_as<soa::Join<aod::McCollisions, aod::McCentFT0Ms>>();
355+
multiplicity = mcColl.centFT0M();
356+
357+
if (!collision.isInAfterAllCuts())
358+
return;
359+
360+
fillHistograms<true>(collision, multiplicity, resotracks);
316361
}
317362
PROCESS_SWITCH(rho770analysis, processMCLight, "Process Event for MC", false);
318363

319-
void processMCTrue(ResoMCCols::iterator const& collision, aod::ResoMCParents const& resoParents)
364+
void processMCTrue(ResoMCCols::iterator const& collision,
365+
aod::ResoCollisionColls const& collisionIndex,
366+
aod::ResoMCParents const& resoParents,
367+
aod::ResoCollisionCandidatesMC const& collisionsMC,
368+
soa::Join<aod::McCollisions, aod::McCentFT0Ms> const&)
320369
{
321-
auto multiplicity = collision.cent();
370+
float multiplicity;
371+
auto linkRow = collisionIndex.iteratorAt(collision.globalIndex());
372+
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision
373+
374+
auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision
375+
376+
if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision
377+
return;
378+
379+
auto mcColl = coll.mcCollision_as<soa::Join<aod::McCollisions, aod::McCentFT0Ms>>();
380+
multiplicity = mcColl.centFT0M();
322381

323382
for (const auto& part : resoParents) { // loop over all pre-filtered MC particles
324-
if (std::abs(part.pdgCode()) != 113)
383+
if (std::abs(part.pdgCode()) != kRho770_0)
325384
continue;
326385
if (!part.producedByGenerator())
327386
continue;
328387
if (part.y() < cfgMinRap || part.y() > cfgMaxRap)
329388
continue;
330-
if (!(std::abs(part.daughterPDG1()) == 211 && std::abs(part.daughterPDG2()) == 211))
389+
if (!(std::abs(part.daughterPDG1()) == kPiPlus && std::abs(part.daughterPDG2()) == kPiPlus))
331390
continue;
332391

333-
TLorentzVector truthpar;
334-
truthpar.SetPxPyPzE(part.px(), part.py(), part.pz(), part.e());
392+
LorentzVectorPxPyPzEVector truthpar;
393+
truthpar = LorentzVectorPxPyPzEVector(part.px(), part.py(), part.pz(), part.e());
335394
auto mass = truthpar.M();
336395

337396
histos.fill(HIST("MCL/hpT_rho770_GEN"), 0, mass, part.pt(), multiplicity);

0 commit comments

Comments
 (0)