[PWGDQ] Modification in dqEfficiency_withAsso task

PWGDQ/Core/HistogramsLibrary.cxx

@@ -889,6 +889,15 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
     hm->AddHistogram(histClass, "Mass", "", false, 500, 0.0, 5.0, VarManager::kMass);
     hm->AddHistogram(histClass, "Eta_Pt", "", false, 40, -2.0, 2.0, VarManager::kEta, 200, 0.0, 20.0, VarManager::kPt);
     hm->AddHistogram(histClass, "Phi_Eta", "#phi vs #eta distribution", false, 200, -5.0, 5.0, VarManager::kEta, 200, -2. * o2::constants::math::PI, 2. * o2::constants::math::PI, VarManager::kPhi);
+    int varspTHE[3] = {VarManager::kMCPt, VarManager::kMCCosThetaHE, VarManager::kMCPhiHE};
+    int varspTCS[3] = {VarManager::kMCPt, VarManager::kMCCosThetaCS, VarManager::kMCPhiCS};
+    int varspTPP[3] = {VarManager::kMCPt, VarManager::kMCCosThetaPP, VarManager::kMCPhiPP};
+    int binspT[3] = {20, 20, 20};
+    double xminpT[3] = {0., -1., -3.14};
+    double xmaxpT[3] = {20., 1., +3.14};
+    hm->AddHistogram(histClass, "Pt_cosThetaHE_phiHE", "", 3, varspTHE, binspT, xminpT, xmaxpT, 0, -1, kFALSE);
+    hm->AddHistogram(histClass, "Pt_cosThetaCS_phiCS", "", 3, varspTCS, binspT, xminpT, xmaxpT, 0, -1, kFALSE);
+    hm->AddHistogram(histClass, "Pt_cosThetaPP_phiPP", "", 3, varspTPP, binspT, xminpT, xmaxpT, 0, -1, kFALSE);
   }
   if (!groupStr.CompareTo("mctruth_quad")) {
     hm->AddHistogram(histClass, "hMass_defaultDileptonMass", "", false, 1000, 3.0, 5.0, VarManager::kQuadDefaultDileptonMass);

PWGDQ/Core/VarManager.cxx

@@ -542,6 +542,18 @@ void VarManager::SetDefaultVarNames()
   fgVariableUnits[kMCVy] = "cm"; // TODO: check the unit
   fgVariableNames[kMCVz] = "MC vz";
   fgVariableUnits[kMCVz] = "cm"; // TODO: check the unit
+  fgVariableNames[kMCCosThetaHE] = "cos#it{#theta}^{MC}_{HE}";
+  fgVariableUnits[kMCCosThetaHE] = "";
+  fgVariableNames[kMCPhiHE] = "#varphi^{MC}_{HE}";
+  fgVariableUnits[kMCPhiHE] = "rad.";
+  fgVariableNames[kMCCosThetaCS] = "cos#it{#theta}^{MC}_{CS}";
+  fgVariableUnits[kMCCosThetaCS] = "";
+  fgVariableNames[kMCPhiCS] = "#varphi^{MC}_{CS}";
+  fgVariableUnits[kMCPhiCS] = "rad.";
+  fgVariableNames[kMCCosThetaPP] = "cos#it{#theta}^{MC}_{PP}";
+  fgVariableUnits[kMCCosThetaPP] = "";
+  fgVariableNames[kMCPhiPP] = "#varphi^{MC}_{PP}";
+  fgVariableUnits[kMCPhiPP] = "rad.";
   fgVariableNames[kCandidateId] = "";
   fgVariableUnits[kCandidateId] = "";
   fgVariableNames[kPairType] = "Pair type";
@@ -1441,6 +1453,12 @@ void VarManager::SetDefaultVarNames()
   fgVarNamesMap["kMCPhi"] = kMCPhi;
   fgVarNamesMap["kMCEta"] = kMCEta;
   fgVarNamesMap["kMCY"] = kMCY;
+  fgVarNamesMap["kMCCosThetaHE"] = kMCCosThetaHE;
+  fgVarNamesMap["kMCCosThetaCS"] = kMCCosThetaCS;
+  fgVarNamesMap["kMCCosThetaPP"] = kMCCosThetaPP;
+  fgVarNamesMap["kMCPhiHE"] = kMCPhiHE;
+  fgVarNamesMap["kMCPhiCS"] = kMCPhiCS;
+  fgVarNamesMap["kMCPhiPP"] = kMCPhiPP;
   fgVarNamesMap["kMCParticleGeneratorId"] = kMCParticleGeneratorId;
   fgVarNamesMap["kNMCParticleVariables"] = kNMCParticleVariables;
   fgVarNamesMap["kMCMotherPdgCode"] = kMCMotherPdgCode;

PWGDQ/Core/VarManager.h

@@ -618,6 +618,14 @@ class VarManager : public TObject
     // MC mother particle variables
     kMCMotherPdgCode,
 
+    // MC pair variables
+    kMCCosThetaHE,
+    kMCCosThetaCS,
+    kMCCosThetaPP,
+    kMCPhiHE,
+    kMCPhiCS,
+    kMCPhiPP,
+
     // Pair variables
     kCandidateId,
     kPairType,
@@ -3242,6 +3250,39 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values)
   values[kMCEta] = v12.Eta();
   values[kMCPhi] = v12.Phi();
   values[kMCY] = -v12.Rapidity();
+  double BeamMomentum = TMath::Sqrt(fgCenterOfMassEnergy * fgCenterOfMassEnergy / 4 - fgMassofCollidingParticle * fgMassofCollidingParticle); // GeV
+  ROOT::Math::PxPyPzEVector Beam1(0., 0., -BeamMomentum, fgCenterOfMassEnergy / 2);
+  ROOT::Math::PxPyPzEVector Beam2(0., 0., BeamMomentum, fgCenterOfMassEnergy / 2);
+
+  // Boost to center of mass frame
+  ROOT::Math::Boost boostv12{v12.BoostToCM()};
+  ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()};
+  ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()};
+  ROOT::Math::XYZVectorF Beam1_CM{(boostv12(Beam1).Vect()).Unit()};
+  ROOT::Math::XYZVectorF Beam2_CM{(boostv12(Beam2).Vect()).Unit()};
+
+  if (fgUsedVars[kMCCosThetaHE] || fgUsedVars[kMCPhiHE]) {
+    ROOT::Math::XYZVectorF zaxis_HE{(v12.Vect()).Unit()};
+    ROOT::Math::XYZVectorF yaxis_HE{(Beam1_CM.Cross(Beam2_CM)).Unit()};
+    ROOT::Math::XYZVectorF xaxis_HE{(yaxis_HE.Cross(zaxis_HE)).Unit()};
+    values[kMCCosThetaHE] = (t1.pdgCode() < 0 ? zaxis_HE.Dot(v1_CM) : zaxis_HE.Dot(v2_CM));
+    values[kMCPhiHE] = (t1.pdgCode() < 0 ? TMath::ATan2(yaxis_HE.Dot(v1_CM), xaxis_HE.Dot(v1_CM)) : TMath::ATan2(yaxis_HE.Dot(v2_CM), xaxis_HE.Dot(v2_CM)));
+  }
+
+  if (fgUsedVars[kMCCosThetaCS] || fgUsedVars[kMCPhiCS]) {
+    ROOT::Math::XYZVectorF zaxis_CS{((Beam1_CM.Unit() - Beam2_CM.Unit()).Unit())};
+    ROOT::Math::XYZVectorF yaxis_CS{(Beam1_CM.Cross(Beam2_CM)).Unit()};
+    ROOT::Math::XYZVectorF xaxis_CS{(yaxis_CS.Cross(zaxis_CS)).Unit()};
+    values[kMCCosThetaCS] = (t1.pdgCode() < 0 ? zaxis_CS.Dot(v1_CM) : zaxis_CS.Dot(v2_CM));
+    values[kMCPhiCS] = (t1.pdgCode() < 0 ? TMath::ATan2(yaxis_CS.Dot(v1_CM), xaxis_CS.Dot(v1_CM)) : TMath::ATan2(yaxis_CS.Dot(v2_CM), xaxis_CS.Dot(v2_CM)));
+  }
+
+  if (fgUsedVars[kMCCosThetaPP] || fgUsedVars[kMCPhiPP]) {
+    ROOT::Math::XYZVectorF zaxis_HE{(v12.Vect()).Unit()};
+    ROOT::Math::XYZVector normalVec = ROOT::Math::XYZVector(v12.Py(), -v12.Px(), 0.f);
+    values[kMCCosThetaPP] = (t1.pdgCode() < 0 ? normalVec.Dot(v1_CM) : normalVec.Dot(v2_CM));
+    values[kMCPhiPP] = (t1.pdgCode() < 0 ? TMath::ATan2((normalVec.Dot(v1_CM)), zaxis_HE.Dot(v1_CM)) : TMath::ATan2((normalVec.Dot(v2_CM)), zaxis_HE.Dot(v2_CM)));
+  }
 }
 
 template <typename T1, typename T2, typename T3>

PWGDQ/Tasks/dqEfficiency_withAssoc.cxx

@@ -1252,7 +1252,7 @@ struct AnalysisSameEventPairing {
     Configurable<std::string> recSignals{"cfgBarrelMCRecSignals", "", "Comma separated list of MC signals (reconstructed)"};
     Configurable<std::string> recSignalsJSON{"cfgMCRecSignalsJSON", "", "Comma separated list of MC signals (reconstructed) via JSON"};
     Configurable<bool> skimSignalOnly{"cfgSkimSignalOnly", false, "Configurable to select only matched candidates"};
-    // Configurable<bool> runMCGenPair{"cfgRunMCGenPair", false, "Do pairing of true MC particles"};
+    Configurable<bool> runMCGenPair{"cfgRunMCGenPair", false, "Do pairing of true MC particles"};
   } fConfigMC;
 
   struct : ConfigurableGroup {
@@ -1981,6 +1981,7 @@ struct AnalysisSameEventPairing {
   // Preslice<ReducedMCTracks> perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId;
   PresliceUnsorted<ReducedMCTracks> perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId;
 
+  template <int TPairType>
   void runMCGen(ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks)
   {
     // loop over mc stack and fill histograms for pure MC truth signals
@@ -2017,29 +2018,27 @@ struct AnalysisSameEventPairing {
     }
 
     if (fHasTwoProngGenMCsignals) {
-      for (auto& event : mcEvents) {
-        auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.globalIndex());
-        groupedMCTracks.bindInternalIndicesTo(&mcTracks);
-        for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) {
-          auto t1_raw = groupedMCTracks.rawIteratorAt(t1.globalIndex());
-          auto t2_raw = groupedMCTracks.rawIteratorAt(t2.globalIndex());
+      for (auto& [t1, t2] : combinations(mcTracks, mcTracks)) {
+        auto t1_raw = mcTracks.rawIteratorAt(t1.globalIndex());
+        auto t2_raw = mcTracks.rawIteratorAt(t2.globalIndex());
+        if (t1_raw.reducedMCeventId() == t2_raw.reducedMCeventId()) {
           for (auto& sig : fGenMCSignals) {
             if (sig->GetNProngs() != 2) { // NOTE: 2-prong signals required here
               continue;
             }
             if (sig->CheckSignal(true, t1_raw, t2_raw)) {
               mcDecision |= (static_cast<uint32_t>(1) << isig);
-              VarManager::FillPairMC<VarManager::kDecayToEE>(t1, t2);
+              VarManager::FillPairMC<TPairType>(t1, t2);
               fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), VarManager::fgValues);
               if (useMiniTree.fConfigMiniTree) {
                 // WARNING! To be checked
-                dileptonMiniTreeGen(mcDecision, event.impactParameter(), t1.pt(), t1.eta(), t1.phi(), t2.pt(), t2.eta(), t2.phi());
+                dileptonMiniTreeGen(mcDecision, -999, t1.pt(), t1.eta(), t1.phi(), t2.pt(), t2.eta(), t2.phi());
               }
             }
             isig++;
-          } // end loop over MC signals
-        } // end loop over pairs
-      } // end loop over events
+          }
+        }
+      }
     }
   } // end runMCGen
 
@@ -2051,9 +2050,9 @@ struct AnalysisSameEventPairing {
     runSameEventPairing<true, VarManager::kDecayToEE, gkEventFillMapWithCov, gkTrackFillMapWithCov>(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks);
     runSameEventPairing<true, VarManager::kDecayToMuMu, gkEventFillMapWithCov, gkMuonFillMapWithCov>(events, muonAssocsPerCollision, muonAssocs, muons, mcEvents, mcTracks);
     // Feature replaced by processMCGen
-    /*if (fConfigMC.runMCGenPair) {
-      runMCGen(mcEvents, mcTracks);
-    }*/
+    if (fConfigMC.runMCGenPair) {
+      runMCGen<VarManager::kDecayToEE>(mcEvents, mcTracks);
+    }
     // runSameEventPairing<true, VarManager::kElectronMuon, gkEventFillMap, gkTrackFillMap>(event, tracks, muons);
   }
 
@@ -2063,9 +2062,9 @@ struct AnalysisSameEventPairing {
   {
     runSameEventPairing<true, VarManager::kDecayToEE, gkEventFillMapWithCov, gkTrackFillMapWithCov>(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks);
     // Feature replaced by processMCGen
-    /*if (fConfigMC.runMCGenPair) {
-      runMCGen(mcEvents, mcTracks);
-    }*/
+    if (fConfigMC.runMCGenPair) {
+      runMCGen<VarManager::kDecayToEE>(mcEvents, mcTracks);
+    }
   }
 
   void processBarrelOnlyWithCollSkimmed(MyEventsVtxCovSelected const& events,
@@ -2074,19 +2073,19 @@ struct AnalysisSameEventPairing {
   {
     runSameEventPairing<true, VarManager::kDecayToEE, gkEventFillMapWithCov, gkTrackFillMapWithCovWithColl>(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks);
     // Feature replaced by processMCGen
-    /*if (fConfigMC.runMCGenPair) {
-      runMCGen(mcEvents, mcTracks);
-    }*/
+    if (fConfigMC.runMCGenPair) {
+      runMCGen<VarManager::kDecayToEE>(mcEvents, mcTracks);
+    }
   }
 
   void processMuonOnlySkimmed(MyEventsVtxCovSelected const& events,
                               soa::Join<aod::ReducedMuonsAssoc, aod::MuonTrackCuts> const& muonAssocs, MyMuonTracksWithCovWithAmbiguities const& muons, ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks)
   {
     runSameEventPairing<true, VarManager::kDecayToMuMu, gkEventFillMapWithCov, gkMuonFillMapWithCov>(events, muonAssocsPerCollision, muonAssocs, muons, mcEvents, mcTracks);
     // Feature replaced by processMCGen
-    /*if (fConfigMC.runMCGenPair) {
-      runMCGen(mcEvents, mcTracks);
-    }*/
+    if (fConfigMC.runMCGenPair) {
+      runMCGen<VarManager::kDecayToMuMu>(mcEvents, mcTracks);
+    }
   }
 
   PresliceUnsorted<ReducedMCTracks> perReducedMcGenEvent = aod::reducedtrackMC::reducedMCeventId;