Skip to content

Commit 5a4ad2d

Browse files
limhyunji999limhyunji999
committed
modify INEL>0 part
1 parent 0ff57fc commit 5a4ad2d

File tree

1 file changed

+94
-49
lines changed

1 file changed

+94
-49
lines changed

‎PWGLF/Tasks/Resonances/rho770analysis.cxx‎

Lines changed: 94 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -12,28 +12,35 @@
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 12/01/2026
16+
17+
#include <Framework/Configurable.h>
18+
#include "Math/Vector4D.h"
19+
20+
#include "DataFormatsParameters/GRPObject.h"
1621

1722
#include "PWGLF/DataModel/LFResonanceTables.h"
23+
#include "PWGLF/DataModel/mcCentrality.h"
24+
#include "PWGLF/Utils/inelGt.h"
1825

26+
#include "Common/Core/RecoDecay.h"
1927
#include "Common/DataModel/Centrality.h"
2028
#include "Common/DataModel/EventSelection.h"
29+
#include "Common/DataModel/PIDResponse.h"
2130

2231
#include "CommonConstants/PhysicsConstants.h"
23-
#include "DataFormatsParameters/GRPObject.h"
2432
#include "Framework/ASoAHelpers.h"
2533
#include "Framework/AnalysisTask.h"
34+
#include "Framework/O2DatabasePDGPlugin.h"
2635
#include "Framework/runDataProcessing.h"
27-
#include <Framework/Configurable.h>
28-
29-
#include "TVector2.h"
30-
#include <TLorentzVector.h>
3136

3237
using namespace o2;
3338
using namespace o2::framework;
3439
using namespace o2::framework::expressions;
3540
using namespace o2::soa;
3641
using namespace o2::constants::physics;
42+
using LorentzVectorPxPyPzMVector = ROOT::Math::PxPyPzMVector;
43+
using LorentzVectorPxPyPzEVector = ROOT::Math::PxPyPzEVector;
3744

3845
struct rho770analysis {
3946
SliceCache cache;
@@ -62,17 +69,18 @@ struct rho770analysis {
6269
// kEtaRange)
6370
Configurable<bool> cfgPVContributor{"cfgPVContributor", true, "PV contributor track selection"}; // PV Contriuibutor
6471
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"};
72+
Configurable<int> cfgTPCcluster{"cfgTPCcluster", 1, "Number of TPC cluster"};
73+
Configurable<bool> cfgUseTPCRefit{"cfgUseTPCRefit", false, "Require TPC Refit"}; //refit is included in global track selection
6774
Configurable<bool> cfgUseITSRefit{"cfgUseITSRefit", false, "Require ITS Refit"};
6875
Configurable<bool> cfgHasTOF{"cfgHasTOF", false, "Require TOF"};
76+
Configurable<int> cfgTPCRows{"cfgTPCRows", 80, "Minimum Number of TPC Crossed Rows "};
6977

7078
// PID
71-
Configurable<double> cMaxTOFnSigmaPion{"cMaxTOFnSigmaPion", 3.0, "TOF nSigma cut for Pion"}; // TOF
72-
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
79+
Configurable<double> cMaxTOFnSigmaPion{"cMaxTOFnSigmaPion", 3.0, "TOF nSigma cut for Pion"}; // TOF
80+
Configurable<double> cMaxTPCnSigmaPion{"cMaxTPCnSigmaPion", 5.0, "TPC nSigma cut for Pion"}; // TPC
81+
Configurable<double> cMaxTPCnSigmaPionnoTOF{"cMaxTPCnSigmaPionnoTOF", 3.0, "TPC nSigma cut for Pion in no TOF case"}; // TPC
7482
Configurable<double> nsigmaCutCombinedPion{"nsigmaCutCombinedPion", 3.0, "Combined nSigma cut for Pion"};
75-
Configurable<int> selectType{"selectType", 0, "PID selection type"};
83+
Configurable<int> selectType{"selectType", 1, "PID selection type"};
7684

7785
// Axis
7886
ConfigurableAxis massAxis{"massAxis", {400, 0.2, 2.2}, "Invariant mass axis"};
@@ -128,7 +136,7 @@ struct rho770analysis {
128136
{
129137
if (std::abs(track.pt()) < cfgMinPt)
130138
return false;
131-
if (std::fabs(track.eta()) > cfgMaxEta)
139+
if (std::abs(track.eta()) > cfgMaxEta)
132140
return false;
133141
if (std::abs(track.dcaXY()) > cfgMaxDCArToPVcut)
134142
return false;
@@ -150,26 +158,28 @@ struct rho770analysis {
150158
return false;
151159
if (cfgGlobalTrack && !track.isGlobalTrack())
152160
return false;
161+
if (track.tpcNClsCrossedRows() < cfgTPCRows)
162+
return false;
153163

154164
return true;
155165
}
156166

157167
template <typename TrackType>
158168
bool selPion(const TrackType track)
159169
{
160-
if (selectType == 0) {
170+
if (selectType == 0) { // TPC or TOF
161171
if (std::fabs(track.tpcNSigmaPi()) >= cMaxTPCnSigmaPion || std::fabs(track.tofNSigmaPi()) >= cMaxTOFnSigmaPion)
162172
return false;
163173
}
164-
if (selectType == 1) {
174+
if (selectType == 1) { // only TPC
165175
if (std::fabs(track.tpcNSigmaPi()) >= cMaxTPCnSigmaPion)
166176
return false;
167177
}
168-
if (selectType == 2) {
178+
if (selectType == 2) { // combining
169179
if (track.tpcNSigmaPi() * track.tpcNSigmaPi() + track.tofNSigmaPi() * track.tofNSigmaPi() >= nsigmaCutCombinedPion * nsigmaCutCombinedPion)
170180
return false;
171181
}
172-
if (selectType == 3) {
182+
if (selectType == 3) { //TPC TOF veto
173183
if (track.hasTOF()) {
174184
if (std::fabs(track.tpcNSigmaPi()) >= cMaxTPCnSigmaPion || std::fabs(track.tofNSigmaPi()) >= cMaxTOFnSigmaPion)
175185
return false;
@@ -208,10 +218,10 @@ struct rho770analysis {
208218
return true;
209219
}
210220

211-
template <bool IsMC, typename CollisionType, typename TracksType>
212-
void fillHistograms(const CollisionType& collision, const TracksType& dTracks)
221+
template <bool IsMC, typename CollisionType, typename CenMult, typename TracksType>
222+
void fillHistograms(const CollisionType& collision, const CenMult& multiplicity, const TracksType& dTracks)
213223
{
214-
TLorentzVector part1, part2, reco;
224+
LorentzVectorPxPyPzMVector part1, part2, reco;
215225
for (const auto& [trk1, trk2] : combinations(CombinationsUpperIndexPolicy(dTracks, dTracks))) {
216226

217227
if (trk1.index() == trk2.index()) {
@@ -232,93 +242,128 @@ struct rho770analysis {
232242
continue;
233243

234244
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);
245+
246+
part1 = LorentzVectorPxPyPzMVector(trk1.px(), trk1.py(), trk1.pz(), massPi);
247+
part2 = LorentzVectorPxPyPzMVector(trk2.px(), trk2.py(), trk2.pz(), massPi);
237248
reco = part1 + part2;
238249

239250
if (reco.Rapidity() > cfgMaxRap || reco.Rapidity() < cfgMinRap)
240251
continue;
241252

242253
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());
254+
histos.fill(HIST("hInvMass_rho770_US"), reco.M(), reco.Pt(), multiplicity);
255+
histos.fill(HIST("hInvMass_K0s_US"), reco.M(), reco.Pt(), multiplicity);
245256

246257
if constexpr (IsMC) {
247258
if (trk1.motherId() != trk2.motherId())
248259
continue;
249260
if (std::abs(trk1.pdgCode()) == 211 && std::abs(trk2.pdgCode()) == 211) {
250261
if (std::abs(trk1.motherPDG()) == 113) {
251-
histos.fill(HIST("MCL/hpT_rho770_REC"), reco.M(), reco.Pt(), collision.cent());
262+
histos.fill(HIST("MCL/hpT_rho770_REC"), reco.M(), reco.Pt(), multiplicity);
252263
} else if (std::abs(trk1.motherPDG()) == 223) {
253-
histos.fill(HIST("MCL/hpT_omega_REC"), reco.M(), reco.Pt(), collision.cent());
264+
histos.fill(HIST("MCL/hpT_omega_REC"), reco.M(), reco.Pt(), multiplicity);
254265
} 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());
266+
histos.fill(HIST("MCL/hpT_K0s_REC"), reco.M(), reco.Pt(), multiplicity);
267+
histos.fill(HIST("MCL/hpT_K0s_pipi_REC"), reco.M(), reco.Pt(), multiplicity);
257268
}
258269
} else if ((std::abs(trk1.pdgCode()) == 211 && std::abs(trk2.pdgCode()) == 321) || (std::abs(trk1.pdgCode()) == 321 && std::abs(trk2.pdgCode()) == 211)) {
259270
if (std::abs(trk1.motherPDG()) == 313) {
260-
histos.fill(HIST("MCL/hpT_Kstar_REC"), reco.M(), reco.Pt(), collision.cent());
271+
histos.fill(HIST("MCL/hpT_Kstar_REC"), reco.M(), reco.Pt(), multiplicity);
261272
}
262273
}
263274
}
264275
} 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());
276+
histos.fill(HIST("hInvMass_rho770_LSpp"), reco.M(), reco.Pt(), multiplicity);
277+
histos.fill(HIST("hInvMass_K0s_LSpp"), reco.M(), reco.Pt(), multiplicity);
267278
} 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());
279+
histos.fill(HIST("hInvMass_rho770_LSmm"), reco.M(), reco.Pt(), multiplicity);
280+
histos.fill(HIST("hInvMass_K0s_LSmm"), reco.M(), reco.Pt(), multiplicity);
270281
}
271282
}
272283

273284
if ((selPion(trk1) && selKaon(trk2)) || (selKaon(trk1) && selPion(trk2))) {
274285
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);
286+
part1 = LorentzVectorPxPyPzMVector(trk1.px(), trk1.py(), trk1.pz(), massPi);
287+
part2 = LorentzVectorPxPyPzMVector(trk2.px(), trk2.py(), trk2.pz(), massKa);
277288
} 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);
289+
part1 = LorentzVectorPxPyPzMVector(trk1.px(), trk1.py(), trk1.pz(), massKa);
290+
part2 = LorentzVectorPxPyPzMVector(trk2.px(), trk2.py(), trk2.pz(), massPi);
280291
}
281292
reco = part1 + part2;
282293

283294
if (reco.Rapidity() > cfgMaxRap || reco.Rapidity() < cfgMinRap)
284295
continue;
285296

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

289300
if constexpr (IsMC) {
290301
if (trk1.motherId() != trk2.motherId())
291302
continue;
292303
if ((std::abs(trk1.pdgCode()) == 211 && std::abs(trk2.pdgCode()) == 321) || (std::abs(trk1.pdgCode()) == 321 && std::abs(trk2.pdgCode()) == 211)) {
293304
if (std::abs(trk1.motherPDG()) == 313) {
294-
histos.fill(HIST("MCL/hpT_Kstar_Kpi_REC"), reco.M(), reco.Pt(), collision.cent());
305+
histos.fill(HIST("MCL/hpT_Kstar_Kpi_REC"), reco.M(), reco.Pt(), multiplicity);
295306
}
296307
}
297308
}
298309
} else if (trk1.sign() > 0 && trk2.sign() > 0) {
299-
histos.fill(HIST("hInvMass_Kstar_LSpp"), reco.M(), reco.Pt(), collision.cent());
310+
histos.fill(HIST("hInvMass_Kstar_LSpp"), reco.M(), reco.Pt(), multiplicity);
300311
} else if (trk1.sign() < 0 && trk2.sign() < 0) {
301-
histos.fill(HIST("hInvMass_Kstar_LSmm"), reco.M(), reco.Pt(), collision.cent());
312+
histos.fill(HIST("hInvMass_Kstar_LSmm"), reco.M(), reco.Pt(), multiplicity);
302313
}
303314
}
304315
}
305316
}
306317

307318
void processData(aod::ResoCollision const& collision, aod::ResoTracks const& resotracks)
308319
{
309-
fillHistograms<false>(collision, resotracks);
320+
auto multiplicity = collision.cent();
321+
fillHistograms<false>(collision, multiplicity, resotracks);
310322
}
311323
PROCESS_SWITCH(rho770analysis, processData, "Process Event for data", true);
312324

313-
void processMCLight(aod::ResoCollision const& collision, soa::Join<aod::ResoTracks, aod::ResoMCTracks> const& resotracks)
325+
void processMCLight(ResoMCCols::iterator const& collision,
326+
aod::ResoCollisionColls const& collisionIndex,
327+
soa::Join<aod::ResoCollisionCandidatesMC, aod::PVMults> const& collisionsMC,
328+
soa::Join<aod::ResoTracks, aod::ResoMCTracks> const& resotracks,
329+
soa::Join<aod::McCollisions, aod::McCentFT0Ms> const&)
314330
{
315-
fillHistograms<true>(collision, resotracks);
331+
float multiplicity;
332+
auto linkRow = collisionIndex.iteratorAt(collision.globalIndex());
333+
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision
334+
335+
auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision
336+
337+
if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision
338+
return;
339+
340+
auto mcColl = coll.mcCollision_as<soa::Join<aod::McCollisions, aod::McCentFT0Ms>>();
341+
multiplicity = mcColl.centFT0M();
342+
343+
if (!collision.isInAfterAllCuts())
344+
return;
345+
346+
fillHistograms<true>(collision, multiplicity, resotracks);
316347
}
317348
PROCESS_SWITCH(rho770analysis, processMCLight, "Process Event for MC", false);
318349

319-
void processMCTrue(ResoMCCols::iterator const& collision, aod::ResoMCParents const& resoParents)
350+
void processMCTrue(ResoMCCols::iterator const& collision,
351+
aod::ResoCollisionColls const& collisionIndex,
352+
aod::ResoMCParents const& resoParents,
353+
aod::ResoCollisionCandidatesMC const& collisionsMC,
354+
soa::Join<aod::McCollisions, aod::McCentFT0Ms> const&)
320355
{
321-
auto multiplicity = collision.cent();
356+
float multiplicity;
357+
auto linkRow = collisionIndex.iteratorAt(collision.globalIndex());
358+
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision
359+
360+
auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision
361+
362+
if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision
363+
return;
364+
365+
auto mcColl = coll.mcCollision_as<soa::Join<aod::McCollisions, aod::McCentFT0Ms>>();
366+
multiplicity = mcColl.centFT0M();
322367

323368
for (const auto& part : resoParents) { // loop over all pre-filtered MC particles
324369
if (std::abs(part.pdgCode()) != 113)
@@ -330,12 +375,12 @@ struct rho770analysis {
330375
if (!(std::abs(part.daughterPDG1()) == 211 && std::abs(part.daughterPDG2()) == 211))
331376
continue;
332377

333-
TLorentzVector truthpar;
334-
truthpar.SetPxPyPzE(part.px(), part.py(), part.pz(), part.e());
378+
LorentzVectorPxPyPzEVector truthpar;
379+
truthpar = LorentzVectorPxPyPzEVector(part.px(), part.py(), part.pz(), part.e());
335380
auto mass = truthpar.M();
336381

337382
histos.fill(HIST("MCL/hpT_rho770_GEN"), 0, mass, part.pt(), multiplicity);
338-
383+
339384
if (collision.isVtxIn10()) {
340385
histos.fill(HIST("MCL/hpT_rho770_GEN"), 1, mass, part.pt(), multiplicity);
341386
}

0 commit comments

Comments
 (0)