[PWGDQ] Add and fix some code for the energy correlator study

PWGDQ/Core/VarManager.cxx

@@ -191,8 +191,8 @@
   // TO Do: add more systems
 
   // set the beam 4-momentum vectors
   float beamAEnergy = energy / 2.0 * sqrt(NumberOfProtonsA * NumberOfProtonsC / NumberOfProtonsC / NumberOfProtonsA); // GeV
   float beamCEnergy = energy / 2.0 * sqrt(NumberOfProtonsC * NumberOfProtonsA / NumberOfProtonsA / NumberOfProtonsC); // GeV
   float beamAMomentum = std::sqrt(beamAEnergy * beamAEnergy - NumberOfNucleonsA * NumberOfNucleonsA * MassProton * MassProton);
   float beamCMomentum = std::sqrt(beamCEnergy * beamCEnergy - NumberOfNucleonsC * NumberOfNucleonsC * MassProton * MassProton);
   fgBeamA.SetPxPyPzE(0, 0, beamAMomentum, beamAEnergy);
@@ -291,7 +291,7 @@
     double mean = calibMeanHist->GetBinContent(binTPCncls, binPin, binEta);
     double sigma = calibSigmaHist->GetBinContent(binTPCncls, binPin, binEta);
     return (nSigmaValue - mean) / sigma; // Return the calibrated nSigma value
   } else if (fgCalibrationType == 2) {
     // get the calibration histograms
     CalibObjects calibMean, calibSigma, calibStatus;
     switch (species) {
@@ -386,14 +386,14 @@
   // Bimodality coefficient = (skewness^2 + 1) / kurtosis
   // return a tuple including the coefficient, mean, RMS, skewness, and kurtosis
   size_t n = data.size();
   if (n < 3) {
     return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0);
   }
   float mean = std::accumulate(data.begin(), data.end(), 0.0) / n;
 
   float m2 = 0.0, m3 = 0.0, m4 = 0.0;
   float diff, diff2;
   for (float x : data) {
     diff = x - mean;
     diff2 = diff * diff;
     m2 += diff2;
@@ -432,7 +432,7 @@
   int nBins = static_cast<int>((max - min) / binWidth);
   std::vector<int> counts(nBins, 0.0);
 
   for (float x : data) {
     if (x < min || x >= max) {
       continue; // skip out-of-range values
     }
@@ -2253,6 +2253,7 @@
   fgVarNamesMap["kMCCosChi"] = kMCCosChi;
   fgVarNamesMap["kMCHadronPt"] = kMCHadronPt;
   fgVarNamesMap["kMCWeight_before"] = kMCWeight_before;
+  fgVarNamesMap["kMCEWeight_before"] = kMCEWeight_before;
   fgVarNamesMap["kMCdeltaeta"] = kMCdeltaeta;
   fgVarNamesMap["kMCHadronPt"] = kMCHadronPt;
   fgVarNamesMap["kMCHadronEta"] = kMCHadronEta;
@@ -2267,6 +2268,7 @@
   fgVarNamesMap["kMCdeltaphi_randomPhi_toward"] = kMCdeltaphi_randomPhi_toward;
   fgVarNamesMap["kMCdeltaphi_randomPhi_away"] = kMCdeltaphi_randomPhi_away;
   fgVarNamesMap["kMCdeltaphi_randomPhi_trans"] = kMCdeltaphi_randomPhi_trans;
+  fgVarNamesMap["kMCHadronpt_randomPhi_trans"] = kMCHadronpt_randomPhi_trans;
   fgVarNamesMap["kMCCosChi_gen"] = kMCCosChi_gen;
   fgVarNamesMap["kMCWeight_gen"] = kMCWeight_gen;
   fgVarNamesMap["kMCdeltaeta_gen"] = kMCdeltaeta_gen;
@@ -2504,6 +2506,7 @@
   fgVarNamesMap["kWeight"] = kWeight;
   fgVarNamesMap["kECWeight"] = kECWeight;
   fgVarNamesMap["kEWeight_before"] = kEWeight_before;
+  fgVarNamesMap["kWeight_before"] = kWeight_before;
   fgVarNamesMap["kPtDau"] = kPtDau;
   fgVarNamesMap["kEtaDau"] = kEtaDau;
   fgVarNamesMap["kPhiDau"] = kPhiDau;
@@ -2516,6 +2519,7 @@
   fgVarNamesMap["kdeltaphi_randomPhi_trans"] = kdeltaphi_randomPhi_trans;
   fgVarNamesMap["kdeltaphi_randomPhi_toward"] = kdeltaphi_randomPhi_toward;
   fgVarNamesMap["kdeltaphi_randomPhi_away"] = kdeltaphi_randomPhi_away;
+  fgVarNamesMap["kPtDau_randomPhi_trans"] = kPtDau_randomPhi_trans;
   fgVarNamesMap["kdileptonmass"] = kdileptonmass;
   fgVarNamesMap["kNCorrelationVariables"] = kNCorrelationVariables;
   fgVarNamesMap["kQuadMass"] = kQuadMass;

PWGDQ/Core/VarManager.h

@@ -692,7 +692,9 @@
     kMCdeltaphi_randomPhi_toward,
     kMCdeltaphi_randomPhi_away,
     kMCdeltaphi_randomPhi_trans,
+    kMCHadronpt_randomPhi_trans,
     kMCWeight_before,
+    kMCEWeight_before,
     kMCCosChi_gen,
     kMCWeight_gen,
     kMCdeltaeta_gen,
@@ -945,6 +947,7 @@
     kPtDau,
     kCosTheta,
     kEWeight_before,
+    kWeight_before,
     kCosChi_randomPhi_trans,
     kCosChi_randomPhi_toward,
     kCosChi_randomPhi_away,
@@ -955,6 +958,7 @@
     kdeltaphi_randomPhi_toward,
     kdeltaphi_randomPhi_away,
     kdileptonmass,
+    kPtDau_randomPhi_trans,
 
     // Dilepton-track-track variables
     kQuadMass,
@@ -1804,9 +1808,9 @@
   }
   if constexpr ((fillMap & MuonCov) > 0 || (fillMap & ReducedMuonCov) > 0) {
     o2::dataformats::GlobalFwdTrack propmuon = PropagateMuon(muontrack, collision);
     double px = propmuon.getP() * sin(M_PI / 2 - atan(mfttrack.tgl())) * cos(mfttrack.phi());
     double py = propmuon.getP() * sin(M_PI / 2 - atan(mfttrack.tgl())) * sin(mfttrack.phi());
     double pz = propmuon.getP() * cos(M_PI / 2 - atan(mfttrack.tgl()));
     double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
     auto mftprop = o2::aod::fwdtrackutils::getTrackParCovFwdShift(mfttrack, fgzShiftFwd);
     values[kX] = mftprop.getX();
@@ -3370,6 +3374,7 @@
   values[kMCAccweight] = Accweight;
   values[kMCCosChi] = CosChi;
   values[kMCWeight_before] = t1.pt() / o2::constants::physics::MassJPsi * Accweight;
+  values[kMCEWeight_before] = t1.e() / o2::constants::physics::MassJPsi * Accweight;
   values[kMCCosTheta] = CosTheta;
   values[kMCdeltaphi] = deltaphi;
   values[kMCdeltaeta] = deltaeta;
@@ -3397,6 +3402,7 @@
     randomPhi_toward = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
     randomPhi_away = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
 
+    values[kMCHadronpt_randomPhi_trans] = v2.pt();
     ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans(v2.pt(), v2.eta(), randomPhi_trans, MassHadron);
     values[kMCCosChi_randomPhi_trans] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_randomPhi_trans);
     values[kMCWeight_randomPhi_trans] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_randomPhi_trans) / v1.M() * Accweight;
@@ -4594,7 +4600,7 @@
         values[kVertexingLxyErr] = (KFPV.GetCovariance(0) + KFGeoTwoProng.GetCovariance(0)) * dxPair2PV * dxPair2PV + (KFPV.GetCovariance(2) + KFGeoTwoProng.GetCovariance(2)) * dyPair2PV * dyPair2PV + 2 * ((KFPV.GetCovariance(1) + KFGeoTwoProng.GetCovariance(1)) * dxPair2PV * dyPair2PV);
         values[kVertexingLzErr] = (KFPV.GetCovariance(5) + KFGeoTwoProng.GetCovariance(5)) * dzPair2PV * dzPair2PV;
         values[kVertexingLxyzErr] = (KFPV.GetCovariance(0) + KFGeoTwoProng.GetCovariance(0)) * dxPair2PV * dxPair2PV + (KFPV.GetCovariance(2) + KFGeoTwoProng.GetCovariance(2)) * dyPair2PV * dyPair2PV + (KFPV.GetCovariance(5) + KFGeoTwoProng.GetCovariance(5)) * dzPair2PV * dzPair2PV + 2 * ((KFPV.GetCovariance(1) + KFGeoTwoProng.GetCovariance(1)) * dxPair2PV * dyPair2PV + (KFPV.GetCovariance(3) + KFGeoTwoProng.GetCovariance(3)) * dxPair2PV * dzPair2PV + (KFPV.GetCovariance(4) + KFGeoTwoProng.GetCovariance(4)) * dyPair2PV * dzPair2PV);
         if (fabs(values[kVertexingLxy]) < 1.e-8f)
           values[kVertexingLxy] = 1.e-8f;
         values[kVertexingLxyErr] = values[kVertexingLxyErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyErr]) / values[kVertexingLxy];
         if (fabs(values[kVertexingLz]) < 1.e-8f)
@@ -5877,7 +5883,8 @@
     values[kWeight] = weight;
     values[kECWeight] = E_boost / v1.M() * weight;
     values[kCosTheta] = LorentzTransformJpsihadroncosChi("costheta", v1, v2);
-    values[kEWeight_before] = v2.Pt() / v1.M() * weight;
+    values[kEWeight_before] = v2.E() / v1.M() * weight;
+    values[kWeight_before] = v2.Pt() / v1.M() * weight;
     values[kPtDau] = v2.pt();
     values[kEtaDau] = v2.eta();
     values[kPhiDau] = RecoDecay::constrainAngle(v2.phi(), -o2::constants::math::PIHalf);
@@ -5901,7 +5908,7 @@
       randomPhi_trans = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
       randomPhi_toward = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
       randomPhi_away = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
-
+      values[kPtDau_randomPhi_trans] = v2.pt();
       ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans(v2.pt(), v2.eta(), randomPhi_trans, o2::constants::physics::MassPionCharged);
       values[kCosChi_randomPhi_trans] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_randomPhi_trans);
       values[kWeight_randomPhi_trans] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_randomPhi_trans) / v1.M() * weight;
@@ -5969,8 +5976,7 @@
       ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_rec(v2_rec.pt(), v2_rec.eta(), randomPhi_trans_rec, o2::constants::physics::MassPionCharged);
       values[kMCCosChi_randomPhi_trans_rec] = LorentzTransformJpsihadroncosChi("coschi", v1_rec, v2_randomPhi_trans_rec);
       values[kMCWeight_randomPhi_trans_rec] = LorentzTransformJpsihadroncosChi("weight_boost", v1_rec, v2_randomPhi_trans_rec) / v1_rec.M() * Effweight_rec;
-      float randomPhi_trans_gen = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
-      ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_gen(v2_gen.pt(), v2_gen.eta(), randomPhi_trans_gen, MassHadron);
+      ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_gen(v2_gen.pt(), v2_gen.eta(), randomPhi_trans_rec, MassHadron);
       values[kMCCosChi_randomPhi_trans_gen] = LorentzTransformJpsihadroncosChi("coschi", v1_gen, v2_randomPhi_trans_gen);
       values[kMCWeight_randomPhi_trans_gen] = LorentzTransformJpsihadroncosChi("weight_boost", v1_gen, v2_randomPhi_trans_gen) / v1_gen.M() * Accweight_gen;
     }

PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx

@@ -104,6 +104,8 @@ struct AnalysisEnergyCorrelator {
     Configurable<std::string> fConfigTrackCuts{"cfgTrackCuts", "electronSelection1_ionut", "Comma separated list of barrel track cuts for electrons"};
     Configurable<std::string> fConfigTrackCutsJSON{"cfgTrackCutsJSON", "", "Additional track cuts in JSON"};
     Configurable<std::string> fConfigAddTrackHistogram{"cfgAddTrackHistogram", "", "Track histograms"};
+    Configurable<std::string> fConfigMCRecTrackSignals{"cfgMCRecTrackSignals", "", "Comma separated list of MC signals (reconstructed)"};
+    Configurable<std::string> fConfigMCRecTrackSignalsJSON{"cfgMCRecTrackSignalsJSON", "", "Additional list of MC signals (reconstructed) via JSON"};
     Configurable<bool> fConfigTrackQA{"cfgTrackQA", false, "If true, fill Track QA histograms"};
   } fConfigTrackOptions;
 
@@ -155,12 +157,14 @@ struct AnalysisEnergyCorrelator {
   std::vector<TString> fTrackCutNames;
   std::vector<TString> fHadronCutNames;
   std::vector<TString> fHistNamesReco;
+  std::vector<TString> fHistNamesMCMatched;
 
   std::map<int, std::vector<TString>> fTrackHistNames;
   std::map<int, std::vector<TString>> fBarrelHistNamesMCmatched;
   std::map<int64_t, bool> fSelMap;
 
-  std::vector<MCSignal*> fRecMCSignals; // MC signals for reconstructed pairs
+  std::vector<MCSignal*> fRecMCTrackSignals; // MC signals for reconstructed tracks
+  std::vector<MCSignal*> fRecMCSignals;      // MC signals for reconstructed pairs
   std::vector<MCSignal*> fGenMCSignals;
   std::vector<MCSignal*> fRecMCTripleSignals; // MC signals for reconstructed triples
 
@@ -175,8 +179,8 @@ struct AnalysisEnergyCorrelator {
 
   TH2F* hAcceptance_rec;
   TH2F* hAcceptance_gen;
-  TH1F* hEfficiency_dilepton;
-  TH1F* hEfficiency_hadron;
+  TH2F* hEfficiency_dilepton;
+  TH2F* hEfficiency_hadron;
   TH1F* hMasswindow;
 
   void init(o2::framework::InitContext& context)
@@ -332,6 +336,28 @@ struct AnalysisEnergyCorrelator {
       }
     }
 
+    // TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function
+    TString sigRecTrackNamesStr = fConfigTrackOptions.fConfigMCRecTrackSignals.value;
+    std::unique_ptr<TObjArray> objRecTrackSigArray(sigRecTrackNamesStr.Tokenize(","));
+    for (int isig = 0; isig < objRecTrackSigArray->GetEntries(); isig++) {
+      MCSignal* sig_RecTrack = o2::aod::dqmcsignals::GetMCSignal(objRecTrackSigArray->At(isig)->GetName());
+      if (sig_RecTrack) {
+        fRecMCTrackSignals.push_back(sig_RecTrack);
+      }
+    }
+
+    // Add the MCSignals from the JSON config
+    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;
+        }
+        fRecMCTrackSignals.push_back(mcIt);
+      }
+    }
+
     VarManager::SetUseVars(AnalysisCut::fgUsedVars);
 
     fHistMan = new HistogramManager("analysisHistos", "", VarManager::kNVars);
@@ -352,6 +378,11 @@ struct AnalysisEnergyCorrelator {
         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());
+        }
       }
       DefineHistograms(fHistMan, histClasses.Data(), fConfigTrackOptions.fConfigAddTrackHistogram.value.data());
     }
@@ -435,8 +466,8 @@ struct AnalysisEnergyCorrelator {
     if (!listAccs) {
       LOG(fatal) << "Problem getting TList object with efficiencies!";
     }
-    hEfficiency_dilepton = static_cast<TH1F*>(listAccs->FindObject("hEfficiency_dilepton"));
-    hEfficiency_hadron = static_cast<TH1F*>(listAccs->FindObject("hEfficiency_hadron"));
+    hEfficiency_dilepton = static_cast<TH2F*>(listAccs->FindObject("hEfficiency_dilepton"));
+    hEfficiency_hadron = static_cast<TH2F*>(listAccs->FindObject("hEfficiency_hadron"));
     hAcceptance_rec = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_rec"));
     hAcceptance_gen = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_gen"));
     hMasswindow = static_cast<TH1F*>(listAccs->FindObject("hMasswindow"));
@@ -479,9 +510,9 @@ struct AnalysisEnergyCorrelator {
       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]);
+      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());
+      float Effhadron = hEfficiency_hadron->Interpolate(hadron.pt(), hadron.eta());
       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
@@ -647,17 +678,32 @@ struct AnalysisEnergyCorrelator {
           fHistMan->FillHistClass("AssocsBarrel_BeforeCuts", VarManager::fgValues);
         }
 
+        uint32_t mcDecision = static_cast<uint32_t>(0);
+        // run MC matching for this pair
+        int isig = 0;
+        mcDecision = 0;
+        for (auto sig = fRecMCTrackSignals.begin(); sig != fRecMCTrackSignals.end(); sig++, isig++) {
+          if (t1.has_mcParticle()) {
+            if ((*sig)->CheckSignal(true, t1.mcParticle())) {
+              mcDecision |= (static_cast<uint32_t>(1) << isig);
+            }
+          }
+        }
         // Apply electron cuts and fill histograms
         int iCut1 = 0;
         for (auto cut1 = fTrackCuts.begin(); cut1 != fTrackCuts.end(); cut1++, iCut1++) {
           if ((*cut1)->IsSelected(VarManager::fgValues)) {
             filter1 |= (static_cast<uint32_t>(1) << iCut1);
             if (fConfigTrackOptions.fConfigTrackQA) {
               fHistMan->FillHistClass(fHistNamesReco[iCut1], VarManager::fgValues);
+              for (size_t isig = 0; isig < fRecMCTrackSignals.size(); isig++) { // loop over MC signals
+                if (mcDecision & (static_cast<uint32_t>(1) << isig)) {
+                  fHistMan->FillHistClass(fHistNamesMCMatched[iCut1 * fRecMCTrackSignals.size() + isig], VarManager::fgValues); // matched signal
+                }
+              }
             }
           }
         }
-
         // Check opposite charge with t2
         for (auto& a2 : groupedAssocs) {
           auto t2 = a2.template track_as<MyBarrelTracksWithCov>();

PWGDQ/Tasks/tableReader_withAssoc.cxx

@@ -3493,8 +3493,8 @@ struct AnalysisDileptonTrack {
 
   TH2F* hAcceptance_rec;
   TH2F* hAcceptance_gen;
-  TH1F* hEfficiency_dilepton;
-  TH1F* hEfficiency_hadron;
+  TH2F* hEfficiency_dilepton;
+  TH2F* hEfficiency_hadron;
   TH1F* hMasswindow;
 
   void init(o2::framework::InitContext& context)
@@ -3760,8 +3760,8 @@ struct AnalysisDileptonTrack {
     if (!listAccs) {
       LOG(fatal) << "Problem getting TList object with efficiencies!";
     }
-    hEfficiency_dilepton = static_cast<TH1F*>(listAccs->FindObject("hEfficiency_dilepton"));
-    hEfficiency_hadron = static_cast<TH1F*>(listAccs->FindObject("hEfficiency_hadron"));
+    hEfficiency_dilepton = static_cast<TH2F*>(listAccs->FindObject("hEfficiency_dilepton"));
+    hEfficiency_hadron = static_cast<TH2F*>(listAccs->FindObject("hEfficiency_hadron"));
     hAcceptance_rec = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_rec"));
     hAcceptance_gen = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_gen"));
     hMasswindow = static_cast<TH1F*>(listAccs->FindObject("hMasswindow"));
@@ -3854,9 +3854,9 @@ struct AnalysisDileptonTrack {
             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());
+            float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt(), dilepton_eta);
             float Masswindow = hMasswindow->Interpolate(dilepton.pt());
-            float Effhadron = hEfficiency_hadron->Interpolate(track.pt());
+            float Effhadron = hEfficiency_hadron->Interpolate(track.pt(), hadron_eta);
             Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow;
           }
           std::vector<float> fTransRange = fConfigTransRange;
@@ -4088,17 +4088,17 @@ struct AnalysisDileptonTrack {
             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());
+            float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt(), dilepton_eta);
             float Masswindow = hMasswindow->Interpolate(dilepton.pt());
-            float Effhadron = hEfficiency_hadron->Interpolate(track.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 {
               Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow; // apply acceptance and efficiency correction for the real pairs
             }
           }
           std::vector<float> fTransRange = fConfigTransRange;
-          VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, fValuesHadron, fTransRange[0], fTransRange[1], fConfigApplyMassEC, fMassBkg->GetRandom(), 1. / Effweight_rec);
+          VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, VarManager::fgValues, fTransRange[0], fTransRange[1], fConfigApplyMassEC, fMassBkg->GetRandom(), 1. / Effweight_rec);
 
           // loop over dilepton leg cuts and track cuts and fill histograms separately for each combination
           for (int icut = 0; icut < fNCuts; icut++) {