Skip to content

Commit 24b9f59

Browse files
alcalivaalcaliva
authored
[PWGLF] replaced custom ITS PID selection with official nSigmaITS (AliceO2Group#9623)
1 parent 7502ae4 commit 24b9f59

File tree

1 file changed

+22
-80
lines changed

1 file changed

+22
-80
lines changed

‎PWGLF/Tasks/Nuspex/nucleiInJets.cxx‎

Lines changed: 22 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,7 @@
4545
#include "Common/DataModel/EventSelection.h"
4646
#include "Common/DataModel/Centrality.h"
4747
#include "Common/DataModel/PIDResponse.h"
48+
#include "Common/DataModel/PIDResponseITS.h"
4849

4950
using namespace std;
5051
using namespace o2;
@@ -95,11 +96,11 @@ struct NucleiInJets {
9596
Configurable<double> zVtx{"zVtx", 10.0, "Maximum zVertex"};
9697
Configurable<int> minNparticlesInJet{"minNparticlesInJet", 2, "Minimum number of particles inside jet"};
9798
Configurable<int> nJetsPerEventMax{"nJetsPerEventMax", 1000, "Maximum number of jets per event"};
98-
Configurable<bool> requireNoOverlap{"requireNoOverlap", false, "require no overlap between jets and UE cones"};
99+
Configurable<bool> requireNoOverlap{"requireNoOverlap", true, "require no overlap between jets and UE cones"};
99100

100101
// Track Parameters
101-
Configurable<double> par0{"par0", 0.004, "par 0"};
102-
Configurable<double> par1{"par1", 0.013, "par 1"};
102+
Configurable<double> par0{"par0", 0.00164, "par 0"};
103+
Configurable<double> par1{"par1", 0.00231, "par 1"};
103104
Configurable<int> minItsNclusters{"minItsNclusters", 5, "minimum number of ITS clusters"};
104105
Configurable<int> minTpcNclusters{"minTpcNclusters", 80, "minimum number of TPC clusters"};
105106
Configurable<int> minTpcNcrossedRows{"minTpcNcrossedRows", 80, "minimum number of TPC crossed pad rows"};
@@ -119,27 +120,15 @@ struct NucleiInJets {
119120
Configurable<bool> requirePvContributor{"requirePvContributor", true, "require that the track is a PV contributor"};
120121
Configurable<bool> setDCAselectionPtDep{"setDCAselectionPtDep", true, "require pt dependent selection"};
121122
Configurable<bool> applyReweighting{"applyReweighting", true, "apply reweighting"};
122-
123-
// Bethe-bloch Parametrization of ITS cluster size
124-
Configurable<double> bbPar0{"bbPar0", 0.00089176700, "Bethe Bloch Par 0"};
125-
Configurable<double> bbPar1{"bbPar1", 33.9651487037, "Bethe Bloch Par 1"};
126-
Configurable<double> bbPar2{"bbPar2", 0.42595677370, "Bethe Bloch Par 2"};
127-
Configurable<double> bbPar3{"bbPar3", 1.39638691440, "Bethe Bloch Par 3"};
128-
Configurable<double> bbPar4{"bbPar4", 7.97312623880, "Bethe Bloch Par 4"};
129-
Configurable<double> bbPar5{"bbPar5", 61.3838254956, "Bethe Bloch Par 5"};
130-
Configurable<double> bbPar6{"bbPar6", 2.30000000000, "Bethe Bloch Par 6"};
131-
Configurable<double> bbPar7{"bbPar7", 0.93827208820, "Bethe Bloch Par 7"};
132-
Configurable<double> bbPar8{"bbPar8", 1.0, "Bethe Bloch Par 8"};
133-
Configurable<double> resolClsSize{"resolClsSize", 0.214, "Resolution of cls size distribution"};
134-
Configurable<double> nSigmaClsSizeMax{"nSigmaClsSizeMax", 2.0, "nSigma cut on cluster size"};
123+
Configurable<bool> applyItsPid{"applyItsPid", true, "apply ITS PID"};
124+
Configurable<double> ptMaxItsPid{"ptMaxItsPid", 1.0, "maximum pt for ITS PID"};
125+
Configurable<double> nSigmaItsMin{"nSigmaItsMin", -2.0, "nSigmaITS min"};
126+
Configurable<double> nSigmaItsMax{"nSigmaItsMax", +2.0, "nSigmaITS max"};
135127
Configurable<std::string> urlToCcdb{"urlToCcdb", "http://alice-ccdb.cern.ch", "url of the personal ccdb"};
136128
Configurable<std::string> pathToFile{"pathToFile", "", "path to file with reweighting"};
137129
Configurable<std::string> histoNameWeightAntipJet{"histoNameWeightAntipJet", "", "reweighting histogram: antip in jet"};
138130
Configurable<std::string> histoNameWeightAntipUe{"histoNameWeightAntipUe", "", "reweighting histogram: antip in ue"};
139131

140-
// Bethe-Bloch
141-
TF1* bbClsSize = nullptr;
142-
143132
TH2F* twoDweightsAntipJet;
144133
TH2F* twoDweightsAntipUe;
145134

@@ -250,17 +239,6 @@ struct NucleiInJets {
250239
registryMC.add("antiproton_eta_pt_pythia", "antiproton_eta_pt_pythia", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}});
251240
registryMC.add("antiproton_eta_pt_jet", "antiproton_eta_pt_jet", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}});
252241
registryMC.add("antiproton_eta_pt_ue", "antiproton_eta_pt_ue", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}});
253-
254-
bbClsSize = new TF1("bbClsSize", betheBloch, 0.1, 10, 9);
255-
bbClsSize->SetParameter(0, bbPar0);
256-
bbClsSize->SetParameter(1, bbPar1);
257-
bbClsSize->SetParameter(2, bbPar2);
258-
bbClsSize->SetParameter(3, bbPar3);
259-
bbClsSize->SetParameter(4, bbPar4);
260-
bbClsSize->SetParameter(5, bbPar5);
261-
bbClsSize->SetParameter(6, bbPar6);
262-
bbClsSize->SetParameter(7, bbPar7);
263-
bbClsSize->SetParameter(8, bbPar8);
264242
}
265243

266244
// Single-Track Selection for Particles inside Jets
@@ -444,36 +422,6 @@ struct NucleiInJets {
444422
return false;
445423
}
446424

447-
double trackInclination(double eta)
448-
{
449-
double lambda(0);
450-
double theta = 2.0 * std::atan(std::exp(-eta));
451-
if (theta <= o2::constants::math::PIHalf)
452-
lambda = o2::constants::math::PIHalf - theta;
453-
if (theta > o2::constants::math::PIHalf)
454-
lambda = theta - o2::constants::math::PIHalf;
455-
return lambda;
456-
}
457-
458-
static double betheBloch(double* x, double* par)
459-
{
460-
// 5 parameters for the bethe bloch from 0 to 4
461-
// 1 parameter for the mip mpar[5]
462-
// 1 parameter for the charge exponent mpar[6]
463-
// 1 parameter for the mass mpar[7]
464-
// 1 parameter for the charge mpar[8]
465-
return par[5] * betheBlochAleph(x[0] / par[7], par[0], par[1], par[2], par[3], par[4]) * std::pow(par[8], par[6]);
466-
}
467-
468-
static double betheBlochAleph(double bg, double kp1, double kp2, double kp3, double kp4, double kp5)
469-
{
470-
double beta = bg / std::sqrt(1.0 + bg * bg);
471-
double aa = std::pow(beta, kp4);
472-
double bb = std::pow(1.0 / bg, kp5);
473-
bb = std::log(kp3 + bb);
474-
return (kp2 - aa - bb) * kp1 / aa;
475-
}
476-
477425
void getReweightingHistograms(o2::framework::Service<o2::ccdb::BasicCCDBManager> const& ccdbObj, TString filepath, TString histname_antip_jet, TString histname_antip_ue)
478426
{
479427
TList* l = ccdbObj->get<TList>(filepath.Data());
@@ -516,6 +464,9 @@ struct NucleiInJets {
516464
// Event Counter: after z-vertex cut
517465
registryData.fill(HIST("number_of_events_data"), 2.5);
518466

467+
// ITS Response
468+
o2::aod::ITSResponse itsResponse;
469+
519470
// List of Tracks
520471
std::vector<TVector3> trk;
521472

@@ -664,7 +615,7 @@ struct NucleiInJets {
664615
for (int j = 0; j < static_cast<int>(jet.size()); j++) { // o2-linter: disable=[const-ref-in-for-loop]
665616
if (isSelected[j] == 0 || i == j)
666617
continue;
667-
if (overlap(jet[i], ue1[j], rJet) || overlap(jet[i], ue2[j], rJet))
618+
if (overlap(jet[i], ue1[j], rJet) || overlap(jet[i], ue2[j], rJet) || overlap(jet[i], jet[j], rJet))
668619
nOverlaps++;
669620
}
670621
}
@@ -701,25 +652,16 @@ struct NucleiInJets {
701652
double dcaxy = track.dcaXY();
702653
double dcaz = track.dcaZ();
703654

704-
// ITS Cluster size
705-
double averageItsClusterSize(0);
706-
int nItsCls(0);
707-
for (int i = 0; i < 7; i++) { // o2-linter: disable=[const-ref-in-for-loop]
708-
int clusterSize = track.itsClsSizeInLayer(i);
709-
averageItsClusterSize += static_cast<double>(clusterSize);
710-
if (clusterSize > 0)
711-
nItsCls++;
712-
}
713-
averageItsClusterSize = averageItsClusterSize / static_cast<double>(nItsCls);
714-
double lambda = trackInclination(track.eta());
715-
double avgClsCosL = averageItsClusterSize * std::cos(lambda);
716-
double nsigma = (avgClsCosL - bbClsSize->Eval(pt)) / (resolClsSize * bbClsSize->Eval(pt));
717-
718-
bool isItsSelected = false;
719-
if (std::fabs(nsigma) < nSigmaClsSizeMax) {
720-
isItsSelected = true;
655+
// Selection on <ITS Cluster size>
656+
bool passedItsPid = false;
657+
if (itsResponse.nSigmaITS<o2::track::PID::Proton>(track) > nSigmaItsMin && itsResponse.nSigmaITS<o2::track::PID::Proton>(track) < nSigmaItsMax) {
658+
passedItsPid = true;
721659
}
722660

661+
bool passedItsPidSelection = true;
662+
if (applyItsPid && pt < ptMaxItsPid && (!passedItsPid))
663+
passedItsPidSelection = false;
664+
723665
TVector3 particleDirection(track.px(), track.py(), track.pz());
724666
double deltaEtaJet = particleDirection.Eta() - jet[i].Eta();
725667
double deltaPhiJet = getDeltaPhi(particleDirection.Phi(), jet[i].Phi());
@@ -753,7 +695,7 @@ struct NucleiInJets {
753695

754696
if (track.sign() < 0) { // only antimatter
755697
// Antiproton
756-
if (isItsSelected) {
698+
if (passedItsPidSelection) {
757699
if (pt < maxPtForNsigmaTpc)
758700
registryData.fill(HIST("antiproton_jet_tpc"), pt, nsigmaTPCPr);
759701
if (pt >= minPtForNsigmaTof && nsigmaTPCPr > minNsigmaTpc && nsigmaTPCPr < maxNsigmaTpc && track.hasTOF())
@@ -785,7 +727,7 @@ struct NucleiInJets {
785727

786728
if (track.sign() < 0) { // only antimatter
787729
// Antiproton
788-
if (isItsSelected) {
730+
if (passedItsPidSelection) {
789731
if (pt < maxPtForNsigmaTpc)
790732
registryData.fill(HIST("antiproton_ue_tpc"), pt, nsigmaTPCPr);
791733
if (pt >= minPtForNsigmaTof && nsigmaTPCPr > minNsigmaTpc && nsigmaTPCPr < maxNsigmaTpc && track.hasTOF())

0 commit comments

Comments
 (0)