Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mhemmer-cern merged 1 commit into from
Feb 8, 2025
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions PWGEM/Dilepton/Utils/MCUtilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,14 @@
}
//_______________________________________________________________________
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()) {
Expand All @@ -230,14 +238,14 @@

int motherid = p.mothersIds()[0]; // first mother index
auto mp_tmp = mcparticles.iteratorAt(motherid);
if (abs(mp_tmp.pdgCode()) < 1e+9 && (std::to_string(abs(mp_tmp.pdgCode()))[std::to_string(abs(mp_tmp.pdgCode())).length() - 2] == '5' && std::to_string(abs(mp_tmp.pdgCode()))[std::to_string(abs(mp_tmp.pdgCode())).length() - 3] == '5') && abs(mp_tmp.pdgCode()) % 2 == 1) {

Check warning on line 241 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
return -999; // reject bottomonia
}

while (motherid > -1) {
if (motherid < mcparticles.size()) { // protect against bad mother indices. why is this needed?
auto mp = mcparticles.iteratorAt(motherid);
if (abs(mp.pdgCode()) < 1e+9 && (std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 3] == '5' || std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 4] == '5')) {

Check warning on line 248 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
return motherid;
}
if (mp.has_mothers()) {
Expand All @@ -263,13 +271,13 @@

int motherid = p.mothersIds()[0]; // first mother index
auto mp_tmp = mcparticles.iteratorAt(motherid);
if (abs(mp_tmp.pdgCode()) < 1e+9 && (std::to_string(abs(mp_tmp.pdgCode()))[std::to_string(abs(mp_tmp.pdgCode())).length() - 2] == '4' && std::to_string(abs(mp_tmp.pdgCode()))[std::to_string(abs(mp_tmp.pdgCode())).length() - 3] == '4') && abs(mp_tmp.pdgCode()) % 2 == 1) {

Check warning on line 274 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
return -999; // reject bottomonia
}
while (motherid > -1) {
if (motherid < mcparticles.size()) { // protect against bad mother indices. why is this needed?
auto mp = mcparticles.iteratorAt(motherid);
if (abs(mp.pdgCode()) < 1e+9 && (std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 3] == '4' || std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 4] == '4')) {

Check warning on line 280 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
return motherid;
}
if (mp.has_mothers()) {
Expand Down Expand Up @@ -384,8 +392,8 @@
if (mid1 == mid2) {
auto common_mp = mcparticles.iteratorAt(mid1);
int mp_pdg = common_mp.pdgCode();
bool is_mp_diquark = (1100 < abs(mp_pdg) && abs(mp_pdg) < 5600) && std::to_string(mp_pdg)[std::to_string(mp_pdg).length() - 2] == '0';

Check warning on line 395 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
if (!is_mp_diquark && abs(mp_pdg) < 1e+9 && (std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 3] == '5' || std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 4] == '5')) {

Check warning on line 396 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
mothers_id1.clear();
mothers_pdg1.clear();
mothers_id2.clear();
Expand All @@ -406,8 +414,8 @@
if (mid1 == mid2) {
auto common_mp = mcparticles.iteratorAt(mid1);
int mp_pdg = common_mp.pdgCode();
bool is_mp_diquark = (1100 < abs(mp_pdg) && abs(mp_pdg) < 5600) && std::to_string(mp_pdg)[std::to_string(mp_pdg).length() - 2] == '0';

Check warning on line 417 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
if (!is_mp_diquark && abs(mp_pdg) < 1e+9 && (std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 3] == '5' || std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 4] == '5')) {

Check warning on line 418 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
is_same_mother_found = true;
}
}
Expand Down Expand Up @@ -455,7 +463,7 @@
} else if (mothersids[1] < mothersids[0]) {
allmothersids.push_back(mothersids[0]);
allmothersids.push_back(mothersids[1]);
} else if ((80 < abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) && abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) < 90) || (100 < abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) && abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) < 110)) { // NOTE: THIS IS GENERATOR DEPENDENT AND WORKS ONLY FOR PYTHIA!

Check warning on line 466 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
for (int i = mothersids[0]; i <= mothersids[1]; i++) {
allmothersids.push_back(i);
}
Expand All @@ -473,7 +481,7 @@
for (int i : allmothersids) {
auto mother = mcParticles.iteratorAt(i);
int mpdg = mother.pdgCode();
if (abs(mpdg) == pdg && mpdg * p.pdgCode() > 0) { // check for quark

Check warning on line 484 in PWGEM/Dilepton/Utils/MCUtilities.h

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
if (quark_id > -1 || next_mother_id > -1) { // we already found a possible candidate in the list of mothers, so now we have (at least) two
// LOG(warning) << "Flavour tracking is ambiguous. Stopping here.";
return -1;
Expand Down
7 changes: 7 additions & 0 deletions PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
17 changes: 16 additions & 1 deletion PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

Expand All @@ -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());
Expand Down
2 changes: 2 additions & 0 deletions PWGEM/PhotonMeson/Utils/NMHistograms.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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|"};
Expand Down
66 changes: 66 additions & 0 deletions PWGJE/Tasks/emcalGammaGammaBcWise.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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)
{
Expand All @@ -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}});
Expand Down Expand Up @@ -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());

Expand All @@ -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());
}
Expand Down
Loading