Skip to content

Commit 607db4d

Browse files
SwatiSaha-1997SwatiSaha-1997
authored
[PWGCF] ITS pid option added and previous PR comments addressed (#9858)
1 parent 9ad3635 commit 607db4d

File tree

1 file changed

+107
-15
lines changed

1 file changed

+107
-15
lines changed

‎PWGCF/EbyEFluctuations/Tasks/netprotonCumulantsMc.cxx‎

Lines changed: 107 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@
3232
#include "Common/DataModel/TrackSelectionTables.h"
3333
#include "Common/DataModel/Centrality.h"
3434
#include "Common/DataModel/PIDResponse.h"
35+
#include "Common/DataModel/PIDResponseITS.h"
3536
#include "Common/DataModel/Multiplicity.h"
3637
#include "Common/Core/trackUtilities.h"
3738
#include "CommonConstants/PhysicsConstants.h"
@@ -133,8 +134,10 @@ struct NetprotonCumulantsMc {
133134
Configurable<float> cfgCutItsChi2NCl{"cfgCutItsChi2NCl", 36.0f, "Maximum ITSchi2NCl"};
134135
Configurable<float> cfgCutDCAxy{"cfgCutDCAxy", 2.0f, "DCAxy range for tracks"};
135136
Configurable<float> cfgCutDCAz{"cfgCutDCAz", 2.0f, "DCAz range for tracks"};
136-
Configurable<int> cfgITScluster{"cfgITScluster", 0, "Number of ITS cluster"};
137-
Configurable<int> cfgTPCcluster{"cfgTPCcluster", 70, "Number of TPC cluster"};
137+
Configurable<int> cfgITScluster{"cfgITScluster", 1, "Minimum Number of ITS cluster"};
138+
Configurable<int> cfgTPCcluster{"cfgTPCcluster", 80, "Minimum Number of TPC cluster"};
139+
Configurable<int> cfgTPCnCrossedRows{"cfgTPCnCrossedRows", 70, "Minimum Number of TPC crossed-rows"};
140+
Configurable<bool> cfgUseItsPid{"cfgUseItsPid", true, "Use ITS nSigma Cut"};
138141

139142
// Calculation of cumulants central/error
140143
Configurable<int> cfgNSubsample{"cfgNSubsample", 10, "Number of subsamples for ERR"};
@@ -149,6 +152,8 @@ struct NetprotonCumulantsMc {
149152
Configurable<bool> cfgLoadEff{"cfgLoadEff", true, "Load efficiency from file"};
150153
Configurable<bool> cfgEvSelkNoSameBunchPileup{"cfgEvSelkNoSameBunchPileup", true, "Pileup removal"};
151154

155+
ConfigurableAxis cfgCentralityBins{"cfgCentralityBins", {90, 0., 90.}, "Centrality/Multiplicity percentile bining"};
156+
152157
// Connect to ccdb
153158
Service<ccdb::BasicCCDBManager> ccdb;
154159
Configurable<int64_t> ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"};
@@ -207,11 +212,13 @@ struct NetprotonCumulantsMc {
207212
AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/#it{c})"};
208213
std::vector<double> etaBinning = {-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};
209214
AxisSpec etaAxis = {etaBinning, "#it{#eta}"};
210-
std::vector<double> centBining = {0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90};
211-
AxisSpec centAxis = {centBining, "Multiplicity percentile from FT0M (%)"};
215+
// std::vector<double> centBining = {0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90};
216+
// AxisSpec centAxis = {centBining, "Multiplicity percentile from FT0M (%)"};
217+
const AxisSpec centAxis{cfgCentralityBins, "Multiplicity percentile from FT0M (%)"};
212218
AxisSpec netprotonAxis = {41, -20.5, 20.5, "net-proton number"};
213219
AxisSpec protonAxis = {21, -0.5, 20.5, "proton number"};
214220
AxisSpec antiprotonAxis = {21, -0.5, 20.5, "antiproton number"};
221+
AxisSpec nSigmaAxis = {200, -5.0, 5.0, "nSigma(Proton)"};
215222

216223
auto noSubsample = static_cast<int>(cfgNSubsample);
217224
float maxSubsample = 1.0 * noSubsample;
@@ -252,6 +259,11 @@ struct NetprotonCumulantsMc {
252259
histos.add("hgen2DEtaVsPtProton", "2D hist of Generated Proton y: eta vs. x: pT", kTH2F, {ptAxis, etaAxis});
253260
histos.add("hgen2DEtaVsPtAntiproton", "2D hist of Generated Anti-proton y: eta vs. x: pT", kTH2F, {ptAxis, etaAxis});
254261

262+
// 2D histograms of nSigma
263+
histos.add("h2DnsigmaTpcVsPt", "2D hist of nSigmaTPC vs. pT", kTH2F, {ptAxis, nSigmaAxis});
264+
histos.add("h2DnsigmaTofVsPt", "2D hist of nSigmaTOF vs. pT", kTH2F, {ptAxis, nSigmaAxis});
265+
histos.add("h2DnsigmaItsVsPt", "2D hist of nSigmaITS vs. pT", kTH2F, {ptAxis, nSigmaAxis});
266+
255267
if (cfgIsCalculateCentral) {
256268
// uncorrected
257269
histos.add("Prof_mu1_netproton", "", {HistType::kTProfile, {centAxis}});
@@ -823,6 +835,48 @@ struct NetprotonCumulantsMc {
823835
template <typename T>
824836
bool selectionPIDold(const T& candidate)
825837
{
838+
if (!candidate.hasTPC())
839+
return false;
840+
841+
//! PID checking as done in Run2 my analysis
842+
//! ----------------------------------------------------------------------
843+
int flag = 0; //! pid check main flag
844+
845+
if (candidate.pt() > 0.2f && candidate.pt() <= cfgCutPtUpperTPC) {
846+
if (std::abs(candidate.tpcNSigmaPr()) < cfgnSigmaCutTPC) {
847+
flag = 1;
848+
}
849+
}
850+
if (candidate.hasTOF() && candidate.pt() > cfgCutPtUpperTPC && candidate.pt() < 5.0f) {
851+
const float combNSigmaPr = std::sqrt(std::pow(candidate.tpcNSigmaPr(), 2.0) + std::pow(candidate.tofNSigmaPr(), 2.0));
852+
const float combNSigmaPi = std::sqrt(std::pow(candidate.tpcNSigmaPi(), 2.0) + std::pow(candidate.tofNSigmaPi(), 2.0));
853+
const float combNSigmaKa = std::sqrt(std::pow(candidate.tpcNSigmaKa(), 2.0) + std::pow(candidate.tofNSigmaKa(), 2.0));
854+
855+
int flag2 = 0;
856+
if (combNSigmaPr < 3.0)
857+
flag2 += 1;
858+
if (combNSigmaPi < 3.0)
859+
flag2 += 1;
860+
if (combNSigmaKa < 3.0)
861+
flag2 += 1;
862+
if (!(flag2 > 1) && !(combNSigmaPr > combNSigmaPi) && !(combNSigmaPr > combNSigmaKa)) {
863+
if (combNSigmaPr < cfgnSigmaCutCombTPCTOF) {
864+
flag = 1;
865+
}
866+
}
867+
}
868+
if (flag == 1)
869+
return true;
870+
else
871+
return false;
872+
}
873+
874+
template <typename T>
875+
bool selectionPIDoldTOFveto(const T& candidate)
876+
{
877+
if (!candidate.hasTPC())
878+
return false;
879+
826880
//! PID checking as done in Run2 my analysis
827881
//! ----------------------------------------------------------------------
828882
int flag = 0; //! pid check main flag
@@ -860,7 +914,7 @@ struct NetprotonCumulantsMc {
860914
}
861915

862916
template <typename T>
863-
bool selectionPIDnew(const T& candidate)
917+
bool selectionPIDnew(const T& candidate) // Victor's BF analysis
864918
{
865919
// electron rejection
866920
if (candidate.tpcNSigmaEl() > -3.0f && candidate.tpcNSigmaEl() < 5.0f && std::abs(candidate.tpcNSigmaPi()) > 3.0f && std::abs(candidate.tpcNSigmaKa()) > 3.0f && std::abs(candidate.tpcNSigmaPr()) > 3.0f) {
@@ -906,12 +960,12 @@ struct NetprotonCumulantsMc {
906960
// Find the pt bin index based on the track's pt value
907961
int binIndex = -1;
908962
// Get the array from the Configurable
909-
auto ptBins = (std::vector<float>)cfgPtBins;
910-
auto effProt = (std::vector<float>)cfgProtonEff;
911-
auto effAntiprot = (std::vector<float>)cfgAntiprotonEff;
963+
// auto ptBins = (std::vector<float>)cfgPtBins;
964+
// auto effProt = (std::vector<float>)cfgProtonEff;
965+
// auto effAntiprot = (std::vector<float>)cfgAntiprotonEff;
912966

913967
for (int i = 0; i < 16; ++i) {
914-
if (candidate.pt() >= ptBins[i] && candidate.pt() < ptBins[i + 1]) {
968+
if (candidate.pt() >= cfgPtBins.value[i] && candidate.pt() < cfgPtBins.value[i + 1]) {
915969
binIndex = i;
916970
break;
917971
}
@@ -921,9 +975,9 @@ struct NetprotonCumulantsMc {
921975
return 0.0; // Default efficiency (0% if outside bins)
922976
}
923977
if (candidate.sign() > 0)
924-
return effProt[binIndex];
978+
return cfgProtonEff.value[binIndex];
925979
if (candidate.sign() < 0)
926-
return effAntiprot[binIndex];
980+
return cfgAntiprotonEff.value[binIndex];
927981
return 0.0;
928982
}
929983
}
@@ -1030,8 +1084,8 @@ struct NetprotonCumulantsMc {
10301084

10311085
if (cfgIsCalculateError) {
10321086

1033-
float l_Random = fRndm->Rndm();
1034-
int sampleIndex = static_cast<int>(cfgNSubsample * l_Random);
1087+
float lRandom = fRndm->Rndm();
1088+
int sampleIndex = static_cast<int>(cfgNSubsample * lRandom);
10351089

10361090
histos.get<TProfile2D>(HIST("GenProf2D_mu1_netproton"))->Fill(cent, sampleIndex, std::pow(netProt, 1.0));
10371091
histos.get<TProfile2D>(HIST("GenProf2D_mu2_netproton"))->Fill(cent, sampleIndex, std::pow(netProt, 2.0));
@@ -1052,6 +1106,8 @@ struct NetprotonCumulantsMc {
10521106

10531107
void processMCRec(MyMCRecCollision const& collision, MyMCTracks const& tracks, aod::McCollisions const&, aod::McParticles const&)
10541108
{
1109+
// auto tracksWithITSPid = soa::Attach<MyMCTracks, aod::pidits::ITSNSigmaEl, aod::pidits::ITSNSigmaPi, aod::pidits::ITSNSigmaKa, aod::pidits::ITSNSigmaPr>(tracks);
1110+
10551111
if (!collision.sel8()) {
10561112
return;
10571113
}
@@ -1071,6 +1127,8 @@ struct NetprotonCumulantsMc {
10711127
std::array<float, 7> fTCP0 = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
10721128
std::array<float, 7> fTCP1 = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
10731129

1130+
o2::aod::ITSResponse itsResponse;
1131+
10741132
// Start of the Monte-Carlo reconstructed tracks
10751133
for (const auto& track : tracks) {
10761134
if (!track.has_mcParticle()) //! check if track has corresponding MC particle
@@ -1086,6 +1144,10 @@ struct NetprotonCumulantsMc {
10861144
if ((particle.pt() < cfgCutPtLower) || (particle.pt() > 5.0f) || (std::abs(particle.eta()) > 0.8f)) {
10871145
continue;
10881146
}
1147+
if (!(track.itsNCls() > cfgITScluster) || !(track.tpcNClsFound() >= cfgTPCcluster) || !(track.tpcNClsCrossedRows() >= cfgTPCnCrossedRows)) {
1148+
continue;
1149+
}
1150+
10891151
if (particle.isPhysicalPrimary()) {
10901152
histos.fill(HIST("hrecPartPtAll"), particle.pt());
10911153
histos.fill(HIST("hrecPtAll"), track.pt());
@@ -1096,12 +1158,24 @@ struct NetprotonCumulantsMc {
10961158

10971159
bool trackSelected = false;
10981160
if (cfgPIDchoice == 0)
1099-
trackSelected = selectionPIDold(track);
1161+
trackSelected = selectionPIDoldTOFveto(track);
11001162
if (cfgPIDchoice == 1)
11011163
trackSelected = selectionPIDnew(track);
1164+
if (cfgPIDchoice == 2)
1165+
trackSelected = selectionPIDold(track);
1166+
1167+
if (cfgUseItsPid) {
1168+
if (std::abs(itsResponse.nSigmaITS<o2::track::PID::Proton>(track)) > 3.0)
1169+
continue;
1170+
}
11021171

11031172
if (trackSelected) {
11041173
recEbyeCompleteCollisions(recCollisions.lastIndex(), particle.pt(), particle.eta(), track.sign());
1174+
// filling nSigma distribution
1175+
histos.fill(HIST("h2DnsigmaTpcVsPt"), track.pt(), track.tpcNSigmaPr());
1176+
histos.fill(HIST("h2DnsigmaTofVsPt"), track.pt(), track.tofNSigmaPr());
1177+
histos.fill(HIST("h2DnsigmaItsVsPt"), track.pt(), itsResponse.nSigmaITS<o2::track::PID::Proton>(track));
1178+
11051179
if (track.sign() > 0) {
11061180
histos.fill(HIST("hrecPartPtProton"), particle.pt()); //! hist for p rec
11071181
histos.fill(HIST("hrecPtProton"), track.pt()); //! hist for p rec
@@ -1959,6 +2033,8 @@ struct NetprotonCumulantsMc {
19592033

19602034
void processDataRec(AodCollisions::iterator const& coll, aod::BCsWithTimestamps const&, AodTracks const& inputTracks)
19612035
{
2036+
// auto inputTracksWithPid = soa::Attach<AodTracks, aod::pidits::ITSNSigmaEl, aod::pidits::ITSNSigmaPi, aod::pidits::ITSNSigmaKa, aod::pidits::ITSNSigmaPr>(inputTracks);
2037+
19622038
if (!coll.sel8()) {
19632039
return;
19642040
}
@@ -1982,6 +2058,8 @@ struct NetprotonCumulantsMc {
19822058
std::array<float, 7> fTCP0 = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
19832059
std::array<float, 7> fTCP1 = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
19842060

2061+
o2::aod::ITSResponse itsResponse;
2062+
19852063
// Start of the Monte-Carlo reconstructed tracks
19862064
for (const auto& track : inputTracks) {
19872065
if (!track.isPVContributor()) //! track check as used in data
@@ -1991,6 +2069,9 @@ struct NetprotonCumulantsMc {
19912069
if ((track.pt() < cfgCutPtLower) || (track.pt() > 5.0f) || (std::abs(track.eta()) > 0.8f)) {
19922070
continue;
19932071
}
2072+
if (!(track.itsNCls() > cfgITScluster) || !(track.tpcNClsFound() >= cfgTPCcluster) || !(track.tpcNClsCrossedRows() >= cfgTPCnCrossedRows)) {
2073+
continue;
2074+
}
19942075

19952076
histos.fill(HIST("hrecPtAll"), track.pt());
19962077
histos.fill(HIST("hrecEtaAll"), track.eta());
@@ -2000,12 +2081,23 @@ struct NetprotonCumulantsMc {
20002081

20012082
bool trackSelected = false;
20022083
if (cfgPIDchoice == 0)
2003-
trackSelected = selectionPIDold(track);
2084+
trackSelected = selectionPIDoldTOFveto(track);
20042085
if (cfgPIDchoice == 1)
20052086
trackSelected = selectionPIDnew(track);
2087+
if (cfgPIDchoice == 2)
2088+
trackSelected = selectionPIDold(track);
2089+
2090+
if (cfgUseItsPid) {
2091+
if (std::abs(itsResponse.nSigmaITS<o2::track::PID::Proton>(track)) > 3.0)
2092+
continue;
2093+
}
20062094

20072095
if (trackSelected) {
20082096
recEbyeCompleteCollisions(recCollisions.lastIndex(), track.pt(), track.eta(), track.sign());
2097+
// filling nSigma distribution
2098+
histos.fill(HIST("h2DnsigmaTpcVsPt"), track.pt(), track.tpcNSigmaPr());
2099+
histos.fill(HIST("h2DnsigmaTofVsPt"), track.pt(), track.tofNSigmaPr());
2100+
histos.fill(HIST("h2DnsigmaItsVsPt"), track.pt(), itsResponse.nSigmaITS<o2::track::PID::Proton>(track));
20092101

20102102
// for protons
20112103
if (track.sign() > 0) {

0 commit comments

Comments
 (0)