Skip to content

PWGEM/PhotonMeson: update MC analysis #2427

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
Apr 18, 2023
Merged

PWGEM/PhotonMeson: update MC analysis #2427

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
31 changes: 18 additions & 13 deletions PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
}
if (TString(histClass) == "Track") {
list->Add(new TH1F("hPt", "pT", 1000, 0.0f, 10));
list->Add(new TH1F("hQoverPt", "Q/pT;Q/p_{T} (GeV/c)^{-1}", 1000, -50, 50));
list->Add(new TH2F("hEtaPhi", "#eta vs. #varphi", 180, 0, TMath::TwoPi(), 40, -2.0f, 2.0f));
list->Add(new TH2F("hDCAxyz", "DCA xy vs. z;DCA_{xy} (cm);DCA_{z} (cm)", 100, -5.0f, 5.0f, 100, -5.0f, 5.0f));
list->Add(new TH1F("hNclsTPC", "number of TPC clusters", 161, -0.5, 160.5));
Expand All @@ -63,12 +64,16 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
list->Add(new TH2F("hAPplot", "AP plot;#alpha;q_{T} (GeV/c)", 200, -1.0f, +1.0f, 250, 0.0f, 0.25f));
list->Add(new TH2F("hGammaPsiPair", "#psi_{pair} for photon conversion;#psi_{pair} (rad.);m_{ee} (GeV/c^{2})", 150, 0, TMath::PiOver2(), 100, 0.0f, 0.1f));
list->Add(new TH2F("hMassGamma", "hMassGamma;R_{xy} (cm);m_{ee} (GeV/c^{2})", 200, 0.0f, 100.0f, 100, 0.0f, 0.1f));
list->Add(new TH2F("hMassGamma_recalc", "recalc. KF hMassGamma;R_{xy} (cm);m_{ee} (GeV/c^{2})", 200, 0.0f, 100.0f, 100, 0.0f, 0.1f));
list->Add(new TH2F("hMassGamma_recalc", "recalc. hMassGamma;R_{xy} (cm);m_{ee} (GeV/c^{2})", 200, 0.0f, 100.0f, 100, 0.0f, 0.1f));
list->Add(new TH2F("hGammaRxy", "conversion point in XY;V_{x} (cm);V_{y} (cm)", 400, -100.0f, 100.0f, 400, -100.0f, 100.0f));
list->Add(new TH2F("hGammaRxy_recalc", "recalc. KF conversion point in XY;V_{x} (cm);V_{y} (cm)", 400, -100.0f, 100.0f, 400, -100.0f, 100.0f));
list->Add(new TH2F("hKFChi2vsR_recalc", "recalc. KF conversion point in XY;R_{xy} (cm);KF chi2/NDF", 250, 0.0f, 250.0f, 500, 0.f, 5000.0f));
list->Add(new TH2F("hKFChi2vsZ_recalc", "recalc. KF conversion point in Z;Z (cm);KF chi2/NDF", 500, -250.0f, 250.0f, 500, 0.f, 5000.0f));
list->Add(new TH2F("hGammaRxy_recalc", "recalc. conversion point in XY;V_{x} (cm);V_{y} (cm)", 400, -100.0f, 100.0f, 400, -100.0f, 100.0f));
list->Add(new TH2F("hKFChi2vsR_recalc", "KF chi2 vs. recalc. conversion point in XY;R_{xy} (cm);KF chi2/NDF", 250, 0.0f, 250.0f, 500, 0.f, 5000.0f));
list->Add(new TH2F("hKFChi2vsZ_recalc", "KF chi2 vs. recalc. conversion point in Z;Z (cm);KF chi2/NDF", 500, -250.0f, 250.0f, 500, 0.f, 5000.0f));
list->Add(new TH1F("hNgamma", "Number of #gamma candidates per collision", 101, -0.5f, 100.5f));
list->Add(new TH1F("hPt_Primary", "pT;p_{T} (GeV/c)", 1000, 0.0f, 10)); // for MC efficiency
list->Add(new TH2F("hEtaPhi_Primary", "#eta vs. #varphi;#varphi (rad.);#eta", 180, 0, TMath::TwoPi(), 40, -2.0f, 2.0f)); // for MC efficiency
list->Add(new TH1F("hPt_FromWD", "pT;p_{T} (GeV/c)", 1000, 0.0f, 10)); // for MC feed down correction
list->Add(new TH2F("hEtaPhi_FromWD", "#eta vs. #varphi;#varphi (rad.);#eta", 180, 0, TMath::TwoPi(), 40, -2.0f, 2.0f)); // for MC feed down correction
}

if (TString(histClass) == "gammagamma_mass_pt") {
Expand All @@ -86,16 +91,16 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
list->Add(new TH1F("hZvtx_after", "vertex z; Zvtx (cm)", 100, -50, +50));

if (TString(subGroup) == "ConversionStudy") {
const float rxy[] = {0, 6, 10, 20, 30, 40, 50, 60, 70, 80, 90};
list->Add(new TH2F("hGammaRxy", "conversion point in XY MC;V_{x} (cm);V_{y} (cm)", 2000, -100.0f, 100.0f, 2000, -100.0f, 100.0f));
list->Add(new TH2F("hGammaRZ", "conversion point in RZ MC;V_{z} (cm);R_{xy} (cm)", 5000, -250.0f, 250.0f, 1000, 0.f, 100.0f));
list->Add(new TH2F("hPhotonRxy", "conversion point in XY MC;V_{x} (cm);V_{y} (cm)", 2000, -100.0f, 100.0f, 2000, -100.0f, 100.0f));
list->Add(new TH2F("hPhotonRZ", "conversion point in RZ MC;V_{z} (cm);R_{xy} (cm)", 5000, -250.0f, 250.0f, 1000, 0.f, 100.0f));
list->Add(new TH2F("hPhotonPhivsRxy", "conversion point of #varphi vs. R_{xy} MC;#varphi (rad.);R_{xy} (cm);N_{e}", 360, 0.0f, TMath::TwoPi(), 900, 0, 90));
}

const int n = sizeof(rxy) / sizeof(rxy[0]);
for (int i = 0; i < n - 1; i++) {
float rmin = rxy[i];
float rmax = rxy[i + 1];
list->Add(new TH1F(Form("hConvPhi_Rxy%d_%dcm", static_cast<int>(rmin), static_cast<int>(rmax)), Form("conversion point of #varphi MC in %d < R_{xy} < %d cm;#varphi (rad.);N_{e}", static_cast<int>(rmin), static_cast<int>(rmax)), 360, 0.0f, TMath::TwoPi()));
}
// Generated, particles
if (TString(subGroup) == "Photon") {
list->Add(new TH1F("hPt_Photon", "photon pT;p_{T} (GeV/c)", 1000, 0.0f, 10));
list->Add(new TH1F("hY_Photon", "photon y;rapidity y", 40, -2.0f, 2.0f));
list->Add(new TH1F("hPhi_Photon", "photon #varphi;#varphi (rad.)", 180, 0, TMath::TwoPi()));
}

////Generated, particles
Expand Down
2 changes: 1 addition & 1 deletion PWGEM/PhotonMeson/DataModel/gammaTables.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(Y, y, //! Particle rapidity
});
} // namespace emmcparticle
// NOTE: This table is nearly identical to the one from Framework (except that it points to the event ID, not the BC id)
// This table contains all MC truth tracks (both barrel and muon)
// This table contains all MC truth tracks (both v0 and calos)
DECLARE_SOA_TABLE_FULL(EMMCParticles, "EMMCParticles", "AOD", "EMMCPARTICLE", //! MC track information (on disk)
o2::soa::Index<>, emmcparticle::EMReducedMCEventId,
mcparticle::PdgCode, mcparticle::StatusCode, mcparticle::Flags,
Expand Down
94 changes: 72 additions & 22 deletions PWGEM/PhotonMeson/TableProducer/createEMReducedMCEvent.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,15 +54,15 @@ struct createEMReducedMCEvent {
hEventCounter->GetXaxis()->SetBinLabel(1, "all");
}

Preslice<aod::McParticles_001> perMcCollision = aod::mcparticle::mcCollisionId;
Preslice<aod::McParticles> perMcCollision = aod::mcparticle::mcCollisionId;
Preslice<aod::V0Photons> perCollision_pcm = aod::v0photon::collisionId;
Preslice<aod::PHOSClusters> perCollision_phos = aod::skimmedcluster::collisionId;
Preslice<aod::SkimEMCClusters> perCollision_emc = aod::skimmedcluster::collisionId;

using MyCollisions = soa::Join<aod::Collisions, aod::Mults, aod::EvSels, aod::McCollisionLabels>;

template <uint8_t system, typename TTracks, typename TPCMs, typename TPCMLegs, typename TPHOSs, typename TEMCs>
void skimmingMC(MyCollisions const& collisions, aod::BCs const&, aod::McCollisions const&, aod::McParticles_001 const& mcTracks, TTracks const&, TPCMs const& v0photons, TPCMLegs const& v0legs, TPHOSs const& phosclusters, TEMCs const& emcclusters)
void skimmingMC(MyCollisions const& collisions, aod::BCs const&, aod::McCollisions const&, aod::McParticles const& mcTracks, TTracks const&, TPCMs const& v0photons, TPCMLegs const& v0legs, TPHOSs const& phosclusters, TEMCs const& emcclusters)
{
// temporary variables used for the indexing of the skimmed MC stack
std::map<uint64_t, int> fNewLabels;
Expand Down Expand Up @@ -132,26 +132,51 @@ struct createEMReducedMCEvent {
auto groupedMcTracks = mcTracks.sliceBy(perMcCollision, mcCollision.globalIndex());

for (auto& mctrack : groupedMcTracks) {
if (mctrack.pt() < 1e-2 || abs(mctrack.y()) > 1.5 || abs(mctrack.vz()) > 250 || sqrt(pow(mctrack.vx(), 2) + pow(mctrack.vy(), 2)) > 500) { // if pT < 10 MeV/c, don't store. Anyway, we don't care such low pT particles.
// if (mctrack.pt() < 1e-2 || abs(mctrack.vz()) > 250 || sqrt(pow(mctrack.vx(), 2) + pow(mctrack.vy(), 2)) > 500) { // if pT < 10 MeV/c, don't store. Anyway, we don't care such low pT particles.
if (mctrack.pt() < 1e-2 || abs(mctrack.y()) > 1.5 || abs(mctrack.vz()) > 250 || sqrt(pow(mctrack.vx(), 2) + pow(mctrack.vy(), 2)) > 260) {
continue;
}

int pdg = mctrack.pdgCode();
if (
abs(pdg) != 11 // electron
&& (abs(pdg) != 22 || !IsPhysicalPrimary(mctrack, mcTracks)) // photon
abs(pdg) != 11 // electron
&& (abs(pdg) != 22) // photon
// light mesons
&& (abs(pdg) != 111 || !IsPhysicalPrimary(mctrack, mcTracks)) // pi0
&& (abs(pdg) != 113 || !IsPhysicalPrimary(mctrack, mcTracks)) // rho(770)
&& (abs(pdg) != 211 || !IsPhysicalPrimary(mctrack, mcTracks)) // changed pion
&& (abs(pdg) != 221 || !IsPhysicalPrimary(mctrack, mcTracks)) // eta
&& (abs(pdg) != 223 || !IsPhysicalPrimary(mctrack, mcTracks)) // omega(782)
&& (abs(pdg) != 331 || !IsPhysicalPrimary(mctrack, mcTracks)) // eta'(958)
&& (abs(pdg) != 333 || !IsPhysicalPrimary(mctrack, mcTracks)) // phi(1020)
&& (abs(pdg) != 111) // pi0
&& (abs(pdg) != 113) // rho(770)
//&& (abs(pdg) != 211 ) // changed pion
&& (abs(pdg) != 221) // eta
&& (abs(pdg) != 223) // omega(782)
&& (abs(pdg) != 331) // eta'(958)
&& (abs(pdg) != 333) // phi(1020)
// strange hadrons
&& (abs(pdg) != 310) // K0S
&& (abs(pdg) != 3122) // Lambda
) {
continue;
}

// float dx = mctrack.vx() - mcCollision.posX();
// float dy = mctrack.vy() - mcCollision.posY();
// float dz = mctrack.vz() - mcCollision.posZ();
// float r3D = sqrt(dx*dx + dy*dy + dz*dz);
// if (
// abs(pdg) != 11 // electron
// && (abs(pdg) != 22 || r3D > 1.0f) // photon
// // light mesons
// && (abs(pdg) != 111 || r3D > 1.0f) // pi0
// && (abs(pdg) != 113 || r3D > 1.0f) // rho(770)
// && (abs(pdg) != 211 || r3D > 1.0f) // changed pion
// && (abs(pdg) != 221 || r3D > 1.0f) // eta
// && (abs(pdg) != 223 || r3D > 1.0f) // omega(782)
// && (abs(pdg) != 331 || r3D > 1.0f) // eta'(958)
// && (abs(pdg) != 333 || r3D > 1.0f) // phi(1020)
// //strange hadrons
// && (abs(pdg) != 310 || r3D > 1.0f) // K0S
// && (abs(pdg) != 3122 || r3D > 1.0f) // Lambda
//) {
// continue;
// }

// LOGF(info,"index = %d , mc track pdg = %d , producedByGenerator = %d , isPhysicalPrimary = %d", mctrack.index(), mctrack.pdgCode(), mctrack.producedByGenerator(), mctrack.isPhysicalPrimary());

// these are used as denominator for efficiency. (i.e. generated information)
Expand All @@ -172,7 +197,10 @@ struct createEMReducedMCEvent {

for (auto& leg : {pos, ele}) { // be carefull of order {pos, ele}!
auto o2track = leg.template track_as<TracksMC>();
auto mctrack = o2track.template mcParticle_as<aod::McParticles_001>();
if (!o2track.has_mcParticle()) {
continue; // If no MC particle is found, skip the track
}
auto mctrack = o2track.template mcParticle_as<aod::McParticles>();

// if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
if (!(fNewLabels.find(mctrack.globalIndex()) != fNewLabels.end())) {
Expand All @@ -196,6 +224,28 @@ struct createEMReducedMCEvent {

std::vector<int> mothers;
if (mctrack.has_mothers()) {
////auto mp = mctrack.template mothers_first_as<TMCs>();
// int motherid = mctrack.mothersIds()[0];//first mother index
// while(motherid > -1) {
// if (motherid < mcTracks.size()) { // protect against bad mother indices. why is this needed?
// auto mp = mcTracks.iteratorAt(motherid);

// if (fNewLabels.find(mp.globalIndex()) != fNewLabels.end()) {
// mothers.push_back(fNewLabels.find(mp.globalIndex())->second);
// }

// if(mp.has_mothers()){
// motherid = mp.mothersIds()[0];//first mother index
// } else {
// motherid = -999;
// }
// }
// else{
// std::cout << "Mother label (" << motherid << ") exceeds the McParticles size (" << mcTracks.size() << ")" << std::endl;
// std::cout << " Check the MC generator" << std::endl;
// }
//}

for (auto& m : mctrack.mothersIds()) {
if (m < mcTracks.size()) { // protect against bad mother indices
if (fNewLabels.find(m) != fNewLabels.end()) {
Expand Down Expand Up @@ -245,40 +295,40 @@ struct createEMReducedMCEvent {
fCounters[1] = 0;
} // end of skimmingMC

void processMC_PCM(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles_001 const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs)
void processMC_PCM(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs)
{
skimmingMC<kPCM>(collisions, bcs, mccollisions, mcTracks, tracks, v0photons, v0legs, nullptr, nullptr);
}
void processMC_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles_001 const& mcTracks, aod::PHOSClusters const& phosclusters)
void processMC_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, aod::PHOSClusters const& phosclusters)
{
skimmingMC<kPHOS>(collisions, bcs, mccollisions, mcTracks, nullptr, nullptr, nullptr, phosclusters, nullptr);
}
void processMC_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles_001 const& mcTracks, aod::SkimEMCClusters const& emcclusters)
void processMC_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, aod::SkimEMCClusters const& emcclusters)
{
skimmingMC<kEMC>(collisions, bcs, mccollisions, mcTracks, nullptr, nullptr, nullptr, nullptr, emcclusters);
}
void processMC_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles_001 const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::PHOSClusters const& phosclusters)
void processMC_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::PHOSClusters const& phosclusters)
{
const uint8_t sysflag = kPCM | kPHOS;
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, tracks, v0photons, v0legs, phosclusters, nullptr);
}
void processMC_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles_001 const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::SkimEMCClusters const& emcclusters)
void processMC_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::SkimEMCClusters const& emcclusters)
{
const uint8_t sysflag = kPCM | kEMC;
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, tracks, v0photons, v0legs, nullptr, emcclusters);
}
void processMC_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles_001 const& mcTracks, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
void processMC_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
{
const uint8_t sysflag = kPHOS | kEMC;
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, nullptr, nullptr, nullptr, phosclusters, emcclusters);
}
void processMC_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles_001 const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
void processMC_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, TracksMC const& tracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
{
const uint8_t sysflag = kPCM | kPHOS | kEMC;
skimmingMC<sysflag>(collisions, bcs, mccollisions, mcTracks, tracks, v0photons, v0legs, phosclusters, emcclusters);
}

void processDummy(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles_001 const& mcTracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs) {}
void processDummy(MyCollisions const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcTracks, aod::V0Photons const& v0photons, aod::V0Legs const& v0legs) {}

PROCESS_SWITCH(createEMReducedMCEvent, processMC_PCM, "create em mc event table for PCM", false);
PROCESS_SWITCH(createEMReducedMCEvent, processMC_PHOS, "create em mc event table for PHOS", false);
Expand Down
Loading