diff --git a/PWGEM/Dilepton/Utils/MCUtilities.h b/PWGEM/Dilepton/Utils/MCUtilities.h index 9bd0797259b..58484f6fe66 100644 --- a/PWGEM/Dilepton/Utils/MCUtilities.h +++ b/PWGEM/Dilepton/Utils/MCUtilities.h @@ -222,6 +222,14 @@ int FindCommonMotherFrom3Prongs(TMCParticle1 const& p1, TMCParticle2 const& p2, } //_______________________________________________________________________ template +int getMotherPDGCode(TMCParticle const& p, TMCParticles const& mcparticles) +{ + int motherid = p.mothersIds()[0]; + auto mother = mcparticles.iteratorAt(motherid); + return (mother.pdgCode()); +} +//_______________________________________________________________________ +template int IsFromBeauty(TMCParticle const& p, TMCParticles const& mcparticles) { if (!p.has_mothers()) { diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h index 03032cf5112..d26dca62045 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h +++ b/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) { diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h index d9d436f4ce9..e39d4165e4e 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h +++ b/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()); diff --git a/PWGEM/PhotonMeson/Utils/NMHistograms.h b/PWGEM/PhotonMeson/Utils/NMHistograms.h index 95abaf92292..a3e62b0bc5a 100644 --- a/PWGEM/PhotonMeson/Utils/NMHistograms.h +++ b/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|"}; diff --git a/PWGJE/Tasks/emcalGammaGammaBcWise.cxx b/PWGJE/Tasks/emcalGammaGammaBcWise.cxx index 7e6843efaab..7007f94bc15 100644 --- a/PWGJE/Tasks/emcalGammaGammaBcWise.cxx +++ b/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 cfgClusterDefinition{"cfgClusterDefinition", 13, "Clusterizer to be selected, e.g. 13 for kV3MostSplitLowSeed"}; Configurable cfgMinClusterEnergy{"cfgMinClusterEnergy", 0.7, "Minimum energy of selected clusters (GeV)"}; Configurable cfgMinM02{"cfgMinM02", 0.1, "Minimum M02 of selected clusters"}; Configurable cfgMaxM02{"cfgMaxM02", 0.7, "Maximum M02 of selected clusters"}; Configurable cfgMinTime{"cfgMinTime", -15, "Minimum time of selected clusters (ns)"}; Configurable cfgMaxTime{"cfgMaxTime", 15, "Maximum time of selected clusters (ns)"}; + Configurable 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 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()); }