Skip to content

Commit

Permalink
added function for antip reweighting (AliceO2Group#7211)
Browse files Browse the repository at this point in the history
  • Loading branch information
@alcaliva
alcaliva authored Aug 8, 2024
1 parent 3f94e78 commit 328dd1a
Showing 1 changed file with 223 additions and 8 deletions.
231 changes: 223 additions & 8 deletions PWGLF/Tasks/Nuspex/nuclei_in_jets.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,9 +188,9 @@ struct nuclei_in_jets {
registryMC.add("antiproton_all_ue", "antiproton_all_ue", HistType::kTH1F, {{100, 0.0, 5.0, "#it{p}_{T} (GeV/#it{c})"}});

// Antiproton Reweighting
registryMC.add("antiproton_incl", "antiproton_incl", HistType::kTH1F, {{100, 0.0, 5.0, "#it{p}_{T} (GeV/#it{c})"}});
registryMC.add("antiproton_jet", "antiproton_jet", HistType::kTH1F, {{100, 0.0, 5.0, "#it{p}_{T} (GeV/#it{c})"}});
registryMC.add("antiproton_ue", "antiproton_ue", HistType::kTH1F, {{100, 0.0, 5.0, "#it{p}_{T} (GeV/#it{c})"}});
registryMC.add("antiproton_eta_pt_pythia", "antiproton_eta_pt_pythia", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {16, -0.8, 0.8, "#it{#eta}"}});
registryMC.add("antiproton_eta_pt_jet", "antiproton_eta_pt_jet", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {16, -0.8, 0.8, "#it{#eta}"}});
registryMC.add("antiproton_eta_pt_ue", "antiproton_eta_pt_ue", HistType::kTH2F, {{100, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {16, -0.8, 0.8, "#it{#eta}"}});
}

// Single-Track Selection for Particles inside Jets
Expand Down Expand Up @@ -1035,29 +1035,244 @@ struct nuclei_in_jets {
float deltaPhi_ue2 = GetDeltaPhi(particle_dir.Phi(), ue_axis2.Phi());
float deltaR_ue2 = sqrt(deltaEta_ue2 * deltaEta_ue2 + deltaPhi_ue2 * deltaPhi_ue2);

// Inclusive antiprotons
registryMC.fill(HIST("antiproton_incl"), track.pt());

if (deltaR_jet < Rmax) {
registryMC.fill(HIST("antiproton_all_jet"), track.pt());
if (particle.isPhysicalPrimary()) {
registryMC.fill(HIST("antiproton_prim_jet"), track.pt());
registryMC.fill(HIST("antiproton_jet"), track.pt());
}
}
if (deltaR_ue1 < Rmax || deltaR_ue2 < Rmax) {
registryMC.fill(HIST("antiproton_all_ue"), track.pt());
if (particle.isPhysicalPrimary()) {
registryMC.fill(HIST("antiproton_prim_ue"), track.pt());
registryMC.fill(HIST("antiproton_ue"), track.pt());
}
}
}
}
}

void processAntipReweighting(o2::aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles)
{
for (const auto& mccollision : mcCollisions) {

// Selection on z_{vertex}
if (abs(mccollision.posZ()) > 10)
continue;

// MC Particles per Collision
auto mcParticles_per_coll = mcParticles.sliceBy(perMCCollision, mccollision.globalIndex());

// Reduced Event
std::vector<int> particle_ID;
std::vector<int> particle_ID_copy;
int leading_ID(0);
double pt_max(0);
int i = -1;

for (auto& particle : mcParticles_per_coll) {

i++;

if (particle.isPhysicalPrimary() && particle.pdgCode() == -2212 && particle.y() > min_y && particle.y() < max_y) {
registryMC.fill(HIST("antiproton_eta_pt_pythia"), particle.pt(), particle.eta());
}

// Select Primary Particles
double dx = particle.vx() - mccollision.posX();
double dy = particle.vy() - mccollision.posY();
double dz = particle.vz() - mccollision.posZ();
double dcaxy = sqrt(dx * dx + dy * dy);
double dcaz = abs(dz);
if (dcaxy > (0.0105 * 0.035 / TMath::Power(particle.pt(), 1.1)))
continue;
if (dcaz > 2.0)
continue;
if (abs(particle.eta()) > 0.8)
continue;
if (particle.pt() < 0.1)
continue;

// PDG Selection
int pdg = abs(particle.pdgCode());
if ((pdg != 11) && (pdg != 211) && (pdg != 321) && (pdg != 2212))
continue;

// Find pt Leading
if (particle.pt() > pt_max) {
leading_ID = i;
pt_max = particle.pt();
}

// Store Array Elements
particle_ID.push_back(i);
particle_ID_copy.push_back(i);
}

// Momentum of the Leading Particle
auto const& leading_track = mcParticles_per_coll.iteratorAt(leading_ID);
TVector3 p_leading(leading_track.px(), leading_track.py(), leading_track.pz());
TVector3 p_jet(leading_track.px(), leading_track.py(), leading_track.pz());

// Skip Events with pt<pt_leading_min
if (pt_max < min_pt_leading)
continue;

// Labels
int exit(0);
int nPartAssociated(0);
int nParticles = static_cast<int>(particle_ID.size());
std::vector<int> jet_particle_ID;
jet_particle_ID.push_back(leading_ID);

// Jet Finder
do {
// Initialization
double distance_jet_min(1e+08);
double distance_bkg_min(1e+08);
int label_jet_particle(0);
int i_jet_particle(0);

for (int i = 0; i < nParticles; i++) {

// Skip Leading Particle & Elements already associated to the Jet
if (particle_ID[i] == leading_ID || particle_ID[i] == -1)
continue;

// Get Particle Momentum
auto stored_track = mcParticles_per_coll.iteratorAt(particle_ID[i]);
TVector3 p_particle(stored_track.px(), stored_track.py(), stored_track.pz());

// Variables
double one_over_pt2_part = 1.0 / (p_particle.Pt() * p_particle.Pt());
double one_over_pt2_lead = 1.0 / (p_jet.Pt() * p_jet.Pt());
double deltaEta = p_particle.Eta() - p_jet.Eta();
double deltaPhi = GetDeltaPhi(p_particle.Phi(), p_jet.Phi());
double min = Minimum(one_over_pt2_part, one_over_pt2_lead);
double Delta2 = deltaEta * deltaEta + deltaPhi * deltaPhi;

// Distances
double distance_jet = min * Delta2 / (Rjet * Rjet);
double distance_bkg = one_over_pt2_part;

// Find Minimum Distance Jet
if (distance_jet < distance_jet_min) {
distance_jet_min = distance_jet;
label_jet_particle = particle_ID[i];
i_jet_particle = i;
}

// Find Minimum Distance Bkg
if (distance_bkg < distance_bkg_min) {
distance_bkg_min = distance_bkg;
}
}

if (distance_jet_min <= distance_bkg_min) {

// Add Particle to Jet
jet_particle_ID.push_back(label_jet_particle);

// Update Momentum of Leading Particle
auto jet_track = mcParticles_per_coll.iteratorAt(label_jet_particle);
TVector3 p_i(jet_track.px(), jet_track.py(), jet_track.pz());
p_jet = p_jet + p_i;

// Remove Element
particle_ID[i_jet_particle] = -1;
nPartAssociated++;
}

if (nPartAssociated >= (nParticles - 1))
exit = 1;
if (distance_jet_min > distance_bkg_min)
exit = 2;
} while (exit == 0);

// Event Counter: Skip Events with jet not fully inside acceptance
if ((TMath::Abs(p_jet.Eta()) + Rmax) > max_eta)
continue;

// Perpendicular Cones for the UE Estimate
TVector3 ue_axis1(0.0, 0.0, 0.0);
TVector3 ue_axis2(0.0, 0.0, 0.0);
get_perpendicular_axis(p_jet, ue_axis1, +1.0);
get_perpendicular_axis(p_jet, ue_axis2, -1.0);

// Protection against delta<0
if (ue_axis1.Mag() == 0 || ue_axis2.Mag() == 0)
continue;

double NchJetPlusUE(0);
double ptJetPlusUE(0);
double ptJet(0);
double ptUE(0);

for (int i = 0; i < nParticles; i++) {

auto track = mcParticles_per_coll.iteratorAt(particle_ID_copy[i]);
TVector3 particle_dir(track.px(), track.py(), track.pz());
double deltaEta_jet = particle_dir.Eta() - p_jet.Eta();
double deltaPhi_jet = GetDeltaPhi(particle_dir.Phi(), p_jet.Phi());
double deltaR_jet = sqrt(deltaEta_jet * deltaEta_jet + deltaPhi_jet * deltaPhi_jet);
double deltaEta_ue1 = particle_dir.Eta() - ue_axis1.Eta();
double deltaPhi_ue1 = GetDeltaPhi(particle_dir.Phi(), ue_axis1.Phi());
double deltaR_ue1 = sqrt(deltaEta_ue1 * deltaEta_ue1 + deltaPhi_ue1 * deltaPhi_ue1);
double deltaEta_ue2 = particle_dir.Eta() - ue_axis2.Eta();
double deltaPhi_ue2 = GetDeltaPhi(particle_dir.Phi(), ue_axis2.Phi());
double deltaR_ue2 = sqrt(deltaEta_ue2 * deltaEta_ue2 + deltaPhi_ue2 * deltaPhi_ue2);

if (deltaR_jet < Rmax) {
NchJetPlusUE++;
ptJetPlusUE = ptJetPlusUE + track.pt();
}
if (deltaR_ue1 < Rmax || deltaR_ue2 < Rmax) {
ptUE = ptUE + track.pt();
}
}
ptJet = ptJetPlusUE - 0.5 * ptUE;

// Skip Events with n. particles in jet less than given value
if (NchJetPlusUE < min_nPartInJet)
continue;

// Skip Events with Jet Pt lower than threshold
if (ptJet < min_jet_pt)
continue;

// Generated Particles
for (auto& particle : mcParticles_per_coll) {

if (!particle.isPhysicalPrimary())
continue;
if (particle.pdgCode() != -2212)
continue;
if (particle.y() < min_y || particle.y() > max_y)
continue;

TVector3 particle_dir(particle.px(), particle.py(), particle.pz());
double deltaEta_jet = particle_dir.Eta() - p_jet.Eta();
double deltaPhi_jet = GetDeltaPhi(particle_dir.Phi(), p_jet.Phi());
double deltaR_jet = sqrt(deltaEta_jet * deltaEta_jet + deltaPhi_jet * deltaPhi_jet);
double deltaEta_ue1 = particle_dir.Eta() - ue_axis1.Eta();
double deltaPhi_ue1 = GetDeltaPhi(particle_dir.Phi(), ue_axis1.Phi());
double deltaR_ue1 = sqrt(deltaEta_ue1 * deltaEta_ue1 + deltaPhi_ue1 * deltaPhi_ue1);
double deltaEta_ue2 = particle_dir.Eta() - ue_axis2.Eta();
double deltaPhi_ue2 = GetDeltaPhi(particle_dir.Phi(), ue_axis2.Phi());
double deltaR_ue2 = sqrt(deltaEta_ue2 * deltaEta_ue2 + deltaPhi_ue2 * deltaPhi_ue2);

if (deltaR_jet < Rmax) {
registryMC.fill(HIST("antiproton_eta_pt_jet"), particle.pt(), particle.eta());
}
if (deltaR_ue1 < Rmax || deltaR_ue2 < Rmax) {
registryMC.fill(HIST("antiproton_eta_pt_ue"), particle.pt(), particle.eta());
}
}
}
}
PROCESS_SWITCH(nuclei_in_jets, processData, "Process Data", true);
PROCESS_SWITCH(nuclei_in_jets, processMC, "process MC", false);
PROCESS_SWITCH(nuclei_in_jets, processSecAntiprotons, "Process sec antip", false);
PROCESS_SWITCH(nuclei_in_jets, processAntipReweighting, "Process antip reweighting", false);
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
Expand Down

0 comments on commit 328dd1a

Please sign in to comment.