[PWGDQ] Add e-mu pairing with MC matching to dqEfficiency_withAssoc

PWGDQ/Core/HistogramsLibrary.cxx

@@ -1971,6 +1971,15 @@
     if (subGroupStr.Contains("correlation-emu")) {
       hm->AddHistogram(histClass, "DeltaPhiPair2_DeltaEtaPair2", "", false, 600, -o2::constants::math::PIHalf, 1.5 * o2::constants::math::PI, VarManager::kDeltaPhiPair2, 350, 1.5, 5.0, VarManager::kDeltaEtaPair2);
       hm->AddHistogram(histClass, "DeltaPhiPair2_Pt", "", false, 600, -o2::constants::math::PIHalf, 1.5 * o2::constants::math::PI, VarManager::kDeltaPhiPair2, 200, 0.0, 20.0, VarManager::kPt);
+      // 4D correlation map: (Delta phi_pair, Delta eta_pair, p_T^{e}, p_T^{mu}). Per-track pT axes use kPt1/kPt2,
+      // which FillPair populates for SE and FillPairME (patched in this PR) populates for ME. Stored as THnSparse so
+      // running with many MC-matched hist classes (track cut x muon cut x signal x QA variant) does not blow the
+      // 1 GB TBufferFile limit during the final ROOT serialization.
+      int varsEmu4D[4] = {VarManager::kDeltaPhiPair2, VarManager::kDeltaEtaPair2, VarManager::kPt1, VarManager::kPt2};
+      int binsEmu4D[4] = {60, 35, 20, 20};
+      double xminEmu4D[4] = {-o2::constants::math::PIHalf, 1.5, 0.0, 0.0};
+      double xmaxEmu4D[4] = {1.5 * o2::constants::math::PI, 5.0, 20.0, 20.0};
+      hm->AddHistogram(histClass, "DeltaPhiPair2_DeltaEtaPair2_PtE_PtMu", "", 4, varsEmu4D, binsEmu4D, xminEmu4D, xmaxEmu4D, nullptr, -1, kTRUE);
     }
     if (subGroupStr.Contains("dielectrons")) {
       if (subGroupStr.Contains("prefilter")) {
@@ -2419,7 +2428,7 @@
     LOG(fatal) << "histClass field should be an array of strings, e.g. [class1, class2]";
     return false;
   }
   for (auto& v : hist->FindMember("histClass")->value.GetArray()) {
     if (!v.IsString()) {
       LOG(fatal) << "histClass field should be an array of strings, e.g. [class1, class2]";
       return false;
@@ -2575,7 +2584,7 @@
     }
   }
   if (isTHn) {
     for (auto& v : hist->FindMember("vars")->value.GetArray()) {
       if (VarManager::fgVarNamesMap.find(v.GetString()) == VarManager::fgVarNamesMap.end()) {
         LOG(fatal) << "Bad variable in vars (" << v.GetString() << ") specified for histogram";
         return false;
@@ -2638,7 +2647,7 @@
 
     // create an array of strings to store the different histogram classes
     std::vector<const char*> histClasses;
     for (auto& v : hist.FindMember("histClass")->value.GetArray()) {
       histClasses.push_back(v.GetString());
     }
     const char* title = hist.FindMember("title")->value.GetString();
@@ -2649,7 +2658,7 @@
 
       int* vars = new int[nDimensions];
       int iDim = 0;
       for (auto& v : hist.FindMember("vars")->value.GetArray()) {
         LOG(debug) << "iDim " << iDim << ": " << v.GetString();
         vars[iDim++] = VarManager::fgVarNamesMap[v.GetString()];
       }
@@ -2663,27 +2672,27 @@
         xmin = new double[nDimensions];
         xmax = new double[nDimensions];
         int iDim = 0;
         for (auto& v : hist.FindMember("nBins")->value.GetArray()) {
           nBins[iDim++] = v.GetInt();
           LOG(debug) << "nBins " << iDim << ": " << nBins[iDim - 1];
         }
         iDim = 0;
         for (auto& v : hist.FindMember("xmin")->value.GetArray()) {
           xmin[iDim++] = v.GetDouble();
           LOG(debug) << "xmin " << iDim << ": " << xmin[iDim - 1];
         }
         iDim = 0;
         for (auto& v : hist.FindMember("xmax")->value.GetArray()) {
           xmax[iDim++] = v.GetDouble();
           LOG(debug) << "xmax " << iDim << ": " << xmax[iDim - 1];
         }
       } else {
         int iDim = 0;
         binLimits = new TArrayD[nDimensions];
         for (auto& v : hist.FindMember("binLimits")->value.GetArray()) {
           double* lims = new double[v.GetArray().Size()];
           int iElem = 0;
           for (auto& lim : v.GetArray()) {
             lims[iElem++] = lim.GetDouble();
           }
           binLimits[iDim++] = TArrayD(v.GetArray().Size(), lims);
@@ -2694,7 +2703,7 @@
       if (hist.HasMember("axLabels")) {
         axLabels = new TString[hist.FindMember("axLabels")->value.GetArray().Size()];
         int iDim = 0;
         for (auto& v : hist.FindMember("axLabels")->value.GetArray()) {
           axLabels[iDim++] = v.GetString();
         }
       }

PWGDQ/Core/MCSignalLibrary.cxx

@@ -459,6 +459,30 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name)
     signal = new MCSignal(name, "Electron-muon pair", {electron, muon}, {-1, -1});
     return signal;
   }
+  if (!nameStr.compare("emuFromOpenHFhadron")) {
+    MCProng electron(2, {11, 902}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    electron.SetSourceBit(0, MCProng::kPhysicalPrimary);
+    MCProng muon(2, {13, 902}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    muon.SetSourceBit(0, MCProng::kPhysicalPrimary);
+    signal = new MCSignal(name, "e and mu each from an open charm or beauty hadron decay", {electron, muon}, {-1, -1});
+    return signal;
+  }
+  if (!nameStr.compare("emuFromOpenCharmHadron")) {
+    MCProng electron(2, {11, 402}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    electron.SetSourceBit(0, MCProng::kPhysicalPrimary);
+    MCProng muon(2, {13, 402}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    muon.SetSourceBit(0, MCProng::kPhysicalPrimary);
+    signal = new MCSignal(name, "e and mu each from an open charm hadron decay", {electron, muon}, {-1, -1});
+    return signal;
+  }
+  if (!nameStr.compare("emuFromOpenBeautyHadron")) {
+    MCProng electron(2, {11, 502}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    electron.SetSourceBit(0, MCProng::kPhysicalPrimary);
+    MCProng muon(2, {13, 502}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    muon.SetSourceBit(0, MCProng::kPhysicalPrimary);
+    signal = new MCSignal(name, "e and mu each from an open beauty hadron decay", {electron, muon}, {-1, -1});
+    return signal;
+  }
   if (!nameStr.compare("dielectronFromPC")) {
     MCProng prong(2, {11, 22}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
     signal = new MCSignal(name, "dielectron from a photon conversion", {prong, prong}, {1, 1});

PWGDQ/Core/VarManager.h

@@ -3998,6 +3998,14 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values)
   values[kPhi] = v12.Phi() > 0 ? v12.Phi() : v12.Phi() + 2. * M_PI;
   values[kRap] = -v12.Rapidity();
 
+  // Per-track quantities so ME histograms can use kPt1/kPt2/kEta1/kEta2/kPhi1/kPhi2 just like SE FillPair does.
+  values[kPt1] = t1.pt();
+  values[kEta1] = t1.eta();
+  values[kPhi1] = t1.phi();
+  values[kPt2] = t2.pt();
+  values[kEta2] = t2.eta();
+  values[kPhi2] = t2.phi();
+
   if (fgUsedVars[kDeltaPhiPair2]) {
     double phipair2ME = v1.Phi() - v2.Phi();
     if (phipair2ME > 3 * TMath::Pi() / 2) {