Skip to content

Commit 3e4433b

Browse files
omvazqueomvazque
authored
[PWGLF] Select MC charged particles (#11850)
1 parent 2d84ca1 commit 3e4433b

File tree

1 file changed

+91
-22
lines changed

1 file changed

+91
-22
lines changed

‎PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx‎

Lines changed: 91 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -69,6 +69,7 @@ struct UccZdc {
6969

7070
static constexpr float kCollEnergy{2.68};
7171
static constexpr float kZero{0.};
72+
static constexpr float kMinCharge{3.f};
7273

7374
// Configurables Event Selection
7475
Configurable<bool> isNoCollInTimeRangeStrict{"isNoCollInTimeRangeStrict", true, "use isNoCollInTimeRangeStrict?"};
@@ -710,6 +711,7 @@ struct UccZdc {
710711

711712
std::vector<float> pTs;
712713
std::vector<float> vecFD;
714+
std::vector<float> vecEff;
713715
std::vector<float> vecOneOverEff;
714716

715717
// Calculates the Nch multiplicity
@@ -747,7 +749,7 @@ struct UccZdc {
747749
// Fill vectors for [pT] measurement
748750
pTs.clear();
749751
vecFD.clear();
750-
vecOneOverEff.clear();
752+
vecEff.clear();
751753
for (const auto& track : tracks) {
752754
// Track Selection
753755
if (!track.isGlobalTrack()) {
@@ -771,7 +773,7 @@ struct UccZdc {
771773
}
772774
if ((effValue > 0.) && (fdValue > 0.)) {
773775
pTs.emplace_back(pt);
774-
vecOneOverEff.emplace_back(1. / effValue);
776+
vecEff.emplace_back(effValue);
775777
vecFD.emplace_back(fdValue);
776778
}
777779
// To calculate event-averaged <pt>
@@ -780,7 +782,7 @@ struct UccZdc {
780782

781783
double p1, p2, p3, p4, w1, w2, w3, w4;
782784
p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0;
783-
getPTpowers(pTs, vecOneOverEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4);
785+
getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4);
784786

785787
// EbE one-particle pT correlation
786788
double oneParCorr{p1 / w1};
@@ -819,6 +821,7 @@ struct UccZdc {
819821

820822
// Preslice<aod::McParticles> perMCCollision = aod::mcparticle::mcCollisionId;
821823
Preslice<TheFilteredSimTracks> perCollision = aod::track::collisionId;
824+
Service<o2::framework::O2DatabasePDG> pdg;
822825
TRandom* randPointer = new TRandom();
823826
void processMCclosure(aod::McCollisions::iterator const& mccollision, soa::SmallGroups<o2::aod::SimCollisions> const& collisions, o2::aod::BCsRun3 const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::McParticles const& mcParticles, TheFilteredSimTracks const& simTracks)
824827
{
@@ -879,9 +882,10 @@ struct UccZdc {
879882
return;
880883
}
881884

882-
std::vector<float> pTs;
883-
std::vector<float> vecFD;
884-
std::vector<float> vecOneOverEff;
885+
std::vector<double> pTs;
886+
std::vector<double> vecFD;
887+
std::vector<double> vecEff;
888+
std::vector<double> vecOneOverEffXFD;
885889
// std::vector<float> wIs;
886890
const auto& groupedTracks{simTracks.sliceBy(perCollision, collision.globalIndex())};
887891

@@ -901,6 +905,8 @@ struct UccZdc {
901905
}
902906

903907
// Calculates the event weight, W_k
908+
const int foundNchBin{efficiency->GetXaxis()->FindBin(nchRaw)};
909+
904910
for (const auto& track : groupedTracks) {
905911
// Track Selection
906912
if (track.eta() < minEta || track.eta() > maxEta) {
@@ -912,32 +918,51 @@ struct UccZdc {
912918
if (!track.isGlobalTrack()) {
913919
continue;
914920
}
921+
if (!track.has_mcParticle()) {
922+
continue;
923+
}
924+
const auto& particle{track.mcParticle()};
915925

916-
float pt{track.pt()};
917-
int foundNchBin{efficiency->GetXaxis()->FindBin(nchRaw)};
918-
int foundPtBin{efficiency->GetYaxis()->FindBin(pt)};
919-
float effValue{1.};
920-
float fdValue{1.};
926+
auto charge{0.};
927+
// Get the MC particle
928+
auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
929+
if (pdgParticle != nullptr) {
930+
charge = pdgParticle->Charge();
931+
} else {
932+
continue;
933+
}
934+
935+
// Is it a charged particle?
936+
if (std::abs(charge) < kMinCharge) {
937+
continue;
938+
}
939+
// Is it a primary particle?
940+
// if (!particle.isPhysicalPrimary()) { continue; }
941+
942+
const double pt{static_cast<double>(track.pt())};
943+
const int foundPtBin{efficiency->GetYaxis()->FindBin(pt)};
944+
double effValue{1.};
945+
double fdValue{1.};
921946

922947
if (applyEff) {
923948
effValue = efficiency->GetBinContent(foundNchBin, foundPtBin);
924949
fdValue = fd->GetBinContent(fd->FindBin(pt));
925950
}
926951
if ((effValue > 0.) && (fdValue > 0.)) {
927952
pTs.emplace_back(pt);
928-
vecOneOverEff.emplace_back(1. / effValue);
953+
vecEff.emplace_back(effValue);
929954
vecFD.emplace_back(fdValue);
955+
vecOneOverEffXFD.emplace_back(fdValue / effValue);
930956
}
931957
}
932-
933-
nchMult = std::accumulate(vecOneOverEff.begin(), vecOneOverEff.end(), 0);
958+
nchMult = std::accumulate(vecOneOverEffXFD.begin(), vecOneOverEffXFD.end(), 0);
934959
if (nchMult < minNchSel) {
935960
return;
936961
}
937962

938963
double p1, p2, p3, p4, w1, w2, w3, w4;
939964
p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0;
940-
getPTpowers(pTs, vecOneOverEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4);
965+
getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4);
941966

942967
const double denTwoParCorr{std::pow(w1, 2.) - w2};
943968
const double numTwoParCorr{std::pow(p1, 2.) - p2};
@@ -971,6 +996,21 @@ struct UccZdc {
971996
if (particle.pt() < minPt || particle.pt() > maxPt) {
972997
continue;
973998
}
999+
1000+
auto charge{0.};
1001+
// Get the MC particle
1002+
auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
1003+
if (pdgParticle != nullptr) {
1004+
charge = pdgParticle->Charge();
1005+
} else {
1006+
continue;
1007+
}
1008+
1009+
// Is it a charged particle?
1010+
if (std::abs(charge) < kMinCharge) {
1011+
continue;
1012+
}
1013+
// Is it a primary particle?
9741014
if (!particle.isPhysicalPrimary()) {
9751015
continue;
9761016
}
@@ -1044,9 +1084,23 @@ struct UccZdc {
10441084
if (!track.has_mcParticle()) {
10451085
continue;
10461086
}
1047-
registry.fill(HIST("Pt_all_ch"), nchRaw, track.pt());
1048-
1087+
// Get the MC particle
10491088
const auto& particle{track.mcParticle()};
1089+
auto charge{0.};
1090+
auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
1091+
if (pdgParticle != nullptr) {
1092+
charge = pdgParticle->Charge();
1093+
} else {
1094+
continue;
1095+
}
1096+
1097+
// Is it a charged particle?
1098+
if (std::abs(charge) < kMinCharge) {
1099+
continue;
1100+
}
1101+
// All charged particles
1102+
registry.fill(HIST("Pt_all_ch"), nchRaw, track.pt());
1103+
// Is it a primary particle?
10501104
if (!particle.isPhysicalPrimary()) {
10511105
continue;
10521106
}
@@ -1075,6 +1129,21 @@ struct UccZdc {
10751129
if (particle.pt() < minPt || particle.pt() > maxPt) {
10761130
continue;
10771131
}
1132+
1133+
auto charge{0.};
1134+
// Get the MC particle
1135+
auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
1136+
if (pdgParticle != nullptr) {
1137+
charge = pdgParticle->Charge();
1138+
} else {
1139+
continue;
1140+
}
1141+
1142+
// Is it a charged particle?
1143+
if (std::abs(charge) < kMinCharge) {
1144+
continue;
1145+
}
1146+
// Is it a primary particle?
10781147
if (!particle.isPhysicalPrimary()) {
10791148
continue;
10801149
}
@@ -1101,14 +1170,14 @@ struct UccZdc {
11011170
PROCESS_SWITCH(UccZdc, processMCclosure, "Process MC closure", false);
11021171

11031172
template <typename T, typename U>
1104-
void getPTpowers(const T& pTs, const T& vecOneOverEff, const T& vecFD, U& pOne, U& wOne, U& pTwo, U& wTwo, U& pThree, U& wThree, U& pFour, U& wFour)
1173+
void getPTpowers(const T& pTs, const T& vecEff, const T& vecFD, U& pOne, U& wOne, U& pTwo, U& wTwo, U& pThree, U& wThree, U& pFour, U& wFour)
11051174
{
11061175
pOne = wOne = pTwo = wTwo = pThree = wThree = pFour = wFour = 0.;
11071176
for (std::size_t i = 0; i < pTs.size(); ++i) {
1108-
const float pTi{pTs.at(i)};
1109-
const float eFFi{vecOneOverEff.at(i)};
1110-
const float fDi{vecFD.at(i)};
1111-
const float wEighti{eFFi * fDi};
1177+
const double pTi{pTs.at(i)};
1178+
const double eFFi{vecEff.at(i)};
1179+
const double fDi{vecFD.at(i)};
1180+
const double wEighti{std::pow(eFFi, -1.) * fDi};
11121181
pOne += wEighti * pTi;
11131182
wOne += wEighti;
11141183
pTwo += std::pow(wEighti * pTi, 2.);

0 commit comments

Comments
 (0)