[PWGJE/EMCal,PWGEM/PhotonMeson] Add true pi0 from same gamma histogram and opening angle cut (#9815)

Co-authored-by: Nicolas Strangmann <nicolas.strangmann@.cern.ch>

Author: nstrangm

PWGEM/Dilepton/Utils/MCUtilities.h

@@ -222,6 +222,14 @@ int FindCommonMotherFrom3Prongs(TMCParticle1 const& p1, TMCParticle2 const& p2,
 }
 //_______________________________________________________________________
 template <typename TMCParticle, typename TMCParticles>
+int getMotherPDGCode(TMCParticle const& p, TMCParticles const& mcparticles)
+{
+  int motherid = p.mothersIds()[0];
+  auto mother = mcparticles.iteratorAt(motherid);
+  return (mother.pdgCode());
+}
+//_______________________________________________________________________
+template <typename TMCParticle, typename TMCParticles>
 int IsFromBeauty(TMCParticle const& p, TMCParticles const& mcparticles)
 {
   if (!p.has_mothers()) {

PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h

@@ -639,6 +639,13 @@ struct Pi0EtaToGammaGamma {
             continue;
           }
 
+          if (pairtype == PairType::kEMCEMC) {
+            float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P()));
+            if (openingAngle < emccuts.minOpenAngle) {
+              continue;
+            }
+          }
+
           fRegistry.fill(HIST("Pair/same/hs"), v12.M(), v12.Pt(), collision.weight());
 
           if constexpr (pairtype == PairType::kEMCEMC) {

PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

@@ -500,7 +500,7 @@ struct Pi0EtaToGammaGammaMC {
           pi0id = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 111, mcparticles);
           etaid = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 221, mcparticles);
 
-          if (pi0id < 0 && etaid < 0) {
+          if (g1mc.globalIndex() != g2mc.globalIndex() && pi0id < 0 && etaid < 0) { // for same gamma no pi0/eta will be found, but we still want to fill the FromSameGamma hist
             continue;
           }
 
@@ -511,6 +511,21 @@ struct Pi0EtaToGammaGammaMC {
             continue;
           }
 
+          if (pairtype == PairType::kEMCEMC) {
+            float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P()));
+            if (openingAngle < emccuts.minOpenAngle) {
+              continue;
+            }
+          }
+
+          if (g1mc.globalIndex() == g2mc.globalIndex()) {
+            if (getMotherPDGCode(g1mc, mcparticles) == 111)
+              fRegistry.fill(HIST("Pair/Pi0/hs_FromSameGamma"), v12.M(), v12.Pt(), collision.weight());
+            else if (getMotherPDGCode(g1mc, mcparticles) == 221)
+              fRegistry.fill(HIST("Pair/Eta/hs_FromSameGamma"), v12.M(), v12.Pt(), collision.weight());
+            continue;
+          }
+
           if (pi0id > 0) {
             auto pi0mc = mcparticles.iteratorAt(pi0id);
             o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, pi0mc, mcparticles, mccollisions, f1fd_k0s_to_pi0, collision.weight());

PWGEM/PhotonMeson/Utils/NMHistograms.h

@@ -48,6 +48,8 @@ void addNMHistograms(HistogramRegistry* fRegistry, bool isMC, const char* pairna
     fRegistry->add("Pair/Pi0/hs_Primary", "rec. true pi0", kTHnSparseD, {axis_mass, axis_pt}, true);
     fRegistry->add("Pair/Pi0/hs_FromWD", "rec. true pi0 from weak decay", kTHnSparseD, {axis_mass, axis_pt}, true);
     fRegistry->add("Pair/Pi0/hs_FromHS", "rec. true pi0 from hadronic shower in material", kTHnSparseD, {axis_mass, axis_pt}, true);
+    fRegistry->add("Pair/Pi0/hs_FromSameGamma", "Two clusters from same gamma that is a pi0 daughter (conversion)", kTHnSparseD, {axis_mass, axis_pt}, true);
+    fRegistry->add("Pair/Eta/hs_FromSameGamma", "Two clusters from same gamma that is a eta daughter (conversion)", kTHnSparseD, {axis_mass, axis_pt}, true);
     fRegistry->add("Pair/Eta/hs_Primary", "rec. true eta", kTHnSparseD, {axis_mass, axis_pt}, true);
 
     const AxisSpec axis_rapidity{{0.0, +0.8, +0.9}, "rapidity |y|"};

PWGJE/Tasks/emcalGammaGammaBcWise.cxx

@@ -25,6 +25,7 @@
 #include "Framework/AnalysisTask.h"
 #include "Framework/HistogramRegistry.h"
 #include "Common/DataModel/EventSelection.h"
+#include "EMCALBase/Geometry.h"
 #include "PWGJE/DataModel/EMCALClusters.h"
 
 using namespace o2;
@@ -49,6 +50,59 @@ struct Photon { // Struct to store photons (unique and ambiguous clusters that p
   float px, py, pz, pt;
 };
 
+/// \brief returns if cluster is too close to edge of EMCal (using rotation background method only for EMCal!)
+bool IsTooCloseToEdge(const int cellID, const int DistanceToBorder = 1, emcal::Geometry* emcalGeom = nullptr)
+{
+  if (DistanceToBorder <= 0) {
+    return false;
+  }
+  if (cellID < 0) {
+    return true;
+  }
+
+  int iBadCell = -1;
+
+  // check distance to border in case the cell is okay
+  auto [iSupMod, iMod, iPhi, iEta] = emcalGeom->GetCellIndex(cellID);
+  auto [irow, icol] = emcalGeom->GetCellPhiEtaIndexInSModule(iSupMod, iMod, iPhi, iEta);
+
+  // Check rows/phi
+  int iRowLast = 24;
+  if (emcalGeom->GetSMType(iSupMod) == o2::emcal::EMCALSMType::EMCAL_HALF) {
+    iRowLast /= 2; // 2/3 sm case
+  } else if (emcalGeom->GetSMType(iSupMod) == o2::emcal::EMCALSMType::EMCAL_THIRD) {
+    iRowLast /= 3; // 1/3 sm case
+  } else if (emcalGeom->GetSMType(iSupMod) == o2::emcal::EMCALSMType::DCAL_EXT) {
+    iRowLast /= 3; // 1/3 sm case
+  }
+
+  if (irow < DistanceToBorder || (iRowLast - irow) <= DistanceToBorder) {
+    iBadCell = 1;
+  }
+
+  if (iBadCell > 0) {
+    return true;
+  }
+  return false;
+}
+
+bool isPhotonAccepted(Photon const& p, emcal::Geometry* emcalGeom = nullptr)
+{
+  int cellID = 0;
+  try {
+    cellID = emcalGeom->GetAbsCellIdFromEtaPhi(p.eta, p.phi);
+    if (IsTooCloseToEdge(cellID, 1, emcalGeom))
+      cellID = -1;
+  } catch (o2::emcal::InvalidPositionException& e) {
+    cellID = -1;
+  }
+
+  if (cellID == -1)
+    return false;
+
+  return true;
+}
+
 struct Meson {
   Meson(Photon p1, Photon p2) : p1(p1), p2(p2)
   {
@@ -65,20 +119,26 @@ struct Meson {
 struct EmcalGammaGammaBcWise {
   HistogramRegistry mHistManager{"EmcalGammaGammaBcWiseHistograms"};
 
+  Configurable<int> cfgClusterDefinition{"cfgClusterDefinition", 13, "Clusterizer to be selected, e.g. 13 for kV3MostSplitLowSeed"};
   Configurable<float> cfgMinClusterEnergy{"cfgMinClusterEnergy", 0.7, "Minimum energy of selected clusters (GeV)"};
   Configurable<float> cfgMinM02{"cfgMinM02", 0.1, "Minimum M02 of selected clusters"};
   Configurable<float> cfgMaxM02{"cfgMaxM02", 0.7, "Maximum M02 of selected clusters"};
   Configurable<float> cfgMinTime{"cfgMinTime", -15, "Minimum time of selected clusters (ns)"};
   Configurable<float> cfgMaxTime{"cfgMaxTime", 15, "Maximum time of selected clusters (ns)"};
+  Configurable<float> cfgMinOpenAngle{"cfgMinOpenAngle", 0.0202, "Minimum opening angle between photons"};
 
+  Filter clusterDefinitionFilter = aod::emcalcluster::definition == cfgClusterDefinition;
   Filter energyFilter = aod::emcalcluster::energy > cfgMinClusterEnergy;
   Filter m02Filter = (aod::emcalcluster::nCells == 1 || (aod::emcalcluster::m02 > cfgMinM02 && aod::emcalcluster::m02 < cfgMaxM02));
   Filter timeFilter = (aod::emcalcluster::time > cfgMinTime && aod::emcalcluster::time < cfgMaxTime);
 
   std::vector<Photon> mPhotons;
 
+  emcal::Geometry* emcalGeom;
+
   void init(InitContext const&)
   {
+    emcalGeom = emcal::Geometry::GetInstanceFromRunNumber(300000);
     mHistManager.add("nBCs", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{3, -0.5, 2.5}});
 
     mHistManager.add("nCollisionsVsClusters", "Number of collisions vs number of clusters;N_{Collisions};N_{Clusters}", HistType::kTH2F, {{4, -0.5, 3.5}, {21, -0.5, 20.5}});
@@ -149,6 +209,8 @@ struct EmcalGammaGammaBcWise {
     for (unsigned int ig1 = 0; ig1 < mPhotons.size(); ++ig1) {
       for (unsigned int ig2 = ig1 + 1; ig2 < mPhotons.size(); ++ig2) {
         // build meson from photons
+        if (mPhotons[ig1].photon.Angle(mPhotons[ig2].photon.Vect()) < cfgMinOpenAngle)
+          continue;
         Meson meson(mPhotons[ig1], mPhotons[ig2]);
         mHistManager.fill(HIST("invMassVsPt"), meson.m(), meson.pT());
 
@@ -171,6 +233,10 @@ struct EmcalGammaGammaBcWise {
         TLorentzVector lvRotationPhoton(mPhotons[ig].px, mPhotons[ig].py, mPhotons[ig].pz, mPhotons[ig].e);
         lvRotationPhoton.Rotate(constants::math::PIHalf, lvRotationPion);
         Photon rotPhoton(lvRotationPhoton.Eta(), lvRotationPhoton.Phi(), lvRotationPhoton.E());
+        if (!isPhotonAccepted(rotPhoton, emcalGeom))
+          continue;
+        if (rotPhoton.photon.Angle(mPhotons[ig3].photon.Vect()) < cfgMinOpenAngle)
+          continue;
         Meson mesonRotated(rotPhoton, mPhotons[ig3]);
         mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated.m(), mesonRotated.pT());
       }