taking off the NUA correction and adding new histos for correlation

JosueMtzGar · JosueMtzGar · commit 8341cc5f0ccf · 2025-10-29T15:19:54.000-06:00
diff --git a/PWGUD/Tasks/upcPhotonuclearAnalysisJMG.cxx b/PWGUD/Tasks/upcPhotonuclearAnalysisJMG.cxx
@@ -149,7 +149,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"};
@@ -198,7 +198,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"}}});
@@ -215,8 +226,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});
@@ -298,11 +307,15 @@ struct UpcPhotonuclearAnalysisJMG {
     // mixedGapSideA.setObject(new CorrelationContainer("mixedEventGapSideA", "mixedEventGapSideA", corrAxis, effAxis, {}));
     // sameGapSideC.setObject(new CorrelationContainer("sameEventGapSideC", "sameEventGapSideC", corrAxis, effAxis, {}));
     // mixedGapSideC.setObject(new CorrelationContainer("mixedEventGapSideC", "mixedEventGapSideC", corrAxis, effAxis, {}));
+
   }
 
   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;
@@ -431,84 +444,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) {
@@ -517,13 +462,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;
         }*/
@@ -535,8 +476,47 @@ struct UpcPhotonuclearAnalysisJMG {
                                     track2.pt(), track1.pt(),
                                     multiplicity,
                                     deltaPhi,
-                                    posZ,
-                                    wPair);
+                                    posZ);
+        if constexpr (std::is_same_v<TTag, SameEventTag>) {
+          if (2 <= multiplicity) {
+            histos.fill(HIST("sameEvent2D"), deltaEta, deltaPhi);
+          }
+          if (2 <= multiplicity && multiplicity <= 10) {
+            histos.fill(HIST("sameEvent_2_10"), deltaEta, deltaPhi);
+          }
+          if (11 <= multiplicity && multiplicity <= 20) {
+            histos.fill(HIST("sameEvent_11_20"), deltaEta, deltaPhi);
+          }
+          if (21 <= multiplicity && multiplicity <= 30) {
+            histos.fill(HIST("sameEvent_21_30"), deltaEta, deltaPhi);
+          }
+          if (31 <= multiplicity && multiplicity <= 40) {
+            histos.fill(HIST("sameEvent_31_40"), deltaEta, deltaPhi);
+          }
+          if (41 <= multiplicity && multiplicity <= 50) {
+            histos.fill(HIST("sameEvent_41_50"), deltaEta, deltaPhi);
+          }
+        }
+        else if constexpr (std::is_same_v<TTag, MixedEventTag>) {
+          if (2 <= multiplicity) {
+            histos.fill(HIST("mixedEvent2D"), deltaEta, deltaPhi);
+          }
+          if (2 <= multiplicity && multiplicity <= 10) {
+            histos.fill(HIST("mixedEvent_2_10"), deltaEta, deltaPhi);
+          }
+          if (11 <= multiplicity && multiplicity <= 20) {
+            histos.fill(HIST("mixedEvent_11_20"), deltaEta, deltaPhi);
+          }
+          if (21 <= multiplicity && multiplicity <= 30) {
+            histos.fill(HIST("mixedEvent_21_30"), deltaEta, deltaPhi);
+          }
+          if (31 <= multiplicity && multiplicity <= 40) {
+            histos.fill(HIST("mixedEvent_31_40"), deltaEta, deltaPhi);
+          }
+          if (41 <= multiplicity && multiplicity <= 50) {
+            histos.fill(HIST("mixedEvent_41_50"), deltaEta, deltaPhi);
+          }
+        }
       }
     }
   }
@@ -563,7 +543,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) {
@@ -717,9 +696,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);
@@ -758,7 +734,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
@@ -909,10 +885,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) {
@@ -921,7 +893,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
