Skip to content

Commit e038e5a

Browse files
pestratmpestratm
committed
Specify photon origins and general improvements
1 parent 2735b91 commit e038e5a

File tree

1 file changed

+102
-61
lines changed

1 file changed

+102
-61
lines changed

‎PWGJE/Tasks/hadronPhotonCorrelation.cxx‎

Lines changed: 102 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818

1919
#include <map>
2020
#include <string>
21+
#include <concepts>
2122
#include <vector>
2223

2324
#include "Framework/ASoA.h"
@@ -67,10 +68,8 @@ struct HadronPhotonCorrelation {
6768

6869
AxisSpec axisPhi = {72, 0., TwoPI, "#phi"}; // Axis for phi distribution
6970
AxisSpec axisDeltaPhi = {72, -PIHalf, 3 * PIHalf, "#Delta #phi"}; // Axis for Delta phi in correlations
70-
AxisSpec axisDeltaPhiDecay = {100, -.5, .5, "#Delta #phi"}; // Axis for Delta phi between neutral hadrons and decay photons
7171
AxisSpec axisEta = {40, -.8, .8, "#eta"}; // Axis for eta distribution
7272
AxisSpec axisDeltaEta = {80, -1.6, 1.6, "#Delta #eta"}; // Axis for Delta eta in correlations
73-
AxisSpec axisDeltaEtaDecay = {100, -.6, .6, "#Delta #eta"}; // Axis for Delta eta between neutral hadrons and decay photons
7473
ConfigurableAxis axisPtTrig = {"axisPtTrig",
7574
{VARIABLE_WIDTH,
7675
3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 9.0f, 11.0f, 15.0f, 20.0f},
@@ -80,8 +79,11 @@ struct HadronPhotonCorrelation {
8079
0.2, 0.5, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 9.0f, 11.0f, 15.0f},
8180
"p_{T, assoc} [GeV]"}; // Axis for associated particle pt distribution
8281
AxisSpec axisDeltaPt = {200, 0., 1.2, "#Delta p_T"}; // Axis for pt ratio between neutral hadrons and decay photons
83-
AxisSpec axisPid = {7, -1.5, 5.5, "pid"}; // Axis for PID of neutral hadrons
82+
AxisSpec axisPid = {9, -3.5, 5.5, "pid"}; // Axis for PID of neutral hadrons
8483
AxisSpec axisMult = {100, 0., 99., "N_{ch}"}; // Axis for mutplipicity
84+
AxisSpec axisAlpha = {100, 0., 1., "alpha"}; // Axis for decay photon pt assymetry
85+
86+
AxisSpec axisDeltaRDecay = {100, 0., 0.8, "#Delta R"}; // Axis for Delta R = sqrt(Delta eta^2 + Delta phi^2) between neutral hadrons and decay photons
8587

8688
float ptMinTrig;
8789
float ptMaxTrig;
@@ -91,31 +93,15 @@ struct HadronPhotonCorrelation {
9193
HistogramRegistry registry{"histogram registry"};
9294

9395
// Particle ids for storing neutral hadrons
94-
std::map<int, int> pidCodes = {
95-
{2212, -1},
96-
{1, -1},
97-
{2, -1},
98-
{3, -1},
99-
{4, -1},
100-
{5, -1},
101-
{6, -1},
102-
{21, -1}, // Protons, quarks, gluons (direct)
103-
{111, 1}, // pi0
104-
{221, 2}, // eta
105-
{223, 3}, // eta'
106-
{331, 4}, // phi
107-
{333, 5}}; // omega
96+
std::map<std::string, int> pidCodes = {
97+
{"pi0", 1}, // pi0
98+
{"eta", 2}, // eta
99+
{"eta'", 3}, // eta'
100+
{"phi", 4}, // phi
101+
{"omega", 5}}; // omega
108102

109103
Service<o2::framework::O2DatabasePDG> pdg;
110104

111-
// Calculate difference between two azimuthal angles, projecting into range [lowerPhi, lowerPhi + pi/2]
112-
float calculateDelta(float phi1, float phi2, float lowerPhi)
113-
{
114-
float dphi = fmod((phi1 - phi2) + 3 * PI, TwoPI) - PI;
115-
dphi = RecoDecay::constrainAngle(dphi, lowerPhi);
116-
return dphi;
117-
}
118-
119105
std::vector<int> eventSelectionBits;
120106
int trackSelection = -1;
121107

@@ -152,8 +138,8 @@ struct HadronPhotonCorrelation {
152138
////Neutral particles
153139
registry.add("generated/neutral/hNeutralCorrelGen", "Generated Trigger-Neutral Hadron Correlation", kTHnSparseF, {axisPtTrig, axisPtAssoc, axisDeltaEta, axisDeltaPhi, axisPid});
154140
registry.add("generated/neutral/hNeutralMultGen", "Generated Neutral Hadron Multiplicity", kTH1F, {axisMult});
155-
registry.add("generated/neutral/hNeutralSpectrumGen", "Generated Neutral Hadron Spectrum", kTHnSparseF, {axisPtAssoc, axisEta, axisPhi, axisPid}); // Particle ID of neutral hadrons
156-
registry.add("generated/neutral/hNeutralDecayGen", "Generated Neutral Hadron-Decay Photon Correlation", kTHnSparseF, {axisPtAssoc, axisDeltaPt, axisDeltaEtaDecay, axisDeltaPhiDecay, axisPid}); // Correlation with decay photons
141+
registry.add("generated/neutral/hNeutralSpectrumGen", "Generated Neutral Hadron Spectrum", kTHnSparseF, {axisPtAssoc, axisEta, axisPhi, axisPid}); // Particle ID of neutral hadrons
142+
registry.add("generated/neutral/hNeutralDecayGen", "Generated Neutral Hadron-Decay Photon Correlation", kTHnSparseF, {axisPtAssoc, axisDeltaPt, axisDeltaRDecay, axisAlpha, axisPid}); // Correlation with decay photons
157143

158144
// Reconstructed histograms
159145
// Triggers
@@ -192,6 +178,7 @@ struct HadronPhotonCorrelation {
192178
bool initTrack(const T& track)
193179
{
194180

181+
// if constexpr (HasHasMcParticle<T>) {
195182
if constexpr (HasHasMcParticle<T>::Value) {
196183
if (!track.has_mcParticle()) {
197184
return false;
@@ -233,6 +220,7 @@ struct HadronPhotonCorrelation {
233220
template <typename T>
234221
bool initTrig(const T& track)
235222
{
223+
// if constexpr (HasHasMcParticle<T>) {
236224
if constexpr (HasHasMcParticle<T>::Value) {
237225
if (!track.has_mcParticle()) {
238226
return false;
@@ -331,6 +319,7 @@ struct HadronPhotonCorrelation {
331319
if (!initTrig(track)) {
332320
continue;
333321
}
322+
334323
registry.fill(HIST("reconstructed/triggers/hTrigSpectrumReco"), track.pt(), track.eta(), track.phi());
335324
nTrigs++;
336325
}
@@ -413,7 +402,7 @@ struct HadronPhotonCorrelation {
413402
if (!initTrig(track)) {
414403
continue;
415404
}
416-
float dphi = calculateDelta(track.phi(), v0.phi(), -PIHalf);
405+
float dphi = RecoDecay::constrainAngle(track.phi() - v0.phi(), -PIHalf);
417406
registry.fill(HIST("reconstructed/photons/hPhotonCorrelReco"), track.pt(), v0.pt(), track.eta() - v0.eta(), dphi);
418407
}
419408
}
@@ -453,7 +442,7 @@ struct HadronPhotonCorrelation {
453442
if (!initTrig(track)) {
454443
continue;
455444
}
456-
float dphi = calculateDelta(track.phi(), v0.phi(), -PIHalf);
445+
float dphi = RecoDecay::constrainAngle(track.phi() - v0.phi(), -PIHalf);
457446
registry.fill(HIST("reconstructed/photons/hPhotonCorrelReco"), track.pt(), v0.pt(), track.eta() - v0.eta(), dphi);
458447
}
459448
}
@@ -462,27 +451,52 @@ struct HadronPhotonCorrelation {
462451
PROCESS_SWITCH(HadronPhotonCorrelation, processPhotonCorrelationsMCReco, "hadron-photon correlation", true);
463452

464453
void processPhotonCorrelationsMCGen(JetMcCollision const&,
465-
JetParticles const& tracks_true)
454+
Join<JetParticles, JMcParticlePIs> const& tracks_true)
466455
{
467456
int nPhotons = 0;
468457
for (const auto& track_assoc : tracks_true) {
469458
if (!initParticle(track_assoc, false)) {
470459
continue;
471460
}
472-
if (std::abs(track_assoc.pdgCode()) != 22) {
461+
if ((PDG_t)std::abs(track_assoc.pdgCode()) != kGamma) {
473462
continue;
474463
}
475-
if (!track_assoc.isPhysicalPrimary() && track_assoc.getGenStatusCode() == -1) {
464+
if (!track_assoc.isPhysicalPrimary() && track_assoc.getGenStatusCode() < 0) {
476465
continue;
477466
}
478467

479468
// Iterate through mother particles until original mother is reached
480-
auto mother = track_assoc.mothers_as<aod::JMcParticles>().at(0);
481-
while (mother.pdgCode() == 22) {
482-
mother = mother.mothers_as<aod::JMcParticles>().at(0);
469+
auto origPhoton = track_assoc;
470+
auto mother = track_assoc.mothers_as<Join<JetParticles, JMcParticlePIs>>().at(0);
471+
while ((PDG_t)mother.pdgCode() == kGamma) {
472+
origPhoton = mother;
473+
mother = mother.mothers_as<Join<JetParticles, JMcParticlePIs>>().at(0);
483474
}
484475

485-
registry.fill(HIST("generated/photons/hPhotonSpectrumGen"), track_assoc.pt(), track_assoc.eta(), track_assoc.phi(), pidCodes[mother.pdgCode()]);
476+
auto pdgMother = pdg->GetParticle(mother.pdgCode());
477+
if (!pdgMother)
478+
continue;
479+
int photonGeneration;
480+
switch (std::abs(origPhoton.getGenStatusCode())) {
481+
case 23: // prompt direct photons
482+
case 33:
483+
photonGeneration = -2;
484+
break;
485+
case 43: // shower photons
486+
case 44:
487+
case 51:
488+
case 52:
489+
photonGeneration = -1;
490+
break;
491+
case 91: // decay photons
492+
photonGeneration = pidCodes[pdgMother->GetName()];
493+
break;
494+
default:
495+
photonGeneration = -3;
496+
break;
497+
}
498+
499+
registry.fill(HIST("generated/photons/hPhotonSpectrumGen"), track_assoc.pt(), track_assoc.eta(), track_assoc.phi(), photonGeneration);
486500

487501
nPhotons++;
488502

@@ -493,8 +507,8 @@ struct HadronPhotonCorrelation {
493507
if (!initParticle(track_assoc)) {
494508
continue;
495509
}
496-
float dphi = calculateDelta(track_trig.phi(), track_assoc.phi(), -PIHalf);
497-
registry.fill(HIST("generated/photons/hPhotonCorrelGen"), track_trig.pt(), track_assoc.pt(), track_trig.eta() - track_assoc.eta(), dphi, pidCodes[mother.pdgCode()]);
510+
float dphi = RecoDecay::constrainAngle(track_trig.phi() - track_assoc.phi(), -PIHalf);
511+
registry.fill(HIST("generated/photons/hPhotonCorrelGen"), track_trig.pt(), track_assoc.pt(), track_trig.eta() - track_assoc.eta(), dphi, photonGeneration);
498512
}
499513
}
500514

@@ -530,11 +544,13 @@ struct HadronPhotonCorrelation {
530544
if (!jetderiveddatautilities::selectTrack(track_trig, trackSelection)) {
531545
continue;
532546
}
533-
547+
if (track_trig == track_assoc) {
548+
continue;
549+
}
534550
if (!initTrig(track_trig)) {
535551
continue;
536552
}
537-
float dphi = calculateDelta(track_trig.phi(), track_assoc.phi(), -PIHalf);
553+
float dphi = RecoDecay::constrainAngle(track_trig.phi() - track_assoc.phi(), -PIHalf);
538554
registry.fill(HIST("reconstructed/hadrons/hadrons/hHadronCorrelReco"), track_trig.pt(), track_assoc.pt(), track_trig.eta() - track_assoc.eta(), dphi);
539555
}
540556
}
@@ -567,7 +583,10 @@ struct HadronPhotonCorrelation {
567583
if (!initTrigParticle(track_trig)) {
568584
continue;
569585
}
570-
float dphi = calculateDelta(track_trig.phi(), track_assoc.phi(), -PIHalf);
586+
if (track_trig == track_assoc) {
587+
continue;
588+
}
589+
float dphi = RecoDecay::constrainAngle(track_trig.phi() - track_assoc.phi(), -PIHalf);
571590

572591
registry.fill(HIST("generated/hadrons/hHadronCorrelGen"), track_trig.pt(), track_assoc.pt(), track_trig.eta() - track_assoc.eta(), dphi);
573592
}
@@ -578,7 +597,7 @@ struct HadronPhotonCorrelation {
578597

579598
void processHadronCorrelationsMCReco(Join<JetCollisions, JMcCollisionLbs>::iterator const& collision_reco,
580599
JetMcCollisions const&,
581-
Join<JTracks, JMcTrackLbs> const& tracks_reco,
600+
Join<JetTracks, JMcTrackLbs> const& tracks_reco,
582601
JetParticles const&)
583602
{
584603
if (!jetderiveddatautilities::selectCollision(collision_reco, eventSelectionBits)) {
@@ -602,8 +621,6 @@ struct HadronPhotonCorrelation {
602621
if (std::abs(particle.pdgCode()) < 100) {
603622
continue;
604623
}
605-
// if(particle.isPhysicalPrimary()){continue;}
606-
// if(std::abs(particle.pdgCode())>1e9){continue;}
607624

608625
registry.fill(HIST("reconstructed/hadrons/hHadronSpectrumReco"), track_assoc.pt(), track_assoc.eta(), track_assoc.phi());
609626

@@ -618,11 +635,13 @@ struct HadronPhotonCorrelation {
618635
if (!jetderiveddatautilities::selectTrack(track_trig, trackSelection)) {
619636
continue;
620637
}
621-
638+
if (track_trig == track_assoc) {
639+
continue;
640+
}
622641
if (!initTrig(track_trig)) {
623642
continue;
624643
}
625-
float dphi = calculateDelta(track_trig.phi(), track_assoc.phi(), -PIHalf);
644+
float dphi = RecoDecay::constrainAngle(track_trig.phi() - track_assoc.phi(), -PIHalf);
626645

627646
registry.fill(HIST("reconstructed/hadrons/hHadronCorrelReco"), track_trig.pt(), track_assoc.pt(), track_trig.eta() - track_assoc.eta(), dphi);
628647
}
@@ -662,11 +681,13 @@ struct HadronPhotonCorrelation {
662681
if (!jetderiveddatautilities::selectTrack(track_trig, trackSelection)) {
663682
continue;
664683
}
665-
684+
if (track_trig == track_assoc) {
685+
continue;
686+
}
666687
if (!initTrig(track_trig)) {
667688
continue;
668689
}
669-
float dphi = calculateDelta(track_trig.phi(), track_assoc.phi(), -PIHalf);
690+
float dphi = RecoDecay::constrainAngle(track_trig.phi() - track_assoc.phi(), -PIHalf);
670691
registry.fill(HIST("reconstructed/charged/hPionCorrelReco"), track_trig.pt(), track_assoc.pt(), track_trig.eta() - track_assoc.eta(), dphi);
671692
}
672693
}
@@ -684,7 +705,7 @@ struct HadronPhotonCorrelation {
684705
if (!initParticle(track_assoc)) {
685706
continue;
686707
}
687-
if (std::abs(track_assoc.pdgCode()) != 211) {
708+
if ((PDG_t)std::abs(track_assoc.pdgCode()) != kPiPlus) {
688709
continue;
689710
}
690711

@@ -695,7 +716,10 @@ struct HadronPhotonCorrelation {
695716
if (!initTrigParticle(track_trig)) {
696717
continue;
697718
}
698-
float dphi = calculateDelta(track_trig.phi(), track_assoc.phi(), -PIHalf);
719+
if (track_trig == track_assoc) {
720+
continue;
721+
}
722+
float dphi = RecoDecay::constrainAngle(track_trig.phi() - track_assoc.phi(), -PIHalf);
699723

700724
registry.fill(HIST("generated/charged/hPionCorrelGen"), track_trig.pt(), track_assoc.pt(), track_trig.eta() - track_assoc.eta(), dphi);
701725
}
@@ -722,7 +746,7 @@ struct HadronPhotonCorrelation {
722746
if (!initTrack(track_assoc)) {
723747
continue;
724748
}
725-
if (std::abs(track_assoc.mcParticle().pdgCode()) != 211) {
749+
if ((PDG_t)std::abs(track_assoc.mcParticle().pdgCode()) != kPiPlus) {
726750
continue;
727751
}
728752

@@ -733,11 +757,14 @@ struct HadronPhotonCorrelation {
733757
if (!jetderiveddatautilities::selectTrack(track_trig, trackSelection)) {
734758
continue;
735759
}
760+
if (track_trig == track_assoc) {
761+
continue;
762+
}
736763

737764
if (!initTrig(track_trig)) {
738765
continue;
739766
}
740-
float dphi = calculateDelta(track_trig.phi(), track_assoc.phi(), -PIHalf);
767+
float dphi = RecoDecay::constrainAngle(track_trig.phi() - track_assoc.phi(), -PIHalf);
741768
registry.fill(HIST("reconstructed/charged/hPionCorrelReco"), track_trig.pt(), track_assoc.pt(), track_trig.eta() - track_assoc.eta(), dphi);
742769
}
743770
}
@@ -766,31 +793,45 @@ struct HadronPhotonCorrelation {
766793
if (pdgParticle->Charge() != 0.) {
767794
continue;
768795
} // remove charged particles
769-
if (track_assoc.pdgCode() < 100 || track_assoc.pdgCode() == 2112) {
796+
if (track_assoc.pdgCode() < 100 || (PDG_t)track_assoc.pdgCode() == kNeutron) {
770797
continue;
771-
} // remove non-hadrons and protons
772-
773-
registry.fill(HIST("generated/neutral/hNeutralSpectrumGen"), track_assoc.pt(), track_assoc.eta(), track_assoc.phi(), pidCodes[track_assoc.pdgCode()]);
798+
} // remove non-hadrons and neutrons
799+
registry.fill(HIST("generated/neutral/hNeutralSpectrumGen"), track_assoc.pt(), track_assoc.eta(), track_assoc.phi(), pidCodes[pdgParticle->GetName()]);
774800
nNeutrals++;
775801

776802
// Get correlations between neutral hadrons and their respective decay photons
777-
for (const auto& daughter : track_assoc.daughters_as<aod::JMcParticles>()) {
778-
if (daughter.pdgCode() != 22)
803+
auto daughters = track_assoc.daughters_as<aod::JMcParticles>();
804+
double alpha = -1;
805+
if (daughters.size() == 2) {
806+
auto daughter = daughters.begin();
807+
double pt1 = daughter.pt();
808+
++daughter;
809+
double pt2 = daughter.pt();
810+
alpha = std::abs((pt1 - pt2) / (pt1 + pt2));
811+
}
812+
813+
for (const auto& daughter : daughters) {
814+
if ((PDG_t)std::abs(daughter.pdgCode()) != kGamma)
779815
continue;
780816
if (!initParticle(daughter, false))
781817
continue;
782818
if (!daughter.isPhysicalPrimary() && daughter.getGenStatusCode() == -1)
783819
continue;
784-
registry.fill(HIST("generated/neutral/hNeutralDecayGen"), track_assoc.pt(), daughter.pt() / track_assoc.pt(), daughter.eta() - track_assoc.eta(), calculateDelta(daughter.phi(), track_assoc.phi(), -PIHalf), pidCodes[track_assoc.pdgCode()]);
820+
double deltaPt = daughter.pt() / track_assoc.pt();
821+
double deltaEta = daughter.eta() - track_assoc.eta();
822+
double deltaPhi = RecoDecay::constrainAngle(daughter.phi() - track_assoc.phi(), -PIHalf);
823+
double deltaR = std::sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi);
824+
825+
registry.fill(HIST("generated/neutral/hNeutralDecayGen"), track_assoc.pt(), deltaPt, deltaR, alpha, pidCodes[pdgParticle->GetName()]);
785826
}
786827

787828
// Get correlations between triggers and neutral hadrons
788829
for (const auto& track_trig : tracks_true) {
789830
if (!initTrigParticle(track_trig)) {
790831
continue;
791832
}
792-
float dphi = calculateDelta(track_assoc.phi(), track_trig.phi(), -PIHalf);
793-
registry.fill(HIST("generated/neutral/hNeutralCorrelGen"), track_trig.pt(), track_assoc.pt(), track_assoc.eta() - track_trig.eta(), dphi, pidCodes[track_assoc.pdgCode()]);
833+
float dphi = RecoDecay::constrainAngle(track_assoc.phi() - track_trig.phi(), -PIHalf);
834+
registry.fill(HIST("generated/neutral/hNeutralCorrelGen"), track_trig.pt(), track_assoc.pt(), track_assoc.eta() - track_trig.eta(), dphi, pidCodes[pdgParticle->GetName()]);
794835
}
795836
}
796837
registry.fill(HIST("generated/neutral/hNeutralMultGen"), nNeutrals);

0 commit comments

Comments
 (0)