Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
219 changes: 109 additions & 110 deletions PWGLF/Tasks/Nuspex/antinucleiInJets.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,7 @@
registryData.add("helium3_jet_tpc", "helium3_jet_tpc", HistType::kTH2F, {{nbins, min * 3, max * 3, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TPC}"}});
registryData.add("helium3_ue_tpc", "helium3_ue_tpc", HistType::kTH2F, {{nbins, min * 3, max * 3, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TPC}"}});

//systematic variations
// systematic variations
registryData.add("antiproton_tpc_syst", "antiproton_tpc_syst", HistType::kTHnSparseF, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TPC}"}, {10, 0, 10, "systematic uncertainty"}});
registryData.add("antiproton_tof_syst", "antiproton_tof_syst", HistType::kTHnSparseF, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TOF}"}, {10, 0, 10, "systematic uncertainty"}});
registryData.add("antideuteron_tpc_syst", "antideuteron_tpc_syst", HistType::kTHnSparseF, {{nbins, min * 2, max * 2, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TPC}"}, {10, 0, 10, "systematic uncertainty"}});
Expand Down Expand Up @@ -785,7 +785,7 @@
if (!particle.isPhysicalPrimary())
continue;

if (particle.pdgCode() == -2212) {

Check warning on line 788 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMC.fill(HIST("antiproton_eta_pt_pythia"), particle.pt(), particle.eta());
}

Expand Down Expand Up @@ -860,17 +860,17 @@
if ((2.0 * track.pt()) > ptMaxItsPidHel)
passedItsPidHel = true;

if (particle.pdgCode() == -2212)

Check warning on line 863 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMC.fill(HIST("antiproton_incl_all"), track.pt());

if (!particle.isPhysicalPrimary())
continue;

if (particle.pdgCode() == -2212)

Check warning on line 869 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registryMC.fill(HIST("antiproton_incl_prim"), track.pt());

// antiprotons
if (particle.pdgCode() == -2212 && passedItsPidProt) {

Check warning on line 873 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
if (nsigmaTPCPr > minNsigmaTpc && nsigmaTPCPr < maxNsigmaTpc) {
registryMC.fill(HIST("antiproton_incl_rec_tpc"), track.pt());
if (track.hasTOF() && nsigmaTOFPr > minNsigmaTof && nsigmaTOFPr < maxNsigmaTof)
Expand All @@ -879,7 +879,7 @@
}

// antideuterons
if (particle.pdgCode() == -1000010020 && passedItsPidDeut) {

Check warning on line 882 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
if (nsigmaTPCDe > minNsigmaTpc && nsigmaTPCDe < maxNsigmaTpc) {
registryMC.fill(HIST("antideuteron_incl_rec_tpc"), track.pt());
if (track.hasTOF() && nsigmaTOFDe > minNsigmaTof && nsigmaTOFDe < maxNsigmaTof)
Expand All @@ -888,7 +888,7 @@
}

// deuterons
if (particle.pdgCode() == 1000010020 && passedItsPidDeut) {

Check warning on line 891 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
if (nsigmaTPCDe > minNsigmaTpc && nsigmaTPCDe < maxNsigmaTpc) {
registryMC.fill(HIST("deuteron_incl_rec_tpc"), track.pt());
if (track.hasTOF() && nsigmaTOFDe > minNsigmaTof && nsigmaTOFDe < maxNsigmaTof)
Expand All @@ -897,14 +897,14 @@
}

// antihelium3
if (particle.pdgCode() == -1000020030 && passedItsPidHel) {

Check warning on line 900 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
if (nsigmaTPCHe > minNsigmaTpc && nsigmaTPCHe < maxNsigmaTpc) {
registryMC.fill(HIST("antihelium3_incl_rec_tpc"), 2.0 * track.pt());
}
}

// helium3
if (particle.pdgCode() == 1000020030 && passedItsPidHel) {

Check warning on line 907 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
if (nsigmaTPCHe > minNsigmaTpc && nsigmaTPCHe < maxNsigmaTpc) {
registryMC.fill(HIST("helium3_incl_rec_tpc"), 2.0 * track.pt());
}
Expand Down Expand Up @@ -993,7 +993,7 @@
if (deltaRUe1 > coneRadius && deltaRUe2 > coneRadius)
continue;

if (particle.pdgCode() != -2212)

Check warning on line 996 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
continue;

registryMC.fill(HIST("antiproton_ue_gen"), particle.pt());
Expand Down Expand Up @@ -1090,7 +1090,7 @@
if (!track.has_mcParticle())
continue;
const auto mcparticle = track.mcParticle();
if (mcparticle.pdgCode() != -2212)

Check warning on line 1093 in PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

Avoid using hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
continue;

// variables
Expand Down Expand Up @@ -1179,134 +1179,133 @@
}
PROCESS_SWITCH(AntinucleiInJets, processJetsMCrec, "process jets MC rec", false);

// Process Systematics
void processSystematicsData(SelectedCollisions::iterator const& collision, FullNucleiTracks const& tracks)
{
const int nSystematics = 10;
int itsNclustersSyst[nSystematics] = {5, 6, 5, 4, 5, 3, 5, 6, 3, 4};
float tpcNcrossedRowsSyst[nSystematics] = {100, 85, 80, 110, 95, 90, 105, 95, 100, 105};
float dcaxySyst[nSystematics] = {0.05, 0.07, 0.10, 0.03, 0.06, 0.15, 0.08, 0.04, 0.09, 0.10};
float dcazSyst[nSystematics] = {0.1, 0.15, 0.3, 0.075, 0.12, 0.18, 0.2, 0.1, 0.15, 0.2};

// Process Systematics
void processSystematicsData(SelectedCollisions::iterator const& collision, FullNucleiTracks const& tracks)
{
const int nSystematics = 10;
int itsNclustersSyst[nSystematics] = {5, 6, 5, 4, 5, 3, 5, 6, 3, 4};
float tpcNcrossedRowsSyst[nSystematics] = {100, 85, 80, 110, 95, 90, 105, 95, 100, 105};
float dcaxySyst[nSystematics] = {0.05, 0.07, 0.10, 0.03, 0.06, 0.15, 0.08, 0.04, 0.09, 0.10};
float dcazSyst[nSystematics] = {0.1, 0.15, 0.3, 0.075, 0.12, 0.18, 0.2, 0.1, 0.15, 0.2};

// event selection
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx)
return;
// event selection
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx)
return;

// loop over reconstructed tracks
int id(-1);
std::vector<fastjet::PseudoJet> fjParticles;
for (auto const& track : tracks) {
id++;
if (!passedTrackSelectionForJetReconstruction(track))
continue;

// 4-momentum representation of a particle
fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged));
fourMomentum.set_user_index(id);
fjParticles.emplace_back(fourMomentum);
}
// loop over reconstructed tracks
int id(-1);
std::vector<fastjet::PseudoJet> fjParticles;
for (auto const& track : tracks) {
id++;
if (!passedTrackSelectionForJetReconstruction(track))
continue;

// reject empty events
if (fjParticles.size() < 1)
return;
// 4-momentum representation of a particle
fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged));
fourMomentum.set_user_index(id);
fjParticles.emplace_back(fourMomentum);
}

// cluster particles using the anti-kt algorithm
fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet);
fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0)); // active_area_explicit_ghosts
fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef);
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
auto [rhoPerp, rhoMPerp] = backgroundSub.estimateRhoPerpCone(fjParticles, jets);
// reject empty events
if (fjParticles.size() < 1)
return;

// loop over reconstructed jets
for (auto& jet : jets) {
// cluster particles using the anti-kt algorithm
fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet);
fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0)); // active_area_explicit_ghosts
fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef);
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
auto [rhoPerp, rhoMPerp] = backgroundSub.estimateRhoPerpCone(fjParticles, jets);

// jet must be fully contained in the acceptance
if ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))
continue;
// loop over reconstructed jets
for (auto& jet : jets) {

// jet pt must be larger than threshold
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jet, rhoPerp, rhoMPerp);
if (getCorrectedPt(jetMinusBkg.pt()) < minJetPt)
continue;
// jet must be fully contained in the acceptance
if ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))
continue;

// get jet constituents
std::vector<fastjet::PseudoJet> jetConstituents = jet.constituents();
o2::aod::ITSResponse itsResponse;
// jet pt must be larger than threshold
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jet, rhoPerp, rhoMPerp);
if (getCorrectedPt(jetMinusBkg.pt()) < minJetPt)
continue;

// loop over jet constituents
for (const auto& particle : jetConstituents) {
for(int i = 0; i < nSystematics; i++) {
// get corresponding track and apply track selection criteria
auto const& track = tracks.iteratorAt(particle.user_index());
// get jet constituents
std::vector<fastjet::PseudoJet> jetConstituents = jet.constituents();
o2::aod::ITSResponse itsResponse;

// variables
double nsigmaTPCPr = track.tpcNSigmaPr();
double nsigmaTOFPr = track.tofNSigmaPr();
double nsigmaTPCDe = track.tpcNSigmaDe();
double nsigmaTOFDe = track.tofNSigmaDe();
double pt = track.pt();
double dcaxy = track.dcaXY();
double dcaz = track.dcaZ();
// loop over jet constituents
for (const auto& particle : jetConstituents) {
for (int i = 0; i < nSystematics; i++) {
// get corresponding track and apply track selection criteria
auto const& track = tracks.iteratorAt(particle.user_index());

if (requirePvContributor && !(track.isPVContributor()))
continue;
if (!track.hasITS())
continue;
if (track.itsNCls() < itsNclustersSyst[i])
continue;
if (!track.hasTPC())
continue;
if (track.tpcNClsCrossedRows() < tpcNcrossedRowsSyst[i])
continue;
if ((static_cast<double>(track.tpcNClsCrossedRows()) / static_cast<double>(track.tpcNClsFindable())) < minTpcNcrossedRowsOverFindable)
continue;
if (track.tpcChi2NCl() > maxChiSquareTpc)
continue;
if (track.itsChi2NCl() > maxChiSquareIts)
continue;
if (track.eta() < minEta || track.eta() > maxEta)
continue;
if (track.pt() < minPt)
continue;
if (std::fabs(dcaxy) > dcaxySyst[i])
continue;
if (std::fabs(dcaz) > dcazSyst[i])
continue;
// variables
double nsigmaTPCPr = track.tpcNSigmaPr();
double nsigmaTOFPr = track.tofNSigmaPr();
double nsigmaTPCDe = track.tpcNSigmaDe();
double nsigmaTOFDe = track.tofNSigmaDe();
double pt = track.pt();
double dcaxy = track.dcaXY();
double dcaz = track.dcaZ();

bool passedItsPidProt(false), passedItsPidDeut(false);
if (itsResponse.nSigmaITS<o2::track::PID::Proton>(track) > nSigmaItsMin && itsResponse.nSigmaITS<o2::track::PID::Proton>(track) < nSigmaItsMax) {
passedItsPidProt = true;
}
if (itsResponse.nSigmaITS<o2::track::PID::Deuteron>(track) > nSigmaItsMin && itsResponse.nSigmaITS<o2::track::PID::Deuteron>(track) < nSigmaItsMax) {
passedItsPidDeut = true;
}
if (!applyItsPid) {
passedItsPidProt = true;
passedItsPidDeut = true;
}
if (pt > ptMaxItsPidProt)
passedItsPidProt = true;
if (pt > ptMaxItsPidDeut)
passedItsPidDeut = true;
if (requirePvContributor && !(track.isPVContributor()))
continue;
if (!track.hasITS())
continue;
if (track.itsNCls() < itsNclustersSyst[i])
continue;
if (!track.hasTPC())
continue;
if (track.tpcNClsCrossedRows() < tpcNcrossedRowsSyst[i])
continue;
if ((static_cast<double>(track.tpcNClsCrossedRows()) / static_cast<double>(track.tpcNClsFindable())) < minTpcNcrossedRowsOverFindable)
continue;
if (track.tpcChi2NCl() > maxChiSquareTpc)
continue;
if (track.itsChi2NCl() > maxChiSquareIts)
continue;
if (track.eta() < minEta || track.eta() > maxEta)
continue;
if (track.pt() < minPt)
continue;
if (std::fabs(dcaxy) > dcaxySyst[i])
continue;
if (std::fabs(dcaz) > dcazSyst[i])
continue;

// antimatter
if (track.sign() < 0) {
if (passedItsPidProt) {
registryData.fill(HIST("antiproton_tpc_syst"), pt, nsigmaTPCPr, i);
if (nsigmaTPCPr > minNsigmaTpc && nsigmaTPCPr < maxNsigmaTpc && track.hasTOF())
registryData.fill(HIST("antiproton_tof_syst"), pt, nsigmaTOFPr, i);
bool passedItsPidProt(false), passedItsPidDeut(false);
if (itsResponse.nSigmaITS<o2::track::PID::Proton>(track) > nSigmaItsMin && itsResponse.nSigmaITS<o2::track::PID::Proton>(track) < nSigmaItsMax) {
passedItsPidProt = true;
}
if (passedItsPidDeut) {
registryData.fill(HIST("antideuteron_tpc_syst"), pt, nsigmaTPCDe, i);
if (nsigmaTPCDe > minNsigmaTpc && nsigmaTPCDe < maxNsigmaTpc && track.hasTOF())
registryData.fill(HIST("antideuteron_tof_syst"), pt, nsigmaTOFDe, i);
if (itsResponse.nSigmaITS<o2::track::PID::Deuteron>(track) > nSigmaItsMin && itsResponse.nSigmaITS<o2::track::PID::Deuteron>(track) < nSigmaItsMax) {
passedItsPidDeut = true;
}
if (!applyItsPid) {
passedItsPidProt = true;
passedItsPidDeut = true;
}
if (pt > ptMaxItsPidProt)
passedItsPidProt = true;
if (pt > ptMaxItsPidDeut)
passedItsPidDeut = true;

// antimatter
if (track.sign() < 0) {
if (passedItsPidProt) {
registryData.fill(HIST("antiproton_tpc_syst"), pt, nsigmaTPCPr, i);
if (nsigmaTPCPr > minNsigmaTpc && nsigmaTPCPr < maxNsigmaTpc && track.hasTOF())
registryData.fill(HIST("antiproton_tof_syst"), pt, nsigmaTOFPr, i);
}
if (passedItsPidDeut) {
registryData.fill(HIST("antideuteron_tpc_syst"), pt, nsigmaTPCDe, i);
if (nsigmaTPCDe > minNsigmaTpc && nsigmaTPCDe < maxNsigmaTpc && track.hasTOF())
registryData.fill(HIST("antideuteron_tof_syst"), pt, nsigmaTOFDe, i);
}
}
}
}
}
}
}
PROCESS_SWITCH(AntinucleiInJets, processSystematicsData, "Process Systematics", false);
};

Expand Down
Loading