Skip to content

Commit 28fd47e

Browse files
dsekihatdsekihat
authored
PWGEM/PhotonMeson: quick fix to 6283 (#6284)
1 parent b36803f commit 28fd47e

2 files changed

Lines changed: 94 additions & 72 deletions

File tree

‎PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGamma.cxx‎

Lines changed: 93 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@ struct Pi0EtaToGammaGamma {
8888
Configurable<int> ndepth{"ndepth", 10, "depth for event mixing"};
8989
ConfigurableAxis ConfVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"};
9090
ConfigurableAxis ConfCentBins{"ConfCentBins", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.f, 999.f}, "Mixing bins - centrality"};
91-
ConfigurableAxis ConfEPBins{"ConfEPBins", {VARIABLE_WIDTH, 0.0f, M_PI / 4, M_PI / 2, M_PI}, "Mixing bins - event plane angle"};
91+
ConfigurableAxis ConfEPBins{"ConfEPBins", {VARIABLE_WIDTH, -M_PI / 2, -M_PI / 4, 0.0f, +M_PI / 4, +M_PI / 2}, "Mixing bins - event plane angle"};
9292

9393
EMEventCut fEMEventCut;
9494
Configurable<float> cfgZvtxMax{"cfgZvtxMax", 10.f, "max. Zvtx"};
@@ -246,6 +246,8 @@ struct Pi0EtaToGammaGamma {
246246

247247
used_photonIds.clear();
248248
used_photonIds.shrink_to_fit();
249+
used_dileptonIds.clear();
250+
used_dileptonIds.shrink_to_fit();
249251
}
250252

251253
void DefineEMEventCut()
@@ -395,44 +397,6 @@ struct Pi0EtaToGammaGamma {
395397
fPHOSCut.SetEnergyRange(phoscuts.cfg_min_Ecluster, 1e+10);
396398
}
397399

398-
template <typename TCollision, typename TTracks, typename TCut>
399-
std::vector<EMTrack> createULSPairs(TCollision const& collision, TTracks const& posTracks, TTracks const& negTracks, TCut const& cut, const float mlepton)
400-
{
401-
std::vector<EMTrack> pairs;
402-
pairs.reserve(posTracks.size() * negTracks.size());
403-
404-
for (auto& [pos, ele] : combinations(CombinationsFullIndexPolicy(posTracks, negTracks))) { // ULS
405-
406-
if (pos.trackId() == ele.trackId()) { // this is protection against pairing identical 2 tracks.
407-
continue;
408-
}
409-
410-
if (dileptoncuts.cfg_pid_scheme == static_cast<int>(DalitzEECut::PIDSchemes::kPIDML)) {
411-
if (!cut.template IsSelectedTrack<true>(pos, collision) || !cut.template IsSelectedTrack<true>(ele, collision)) {
412-
continue;
413-
}
414-
} else { // cut-based
415-
if (!cut.template IsSelectedTrack<false>(pos, collision) || !cut.template IsSelectedTrack<false>(ele, collision)) {
416-
continue;
417-
}
418-
}
419-
420-
if (!cut.IsSelectedPair(pos, ele, collision.bz())) {
421-
continue;
422-
}
423-
424-
ROOT::Math::PtEtaPhiMVector v1(pos.pt(), pos.eta(), pos.phi(), mlepton);
425-
ROOT::Math::PtEtaPhiMVector v2(ele.pt(), ele.eta(), ele.phi(), mlepton);
426-
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
427-
float dca_pos_3d = pos.dca3DinSigma();
428-
float dca_ele_3d = ele.dca3DinSigma();
429-
float dca_ee_3d = std::sqrt((dca_pos_3d * dca_pos_3d + dca_ele_3d * dca_ele_3d) / 2.);
430-
pairs.emplace_back(EMTrack(collision.globalIndex(), ndilepton, v12.Pt(), v12.Eta(), v12.Phi(), v12.M(), 0, dca_ee_3d));
431-
ndilepton++;
432-
}
433-
return pairs;
434-
}
435-
436400
/// \brief Calculate background (using rotation background method only for EMCal!)
437401
template <typename TPhotons>
438402
void RotationBackground(const ROOT::Math::PtEtaPhiMVector& meson, ROOT::Math::PtEtaPhiMVector photon1, ROOT::Math::PtEtaPhiMVector photon2, TPhotons const& photons_coll, unsigned int ig1, unsigned int ig2, EMCPhotonCut const& cut, aod::SkimEMCMTs const&)
@@ -490,6 +454,7 @@ struct Pi0EtaToGammaGamma {
490454
std::map<std::tuple<int, int, int>, std::vector<std::pair<int, int64_t>>> map_mix_bins; // zvtx, centrality, event plane bins -> pair<df index, vector of collisionId>
491455
std::map<std::pair<int, int64_t>, std::vector<EMTrack>> map_photons_to_collision; // pair<df index, collisionId> -> vector of track
492456
std::vector<std::pair<int, int>> used_photonIds; // <ndf, trackId>
457+
std::vector<std::tuple<int, int, int>> used_dileptonIds; // <ndf, trackId>
493458

494459
template <typename TCollisions, typename TPhotons1, typename TPhotons2, typename TSubInfos1, typename TSubInfos2, typename TPreslice1, typename TPreslice2, typename TCut1, typename TCut2, typename TTracksMatchedWithEMC, typename TTracksMatchedWithPHOS>
495460
void runPairing(TCollisions const& collisions,
@@ -585,66 +550,125 @@ struct Pi0EtaToGammaGamma {
585550
auto photons1_per_collision = photons1.sliceBy(perCollision1, collision.globalIndex());
586551
auto positrons_per_collision = positrons->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache);
587552
auto electrons_per_collision = electrons->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache);
588-
auto photons2_per_collision = createULSPairs(collision, positrons_per_collision, electrons_per_collision, cut2, o2::constants::physics::MassElectron);
589553

590554
for (auto& g1 : photons1_per_collision) {
591-
for (auto& g2 : photons2_per_collision) {
592-
if (!cut1.template IsSelected<TSubInfos1>(g1)) {
555+
if (!cut1.template IsSelected<TSubInfos1>(g1)) {
556+
continue;
557+
}
558+
auto pos1 = g1.template posTrack_as<TSubInfos1>();
559+
auto ele1 = g1.template negTrack_as<TSubInfos1>();
560+
ROOT::Math::PtEtaPhiMVector v_gamma(g1.pt(), g1.eta(), g1.phi(), 0.);
561+
562+
for (auto& [pos2, ele2] : combinations(CombinationsFullIndexPolicy(positrons_per_collision, electrons_per_collision))) {
563+
564+
if (pos2.trackId() == ele2.trackId()) { // this is protection against pairing identical 2 tracks.
593565
continue;
594566
}
595-
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
596-
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), g2.mass());
597-
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
598-
if (abs(v12.Rapidity()) > maxY) {
567+
if (pos1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) {
599568
continue;
600569
}
601-
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillPairInfo<0, pairtype>(&fRegistry, collision, v12, cfgDoFlow);
570+
571+
if (dileptoncuts.cfg_pid_scheme == static_cast<int>(DalitzEECut::PIDSchemes::kPIDML)) {
572+
if (!cut2.template IsSelectedTrack<true>(pos2, collision) || !cut2.template IsSelectedTrack<true>(ele2, collision)) {
573+
continue;
574+
}
575+
} else { // cut-based
576+
if (!cut2.template IsSelectedTrack<false>(pos2, collision) || !cut2.template IsSelectedTrack<false>(ele2, collision)) {
577+
continue;
578+
}
579+
}
580+
581+
if (!cut2.IsSelectedPair(pos2, ele2, collision.bz())) {
582+
continue;
583+
}
584+
585+
ROOT::Math::PtEtaPhiMVector v_pos(pos2.pt(), pos2.eta(), pos2.phi(), o2::constants::physics::MassElectron);
586+
ROOT::Math::PtEtaPhiMVector v_ele(ele2.pt(), ele2.eta(), ele2.phi(), o2::constants::physics::MassElectron);
587+
ROOT::Math::PtEtaPhiMVector v_ee = v_pos + v_ele;
588+
ROOT::Math::PtEtaPhiMVector veeg = v_gamma + v_pos + v_ele;
589+
if (abs(veeg.Rapidity()) > maxY) {
590+
continue;
591+
}
592+
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillPairInfo<0, pairtype>(&fRegistry, collision, veeg, cfgDoFlow);
593+
float dca_pos_3d = pos2.dca3DinSigma();
594+
float dca_ele_3d = ele2.dca3DinSigma();
595+
float dca_ee_3d = std::sqrt((dca_pos_3d * dca_pos_3d + dca_ele_3d * dca_ele_3d) / 2.);
602596

603597
std::pair<int, int> pair_tmp_id1 = std::make_pair(ndf, g1.globalIndex());
604-
std::pair<int, int> pair_tmp_id2 = std::make_pair(ndf, g2.trackId());
598+
std::tuple<int, int, int> tuple_tmp_id2 = std::make_tuple(ndf, pos2.trackId(), ele2.trackId());
605599
if (std::find(used_photonIds.begin(), used_photonIds.end(), pair_tmp_id1) == used_photonIds.end()) {
606600
emh1->AddTrackToEventPool(key_df_collision, EMTrack(collision.globalIndex(), g1.globalIndex(), g1.pt(), g1.eta(), g1.phi(), 0, 0, 0));
607601
used_photonIds.emplace_back(pair_tmp_id1);
608602
}
609-
if (std::find(used_photonIds.begin(), used_photonIds.end(), pair_tmp_id2) == used_photonIds.end()) {
610-
emh2->AddTrackToEventPool(key_df_collision, EMTrack(collision.globalIndex(), g2.trackId(), g2.pt(), g2.eta(), g2.phi(), g2.mass(), g2.sign(), g2.dca3DinSigma()));
611-
used_photonIds.emplace_back(pair_tmp_id2);
603+
if (std::find(used_dileptonIds.begin(), used_dileptonIds.end(), tuple_tmp_id2) == used_dileptonIds.end()) {
604+
emh2->AddTrackToEventPool(key_df_collision, EMTrack(collision.globalIndex(), -1, v_ee.Pt(), v_ee.Eta(), v_ee.Phi(), v_ee.M(), 0, dca_ee_3d));
605+
used_dileptonIds.emplace_back(tuple_tmp_id2);
612606
}
613607
ndiphoton++;
614-
} // end of g2 loop
608+
} // end of dielectron loop
615609
} // end of g1 loop
616610
} else if constexpr (pairtype == PairType::kPCMDalitzMuMu) {
617611
auto photons1_per_collision = photons1.sliceBy(perCollision1, collision.globalIndex());
618612
auto muons_pos_per_collision = muons_pos->sliceByCached(o2::aod::emprimarymuon::emeventId, collision.globalIndex(), cache);
619613
auto muons_neg_per_collision = muons_neg->sliceByCached(o2::aod::emprimarymuon::emeventId, collision.globalIndex(), cache);
620-
auto photons2_per_collision = createULSPairs(collision, muons_pos_per_collision, muons_neg_per_collision, cut2, o2::constants::physics::MassMuon);
621614

622615
for (auto& g1 : photons1_per_collision) {
623-
for (auto& g2 : photons2_per_collision) {
624-
if (!cut1.template IsSelected<TSubInfos1>(g1)) {
616+
if (!cut1.template IsSelected<TSubInfos1>(g1)) {
617+
continue;
618+
}
619+
auto pos1 = g1.template posTrack_as<TSubInfos1>();
620+
auto ele1 = g1.template negTrack_as<TSubInfos1>();
621+
ROOT::Math::PtEtaPhiMVector v_gamma(g1.pt(), g1.eta(), g1.phi(), 0.);
622+
623+
for (auto& [muplus, muminus] : combinations(CombinationsFullIndexPolicy(muons_pos_per_collision, muons_neg_per_collision))) {
624+
625+
if (muplus.trackId() == muminus.trackId()) { // this is protection against pairing identical 2 tracks.
625626
continue;
626627
}
627-
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
628-
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), g2.mass());
629-
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
630-
if (abs(v12.Rapidity()) > maxY) {
628+
if (pos1.trackId() == muplus.trackId() || ele1.trackId() == muminus.trackId()) {
631629
continue;
632630
}
633-
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillPairInfo<0, pairtype>(&fRegistry, collision, v12, cfgDoFlow);
631+
632+
if (dileptoncuts.cfg_pid_scheme == static_cast<int>(DalitzEECut::PIDSchemes::kPIDML)) {
633+
if (!cut2.template IsSelectedTrack<true>(muplus, collision) || !cut2.template IsSelectedTrack<true>(muminus, collision)) {
634+
continue;
635+
}
636+
} else { // cut-based
637+
if (!cut2.template IsSelectedTrack<false>(muplus, collision) || !cut2.template IsSelectedTrack<false>(muminus, collision)) {
638+
continue;
639+
}
640+
}
641+
642+
if (!cut2.IsSelectedPair(muplus, muminus, collision.bz())) {
643+
continue;
644+
}
645+
646+
ROOT::Math::PtEtaPhiMVector v_pos(muplus.pt(), muplus.eta(), muplus.phi(), o2::constants::physics::MassElectron);
647+
ROOT::Math::PtEtaPhiMVector v_ele(muminus.pt(), muminus.eta(), muminus.phi(), o2::constants::physics::MassElectron);
648+
ROOT::Math::PtEtaPhiMVector v_ee = v_pos + v_ele;
649+
ROOT::Math::PtEtaPhiMVector veeg = v_gamma + v_pos + v_ele;
650+
if (abs(veeg.Rapidity()) > maxY) {
651+
continue;
652+
}
653+
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillPairInfo<0, pairtype>(&fRegistry, collision, veeg, cfgDoFlow);
654+
float dca_pos_3d = muplus.dca3DinSigma();
655+
float dca_ele_3d = muminus.dca3DinSigma();
656+
float dca_ee_3d = std::sqrt((dca_pos_3d * dca_pos_3d + dca_ele_3d * dca_ele_3d) / 2.);
634657

635658
std::pair<int, int> pair_tmp_id1 = std::make_pair(ndf, g1.globalIndex());
636-
std::pair<int, int> pair_tmp_id2 = std::make_pair(ndf, g2.trackId());
659+
std::tuple<int, int, int> tuple_tmp_id2 = std::make_tuple(ndf, muplus.trackId(), muminus.trackId());
637660
if (std::find(used_photonIds.begin(), used_photonIds.end(), pair_tmp_id1) == used_photonIds.end()) {
638661
emh1->AddTrackToEventPool(key_df_collision, EMTrack(collision.globalIndex(), g1.globalIndex(), g1.pt(), g1.eta(), g1.phi(), 0, 0, 0));
639662
used_photonIds.emplace_back(pair_tmp_id1);
640663
}
641-
if (std::find(used_photonIds.begin(), used_photonIds.end(), pair_tmp_id2) == used_photonIds.end()) {
642-
emh2->AddTrackToEventPool(key_df_collision, EMTrack(collision.globalIndex(), g2.trackId(), g2.pt(), g2.eta(), g2.phi(), g2.mass(), g2.sign(), g2.dca3DinSigma()));
643-
used_photonIds.emplace_back(pair_tmp_id2);
664+
if (std::find(used_dileptonIds.begin(), used_dileptonIds.end(), tuple_tmp_id2) == used_dileptonIds.end()) {
665+
emh2->AddTrackToEventPool(key_df_collision, EMTrack(collision.globalIndex(), -1, v_ee.Pt(), v_ee.Eta(), v_ee.Phi(), v_ee.M(), 0, dca_ee_3d));
666+
used_dileptonIds.emplace_back(tuple_tmp_id2);
644667
}
645668
ndiphoton++;
646-
} // end of g2 loop
647-
} // end of g1 loop
669+
} // end of dimuon loop
670+
} // end of g1 loop
671+
648672
} else { // PCM-EMC, PCM-PHOS. Nightmare. don't run these pairs.
649673
auto photons1_per_collision = photons1.sliceBy(perCollision1, collision.globalIndex());
650674
auto photons2_per_collision = photons2.sliceBy(perCollision2, collision.globalIndex());
@@ -714,7 +738,8 @@ struct Pi0EtaToGammaGamma {
714738
o2::aod::pwgem::photonmeson::utils::nmhistogram::fillPairInfo<1, pairtype>(&fRegistry, collision, v12, cfgDoFlow);
715739
}
716740
}
717-
} // end of loop over mixed event pool
741+
} // end of loop over mixed event pool
742+
718743
} else { //[photon1 from event1, photon2 from event2] and [photon1 from event2, photon2 from event1]
719744

720745
for (auto& mix_dfId_collisionId : collisionIds2_in_mixing_pool) {
@@ -732,11 +757,10 @@ struct Pi0EtaToGammaGamma {
732757
for (auto& g2 : photons2_from_event_pool) {
733758
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
734759
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
735-
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
736-
737760
if constexpr (pairtype == PairType::kPCMDalitzEE || pairtype == PairType::kPCMDalitzMuMu) { //[photon from event1, dilepton from event2] and [photon from event2, dilepton from event1]
738761
v2.SetM(g2.mass());
739762
}
763+
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
740764
if (abs(v12.Rapidity()) > maxY) {
741765
continue;
742766
}
@@ -785,10 +809,8 @@ struct Pi0EtaToGammaGamma {
785809
using FilteredMyCollisions = soa::Filtered<MyCollisions>;
786810

787811
int ndf = 0;
788-
int ndilepton = 0;
789812
void processAnalysis(FilteredMyCollisions const& collisions, Types const&... args)
790813
{
791-
ndilepton = 0;
792814
// LOGF(info, "ndf = %d", ndf);
793815
if constexpr (pairtype == PairType::kPCMPCM) {
794816
auto v0photons = std::get<0>(std::tie(args...));

‎PWGEM/PhotonMeson/Utils/NMHistograms.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@ void addNMHistograms(HistogramRegistry* fRegistry, bool doFlow, bool isMC, const
4040
const AxisSpec axis_pt{ptbins, Form("p_{T,%s} (GeV/c)", pairname)};
4141

4242
const AxisSpec axis_mass{400, 0, 0.8, Form("m_{%s} (GeV/c^{2})", pairname)};
43-
const AxisSpec sp2_mass{50, -5, 5, Form("u_{2}^{%s} #upoint Q_{2}", pairname)};
43+
const AxisSpec sp2_mass{100, -5, 5, Form("u_{2}^{%s} #upoint Q_{2}", pairname)};
4444

4545
if (isMC) {
4646
fRegistry->add("Generated/Pi0/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{2000, 0.0f, 20}});

0 commit comments

Comments
 (0)