Skip to content

Commit

Permalink
PWGEM/PhotonMeson: update dilepton QC MC (AliceO2Group#3906)
Browse files Browse the repository at this point in the history
  • Loading branch information
@dsekihat
dsekihat authored Nov 20, 2023
1 parent 496e128 commit 9e00cf5
Show file tree
Hide file tree
Showing 5 changed files with 170 additions and 20 deletions.
16 changes: 9 additions & 7 deletions PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,11 +202,6 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
list->Add(new TH1F("hNpair_lspp", "Number of LS++ pairs per collision", 101, -0.5f, 100.5f));
list->Add(new TH1F("hNpair_lsmm", "Number of LS-- pairs per collision", 101, -0.5f, 100.5f));

if (TString(histClass) == "Generated") {
TH2F* hMvsPt = new TH2F("hMvsPt", "m_{ll} vs. p_{T,ll};m_{ll} (GeV/c^{2});p_{T,ll}", 110, 0, 1.1f, 100, 0, 10.f);
hMvsPt->Sumw2();
list->Add(hMvsPt);
}
} // end of Dalitz
if (TString(histClass) == "Track") {
float maxP = 10.f;
Expand Down Expand Up @@ -387,14 +382,12 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
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(), 200, 0, 200));
}

// 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
if (TString(subGroup) == "Pi0Eta") {
static constexpr std::string_view parnames[2] = {"Pi0", "Eta"};
for (int i = 0; i < 2; i++) {
Expand All @@ -403,6 +396,15 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
list->Add(new TH1F(Form("hPhi_%s", parnames[i].data()), Form("%s #varphi;#varphi (rad.)", parnames[i].data()), 180, 0, TMath::TwoPi()));
}
}
if (TString(subGroup) == "dielectron") {
TH2F* hMvsPt = new TH2F("hMvsPt", "m_{ee} vs. p_{T,ee};m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c)", 110, 0, 1.1f, 100, 0, 10.f);
hMvsPt->Sumw2();
list->Add(hMvsPt);
} else if (TString(subGroup) == "dimuon") {
TH2F* hMvsPt = new TH2F("hMvsPt", "m_{#mu#mu} vs. p_{T,#mu#mu};m_{#mu#mu} (GeV/c^{2});p_{T,#mu#mu} (GeV/c)", 90, 0.2, 1.1f, 20, 0, 2.f);
hMvsPt->Sumw2();
list->Add(hMvsPt);
}
}

const int nmgg04 = 201;
Expand Down
7 changes: 2 additions & 5 deletions PWGEM/PhotonMeson/TableProducer/associateMCinfo.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ struct AssociateMCInfo {
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)) > max_rxy_gen) {
if (mctrack.pt() < 1e-3 || abs(mctrack.y()) > 1.5 || abs(mctrack.vz()) > 250 || sqrt(pow(mctrack.vx(), 2) + pow(mctrack.vy(), 2)) > max_rxy_gen) {
continue;
}
int pdg = mctrack.pdgCode();
Expand Down Expand Up @@ -143,9 +143,6 @@ struct AssociateMCInfo {
auto ele = v0.template negTrack_as<aod::V0Legs>();
auto pos = v0.template posTrack_as<aod::V0Legs>();

// auto o2track_ele = ele.template track_as<TracksMC>();
// auto o2track_pos = pos.template track_as<TracksMC>();

auto o2track_ele = o2tracks.iteratorAt(pos.trackId());
auto o2track_pos = o2tracks.iteratorAt(ele.trackId());

Expand Down Expand Up @@ -192,7 +189,7 @@ struct AssociateMCInfo {
} // end of em primary track loop
}
if constexpr (static_cast<bool>(system & kDalitzMuMu)) {
// for dalitz ee
// for dalitz mumu
auto emprimarymuons_coll = emprimarymuons.sliceBy(perCollision_mu, collision.globalIndex());
for (auto& emprimarymuon : emprimarymuons_coll) {
auto o2track = o2tracks.iteratorAt(emprimarymuon.trackId());
Expand Down
33 changes: 25 additions & 8 deletions PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -347,13 +347,6 @@ struct Pi0EtaToGammaGammaMC {
int photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles);
int photonid2 = FindCommonMotherFrom2Prongs(pos2mc, ele2mc, -11, 11, 22, mcparticles);

// if (photonid1 < 0) { // check swap, true electron is reconstructed as positron and vice versa.
// photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, 11, -11, 22, mcparticles);
// }
// if (photonid2 < 0) { // check swap, true electron is reconstructed as positron and vice versa.
// photonid2 = FindCommonMotherFrom2Prongs(pos2mc, ele2mc, 11, -11, 22, mcparticles);
// }

// LOGF(info,"photonid1 = %d , photonid2 = %d", photonid1, photonid2);
if (photonid1 < 0 || photonid2 < 0) {
continue;
Expand Down Expand Up @@ -435,14 +428,38 @@ struct Pi0EtaToGammaGammaMC {
auto g1mc = mcparticles.iteratorAt(photonid1);
pi0id = FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -11, 11, 111, mcparticles);
etaid = FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -11, 11, 221, mcparticles);

} else if constexpr (pairtype == PairType::kPCMDalitzMuMu) { // check 4 legs
auto pos1 = g1.template posTrack_as<MyMCV0Legs>();
auto ele1 = g1.template negTrack_as<MyMCV0Legs>();
auto pos2 = g2.template posTrack_as<MyMCMuons>();
auto ele2 = g2.template negTrack_as<MyMCMuons>();
if (pos1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) {
continue;
}

auto pos1mc = pos1.template emmcparticle_as<aod::EMMCParticles>();
auto ele1mc = ele1.template emmcparticle_as<aod::EMMCParticles>();
auto pos2mc = pos2.template emmcparticle_as<aod::EMMCParticles>();
auto ele2mc = ele2.template emmcparticle_as<aod::EMMCParticles>();
// LOGF(info,"pos1mc.globalIndex() = %d , ele1mc.globalIndex() = %d , pos2mc.globalIndex() = %d , ele2mc.globalIndex() = %d", pos1mc.globalIndex(), ele1mc.globalIndex(), pos2mc.globalIndex(), ele2mc.globalIndex());

int photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); // real photon
if (photonid1 < 0) {
continue;
}
auto g1mc = mcparticles.iteratorAt(photonid1);
pi0id = -1;
etaid = FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -13, 13, 221, mcparticles);
}

if (pi0id < 0 && etaid < 0) {
continue;
}

ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
if constexpr (pairtype == PairType::kPCMDalitzEE) {
if constexpr (pairtype == PairType::kPCMDalitzEE || pairtype == PairType::kPCMDalitzMuMu) {
v2.SetM(g2.mass());
}
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
Expand Down
66 changes: 66 additions & 0 deletions PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ struct DalitzEEQCMC {
OutputObj<THashList> fOutputEvent{"Event"};
OutputObj<THashList> fOutputTrack{"Track"};
OutputObj<THashList> fOutputDalitzEE{"DalitzEE"};
OutputObj<THashList> fOutputGen{"Generated"};
THashList* fMainList = new THashList();

void addhistograms()
Expand All @@ -71,6 +72,10 @@ struct DalitzEEQCMC {
o2::aod::emphotonhistograms::AddHistClass(fMainList, "DalitzEE");
THashList* list_dalitzee = reinterpret_cast<THashList*>(fMainList->FindObject("DalitzEE"));

o2::aod::emphotonhistograms::AddHistClass(fMainList, "Generated");
THashList* list_gen = reinterpret_cast<THashList*>(fMainList->FindObject("Generated"));
o2::aod::emphotonhistograms::DefineHistograms(list_gen, "Generated", "dielectron");

for (const auto& cut : fDalitzEECuts) {
const char* cutname = cut.GetName();
o2::aod::emphotonhistograms::AddHistClass(list_track, cutname);
Expand Down Expand Up @@ -114,6 +119,7 @@ struct DalitzEEQCMC {
fOutputEvent.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Event")));
fOutputTrack.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Track")));
fOutputDalitzEE.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("DalitzEE")));
fOutputGen.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Generated")));
}

template <typename TTrack, typename TMCTracks>
Expand Down Expand Up @@ -220,6 +226,66 @@ struct DalitzEEQCMC {
} // end of process
PROCESS_SWITCH(DalitzEEQCMC, processQCMC, "run Dalitz QC", true);

Configurable<float> min_mcPt{"min_mcPt", 0.05, "min. MC pT"};
Configurable<float> max_mcPt{"max_mcPt", 1e+10, "max. MC pT"};
Configurable<float> max_mcEta{"max_mcEta", 0.9, "max. MC eta"};
Partition<aod::EMMCParticles> posTracks = o2::aod::mcparticle::pdgCode == -11; // e+
Partition<aod::EMMCParticles> negTracks = o2::aod::mcparticle::pdgCode == +11; // e-
PresliceUnsorted<aod::EMMCParticles> perMcCollision = aod::emmcparticle::emreducedmceventId;
void processGen(MyCollisions const& collisions, aod::EMReducedMCEvents const&, aod::EMMCParticles const& mcparticles)
{
// 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& collision : collisions) {
auto mccollision = collision.emreducedmcevent();
// LOGF(info, "mccollision.globalIndex() = %d", mccollision.globalIndex());
// auto mctracks_coll = mcparticles.sliceBy(perMcCollision, mccollision.globalIndex());

reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(1.0);
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hZvtx_before"))->Fill(mccollision.posZ());
if (!collision.sel8()) {
continue;
}
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(2.0);

if (collision.numContrib() < 0.5) {
continue;
}
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(3.0);

if (abs(collision.posZ()) > 10.0) {
continue;
}
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(4.0);
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hZvtx_after"))->Fill(mccollision.posZ());
auto posTracks_per_coll = posTracks->sliceByCachedUnsorted(o2::aod::emmcparticle::emreducedmceventId, mccollision.globalIndex(), cache);
auto negTracks_per_coll = negTracks->sliceByCachedUnsorted(o2::aod::emmcparticle::emreducedmceventId, mccollision.globalIndex(), cache);
for (auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) {
// LOGF(info, "pdg1 = %d, pdg2 = %d", t1.pdgCode(), t2.pdgCode());

if (t1.pt() < min_mcPt || max_mcPt < t1.pt() || abs(t1.eta()) > max_mcEta) {
continue;
}
if (t2.pt() < min_mcPt || max_mcPt < t2.pt() || abs(t2.eta()) > max_mcEta) {
continue;
}

int mother_id = FindLF(t1, t2, mcparticles);
if (mother_id < 0) {
continue;
}
auto mcmother = mcparticles.iteratorAt(mother_id);
if (IsPhysicalPrimary(mcmother.emreducedmcevent(), mcmother, mcparticles)) {
ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron);
ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron);
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
reinterpret_cast<TH2F*>(fMainList->FindObject("Generated")->FindObject("hMvsPt"))->Fill(v12.M(), v12.Pt());
} // end of LF
} // end of true ULS pair loop
} // end of collision loop
}
PROCESS_SWITCH(DalitzEEQCMC, processGen, "run genrated info", true);

void processDummy(MyCollisions const& collisions) {}
PROCESS_SWITCH(DalitzEEQCMC, processDummy, "Dummy function", false);
};
Expand Down
68 changes: 68 additions & 0 deletions PWGEM/PhotonMeson/Tasks/dalitzMuMuQCMC.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ struct DalitzMuMuQCMC {
OutputObj<THashList> fOutputEvent{"Event"};
OutputObj<THashList> fOutputTrack{"Track"};
OutputObj<THashList> fOutputDalitzMuMu{"DalitzMuMu"};
OutputObj<THashList> fOutputGen{"Generated"};
THashList* fMainList = new THashList();

void addhistograms()
Expand All @@ -71,6 +72,10 @@ struct DalitzMuMuQCMC {
o2::aod::emphotonhistograms::AddHistClass(fMainList, "DalitzMuMu");
THashList* list_dalitzmumu = reinterpret_cast<THashList*>(fMainList->FindObject("DalitzMuMu"));

o2::aod::emphotonhistograms::AddHistClass(fMainList, "Generated");
THashList* list_gen = reinterpret_cast<THashList*>(fMainList->FindObject("Generated"));
o2::aod::emphotonhistograms::DefineHistograms(list_gen, "Generated", "dimuon");

for (const auto& cut : fDalitzMuMuCuts) {
const char* cutname = cut.GetName();
o2::aod::emphotonhistograms::AddHistClass(list_track, cutname);
Expand Down Expand Up @@ -114,6 +119,7 @@ struct DalitzMuMuQCMC {
fOutputEvent.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Event")));
fOutputTrack.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Track")));
fOutputDalitzMuMu.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("DalitzMuMu")));
fOutputGen.setObject(reinterpret_cast<THashList*>(fMainList->FindObject("Generated")));
}

template <typename TTrack, typename TMCTracks>
Expand Down Expand Up @@ -177,6 +183,8 @@ struct DalitzMuMuQCMC {
auto posmc = pos.template emmcparticle_as<aod::EMMCParticles>();
auto elemc = ele.template emmcparticle_as<aod::EMMCParticles>();

// LOGF(info, "posmc.pdgCode() = %d, , posmc.has_mothers()() = %d | elemc.pdgCode() = %d, elemc.has_mothers() = %d", posmc.pdgCode(), posmc.has_mothers(), elemc.pdgCode(), elemc.has_mothers());

int mother_id = FindLF(posmc, elemc, mcparticles);
if (mother_id < 0) {
continue;
Expand Down Expand Up @@ -209,6 +217,66 @@ struct DalitzMuMuQCMC {
} // end of process
PROCESS_SWITCH(DalitzMuMuQCMC, processQCMC, "run Dalitz QC", true);

Configurable<float> min_mcPt{"min_mcPt", 0.05, "min. MC pT"};
Configurable<float> max_mcPt{"max_mcPt", 0.6, "max. MC pT"};
Configurable<float> max_mcEta{"max_mcEta", 0.9, "max. MC eta"};
Partition<aod::EMMCParticles> posTracks = o2::aod::mcparticle::pdgCode == -13; // mu+
Partition<aod::EMMCParticles> negTracks = o2::aod::mcparticle::pdgCode == +13; // mu-
PresliceUnsorted<aod::EMMCParticles> perMcCollision = aod::emmcparticle::emreducedmceventId;
void processGen(MyCollisions const& collisions, aod::EMReducedMCEvents const&, aod::EMMCParticles const& mcparticles)
{
// 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& collision : collisions) {
auto mccollision = collision.emreducedmcevent();
// LOGF(info, "mccollision.globalIndex() = %d", mccollision.globalIndex());
// auto mctracks_coll = mcparticles.sliceBy(perMcCollision, mccollision.globalIndex());

reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(1.0);
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hZvtx_before"))->Fill(mccollision.posZ());
if (!collision.sel8()) {
continue;
}
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(2.0);

if (collision.numContrib() < 0.5) {
continue;
}
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(3.0);

if (abs(collision.posZ()) > 10.0) {
continue;
}
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(4.0);
reinterpret_cast<TH1F*>(fMainList->FindObject("Generated")->FindObject("hZvtx_after"))->Fill(mccollision.posZ());
auto posTracks_per_coll = posTracks->sliceByCachedUnsorted(o2::aod::emmcparticle::emreducedmceventId, mccollision.globalIndex(), cache);
auto negTracks_per_coll = negTracks->sliceByCachedUnsorted(o2::aod::emmcparticle::emreducedmceventId, mccollision.globalIndex(), cache);
for (auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) {
// LOGF(info, "pdg1 = %d, pdg2 = %d", t1.pdgCode(), t2.pdgCode());

if (t1.pt() < min_mcPt || max_mcPt < t1.pt() || abs(t1.eta()) > max_mcEta) {
continue;
}
if (t2.pt() < min_mcPt || max_mcPt < t2.pt() || abs(t2.eta()) > max_mcEta) {
continue;
}

int mother_id = FindLF(t1, t2, mcparticles);
if (mother_id < 0) {
continue;
}
auto mcmother = mcparticles.iteratorAt(mother_id);
if (IsPhysicalPrimary(mcmother.emreducedmcevent(), mcmother, mcparticles)) {
ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassMuon);
ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassMuon);
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
reinterpret_cast<TH2F*>(fMainList->FindObject("Generated")->FindObject("hMvsPt"))->Fill(v12.M(), v12.Pt());
} // end of LF
} // end of true ULS pair loop
} // end of collision loop
}
PROCESS_SWITCH(DalitzMuMuQCMC, processGen, "run genrated info", true);

void processDummy(MyCollisions const& collisions) {}
PROCESS_SWITCH(DalitzMuMuQCMC, processDummy, "Dummy function", false);
};
Expand Down

0 comments on commit 9e00cf5

Please sign in to comment.