diff --git a/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx b/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx index 596270ca1f3..50174ab9e08 100644 --- a/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx +++ b/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx @@ -204,19 +204,19 @@ void o2::aod::pwgem::photon::histogram::DefineHistograms(THashList* list, const } if (TString(histClass).Contains("EE")) { - const int nm = 167; + const int nm = 145; double mee[nm] = {0.f}; for (int i = 0; i < 110; i++) { mee[i] = 0.01 * (i - 0) + 0.0; // every 0.01 GeV/c2 up to 1.1 GeV/c2 } - for (int i = 110; i < 128; i++) { - mee[i] = 0.1 * (i - 110) + 1.1; // every 0.1 GeV/c2 from 1.1 to 2.9 GeV/c2 + for (int i = 110; i < 126; i++) { + mee[i] = 0.1 * (i - 110) + 1.1; // every 0.1 GeV/c2 from 1.1 to 2.7 GeV/c2 } - for (int i = 128; i < 158; i++) { - mee[i] = 0.01 * (i - 128) + 2.9; // every 0.01 GeV/c2 from 2.9 to 3.2 GeV/c2 + for (int i = 126; i < 136; i++) { + mee[i] = 0.05 * (i - 126) + 2.7; // every 0.05 GeV/c2 from 2.7 to 3.2 GeV/c2 } - for (int i = 158; i < nm; i++) { - mee[i] = 0.1 * (i - 158) + 3.2; // every 0.01 GeV/c2 from 3.2 to 3.5 GeV/c2 + for (int i = 136; i < nm; i++) { + mee[i] = 0.1 * (i - 136) + 3.2; // every 0.1 GeV/c2 from 3.2 to 4.0 GeV/c2 } const int npt = 61; @@ -228,39 +228,45 @@ void o2::aod::pwgem::photon::histogram::DefineHistograms(THashList* list, const pt[i] = 0.5 * (i - 50) + 5.0; } - const int ndim = 4; // m, pt, dca, phiv - const int nbins[ndim] = {nm - 1, npt - 1, ndca - 1, 18}; - const double xmin[ndim] = {0.0, 0.0, 0.0, 0.0}; - const double xmax[ndim] = {4.0, 10.0, 5.0, M_PI}; + const int ndim = 3; // m, pt, dca + const int nbins[ndim] = {nm - 1, npt - 1, ndca - 1}; + const double xmin[ndim] = {0.0, 0.0, 0.0}; + const double xmax[ndim] = {4.0, 10.0, 5.0}; - hs_dilepton_uls_same = new THnSparseF("hs_dilepton_uls_same", "hs_dilepton_uls;m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c);DCA_{ee}^{3D} (#sigma);#varphi_{V} (rad.);", ndim, nbins, xmin, xmax); + hs_dilepton_uls_same = new THnSparseF("hs_dilepton_uls_same", "hs_dilepton_uls;m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c);DCA_{ee}^{3D} (#sigma);", ndim, nbins, xmin, xmax); hs_dilepton_uls_same->SetBinEdges(0, mee); hs_dilepton_uls_same->SetBinEdges(1, pt); hs_dilepton_uls_same->SetBinEdges(2, dca); hs_dilepton_uls_same->Sumw2(); list->Add(hs_dilepton_uls_same); - hs_dilepton_lspp_same = new THnSparseF("hs_dilepton_lspp_same", "hs_dilepton_lspp;m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c);DCA_{ee}^{3D} (#sigma);#varphi_{V} (rad.);", ndim, nbins, xmin, xmax); + hs_dilepton_lspp_same = new THnSparseF("hs_dilepton_lspp_same", "hs_dilepton_lspp;m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c);DCA_{ee}^{3D} (#sigma);", ndim, nbins, xmin, xmax); hs_dilepton_lspp_same->SetBinEdges(0, mee); hs_dilepton_lspp_same->SetBinEdges(1, pt); hs_dilepton_lspp_same->SetBinEdges(2, dca); hs_dilepton_lspp_same->Sumw2(); list->Add(hs_dilepton_lspp_same); - hs_dilepton_lsmm_same = new THnSparseF("hs_dilepton_lsmm_same", "hs_dilepton_lsmm;m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c);DCA_{ee}^{3D} (#sigma);#varphi_{V} (rad.);", ndim, nbins, xmin, xmax); + hs_dilepton_lsmm_same = new THnSparseF("hs_dilepton_lsmm_same", "hs_dilepton_lsmm;m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c);DCA_{ee}^{3D} (#sigma);", ndim, nbins, xmin, xmax); hs_dilepton_lsmm_same->SetBinEdges(0, mee); hs_dilepton_lsmm_same->SetBinEdges(1, pt); hs_dilepton_lsmm_same->SetBinEdges(2, dca); hs_dilepton_lsmm_same->Sumw2(); list->Add(hs_dilepton_lsmm_same); + TH2F* hMvsPhiV_uls_same = new TH2F("hMvsPhiV_uls_same", "m_{ee} vs. #varphi_{V};#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", 72, 0, M_PI, 100, 0.0f, 0.1f); + hMvsPhiV_uls_same->Sumw2(); + list->Add(hMvsPhiV_uls_same); + list->Add(reinterpret_cast(hMvsPhiV_uls_same->Clone("hMvsPhiV_lspp_same"))); + list->Add(reinterpret_cast(hMvsPhiV_uls_same->Clone("hMvsPhiV_lsmm_same"))); + if (TString(subGroup).Contains("mix")) { - THnSparseF* hs_dilepton_uls_mix = reinterpret_cast(hs_dilepton_uls_same->Clone("hs_dilepton_uls_mix")); - THnSparseF* hs_dilepton_lspp_mix = reinterpret_cast(hs_dilepton_lspp_same->Clone("hs_dilepton_lspp_mix")); - THnSparseF* hs_dilepton_lsmm_mix = reinterpret_cast(hs_dilepton_lsmm_same->Clone("hs_dilepton_lsmm_mix")); - list->Add(hs_dilepton_uls_mix); - list->Add(hs_dilepton_lspp_mix); - list->Add(hs_dilepton_lsmm_mix); + list->Add(reinterpret_cast(hs_dilepton_uls_same->Clone("hs_dilepton_uls_mix"))); + list->Add(reinterpret_cast(hs_dilepton_lspp_same->Clone("hs_dilepton_lspp_mix"))); + list->Add(reinterpret_cast(hs_dilepton_lsmm_same->Clone("hs_dilepton_lsmm_mix"))); + list->Add(reinterpret_cast(hMvsPhiV_uls_same->Clone("hMvsPhiV_uls_mix"))); + list->Add(reinterpret_cast(hMvsPhiV_uls_same->Clone("hMvsPhiV_lspp_mix"))); + list->Add(reinterpret_cast(hMvsPhiV_uls_same->Clone("hMvsPhiV_lsmm_mix"))); } if (TString(subGroup).Contains("dca")) { @@ -312,33 +318,30 @@ void o2::aod::pwgem::photon::histogram::DefineHistograms(THashList* list, const } } // end of mc } else if (TString(histClass).Contains("MuMu")) { - const int ndim = 4; // m, pt, dca, phiv - const int nbins[ndim] = {90, 20, ndca - 1, 1}; - const double xmin[ndim] = {0.2, 0.0, 0.0, 0.0}; - const double xmax[ndim] = {1.1, 2.0, 5.0, 3.2}; + const int ndim = 3; // m, pt, dca + const int nbins[ndim] = {90, 20, ndca - 1}; + const double xmin[ndim] = {0.2, 0.0, 0.0}; + const double xmax[ndim] = {1.1, 2.0, 5.0}; - hs_dilepton_uls_same = new THnSparseF("hs_dilepton_uls_same", "hs_dilepton_uls;m_{#mu#mu} (GeV/c^{2});p_{T,#mu#mu} (GeV/c);DCA_{#mu#mu}^{3D} (#sigma);#varphi_{V} (rad.);", ndim, nbins, xmin, xmax); + hs_dilepton_uls_same = new THnSparseF("hs_dilepton_uls_same", "hs_dilepton_uls;m_{#mu#mu} (GeV/c^{2});p_{T,#mu#mu} (GeV/c);DCA_{#mu#mu}^{3D} (#sigma);", ndim, nbins, xmin, xmax); hs_dilepton_uls_same->Sumw2(); hs_dilepton_uls_same->SetBinEdges(2, dca); list->Add(hs_dilepton_uls_same); - hs_dilepton_lspp_same = new THnSparseF("hs_dilepton_lspp_same", "hs_dilepton_lspp;m_{#mu#mu} (GeV/c^{2});p_{T,#mu#mu} (GeV/c);DCA_{#mu#mu}^{3D} (#sigma);#varphi_{V} (rad.);", ndim, nbins, xmin, xmax); + hs_dilepton_lspp_same = new THnSparseF("hs_dilepton_lspp_same", "hs_dilepton_lspp;m_{#mu#mu} (GeV/c^{2});p_{T,#mu#mu} (GeV/c);DCA_{#mu#mu}^{3D} (#sigma);", ndim, nbins, xmin, xmax); hs_dilepton_lspp_same->Sumw2(); hs_dilepton_lspp_same->SetBinEdges(2, dca); list->Add(hs_dilepton_lspp_same); - hs_dilepton_lsmm_same = new THnSparseF("hs_dilepton_lsmm_same", "hs_dilepton_lsmm;m_{#mu#mu} (GeV/c^{2});p_{T,#mu#mu} (GeV/c);DCA_{#mu#mu}^{3D} (#sigma);#varphi_{V} (rad.);", ndim, nbins, xmin, xmax); + hs_dilepton_lsmm_same = new THnSparseF("hs_dilepton_lsmm_same", "hs_dilepton_lsmm;m_{#mu#mu} (GeV/c^{2});p_{T,#mu#mu} (GeV/c);DCA_{#mu#mu}^{3D} (#sigma);", ndim, nbins, xmin, xmax); hs_dilepton_lsmm_same->Sumw2(); hs_dilepton_lsmm_same->SetBinEdges(2, dca); list->Add(hs_dilepton_lsmm_same); if (TString(subGroup).Contains("mix")) { - THnSparseF* hs_dilepton_uls_mix = reinterpret_cast(hs_dilepton_uls_same->Clone("hs_dilepton_uls_mix")); - THnSparseF* hs_dilepton_lspp_mix = reinterpret_cast(hs_dilepton_lspp_same->Clone("hs_dilepton_lspp_mix")); - THnSparseF* hs_dilepton_lsmm_mix = reinterpret_cast(hs_dilepton_lsmm_same->Clone("hs_dilepton_lsmm_mix")); - list->Add(hs_dilepton_uls_mix); - list->Add(hs_dilepton_lspp_mix); - list->Add(hs_dilepton_lsmm_mix); + list->Add(reinterpret_cast(hs_dilepton_uls_same->Clone("hs_dilepton_uls_mix"))); + list->Add(reinterpret_cast(hs_dilepton_lspp_same->Clone("hs_dilepton_lspp_mix"))); + list->Add(reinterpret_cast(hs_dilepton_lsmm_same->Clone("hs_dilepton_lsmm_mix"))); } } else { LOGF(info, "EE or MuMu are supported."); diff --git a/PWGEM/PhotonMeson/Tasks/dalitzEEQC.cxx b/PWGEM/PhotonMeson/Tasks/dalitzEEQC.cxx index f7b3f26b1dc..c7f03892676 100644 --- a/PWGEM/PhotonMeson/Tasks/dalitzEEQC.cxx +++ b/PWGEM/PhotonMeson/Tasks/dalitzEEQC.cxx @@ -146,7 +146,7 @@ struct DalitzEEQC { THashList* list_ev_after = static_cast(fMainList->FindObject("Event")->FindObject(event_types[1].data())); THashList* list_dalitzee = static_cast(fMainList->FindObject("DalitzEE")); THashList* list_track = static_cast(fMainList->FindObject("Track")); - double values[4] = {0, 0, 0, 0}; + double values[3] = {0, 0, 0}; double values_single[4] = {0, 0, 0, 0}; float dca_pos_3d = 999.f, dca_ele_3d = 999.f, dca_ee_3d = 999.f; @@ -195,8 +195,8 @@ struct DalitzEEQC { values[0] = uls_pair.mass(); values[1] = uls_pair.pt(); values[2] = dca_ee_3d; - values[3] = uls_pair.phiv(); reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_uls_same"))->Fill(values); + reinterpret_cast(list_dalitzee_cut->FindObject("hMvsPhiV_uls_same"))->Fill(uls_pair.phiv(), uls_pair.mass()); nuls++; for (auto& track : {pos, ele}) { if (std::find(used_trackIds.begin(), used_trackIds.end(), track.globalIndex()) == used_trackIds.end()) { @@ -226,8 +226,8 @@ struct DalitzEEQC { values[0] = lspp_pair.mass(); values[1] = lspp_pair.pt(); values[2] = dca_ee_3d; - values[3] = lspp_pair.phiv(); reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_lspp_same"))->Fill(values); + reinterpret_cast(list_dalitzee_cut->FindObject("hMvsPhiV_lspp_same"))->Fill(lspp_pair.phiv(), lspp_pair.mass()); nlspp++; } } // end of lspp pair loop @@ -252,8 +252,8 @@ struct DalitzEEQC { values[0] = lsmm_pair.mass(); values[1] = lsmm_pair.pt(); values[2] = dca_ee_3d; - values[3] = lsmm_pair.phiv(); reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_lsmm_same"))->Fill(values); + reinterpret_cast(list_dalitzee_cut->FindObject("hMvsPhiV_lsmm_same"))->Fill(lsmm_pair.phiv(), lsmm_pair.mass()); nlsmm++; } } // end of lsmm pair loop @@ -286,7 +286,7 @@ struct DalitzEEQC { void MixedEventPairing(TEvents const& collisions, TTracks const& tracks, TMixedBinning const& colBinning) { THashList* list_dalitzee = static_cast(fMainList->FindObject("DalitzEE")); - double values[4] = {0, 0, 0, 0}; + double values[3] = {0, 0, 0}; ROOT::Math::PtEtaPhiMVector v1, v2, v12; float phiv = 0; double values_single[4] = {0, 0, 0, 0}; @@ -317,6 +317,13 @@ struct DalitzEEQC { v2 = ROOT::Math::PtEtaPhiMVector(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); v12 = v1 + v2; phiv = getPhivPair(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), t1.sign(), t2.sign(), collision1.bz()); + // if (!std::isfinite(phiv)) { + // LOGF(info, "t1.px() = %f, t1.py() = %f, t1.pz() = %f, t2.px() = %f, t2.py() = %f, t2.pz() = %f, phiv = %f", t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), phiv); + // } + + if (t1.trackId() == t2.trackId()) { // this is protection against pairing identical 2 tracks. This happens, when TTCA is used. TTCA can assign a track to several possible collisions. + continue; + } dca_pos_3d = t1.dca3DinSigma(); dca_ele_3d = t2.dca3DinSigma(); @@ -330,15 +337,17 @@ struct DalitzEEQC { values[0] = v12.M(); values[1] = v12.Pt(); values[2] = dca_ee_3d; - values[3] = phiv; if (cut.IsSelectedTrack(t1) && cut.IsSelectedTrack(t2) && cut.IsSelectedPair(v12.M(), dca_ee_3d, phiv)) { if (t1.sign() * t2.sign() < 0) { reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_uls_mix"))->Fill(values); + reinterpret_cast(list_dalitzee_cut->FindObject("hMvsPhiV_uls_mix"))->Fill(phiv, v12.M()); } else if (t1.sign() > 0 && t2.sign() > 0) { reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_lspp_mix"))->Fill(values); + reinterpret_cast(list_dalitzee_cut->FindObject("hMvsPhiV_lspp_mix"))->Fill(phiv, v12.M()); } else if (t1.sign() < 0 && t2.sign() < 0) { reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_lsmm_mix"))->Fill(values); + reinterpret_cast(list_dalitzee_cut->FindObject("hMvsPhiV_lsmm_mix"))->Fill(phiv, v12.M()); } else { LOGF(info, "This should not happen."); } diff --git a/PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx b/PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx index c8f7bcf54cd..529fe4031fe 100644 --- a/PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx @@ -170,7 +170,7 @@ struct DalitzEEQCMC { THashList* list_ev_after = static_cast(fMainList->FindObject("Event")->FindObject(event_types[1].data())); THashList* list_dalitzee = static_cast(fMainList->FindObject("DalitzEE")); THashList* list_track = static_cast(fMainList->FindObject("Track")); - double values[4] = {0, 0, 0, 0}; + double values[3] = {0, 0, 0}; double values_single[4] = {0, 0, 0, 0}; float dca_pos_3d = 999.f, dca_ele_3d = 999.f, dca_ee_3d = 999.f; @@ -227,7 +227,6 @@ struct DalitzEEQCMC { values[0] = uls_pair.mass(); values[1] = uls_pair.pt(); values[2] = dca_ee_3d; - values[3] = uls_pair.phiv(); if (mother_id > -1) { auto mcmother = mcparticles.iteratorAt(mother_id); @@ -354,7 +353,6 @@ struct DalitzEEQCMC { values[0] = lspp_pair.mass(); values[1] = lspp_pair.pt(); values[2] = dca_ee_3d; - values[3] = lspp_pair.phiv(); if ((posmc.isPhysicalPrimary() || posmc.producedByGenerator()) && (elemc.isPhysicalPrimary() || elemc.producedByGenerator())) { switch (hfee_type) { @@ -418,7 +416,6 @@ struct DalitzEEQCMC { values[0] = lsmm_pair.mass(); values[1] = lsmm_pair.pt(); values[2] = dca_ee_3d; - values[3] = lsmm_pair.phiv(); if ((posmc.isPhysicalPrimary() || posmc.producedByGenerator()) && (elemc.isPhysicalPrimary() || elemc.producedByGenerator())) { switch (hfee_type) { diff --git a/PWGEM/PhotonMeson/Tasks/dalitzMuMuQC.cxx b/PWGEM/PhotonMeson/Tasks/dalitzMuMuQC.cxx index b41b89d319f..00d4cffbff4 100644 --- a/PWGEM/PhotonMeson/Tasks/dalitzMuMuQC.cxx +++ b/PWGEM/PhotonMeson/Tasks/dalitzMuMuQC.cxx @@ -156,7 +156,7 @@ struct DalitzMuMuQC { THashList* list_ev_after = static_cast(fMainList->FindObject("Event")->FindObject(event_types[1].data())); THashList* list_dalitzmumu = static_cast(fMainList->FindObject("DalitzMuMu")); THashList* list_track = static_cast(fMainList->FindObject("Track")); - double values[4] = {0, 0, 0, 0}; + double values[3] = {0, 0, 0}; float dca_pos_3d = 999.f, dca_ele_3d = 999.f, dca_ee_3d = 999.f; for (auto& collision : grouped_collisions) { @@ -195,7 +195,6 @@ struct DalitzMuMuQC { values[0] = uls_pair.mass(); values[1] = uls_pair.pt(); values[2] = dca_ee_3d; - values[3] = uls_pair.phiv(); reinterpret_cast(list_dalitzmumu_cut->FindObject("hs_dilepton_uls_same"))->Fill(values); nuls++; for (auto& track : {pos, ele}) { @@ -219,7 +218,6 @@ struct DalitzMuMuQC { values[0] = lspp_pair.mass(); values[1] = lspp_pair.pt(); values[2] = dca_ee_3d; - values[3] = lspp_pair.phiv(); reinterpret_cast(list_dalitzmumu_cut->FindObject("hs_dilepton_lspp_same"))->Fill(values); nlspp++; } @@ -237,7 +235,6 @@ struct DalitzMuMuQC { values[0] = lsmm_pair.mass(); values[1] = lsmm_pair.pt(); values[2] = dca_ee_3d; - values[3] = lsmm_pair.phiv(); reinterpret_cast(list_dalitzmumu_cut->FindObject("hs_dilepton_lsmm_same"))->Fill(values); nlsmm++; } @@ -271,7 +268,7 @@ struct DalitzMuMuQC { void MixedEventPairing(TEvents const& collisions, TTracks const& tracks, TMixedBinning const& colBinning) { THashList* list_dalitzmumu = static_cast(fMainList->FindObject("DalitzMuMu")); - double values[4] = {0, 0, 0, 0}; + double values[3] = {0, 0, 0}; ROOT::Math::PtEtaPhiMVector v1, v2, v12; float phiv = 0; float dca_pos_3d = 999.f, dca_ele_3d = 999.f, dca_ee_3d = 999.f; @@ -297,6 +294,10 @@ struct DalitzMuMuQC { for (auto& cut : fDalitzMuMuCuts) { THashList* list_dalitzmumu_cut = static_cast(list_dalitzmumu->FindObject(cut.GetName())); for (auto& [t1, t2] : combinations(soa::CombinationsFullIndexPolicy(tracks_coll1, tracks_coll2))) { + if (t1.trackId() == t2.trackId()) { // this is protection against pairing identical 2 tracks. This happens, when TTCA is used. TTCA can assign a track to several possible collisions. + continue; + } + v1 = ROOT::Math::PtEtaPhiMVector(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassMuon); v2 = ROOT::Math::PtEtaPhiMVector(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassMuon); v12 = v1 + v2; @@ -308,7 +309,6 @@ struct DalitzMuMuQC { values[0] = v12.M(); values[1] = v12.Pt(); values[2] = dca_ee_3d; - values[3] = phiv; if (cut.IsSelectedTrack(t1) && cut.IsSelectedTrack(t2) && cut.IsSelectedPair(v12.M(), dca_ee_3d, phiv)) { if (t1.sign() * t2.sign() < 0) { diff --git a/PWGEM/PhotonMeson/Tasks/dalitzMuMuQCMC.cxx b/PWGEM/PhotonMeson/Tasks/dalitzMuMuQCMC.cxx index ae44a4fc2da..b45dbb6e6f2 100644 --- a/PWGEM/PhotonMeson/Tasks/dalitzMuMuQCMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/dalitzMuMuQCMC.cxx @@ -161,7 +161,7 @@ struct DalitzMuMuQCMC { THashList* list_ev_after = static_cast(fMainList->FindObject("Event")->FindObject(event_types[1].data())); THashList* list_dalitzmumu = static_cast(fMainList->FindObject("DalitzMuMu")); THashList* list_track = static_cast(fMainList->FindObject("Track")); - double values[4] = {0, 0, 0, 0}; + double values[3] = {0, 0, 0}; float dca_pos_3d = 999.f, dca_ele_3d = 999.f, dca_ee_3d = 999.f; for (auto& collision : grouped_collisions) { @@ -211,7 +211,6 @@ struct DalitzMuMuQCMC { values[0] = uls_pair.mass(); values[1] = uls_pair.pt(); values[2] = dca_ee_3d; - values[3] = uls_pair.phiv(); reinterpret_cast(list_dalitzmumu_cut->FindObject("hs_dilepton_uls_same"))->Fill(values); nuls++; diff --git a/PWGEM/PhotonMeson/Utils/PCMUtilities.h b/PWGEM/PhotonMeson/Utils/PCMUtilities.h index 389eb79ed0f..ecf23c85a6b 100644 --- a/PWGEM/PhotonMeson/Utils/PCMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/PCMUtilities.h @@ -165,8 +165,8 @@ float getPhivPair(float pxpos, float pypos, float pzpos, float pxneg, float pyne float uy = (pypos + pyneg) / RecoDecay::sqrtSumOfSquares(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); float uz = (pzpos + pzneg) / RecoDecay::sqrtSumOfSquares(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); - float ax = uy / TMath::Sqrt(ux * ux + uy * uy); - float ay = -ux / TMath::Sqrt(ux * ux + uy * uy); + float ax = uy / std::sqrt(ux * ux + uy * uy); + float ay = -ux / std::sqrt(ux * ux + uy * uy); // The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz) float wx = uy * vz - uz * vy; @@ -175,13 +175,16 @@ float getPhivPair(float pxpos, float pypos, float pzpos, float pxneg, float pyne // The angle between them should be small if the pair is conversion. This function then returns values close to pi! auto clipToPM1 = [](float x) { return x < -1.f ? -1.f : (x > 1.f ? 1.f : x); }; - // if(!std::isfinite(std::acos(wx * ax + wy * ay))){ + // if (!std::isfinite(std::acos(wx * ax + wy * ay))) { + // LOGF(info, "pxpos = %f, pypos = %f, pzpos = %f", pxpos, pypos, pzpos); + // LOGF(info, "pxneg = %f, pyneg = %f, pzneg = %f", pxneg, pyneg, pzneg); // LOGF(info, "pos_x_ele[0] = %f, pos_x_ele[1] = %f, pos_x_ele[2] = %f", pos_x_ele[0], pos_x_ele[1], pos_x_ele[2]); // LOGF(info, "ux = %f, uy = %f, uz = %f", ux, uy, uz); // LOGF(info, "ax = %f, ay = %f", ax, ay); // LOGF(info, "wx = %f, wy = %f", wx, wy); // LOGF(info, "wx * ax + wy * ay = %f", wx * ax + wy * ay); - // LOGF(info, "std::acos(wx * ax + wy * ay) = %f", std::acos(clipToPM1(wx * ax + wy * ay))); + // LOGF(info, "std::acos(wx * ax + wy * ay) = %f", std::acos(wx * ax + wy * ay)); + // LOGF(info, "std::acos(clipToPM1(wx * ax + wy * ay)) = %f", std::acos(clipToPM1(wx * ax + wy * ay))); // } return std::acos(clipToPM1(wx * ax + wy * ay)); // phiv in [0,pi] //cosPhiV = wx * ax + wy * ay;