Skip to content

[PWGUD] Adding new histos for angular correlation studies #13596

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rolavick merged 3 commits into from
Oct 30, 2025
Merged

[PWGUD] Adding new histos for angular correlation studies #13596

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8341cc5
taking off the NUA correction and adding new histos for correlation
JosueMtzGar Oct 29, 2025
0bb3ffa
taking off the NUA correction and adding new histos for correlation
JosueMtzGar Oct 29, 2025
055c5cb
resolving o2 linter issues
JosueMtzGar Oct 29, 2025
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
162 changes: 72 additions & 90 deletions PWGUD/Tasks/upcPhotonuclearAnalysisJMG.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,16 @@ struct UpcPhotonuclearAnalysisJMG {
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};

// Declare configurables on events/collisions
Configurable<int> minMultiplicity{"minMultiplicity", 2, {"Range on multiplicity"}};
Configurable<int> range1Max{"range1Max", 10, {"Range on multiplicity"}};
Configurable<int> range2Min{"range2Min", 11, {"Range on multiplicity"}};
Configurable<int> range2Max{"range2Max", 20, {"Range on multiplicity"}};
Configurable<int> range3Min{"range3Min", 21, {"Range on multiplicity"}};
Configurable<int> range3Max{"range3Max", 30, {"Range on multiplicity"}};
Configurable<int> range4Min{"range4Min", 31, {"Range on multiplicity"}};
Configurable<int> range4Max{"range4Max", 40, {"Range on multiplicity"}};
Configurable<int> range5Min{"range5Min", 41, {"Range on multiplicity"}};
Configurable<int> range5Max{"range5Max", 50, {"Range on multiplicity"}};
Configurable<int> nEventsMixed{"nEventsMixed", 3, {"Events to be Mixed"}};
Configurable<int> factorEventsMixed{"factorEventsMixed", 100, {"factorEventsMixed to events mixed"}};
Configurable<float> myZVtxCut{"myZVtxCut", 10., {"My collision cut"}};
Expand Down Expand Up @@ -149,7 +159,7 @@ struct UpcPhotonuclearAnalysisJMG {
"Pair cuts on various particles"};
Configurable<float> cfgTwoTrackCut{"cfgTwoTrackCut", -1, {"Two track cut"}};
ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for histograms"};
ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, kThreeHalfPi}, "delta phi axis for histograms"};
ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {32, -PIHalf, kThreeHalfPi}, "delta phi axis for histograms"};
ConfigurableAxis axisDeltaEta{"axisDeltaEta", {32, -1.6, 1.6}, "delta eta axis for histograms"};
ConfigurableAxis axisPtTrigger{"axisPtTrigger", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt trigger axis for histograms"};
ConfigurableAxis axisPtAssoc{"axisPtAssoc", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0}, "pt associated axis for histograms"};
Expand Down Expand Up @@ -198,7 +208,18 @@ struct UpcPhotonuclearAnalysisJMG {
histos.add("yields", "multiplicity vs pT vs eta", {HistType::kTH3F, {{100, 0, 100, "multiplicity"}, {40, 0, 20, "p_{T}"}, {100, -2, 2, "#eta"}}});
histos.add("etaphi", "multiplicity vs eta vs phi", {HistType::kTH3F, {{100, 0, 100, "multiplicity"}, {100, -2, 2, "#eta"}, {64, 0., TwoPI, "#varphi"}}});
histos.add("etaphiVtx", "vertex Z vs eta vs phi", {HistType::kTH3F, {{20, -10., 10., "vertex Z"}, {32, -0.8, 0.8, "#eta"}, {64, 0., TwoPI, "#varphi"}}});
histos.add("weightNUA", "weight per bin", {HistType::kTH3F, {{20, -10., 10., "vertex Z"}, {32, -0.8, 0.8, "#eta"}, {64, 0., TwoPI, "#varphi"}}});
histos.add("sameEvent2D", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("sameEvent_2_10", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("sameEvent_11_20", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("sameEvent_21_30", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("sameEvent_31_40", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("sameEvent_41_50", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("mixedEvent2D", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("mixedEvent_2_10", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("mixedEvent_11_20", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("mixedEvent_21_30", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("mixedEvent_31_40", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});
histos.add("mixedEvent_41_50", "#Delta #eta vs #Delta #phi", {HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}});

const int maxMixBin = axisMultiplicity->size() * axisVertex->size();
histos.add("eventcount", "bin", {HistType::kTH1F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}}});
Expand All @@ -215,8 +236,6 @@ struct UpcPhotonuclearAnalysisJMG {
histos.add("Events/hCountCollisions", "0 total - 1 side A - 2 side C - 3 both side; Number of analysed collision; counts", kTH1F, {axisCollision});
histos.add("Events/hCountCollisionsMixed", "0 total - 1 side A - 2 side C - 3 both side; Number of analysed collision; counts", kTH1F, {axisCollision});
histos.add("Tracks/hTracksAfterCuts", " ; ; counts", kTH1F, {axisCountTracks});
histos.add("Tracks/hTrackPhiBeforeCorr", "#it{#phi} distribution before NUA correction; #it{#phi}; counts", kTH1F, {axisPhi});
histos.add("Tracks/hTrackPhiAfterCorr", "#it{#phi} distribution after NUA correction; #it{#phi}; counts", kTH1F, {axisPhi});

// histos to selection gap in side A
histos.add("Tracks/SGsideA/hTrackPt", "#it{p_{T}} distribution; #it{p_{T}}; counts", kTH1F, {axisPt});
Expand Down Expand Up @@ -303,6 +322,11 @@ struct UpcPhotonuclearAnalysisJMG {
std::vector<double> vtxBinsEdges{VARIABLE_WIDTH, -10.0f, -7.0f, -5.0f, -2.5f, 0.0f, 2.5f, 5.0f, 7.0f, 10.0f};
std::vector<double> gapSideBinsEdges{VARIABLE_WIDTH, -0.5, 0.5, 1.5};

struct SameEventTag {
};
struct MixedEventTag {
};

SliceCache cache;
// int countEvents = 0;
// int countGapA = 0;
Expand Down Expand Up @@ -431,84 +455,16 @@ struct UpcPhotonuclearAnalysisJMG {
return true;
}

void makeNUAWeights(std::shared_ptr<TH3> histoRaw3D)
{
const int nPhi = histoRaw3D->GetZaxis()->GetNbins();
const int nEta = histoRaw3D->GetYaxis()->GetNbins();
const int nVz = histoRaw3D->GetXaxis()->GetNbins();

for (int jEtha = 1; jEtha <= nEta; ++jEtha) {
for (int iVtxZ = 1; iVtxZ <= nVz; ++iVtxZ) {
// average on phi to (eta_jEtha, vz_iVtxZ)
double sum = 0.0;
double nMax = 0.0;
int count = 0;
for (int kPhi = 1; kPhi <= nPhi; ++kPhi) {
double nEntry = histoRaw3D->GetBinContent(iVtxZ, jEtha, kPhi);
sum += nEntry;
count += 1.0;
if (nEntry > nMax) {
nMax = nEntry;
}
}
double nMean;
if (useNMax) {
nMean = nMax;
} else {
nMean = (count > 0) ? sum / count : 0.0;
}

for (int kPhi = 1; kPhi <= nPhi; ++kPhi) {
double nEntry = histoRaw3D->GetBinContent(iVtxZ, jEtha, kPhi);
double w;
if (useEpsilon) {
if (nMean > 0) {
w = nMean / std::max(nEntry, static_cast<double>(myEpsilonToWeight));
} else {
w = 1.0;
}
} else {
if (nMean > 0) {
w = nMean / nEntry;
} else {
w = 1.0;
}
}
if (w < myWeightMin)
w = myWeightMin;
if (w > myWeightMax)
w = myWeightMax;
if (auto histoWeightNUA = histos.get<TH3>(HIST("weightNUA"))) {
histoWeightNUA->SetBinContent(iVtxZ, jEtha, kPhi, w);
}
}
}
}
}

float getNUAWeight(float vz, float eta, float phi)
{
auto hWeight = histos.get<TH3>(HIST("weightNUA"));
phi = RecoDecay::constrainAngle(phi, 0.f);
int iPhi = hWeight->GetZaxis()->FindBin(phi);
int iEta = hWeight->GetYaxis()->FindBin(eta);
int iVz = hWeight->GetXaxis()->FindBin(vz);
return hWeight->GetBinContent(iVz, iEta, iPhi);
}

template <typename TTarget, typename TTracks>
void fillCorrelationsUD(TTarget target, const TTracks tracks1, const TTracks tracks2, float multiplicity, float posZ)
template <typename TTarget, typename TTracks, typename TTag>
void fillCorrelationsUD(TTarget target, const TTracks& tracks1, const TTracks& tracks2, float multiplicity, float posZ, TTag)
{
// multiplicity = tracks1.size();
for (const auto& track1 : tracks1) {
if (isTrackCut(track1) == false) {
return;
}
// weight NUA for track1
float phi1 = phi(track1.px(), track1.py());
phi1 = RecoDecay::constrainAngle(phi1, 0.f);
float eta1 = eta(track1.px(), track1.py(), track1.pz());
float w1 = getNUAWeight(posZ, eta1, phi1);
target->getTriggerHist()->Fill(CorrelationContainer::kCFStepReconstructed, track1.pt(), multiplicity, posZ, 1.0);
for (const auto& track2 : tracks2) {
if (track1 == track2) {
Expand All @@ -517,13 +473,9 @@ struct UpcPhotonuclearAnalysisJMG {
if (isTrackCut(track2) == false) {
return;
}
// weight NUA for track 2
float phi2 = phi(track2.px(), track2.py());
phi2 = RecoDecay::constrainAngle(phi2, 0.f);
float eta2 = eta(track2.px(), track2.py(), track2.pz());
float w2 = getNUAWeight(posZ, eta2, phi2);
// total weight
float wPair = w1 * w2;
/*if (doPairCuts && mPairCuts.conversionCuts(track1, track2)) {
continue;
}*/
Expand All @@ -535,8 +487,46 @@ struct UpcPhotonuclearAnalysisJMG {
track2.pt(), track1.pt(),
multiplicity,
deltaPhi,
posZ,
wPair);
posZ);
if constexpr (std::is_same_v<TTag, SameEventTag>) {
if (minMultiplicity <= multiplicity) {
histos.fill(HIST("sameEvent2D"), deltaEta, deltaPhi);
}
if (minMultiplicity <= multiplicity && multiplicity <= range1Max) {
histos.fill(HIST("sameEvent_2_10"), deltaEta, deltaPhi);
}
if (range2Min <= multiplicity && multiplicity <= range2Max) {
histos.fill(HIST("sameEvent_11_20"), deltaEta, deltaPhi);
}
if (range3Min <= multiplicity && multiplicity <= range3Max) {
histos.fill(HIST("sameEvent_21_30"), deltaEta, deltaPhi);
}
if (range4Min <= multiplicity && multiplicity <= range4Max) {
histos.fill(HIST("sameEvent_31_40"), deltaEta, deltaPhi);
}
if (range5Min <= multiplicity && multiplicity <= range5Max) {
histos.fill(HIST("sameEvent_41_50"), deltaEta, deltaPhi);
}
} else if constexpr (std::is_same_v<TTag, MixedEventTag>) {
if (minMultiplicity <= multiplicity) {
histos.fill(HIST("mixedEvent2D"), deltaEta, deltaPhi);
}
if (minMultiplicity <= multiplicity && multiplicity <= range1Max) {
histos.fill(HIST("mixedEvent_2_10"), deltaEta, deltaPhi);
}
if (range2Min <= multiplicity && multiplicity <= range2Max) {
histos.fill(HIST("mixedEvent_11_20"), deltaEta, deltaPhi);
}
if (range3Min <= multiplicity && multiplicity <= range3Max) {
histos.fill(HIST("mixedEvent_21_30"), deltaEta, deltaPhi);
}
if (range4Min <= multiplicity && multiplicity <= range4Max) {
histos.fill(HIST("mixedEvent_31_40"), deltaEta, deltaPhi);
}
if (range5Min <= multiplicity && multiplicity <= range5Max) {
histos.fill(HIST("mixedEvent_41_50"), deltaEta, deltaPhi);
}
}
}
}
}
Expand All @@ -563,7 +553,6 @@ struct UpcPhotonuclearAnalysisJMG {
}
float phiVal = RecoDecay::constrainAngle(phi(track.px(), track.py()), 0.f);
histos.fill(HIST("etaphiVtx"), reconstructedCollision.posZ(), eta(track.px(), track.py(), track.pz()), phiVal);
histos.fill(HIST("Tracks/hTrackPhiBeforeCorr"), phiVal);
}

switch (sgSide) {
Expand Down Expand Up @@ -717,9 +706,6 @@ struct UpcPhotonuclearAnalysisJMG {
// maxCountGapC = histEventCount->GetBinContent(binC) * factorEventsMixed;
// }

auto histoEthaPhiVtxZ = histos.get<TH3>(HIST("etaphiVtx"));
makeNUAWeights(histoEthaPhiVtxZ);

BinningType bindingOnVtx{{vtxBinsEdges, gapSideBinsEdges}, true};
// BinningType bindingOnVtx{{vtxBinsEdges}, true};
auto tracksTuple = std::make_tuple(reconstructedTracks);
Expand Down Expand Up @@ -758,7 +744,7 @@ struct UpcPhotonuclearAnalysisJMG {
histos.fill(HIST("Events/hCountCollisionsMixed"), 2);
// histos.fill(HIST("eventcount"), bindingOnVtx.getBin({collision1.posZ()}));
histos.fill(HIST("eventcount"), bindingOnVtx.getBin({collision1.posZ(), collision1.gapSide()}));
fillCorrelationsUD(mixed, tracks1, tracks2, multiplicity, collision1.posZ());
fillCorrelationsUD(mixed, tracks1, tracks2, multiplicity, collision1.posZ(), MixedEventTag{});
// LOGF(info, "Filling mixed events");

// if (collision1.gapSide() == 0 && collision2.gapSide() == 0) { gap on side A
Expand Down Expand Up @@ -909,10 +895,6 @@ struct UpcPhotonuclearAnalysisJMG {
continue;
}
++multiplicity;

float weightNUA = getNUAWeight(reconstructedCollision.posZ(), eta(track.px(), track.py(), track.pz()), phi(track.px(), track.py()));
float phiVal = RecoDecay::constrainAngle(phi(track.px(), track.py()), 0.f);
histos.fill(HIST("Tracks/hTrackPhiAfterCorr"), phiVal, weightNUA);
}
// multiplicity = reconstructedTracks.size();
if (fillCollisionUD(same, multiplicity) == false) {
Expand All @@ -921,7 +903,7 @@ struct UpcPhotonuclearAnalysisJMG {
// LOGF(debug, "Filling same events");
histos.fill(HIST("eventcount"), -2);
fillQAUD(reconstructedTracks, multiplicity);
fillCorrelationsUD(same, reconstructedTracks, reconstructedTracks, multiplicity, reconstructedCollision.posZ());
fillCorrelationsUD(same, reconstructedTracks, reconstructedTracks, multiplicity, reconstructedCollision.posZ(), SameEventTag{});

/*switch (sgSide) {
case 0: // gap for side A
Expand Down
Loading