Skip to content

Commit 69ad281

Browse files
lorber98lorber98
committed
Add process function for MC model predictions
1 parent b6c8e03 commit 69ad281

1 file changed

Lines changed: 290 additions & 3 deletions

File tree

‎PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx‎

Lines changed: 290 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -334,7 +334,10 @@ struct StrangenessInJetsIons {
334334
}
335335

336336
// Histograms for mc generated
337-
if (doprocessMCgenerated) {
337+
if (doprocessMCgenerated || doprocessMCmodels) {
338+
if (doprocessMCgenerated && doprocessMCmodels) {
339+
LOG(fatal) << "processMCgenerated and processMCmodels cannot be selected at the same time. Exit." << endl;
340+
}
338341

339342
// Event counter
340343
registryMC.add("number_of_events_mc_gen", "number of gen events in mc", HistType::kTH1D, {{10, 0, 10, "Event Cuts"}});
@@ -1952,7 +1955,6 @@ struct StrangenessInJetsIons {
19521955
if (fjParticles.size() < 1)
19531956
continue;
19541957
registryMC.fill(HIST("number_of_events_mc_gen"), 2.5);
1955-
19561958

19571959
// Cluster MC particles into jets using anti-kt algorithm
19581960
fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef);
@@ -2237,7 +2239,7 @@ struct StrangenessInJetsIons {
22372239
const auto& mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, mcCollision.globalIndex());
22382240

22392241
FillMBEventHistoMCREC(collision, mcParticles, v0sPerColl,
2240-
cascPerColl, tracksPerColl, multiplicity);
2242+
cascPerColl, tracksPerColl, multiplicity);
22412243

22422244
// Loop over reconstructed tracks
22432245
for (auto const& track : tracksPerColl) {
@@ -2570,6 +2572,291 @@ struct StrangenessInJetsIons {
25702572
}
25712573
}
25722574
PROCESS_SWITCH(StrangenessInJetsIons, processMCreconstructed, "process reconstructed events", false);
2575+
2576+
void processMCmodels(soa::Join<aod::McCollisions, aod::McCentFT0Ms, aod::McCentFT0Cs> const& collisions,
2577+
aod::McParticles const& mcParticles)
2578+
{
2579+
// Define per-event particle containers
2580+
std::vector<fastjet::PseudoJet> fjParticles;
2581+
std::vector<TVector3> strHadronMomentum;
2582+
std::vector<int> pdg;
2583+
2584+
// Jet and area definitions
2585+
fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet);
2586+
fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0));
2587+
2588+
// Loop over MC collisions
2589+
for (const auto& collision : collisions) {
2590+
// Clear containers at the start of the event loop
2591+
fjParticles.clear();
2592+
strHadronMomentum.clear();
2593+
pdg.clear();
2594+
2595+
// Fill event counter before any selection
2596+
registryMC.fill(HIST("number_of_events_mc_gen"), 0.5);
2597+
2598+
// Require vertex position within the allowed z range
2599+
// if (std::fabs(collision.posZ()) > zVtx)
2600+
// continue;
2601+
2602+
// Fill event counter after selection on z-vertex
2603+
registryMC.fill(HIST("number_of_events_mc_gen"), 1.5);
2604+
2605+
// Multiplicity of generated event retrived from corresponding MC RECO collision
2606+
// float genMultiplicity = multiplicity;
2607+
2608+
// Multiplicity of generated event
2609+
float genMultiplicity;
2610+
if (centrEstimator == 0) {
2611+
genMultiplicity = collision.centFT0C();
2612+
} else {
2613+
genMultiplicity = collision.centFT0M();
2614+
}
2615+
registryMC.fill(HIST("number_of_events_vsmultiplicity_gen_MB"), genMultiplicity);
2616+
2617+
// MC particles per collision
2618+
auto mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, collision.globalIndex());
2619+
2620+
// Loop over all MC particles and select physical primaries within acceptance
2621+
for (const auto& particle : mcParticlesPerColl) {
2622+
// Store properties of strange hadrons
2623+
int pdgAbs = std::abs(particle.pdgCode());
2624+
if (particle.isPhysicalPrimary() && (pdgAbs == kK0Short || pdgAbs == kLambda0 || pdgAbs == kXiMinus || pdgAbs == kOmegaMinus)) { // TODO: add protons, kaons, pions
2625+
pdg.emplace_back(particle.pdgCode());
2626+
strHadronMomentum.emplace_back(particle.px(), particle.py(), particle.pz());
2627+
}
2628+
2629+
double minPtParticle = 0.1f;
2630+
if (particle.eta() < configTracks.etaMin || particle.eta() > configTracks.etaMax || particle.pt() < minPtParticle)
2631+
continue;
2632+
2633+
// false --> do not fill histograms with reco events
2634+
FillMBEventHistoMCGEN(particle, genMultiplicity, false);
2635+
2636+
// Select physical primary particles or HF decay products
2637+
if (!isPhysicalPrimaryOrFromHF(particle, mcParticles))
2638+
continue;
2639+
2640+
// Build 4-momentum assuming charged pion mass
2641+
static constexpr float kMassPionChargedSquared = o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged;
2642+
const double energy = std::sqrt(particle.p() * particle.p() + kMassPionChargedSquared);
2643+
fastjet::PseudoJet fourMomentum(particle.px(), particle.py(), particle.pz(), energy);
2644+
fourMomentum.set_user_index(particle.pdgCode());
2645+
fjParticles.emplace_back(fourMomentum);
2646+
}
2647+
2648+
// Skip events with no particles
2649+
if (fjParticles.size() < 1)
2650+
continue;
2651+
registryMC.fill(HIST("number_of_events_mc_gen"), 2.5);
2652+
2653+
// Cluster MC particles into jets using anti-kt algorithm
2654+
fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef);
2655+
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
2656+
2657+
// Estimate background energy density (rho) in perpendicular cone
2658+
auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet);
2659+
2660+
// Loop over clustered jets
2661+
int countSelJet = 0; // number of selected jets
2662+
for (const auto& jet : jets) {
2663+
2664+
// Jet must be fully contained in acceptance
2665+
if ((std::fabs(jet.eta()) + rJet) > (configTracks.etaMax - deltaEtaEdge))
2666+
continue;
2667+
2668+
// Subtract background energy from jet
2669+
auto jetForSub = jet;
2670+
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp);
2671+
2672+
// Apply jet pT threshold
2673+
if (jetMinusBkg.pt() < minJetPt)
2674+
continue;
2675+
countSelJet++;
2676+
registryMC.fill(HIST("number_of_events_mc_gen"), 3.5);
2677+
registryMC.fill(HIST("number_of_events_vsmultiplicity_gen"), genMultiplicity);
2678+
2679+
// Fill jet counter
2680+
registryMC.fill(HIST("n_jets_vs_mult_pT_mc_gen"), genMultiplicity, jetMinusBkg.pt());
2681+
2682+
// Set up two perpendicular cone axes for underlying event estimation
2683+
TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
2684+
double coneRadius = std::sqrt(jet.area() / PI);
2685+
TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);
2686+
getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2);
2687+
if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) {
2688+
continue;
2689+
}
2690+
2691+
// Loop over strange hadrons
2692+
int index = -1;
2693+
for (const auto& hadron : strHadronMomentum) {
2694+
// Particle index
2695+
index++;
2696+
2697+
// Compute distance of particles from jet and UE axes
2698+
double deltaEtaJet = hadron.Eta() - jetAxis.Eta();
2699+
double deltaPhiJet = getDeltaPhi(hadron.Phi(), jetAxis.Phi());
2700+
double deltaRJet = std::sqrt(deltaEtaJet * deltaEtaJet + deltaPhiJet * deltaPhiJet);
2701+
double deltaEtaUe1 = hadron.Eta() - ueAxis1.Eta();
2702+
double deltaPhiUe1 = getDeltaPhi(hadron.Phi(), ueAxis1.Phi());
2703+
double deltaRUe1 = std::sqrt(deltaEtaUe1 * deltaEtaUe1 + deltaPhiUe1 * deltaPhiUe1);
2704+
double deltaEtaUe2 = hadron.Eta() - ueAxis2.Eta();
2705+
double deltaPhiUe2 = getDeltaPhi(hadron.Phi(), ueAxis2.Phi());
2706+
double deltaRUe2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2);
2707+
2708+
// Select particles inside jet
2709+
if (deltaRJet < coneRadius) {
2710+
switch (pdg[index]) {
2711+
case kK0Short:
2712+
if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) {
2713+
registryMC.fill(HIST("K0s_generated_jet"), genMultiplicity, hadron.Pt());
2714+
}
2715+
break;
2716+
case kLambda0:
2717+
if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) {
2718+
registryMC.fill(HIST("Lambda_generated_jet"), genMultiplicity, hadron.Pt());
2719+
}
2720+
break;
2721+
case kLambda0Bar:
2722+
if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) {
2723+
registryMC.fill(HIST("AntiLambda_generated_jet"), genMultiplicity, hadron.Pt());
2724+
}
2725+
break;
2726+
case kXiMinus:
2727+
if (particleOfInterestDict[ParticleOfInterest::kCascades]) {
2728+
registryMC.fill(HIST("XiNeg_generated_jet"), genMultiplicity, hadron.Pt());
2729+
}
2730+
break;
2731+
case kXiPlusBar:
2732+
if (particleOfInterestDict[ParticleOfInterest::kCascades]) {
2733+
registryMC.fill(HIST("XiPos_generated_jet"), genMultiplicity, hadron.Pt());
2734+
}
2735+
break;
2736+
case kOmegaMinus:
2737+
if (particleOfInterestDict[ParticleOfInterest::kCascades]) {
2738+
registryMC.fill(HIST("OmegaNeg_generated_jet"), genMultiplicity, hadron.Pt());
2739+
}
2740+
break;
2741+
case kOmegaPlusBar:
2742+
if (particleOfInterestDict[ParticleOfInterest::kCascades]) {
2743+
registryMC.fill(HIST("OmegaPos_generated_jet"), genMultiplicity, hadron.Pt());
2744+
}
2745+
break;
2746+
case kPiPlus:
2747+
if (particleOfInterestDict[ParticleOfInterest::kPions]) {
2748+
registryMC.fill(HIST("PionPos_generated_jet"), genMultiplicity, hadron.Pt());
2749+
}
2750+
break;
2751+
case kKPlus:
2752+
if (particleOfInterestDict[ParticleOfInterest::kKaons]) {
2753+
registryMC.fill(HIST("KaonPos_generated_jet"), genMultiplicity, hadron.Pt());
2754+
}
2755+
break;
2756+
case kProton:
2757+
if (particleOfInterestDict[ParticleOfInterest::kProtons]) {
2758+
registryMC.fill(HIST("ProtonPos_generated_jet"), genMultiplicity, hadron.Pt());
2759+
}
2760+
break;
2761+
case kPiMinus:
2762+
if (particleOfInterestDict[ParticleOfInterest::kPions]) {
2763+
registryMC.fill(HIST("PionNeg_generated_jet"), genMultiplicity, hadron.Pt());
2764+
}
2765+
break;
2766+
case kKMinus:
2767+
if (particleOfInterestDict[ParticleOfInterest::kKaons]) {
2768+
registryMC.fill(HIST("KaonNeg_generated_jet"), genMultiplicity, hadron.Pt());
2769+
}
2770+
break;
2771+
case kProtonBar:
2772+
if (particleOfInterestDict[ParticleOfInterest::kProtons]) {
2773+
registryMC.fill(HIST("ProtonNeg_generated_jet"), genMultiplicity, hadron.Pt());
2774+
}
2775+
break;
2776+
default:
2777+
break;
2778+
}
2779+
}
2780+
2781+
// Select particles inside UE cones
2782+
if (deltaRUe1 < coneRadius || deltaRUe2 < coneRadius) {
2783+
switch (pdg[index]) {
2784+
case kK0Short:
2785+
if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) {
2786+
registryMC.fill(HIST("K0s_generated_ue"), genMultiplicity, hadron.Pt());
2787+
}
2788+
break;
2789+
case kLambda0:
2790+
if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) {
2791+
registryMC.fill(HIST("Lambda_generated_ue"), genMultiplicity, hadron.Pt());
2792+
}
2793+
break;
2794+
case kLambda0Bar:
2795+
if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) {
2796+
registryMC.fill(HIST("AntiLambda_generated_ue"), genMultiplicity, hadron.Pt());
2797+
}
2798+
break;
2799+
case kXiMinus:
2800+
if (particleOfInterestDict[ParticleOfInterest::kCascades]) {
2801+
registryMC.fill(HIST("XiNeg_generated_ue"), genMultiplicity, hadron.Pt());
2802+
}
2803+
break;
2804+
case kXiPlusBar:
2805+
if (particleOfInterestDict[ParticleOfInterest::kCascades]) {
2806+
registryMC.fill(HIST("XiPos_generated_ue"), genMultiplicity, hadron.Pt());
2807+
}
2808+
break;
2809+
case kOmegaMinus:
2810+
if (particleOfInterestDict[ParticleOfInterest::kCascades]) {
2811+
registryMC.fill(HIST("OmegaNeg_generated_ue"), genMultiplicity, hadron.Pt());
2812+
}
2813+
break;
2814+
case kOmegaPlusBar:
2815+
if (particleOfInterestDict[ParticleOfInterest::kCascades]) {
2816+
registryMC.fill(HIST("OmegaPos_generated_ue"), genMultiplicity, hadron.Pt());
2817+
}
2818+
break;
2819+
case kPiPlus:
2820+
if (particleOfInterestDict[ParticleOfInterest::kPions]) {
2821+
registryMC.fill(HIST("PionPos_generated_ue"), genMultiplicity, hadron.Pt());
2822+
}
2823+
break;
2824+
case kKPlus:
2825+
if (particleOfInterestDict[ParticleOfInterest::kKaons]) {
2826+
registryMC.fill(HIST("KaonPos_generated_ue"), genMultiplicity, hadron.Pt());
2827+
}
2828+
break;
2829+
case kProton:
2830+
if (particleOfInterestDict[ParticleOfInterest::kProtons]) {
2831+
registryMC.fill(HIST("ProtonPos_generated_ue"), genMultiplicity, hadron.Pt());
2832+
}
2833+
break;
2834+
case kPiMinus:
2835+
if (particleOfInterestDict[ParticleOfInterest::kPions]) {
2836+
registryMC.fill(HIST("PionNeg_generated_ue"), genMultiplicity, hadron.Pt());
2837+
}
2838+
break;
2839+
case kKMinus:
2840+
if (particleOfInterestDict[ParticleOfInterest::kKaons]) {
2841+
registryMC.fill(HIST("KaonNeg_generated_ue"), genMultiplicity, hadron.Pt());
2842+
}
2843+
break;
2844+
case kProtonBar:
2845+
if (particleOfInterestDict[ParticleOfInterest::kProtons]) {
2846+
registryMC.fill(HIST("ProtonNeg_generated_ue"), genMultiplicity, hadron.Pt());
2847+
}
2848+
break;
2849+
default:
2850+
break;
2851+
}
2852+
}
2853+
}
2854+
}
2855+
// Fill jet counter
2856+
registryMC.fill(HIST("n_jets_vs_mult_mc_gen"), genMultiplicity, countSelJet);
2857+
}
2858+
}
2859+
PROCESS_SWITCH(StrangenessInJetsIons, processMCmodels, "Process MC generated events (with models)", false);
25732860
};
25742861

25752862
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)