Skip to content
Open
Changes from all commits
Commits
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
91 changes: 62 additions & 29 deletions PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@
#include <Math/Vector4Dfwd.h>
#include <TF1.h>
#include <TH2.h>
#include <TPDGCode.h>
#include <TString.h>

#include <cmath>
Expand All @@ -74,24 +75,25 @@

template <o2::aod::pwgem::photonmeson::photonpair::PairType pairtype, o2::soa::is_table... Types>
struct Pi0EtaToGammaGammaMC {
o2::framework::Configurable<std::string> ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
o2::framework::Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
o2::framework::Configurable<std::string> grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"};
o2::framework::Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
o2::framework::Configurable<bool> skipGRPOquery{"skipGRPOquery", true, "skip grpo query"};
o2::framework::Configurable<float> d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"};
o2::framework::Configurable<float> dBzInput{"dBzInput", -999, "bz field in kG, -999 is automatic"};

o2::framework::Configurable<int> cfgQvecEstimator{"cfgQvecEstimator", 0, "FT0M:0, FT0A:1, FT0C:2"};
o2::framework::Configurable<int> cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"};
o2::framework::Configurable<float> cfgCentMin{"cfgCentMin", -1, "min. centrality"};
o2::framework::Configurable<float> cfgCentMax{"cfgCentMax", 999, "max. centrality"};
o2::framework::Configurable<float> maxY_rec{"maxY_rec", 0.9, "maximum rapidity for reconstructed particles"};
o2::framework::Configurable<std::string> fd_k0s_to_pi0{"fd_k0s_pi0", "1.0", "feed down correction to pi0"};
o2::framework::Configurable<float> maxYRec{"maxYRec", 0.9, "maximum rapidity for reconstructed particles"};
o2::framework::Configurable<std::string> fdK0sToPi0{"fdK0sPi0", "1.0", "feed down correction to pi0"};

Check failure on line 89 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<bool> cfgRequireTrueAssociation{"cfgRequireTrueAssociation", false, "flag to require true mc collision association"};

o2::framework::Configurable<int> cfgAlphaMesonCut{"cfgAlphaMesonCut", 0, "flag for photon energy asymmetry distribution cut: 0: no cut, 1: cut specific value, 2: cut depending on pT"};
o2::framework::Configurable<float> cfgAlphaMeson{"cfgAlphaMeson", 0.65, "photon energy asymmetry distribution parameter for specific value cut"};
o2::framework::Configurable<float> cfgAlphaMesonA{"cfgAlphaMesonA", 0.65, "photon energy asymmetry distribution parameter A for pT dependent cut (A * tanh(B*pT))"};
o2::framework::Configurable<float> cfgAlphaMesonB{"cfgAlphaMesonB", 1.2, "photon energy asymmetry distribution parameter B for pT dependent cut (A * tanh(B*pT))"};
o2::framework::Configurable<bool> cfgGGContaCheck{"cfgGGContaCheck", false, "check gamma gamma contamination of dalitz"};

EMPhotonEventCut fEMEventCut;
struct : o2::framework::ConfigurableGroup {
Expand All @@ -117,15 +119,15 @@
V0PhotonCut fV0PhotonCut;
struct : o2::framework::ConfigurableGroup {
std::string prefix = "pcmcut_group";
o2::framework::Configurable<bool> cfg_require_v0_with_itstpc{"cfg_require_v0_with_itstpc", false, "flag to select V0s with ITS-TPC matched tracks"};

Check failure on line 122 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<bool> cfg_require_v0_with_itsonly{"cfg_require_v0_with_itsonly", false, "flag to select V0s with ITSonly tracks"};

Check failure on line 123 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<bool> cfg_require_v0_with_tpconly{"cfg_require_v0_with_tpconly", false, "flag to select V0s with TPConly tracks"};

Check failure on line 124 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<float> cfg_min_pt_v0{"cfg_min_pt_v0", 0.1, "min pT for v0 photons at PV"};

Check failure on line 125 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<float> cfg_max_pt_v0{"cfg_max_pt_v0", 1e+10, "max pT for v0 photons at PV"};

Check failure on line 126 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<float> cfg_min_eta_v0{"cfg_min_eta_v0", -0.8, "min eta for v0 photons at PV"};

Check failure on line 127 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<float> cfg_max_eta_v0{"cfg_max_eta_v0", 0.8, "max eta for v0 photons at PV"};

Check failure on line 128 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<float> cfg_min_v0radius{"cfg_min_v0radius", 4.0, "min v0 radius"};

Check failure on line 129 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<float> cfg_max_v0radius{"cfg_max_v0radius", 90.0, "max v0 radius"};

Check failure on line 130 in PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
o2::framework::Configurable<float> cfg_max_alpha_ap{"cfg_max_alpha_ap", 0.95, "max alpha for AP cut"};
o2::framework::Configurable<float> cfg_max_qt_ap{"cfg_max_qt_ap", 0.01, "max qT for AP cut"};
o2::framework::Configurable<float> cfg_min_cospa{"cfg_min_cospa", 0.999, "min V0 CosPA"};
Expand Down Expand Up @@ -224,7 +226,7 @@
o2::framework::Configurable<float> cfg_min_Ecluster{"cfg_min_Ecluster", 0.3, "Minimum cluster energy for PHOS in GeV"};
} phoscuts;

TF1* f1fd_k0s_to_pi0;
TF1* f1fdK0sToPi0;
o2::framework::HistogramRegistry fRegistry{"output", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject, false, false};
// static constexpr std::string_view event_types[2] = {"before/", "after/"};
// static constexpr std::string_view event_pair_types[2] = {"same/", "mix/"};
Expand All @@ -249,14 +251,20 @@
DefineEMCCut();
DefinePHOSCut();

f1fd_k0s_to_pi0 = new TF1("f1fd_k0s_to_pi0", TString(fd_k0s_to_pi0), 0.f, 100.f);
f1fdK0sToPi0 = new TF1("f1fdK0sToPi0", TString(fdK0sToPi0), 0.f, 100.f);

fRegistry.add("Event/hNrecPerMCCollision", "Nrec per mc collision;N_{rec} collisions per MC collision", o2::framework::HistType::kTH1F, {{21, -0.5f, 20.5f}}, false);
if (cfgGGContaCheck) {
fRegistry.add("Event/hNGGContamEta", "Number of Eta from etaToGammaGamma; p_{T, #eta} (GeV/#it{c}); N", o2::framework::HistType::kTH1F, {{40, -0.5f, 20.5f}}, false);
fRegistry.add("Event/hNGGContamPion", "Number of Pion from etaToGammaGamma; p_{T, #pi} (GeV/#it{c}); N", o2::framework::HistType::kTH1F, {{40, -0.5f, 20.5f}}, false);
}
fRegistry.add("Event/hNDalitzEtaPt", "Number of DalitzEta; p_{T, #eta} (GeV/#it{c}); N", o2::framework::HistType::kTH1F, {{40, -0.5f, 20.5f}}, false);
fRegistry.add("Event/hNDalitzPionPt", "Number of DalitzPion; p_{T, #pi} (GeV/#it{c}) ; N", o2::framework::HistType::kTH1F, {{40, -0.5f, 20.5f}}, false);

mRunNumber = 0;
d_bz = 0;

ccdb->setURL(ccdburl);
ccdb->setURL(ccdbUrl);
ccdb->setCaching(true);
ccdb->setLocalObjectValidityChecking();
ccdb->setFatalWhenNull(false);
Expand All @@ -270,10 +278,12 @@
}

// In case override, don't proceed, please - no CCDB access required
if (d_bz_input > -990) {
d_bz = d_bz_input;
constexpr float bzInput = -990.0f;
if (dBzInput > bzInput) {
d_bz = dBzInput;
o2::parameters::GRPMagField grpmag;
if (std::fabs(d_bz) > 1e-5) {
float bzThreshold = 1e-5;
if (std::fabs(d_bz) > bzThreshold) {
grpmag.setL3Current(30000.f / (d_bz / 5.0f));
}
mRunNumber = collision.runNumber();
Expand Down Expand Up @@ -304,8 +314,8 @@

~Pi0EtaToGammaGammaMC()
{
delete f1fd_k0s_to_pi0;
f1fd_k0s_to_pi0 = 0x0;
delete f1fdK0sToPi0;
f1fdK0sToPi0 = 0x0;
}

void DefineEMEventCut()
Expand Down Expand Up @@ -359,7 +369,7 @@
fV0PhotonCut.SetLoadMlModelsFromCCDB(pcmcuts.cfg_load_ml_models_from_ccdb);
fV0PhotonCut.SetNClassesMl(pcmcuts.cfg_nclasses_ml);
fV0PhotonCut.SetMlTimestampCCDB(pcmcuts.cfg_timestamp_ccdb);
fV0PhotonCut.SetCcdbUrl(ccdburl);
fV0PhotonCut.SetCcdbUrl(ccdbUrl);
CentType mCentralityTypeMlEnum;
mCentralityTypeMlEnum = static_cast<CentType>(cfgCentEstimator.value);
fV0PhotonCut.SetCentralityTypeMl(mCentralityTypeMlEnum);
Expand Down Expand Up @@ -572,18 +582,19 @@
TMCCollisions const& mccollisions, TMCParticles const& mcparticles,
TLegs const& /*legs*/ = nullptr, TMatchedTracks const& matchedTracks = nullptr, TMatchedSecondaries const& matchedSecondaries = nullptr)
{
for (auto& collision : collisions) {
for (auto const& collision : collisions) {
initCCDB(collision);
if ((pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPHOS) && !collision.alias_bit(triggerAliases::kTVXinPHOS)) {
continue;
}

float weight = 1.f;
float weightThreshold = 1e-10;
if constexpr (std::is_same_v<std::decay_t<TCollisions>, o2::soa::Filtered<o2::soa::Join<o2::soa::Join<o2::aod::PMEvents, o2::aod::EMEventsAlias, o2::aod::EMEventsMult_000, o2::aod::EMEventsCent_000, o2::aod::EMMCEventLabels>, o2::aod::EMEventsWeight>>>) {
weight = collision.weight();
}

if (eventcuts.onlyKeepWeightedEvents && std::fabs(weight - 1.0) < 1e-10) {
if (eventcuts.onlyKeepWeightedEvents && std::fabs(weight - 1.0) < weightThreshold) {
continue;
}

Expand Down Expand Up @@ -680,7 +691,7 @@
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
if (std::fabs(v12.Rapidity()) > maxY_rec) {
if (std::fabs(v12.Rapidity()) > maxYRec) {
continue;
}

Expand Down Expand Up @@ -720,9 +731,9 @@
}

if (g1mc.globalIndex() == g2mc.globalIndex()) {
if (o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(g1mc, mcparticles) == 111)
if (o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(g1mc, mcparticles) == PDG_t::kPi0)
fRegistry.fill(HIST("Pair/Pi0/hs_FromSameGamma"), v12.M(), v12.Pt(), wpair);
else if (o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(g1mc, mcparticles) == 221)
else if (o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(g1mc, mcparticles) == o2::constants::physics::Pdg::kEta)
fRegistry.fill(HIST("Pair/Eta/hs_FromSameGamma"), v12.M(), v12.Pt(), wpair);
continue;
}
Expand All @@ -732,13 +743,13 @@
if (cfgRequireTrueAssociation && (pi0mc.emmceventId() != collision.emmceventId())) {
continue;
}
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, pi0mc, mcparticles, mccollisions, f1fd_k0s_to_pi0, wpair);
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, pi0mc, mcparticles, mccollisions, f1fdK0sToPi0, wpair);
} else if (etaid > 0) {
auto etamc = mcparticles.iteratorAt(etaid);
if (cfgRequireTrueAssociation && (etamc.emmceventId() != collision.emmceventId())) {
continue;
}
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, etamc, mcparticles, mccollisions, f1fd_k0s_to_pi0, wpair);
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, etamc, mcparticles, mccollisions, f1fdK0sToPi0, wpair);
}
} // end of pairing loop
} else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) {
Expand Down Expand Up @@ -792,6 +803,25 @@

auto pos2mc = mcparticles.iteratorAt(pos2.emmcparticleId());
auto ele2mc = mcparticles.iteratorAt(ele2.emmcparticleId());
if (cfgGGContaCheck) {
photonid2 = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom2Prongs(pos2mc, ele2mc, -11, 11, 22, mcparticles); // check possible contamination
if (photonid2 > 0) {
auto photon2 = mcparticles.iteratorAt(photonid2);
int photon2pdg = photon2.pdgCode();
int photon2mothid = photon2.mothersIds()[0];
auto photon2moth = mcparticles.iteratorAt(photon2mothid);
if (photon2pdg == PDG_t::kGamma && (o2::aod::pwgem::photonmeson::utils::mcutil::isGammaGammaDecay(photon2moth, mcparticles))) {
int mothID = o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(photon2, mcparticles);
if (mothID == o2::constants::physics::Pdg::kEta) {
fRegistry.fill(HIST("Event/hNGGContamEta"), photon2moth.pt());
}
if (mothID == PDG_t::kPi0) {
fRegistry.fill(HIST("Event/hNGGContamPion"), photon2moth.pt());
}
}
}
}

pi0id = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -11, 11, 111, mcparticles);
etaid = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -11, 11, 221, mcparticles);
if (pi0id < 0 && etaid < 0) {
Expand All @@ -800,21 +830,23 @@
ROOT::Math::PtEtaPhiMVector v_pos(pos2.pt(), pos2.eta(), pos2.phi(), o2::constants::physics::MassElectron);
ROOT::Math::PtEtaPhiMVector v_ele(ele2.pt(), ele2.eta(), ele2.phi(), o2::constants::physics::MassElectron);
ROOT::Math::PtEtaPhiMVector veeg = v_gamma + v_pos + v_ele;
if (std::fabs(veeg.Rapidity()) > maxY_rec) {
if (std::fabs(veeg.Rapidity()) > maxYRec) {
continue;
}
if (pi0id > 0) {
auto pi0mc = mcparticles.iteratorAt(pi0id);
fRegistry.fill(HIST("Event/hNDalitzPionPt"), pi0mc.pt());
if (cfgRequireTrueAssociation && (pi0mc.emmceventId() != collision.emmceventId())) {
continue;
}
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, veeg, pi0mc, mcparticles, mccollisions, f1fd_k0s_to_pi0, weight);
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, veeg, pi0mc, mcparticles, mccollisions, f1fdK0sToPi0, weight);
} else if (etaid > 0) {
auto etamc = mcparticles.iteratorAt(etaid);
fRegistry.fill(HIST("Event/hNDalitzEtaPt"), etamc.pt());
if (cfgRequireTrueAssociation && (etamc.emmceventId() != collision.emmceventId())) {
continue;
}
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, veeg, etamc, mcparticles, mccollisions, f1fd_k0s_to_pi0, weight);
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, veeg, etamc, mcparticles, mccollisions, f1fdK0sToPi0, weight);
}
} // end of dielectron loop
} // end of pcm loop
Expand Down Expand Up @@ -853,15 +885,15 @@
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
if (std::fabs(v12.Rapidity()) > maxY_rec) {
if (std::fabs(v12.Rapidity()) > maxYRec) {
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, weight);
// o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, pi0mc, mcparticles, mccollisions, f1fdK0sToPi0, weight);
// } else if (etaid > 0) {
// auto etamc = mcparticles.iteratorAt(etaid);
// o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, etamc, mcparticles, mccollisions, f1fd_k0s_to_pi0, weight);
// o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, etamc, mcparticles, mccollisions, f1fdK0sToPi0, weight);
// }
} // end of pairing loop
} // end of pairing in same event
Expand Down Expand Up @@ -900,23 +932,24 @@
// loop over mc stack and fill histograms for pure MC truth signals
// all MC tracks which belong to the MC event corresponding to the current reconstructed event

for (auto& mccollision : mccollisions) {
for (auto const& mccollision : mccollisions) {
auto collision_per_mccoll = collisions.sliceBy(rec_perMcCollision, mccollision.globalIndex());
int nrec_per_mc = collision_per_mccoll.size();
fRegistry.fill(HIST("Event/hNrecPerMCCollision"), nrec_per_mc);
}

for (auto& collision : collisions) {
for (auto const& collision : collisions) {
if ((pairtype == o2::aod::pwgem::photonmeson::photonpair::kPHOSPHOS || pairtype == o2::aod::pwgem::photonmeson::photonpair::kPCMPHOS) && !collision.alias_bit(triggerAliases::kTVXinPHOS)) {
continue; // I don't know why this is necessary in simulation.
}

float weight = 1.f;
float weightTreshhold = 1e-10;
if constexpr (std::is_same_v<std::decay_t<TCollisions>, o2::soa::Filtered<o2::soa::Join<o2::soa::Join<o2::aod::PMEvents, o2::aod::EMEventsAlias, o2::aod::EMEventsMult_000, o2::aod::EMEventsCent_000, o2::aod::EMMCEventLabels>, o2::aod::EMEventsWeight>>>) {
weight = collision.weight();
}

if (eventcuts.onlyKeepWeightedEvents && std::fabs(weight - 1.0) < 1e-10) {
if (eventcuts.onlyKeepWeightedEvents && std::fabs(weight - 1.0) < weightTreshhold) {
continue;
}

Expand Down
Loading