[PWGDQ] Change some code for the energy correlator study

PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx

@@ -215,7 +215,7 @@
     TString eventCutJSONStr = fConfigEventOptions.fConfigEventCutsJSON.value;
     if (eventCutJSONStr != "") {
       std::vector<AnalysisCut*> jsonCuts = dqcuts::GetCutsFromJSON(eventCutJSONStr.Data());
       for (auto& cutIt : jsonCuts) {
         fEventCut->AddCut(cutIt);
       }
     }
@@ -233,7 +233,7 @@
     TString trackCutsJSON = fConfigTrackOptions.fConfigTrackCutsJSON.value;
     if (trackCutsJSON != "") {
       std::vector<AnalysisCut*> addTrackCuts = dqcuts::GetCutsFromJSON(trackCutsJSON.Data());
       for (auto& t : addTrackCuts) {
         fTrackCuts.push_back(reinterpret_cast<AnalysisCompositeCut*>(t));
         fTrackCutNames.push_back(t->GetName());
         trackCutStr += Form(",%s", t->GetName());
@@ -273,7 +273,7 @@
     TString addMCSignalsStr = fConfigPairOptions.recSignalsJSON.value;
     if (addMCSignalsStr != "") {
       std::vector<MCSignal*> addMCSignals = dqmcsignals::GetMCSignalsFromJSON(addMCSignalsStr.Data());
       for (auto& mcIt : addMCSignals) {
         if (mcIt->GetNProngs() != 2) { // NOTE: only 2 prong signals
           continue;
         }
@@ -285,7 +285,7 @@
     TString addMCTripleSignalsStr = fConfigDileptonHadronOptions.fConfigMCRecSignalsJSON.value;
     if (addMCTripleSignalsStr != "") {
       std::vector<MCSignal*> addMCTripleSignals = dqmcsignals::GetMCSignalsFromJSON(addMCTripleSignalsStr.Data());
       for (auto& mcIt : addMCTripleSignals) {
         if (mcIt->GetNProngs() != 3) {
           LOG(fatal) << "Signal at reconstructed level requested (" << mcIt->GetName() << ") " << "does not have 3 prongs! Fix it";
         }
@@ -306,7 +306,7 @@
     TString hadronCutsJSON = fConfigDileptonHadronOptions.fConfigHadronCutsJSON.value;
     if (hadronCutsJSON != "") {
       std::vector<AnalysisCut*> addHadronCuts = dqcuts::GetCutsFromJSON(hadronCutsJSON.Data());
       for (auto& t : addHadronCuts) {
         fHadronCuts.push_back(reinterpret_cast<AnalysisCompositeCut*>(t));
         fHadronCutNames.push_back(t->GetName());
         hadronCutStr += Form(",%s", t->GetName());
@@ -328,7 +328,7 @@
     TString addMCSignalsGenStr = fConfigDileptonHadronOptions.fConfigMCGenSignalsJSON.value;
     if (addMCSignalsGenStr != "") {
       std::vector<MCSignal*> addMCSignals = dqmcsignals::GetMCSignalsFromJSON(addMCSignalsGenStr.Data());
       for (auto& mcIt : addMCSignals) {
         if (mcIt->GetNProngs() > 2) { // NOTE: only 2 prong signals
           continue;
         }
@@ -350,7 +350,7 @@
     TString addRecTrackSignalsGenStr = fConfigTrackOptions.fConfigMCRecTrackSignalsJSON.value;
     if (addRecTrackSignalsGenStr != "") {
       std::vector<MCSignal*> addMCRecTrackSignals = dqmcsignals::GetMCSignalsFromJSON(addRecTrackSignalsGenStr.Data());
       for (auto& mcIt : addMCRecTrackSignals) {
         if (mcIt->GetNProngs() > 2) { // NOTE: only 2 prong signals
           continue;
         }
@@ -374,11 +374,11 @@
       TString histClasses = "AssocsBarrel_BeforeCuts;";
       // Configure histogram classes for each track cut;
       // Add histogram classes for each track cut and for each requested MC signal (reconstructed tracks with MC truth)
       for (auto& cut : fTrackCuts) {
         TString nameStr = Form("AssocsBarrel_%s", cut->GetName());
         fHistNamesReco.push_back(nameStr);
         histClasses += Form("%s;", nameStr.Data());
         for (auto& sig : fRecMCTrackSignals) {
           TString nameStr2 = Form("AssocsBarrelMatched_%s_%s", cut->GetName(), sig->GetName());
           fHistNamesMCMatched.push_back(nameStr2);
           histClasses += Form("%s;", nameStr2.Data());
@@ -401,7 +401,7 @@
           Form("PairsBarrelSEPM_%s", objArray->At(icut)->GetName()),
           Form("PairsBarrelSEPP_%s", objArray->At(icut)->GetName()),
           Form("PairsBarrelSEMM_%s", objArray->At(icut)->GetName())};
         for (auto& n : names) {
           histNames += Form("%s;", n.Data());
         }
         fTrackHistNames[icut] = names;
@@ -475,6 +475,23 @@
       LOG(fatal) << "Problem getting histograms from the TList object with efficiencies!";
     }
   }
+
+  float GetSafeInterpolationWeight(TH2* hEff, float x, float y)
+  {
+    if (!hEff)
+      return 1.0;
+    float minX = hEff->GetXaxis()->GetBinCenter(1);
+    float maxX = hEff->GetXaxis()->GetBinCenter(hEff->GetXaxis()->GetNbins());
+
+    float minY = hEff->GetYaxis()->GetBinCenter(1);
+    float maxY = hEff->GetYaxis()->GetBinCenter(hEff->GetYaxis()->GetNbins());
+
+    float safeX = std::max(minX, std::min(x, maxX));
+    float safeY = std::max(minY, std::min(y, maxY));
+
+    return hEff->Interpolate(safeX, safeY);
+  }
+
   template <bool MixedEvent, uint32_t TTrackFillMap, typename TTrack1, typename TTrack2, typename THadron, typename TEvent>
   void runDileptonHadron(TTrack1 const& track1, TTrack2 const& track2, int iEleCut,
                          THadron const& hadron, TEvent const& event, aod::McParticles const& /*mcParticles*/)
@@ -503,16 +520,18 @@
     float Effweight_rec = 1.0f;
     float Accweight_gen = 1.0f;
     if (fConfigDileptonHadronOptions.fConfigApplyEfficiency) {
+      float dilepton_pt = VarManager::fgValues[VarManager::kPt];
       float dilepton_eta = VarManager::fgValues[VarManager::kEta];
       float dilepton_phi = VarManager::fgValues[VarManager::kPhi];
+      float dilepton_rap = VarManager::fgValues[VarManager::kRap];
       float hadron_eta = hadron.eta();
       float hadron_phi = hadron.phi();
       float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI);
-      Effweight_rec = hAcceptance_rec->Interpolate(dilepton_eta - hadron_eta, deltaphi);
-      Accweight_gen = hAcceptance_gen->Interpolate(dilepton_eta - hadron_eta, deltaphi);
-      float Effdilepton = hEfficiency_dilepton->Interpolate(VarManager::fgValues[VarManager::kPt], dilepton_eta);
-      float Masswindow = hMasswindow->Interpolate(VarManager::fgValues[VarManager::kPt]);
-      float Effhadron = hEfficiency_hadron->Interpolate(hadron.pt(), hadron.eta());
+      float Effweight_rec = GetSafeInterpolationWeight(hAcceptance_rec, dilepton_eta - hadron_eta, deltaphi);
+      float Accweight_gen = GetSafeInterpolationWeight(hAcceptance_gen, dilepton_eta - hadron_eta, deltaphi);
+      float Effdilepton = GetSafeInterpolationWeight(hEfficiency_dilepton, dilepton_rap, dilepton_pt);
+      float Effhadron = GetSafeInterpolationWeight(hEfficiency_hadron, hadron.pt(), hadron_eta);
+      float Masswindow = hMasswindow->Interpolate(dilepton_pt);
       Accweight_gen = Accweight_gen * Effdilepton * Effhadron;
       if (fConfigDileptonHadronOptions.fConfigApplyEfficiencyME) {
         Effweight_rec = Effdilepton * Effhadron * Masswindow; // for the moment, apply the efficiency correction also for the mixed event pairs, but this can be changed in case we want to apply it only for the same event pairs

PWGDQ/Tasks/tableReader_withAssoc.cxx

@@ -3770,6 +3770,22 @@ struct AnalysisDileptonTrack {
     }
   }
 
+  float GetSafeInterpolationWeight(TH2* hEff, float x, float y)
+  {
+    if (!hEff)
+      return 1.0;
+    float minX = hEff->GetXaxis()->GetBinCenter(1);
+    float maxX = hEff->GetXaxis()->GetBinCenter(hEff->GetXaxis()->GetNbins());
+
+    float minY = hEff->GetYaxis()->GetBinCenter(1);
+    float maxY = hEff->GetYaxis()->GetBinCenter(hEff->GetYaxis()->GetNbins());
+
+    float safeX = std::max(minX, std::min(x, maxX));
+    float safeY = std::max(minY, std::min(y, maxY));
+
+    return hEff->Interpolate(safeX, safeY);
+  }
+
   // Template function to run pair - hadron combinations
   template <int TCandidateType, uint32_t TEventFillMap, uint32_t TTrackFillMap, typename TEvent, typename TTracks, typename TTrackAssocs, typename TDileptons>
   void runDileptonHadron(TEvent const& event, TTrackAssocs const& assocs, TTracks const& tracks, TDileptons const& dileptons)
@@ -3853,10 +3869,10 @@ struct AnalysisDileptonTrack {
             float hadron_eta = track.eta();
             float hadron_phi = track.phi();
             float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI);
-            Effweight_rec = hAcceptance_rec->Interpolate(dilepton_eta - hadron_eta, deltaphi);
-            float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt(), dilepton_eta);
+            Effweight_rec = GetSafeInterpolationWeight(hAcceptance_rec, dilepton_eta - hadron_eta, deltaphi);
+            float Effdilepton = GetSafeInterpolationWeight(hEfficiency_dilepton, dilepton.rap(), dilepton.pt());
+            float Effhadron = GetSafeInterpolationWeight(hEfficiency_hadron, track.eta(), track.pt());
             float Masswindow = hMasswindow->Interpolate(dilepton.pt());
-            float Effhadron = hEfficiency_hadron->Interpolate(track.pt(), hadron_eta);
             Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow;
           }
           std::vector<float> fTransRange = fConfigTransRange;
@@ -4087,10 +4103,10 @@ struct AnalysisDileptonTrack {
             float hadron_eta = track.eta();
             float hadron_phi = track.phi();
             float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI);
-            Effweight_rec = hAcceptance_rec->Interpolate(dilepton_eta - hadron_eta, deltaphi);
-            float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt(), dilepton_eta);
+            Effweight_rec = GetSafeInterpolationWeight(hAcceptance_rec, dilepton_eta - hadron_eta, deltaphi);
+            float Effdilepton = GetSafeInterpolationWeight(hEfficiency_dilepton, dilepton.rap(), dilepton.pt());
+            float Effhadron = GetSafeInterpolationWeight(hEfficiency_hadron, track.eta(), track.pt());
             float Masswindow = hMasswindow->Interpolate(dilepton.pt());
-            float Effhadron = hEfficiency_hadron->Interpolate(track.pt(), hadron_eta);
             if (fConfigApplyEfficiencyME) {
               Effweight_rec = Effdilepton * Effhadron * Masswindow; // for the moment, apply the efficiency correction also for the mixed event pairs, but this can be changed in case we want to apply it only for the same event pairs
             } else {