Skip to content

Commit a3e82c9

Browse files
dsekihatdsekihat
authored
[PWGEM/Dilepton] update track propagation (#8447)
1 parent a395f40 commit a3e82c9

File tree

2 files changed

+132
-20
lines changed

2 files changed

+132
-20
lines changed

‎PWGEM/Dilepton/Core/Dilepton.h‎

Lines changed: 68 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,9 @@
3939
#include "CommonConstants/LHCConstants.h"
4040
#include "DataFormatsParameters/GRPLHCIFData.h"
4141
#include "DataFormatsParameters/GRPECSObject.h"
42+
#include "MathUtils/Utils.h"
4243

44+
#include "DetectorsBase/Propagator.h"
4345
#include "DetectorsBase/GeometryManager.h"
4446
#include "DataFormatsParameters/GRPObject.h"
4547
#include "DataFormatsParameters/GRPMagField.h"
@@ -172,9 +174,13 @@ struct Dilepton {
172174
Configurable<bool> cfg_apply_detadphi{"cfg_apply_detadphi", false, "flag to apply deta-dphi elliptic cut"};
173175
Configurable<float> cfg_min_deta{"cfg_min_deta", 0.02, "min deta between 2 electrons (elliptic cut)"};
174176
Configurable<float> cfg_min_dphi{"cfg_min_dphi", 0.2, "min dphi between 2 electrons (elliptic cut)"};
177+
Configurable<bool> cfg_apply_detadphi_geom{"cfg_apply_detadphi_geom", false, "flag to apply generator deta-dphi elliptic cut"};
178+
Configurable<float> cfg_min_deta_geom{"cfg_min_deta_geom", 0.02, "geometrical min deta between 2 electrons (elliptic cut)"};
179+
Configurable<float> cfg_min_dphi_geom{"cfg_min_dphi_geom", 0.2, "geometrical min dphi between 2 electrons (elliptic cut)"};
175180
Configurable<float> cfg_min_opang{"cfg_min_opang", 0.0, "min opening angle"};
176181
Configurable<float> cfg_max_opang{"cfg_max_opang", 6.4, "max opening angle"};
177182
Configurable<bool> cfg_require_diff_sides{"cfg_require_diff_sides", false, "flag to require 2 tracks are from different sides."};
183+
Configurable<float> cfg_x_to_go{"cfg_x_to_go", -1, "x (cm) to be propagated in local coordinate"};
178184

179185
Configurable<float> cfg_min_pt_track{"cfg_min_pt_track", 0.2, "min pT for single track"};
180186
Configurable<float> cfg_min_eta_track{"cfg_min_eta_track", -0.8, "min eta for single track"};
@@ -264,7 +270,7 @@ struct Dilepton {
264270
float d_bz;
265271
// o2::vertexing::DCAFitterN<2> fitter;
266272
// o2::vertexing::FwdDCAFitterN<2> fwdfitter;
267-
// o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
273+
o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
268274

269275
HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false};
270276
static constexpr std::string_view event_cut_types[2] = {"before/", "after/"};
@@ -419,6 +425,7 @@ struct Dilepton {
419425
if (fabs(d_bz) > 1e-5) {
420426
grpmag.setL3Current(30000.f / (d_bz / 5.0f));
421427
}
428+
o2::base::Propagator::initFieldFromGRP(&grpmag);
422429
mRunNumber = collision.runNumber();
423430
// fitter.setBz(d_bz);
424431
// fwdfitter.setBz(d_bz);
@@ -431,6 +438,7 @@ struct Dilepton {
431438
if (!skipGRPOquery)
432439
grpo = ccdb->getForTimeStamp<o2::parameters::GRPObject>(grpPath, run3grp_timestamp);
433440
if (grpo) {
441+
o2::base::Propagator::initFieldFromGRP(grpo);
434442
// Fetch magnetic field from ccdb for current collision
435443
d_bz = grpo->getNominalL3Field();
436444
LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
@@ -439,6 +447,7 @@ struct Dilepton {
439447
if (!grpmag) {
440448
LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp;
441449
}
450+
o2::base::Propagator::initFieldFromGRP(grpmag);
442451
// Fetch magnetic field from ccdb for current collision
443452
d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f);
444453
LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
@@ -515,8 +524,8 @@ struct Dilepton {
515524

516525
if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kQC)) {
517526
fRegistry.add("Pair/same/uls/hs", "dilepton", kTHnSparseD, {axis_mass, axis_pt, axis_dca}, true);
518-
fRegistry.add("Pair/same/uls/hsDeltaP", "difference of p between 2 tracks;|p_{T,1} - p_{T,2}|/|p_{T,1} + p_{T,2}|;#Delta#eta;#Delta#varphi (rad.);", kTHnSparseD, {{20, 0, 1}, {100, -0.5, +0.5}, {100, -0.5, 0.5}}, true);
519-
527+
fRegistry.add("Pair/same/uls/hsDeltaP", "difference of p between 2 tracks;|p_{T,1} - p_{T,2}|/|p_{T,1} + p_{T,2}|;#Delta#eta;#Delta#varphi (rad.);", kTHnSparseD, {{20, 0, 1}, {100, -0.5, +0.5}, {180, -M_PI, M_PI}}, true);
528+
fRegistry.add("Pair/same/uls/hGeomDeltaEtaDeltaPhi", Form("difference in #eta-#varphi plane between 2 tracks at X = %2.1f cm;#Delta#varphi (rad.);#Delta#eta;", dielectroncuts.cfg_x_to_go.value), kTH2D, {{180, -M_PI, M_PI}, {100, -0.5, +0.5}}, true);
520529
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
521530
fRegistry.add("Pair/same/uls/hMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", kTH2D, {{90, 0, M_PI}, {100, 0.0f, 1.0f}}, true); // phiv is only for dielectron
522531
fRegistry.add("Pair/same/uls/hMvsOpAng", "m_{ee} vs. angle between 2 tracks;#omega (rad.);m_{ee} (GeV/c^{2})", kTH2D, {{100, 0, 2.0}, {20, 0.0f, 3.2}}, true);
@@ -807,10 +816,36 @@ struct Dilepton {
807816
}
808817
}
809818

819+
float deta_geom = 999.f;
820+
float dphi_geom = 999.f;
810821
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
811822
if (!cut.template IsSelectedPair(t1, t2, d_bz)) {
812823
return false;
813824
}
825+
if (dielectroncuts.cfg_x_to_go) {
826+
auto track_par_cov1 = getTrackParCov(t1);
827+
track_par_cov1.setPID(o2::track::PID::Electron);
828+
o2::base::Propagator::Instance()->propagateToX(track_par_cov1, dielectroncuts.cfg_x_to_go, d_bz, o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, matCorr);
829+
auto xyz1 = track_par_cov1.getXYZGlo();
830+
float eta1 = RecoDecay::eta(std::array{xyz1.X(), xyz1.Y(), xyz1.Z()});
831+
float phi1 = RecoDecay::phi(std::array{xyz1.X(), xyz1.Y()});
832+
o2::math_utils::bringTo02Pi(phi1);
833+
834+
auto track_par_cov2 = getTrackParCov(t2);
835+
track_par_cov2.setPID(o2::track::PID::Electron);
836+
o2::base::Propagator::Instance()->propagateToX(track_par_cov2, dielectroncuts.cfg_x_to_go, d_bz, o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, matCorr);
837+
auto xyz2 = track_par_cov2.getXYZGlo();
838+
float eta2 = RecoDecay::eta(std::array{xyz2.X(), xyz2.Y(), xyz2.Z()});
839+
float phi2 = RecoDecay::phi(std::array{xyz2.X(), xyz2.Y()});
840+
o2::math_utils::bringTo02Pi(phi2);
841+
842+
deta_geom = t1.sign() * t1.pt() > t2.sign() * t2.pt() ? eta1 - eta2 : eta2 - eta1;
843+
dphi_geom = t1.sign() * t1.pt() > t2.sign() * t2.pt() ? phi1 - phi2 : phi2 - phi1;
844+
o2::math_utils::bringToPMPi(dphi_geom);
845+
if (dielectroncuts.cfg_apply_detadphi_geom && std::pow(deta_geom / dielectroncuts.cfg_min_deta_geom, 2) + std::pow(dphi_geom / dielectroncuts.cfg_min_dphi_geom, 2) < 1.f) {
846+
return false;
847+
}
848+
}
814849
} else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) {
815850
if (!cut.template IsSelectedPair(t1, t2)) {
816851
return false;
@@ -852,26 +887,30 @@ struct Dilepton {
852887
float deta = t1.sign() * v1.Pt() > t2.sign() * v2.Pt() ? v1.Eta() - v2.Eta() : v2.Eta() - v1.Eta();
853888
float dphi = t1.sign() * v1.Pt() > t2.sign() * v2.Pt() ? v1.Phi() - v2.Phi() : v2.Phi() - v1.Phi();
854889
o2::math_utils::bringToPMPi(dphi);
890+
855891
float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), t1.sign(), t2.sign(), d_bz);
856892
float opAng = o2::aod::pwgem::dilepton::utils::pairutil::getOpeningAngle(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz());
857893

858894
if (t1.sign() * t2.sign() < 0) { // ULS
859895
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, weight);
860896
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hsDeltaP"), dpt, deta, dphi, weight);
897+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hGeomDeltaEtaDeltaPhi"), dphi_geom, deta_geom, weight);
861898
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
862899
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hMvsPhiV"), phiv, v12.M(), weight);
863900
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hMvsOpAng"), opAng, v12.M(), weight);
864901
}
865902
} else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
866903
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, weight);
867904
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hsDeltaP"), dpt, deta, dphi, weight);
905+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hGeomDeltaEtaDeltaPhi"), dphi_geom, deta_geom, weight);
868906
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
869907
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hMvsPhiV"), phiv, v12.M(), weight);
870908
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hMvsOpAng"), opAng, v12.M(), weight);
871909
}
872910
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
873911
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, weight);
874912
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hsDeltaP"), dpt, deta, dphi, weight);
913+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hGeomDeltaEtaDeltaPhi"), dphi_geom, deta_geom, weight);
875914
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
876915
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hMvsPhiV"), phiv, v12.M(), weight);
877916
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hMvsOpAng"), opAng, v12.M(), weight);
@@ -895,7 +934,6 @@ struct Dilepton {
895934
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
896935
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, aco, asym, abs(dphi_e_ee), abs(cos_thetaCS), weight);
897936
}
898-
899937
} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kFlowV2) || cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kFlowV3)) {
900938
std::array<float, 2> q2ft0m = {collision.q2xft0m(), collision.q2yft0m()};
901939
std::array<float, 2> q2ft0a = {collision.q2xft0a(), collision.q2yft0a()};
@@ -923,13 +961,10 @@ struct Dilepton {
923961
float sp = RecoDecay::dotProd(std::array<float, 2>{static_cast<float>(std::cos(nmod * v12.Phi())), static_cast<float>(std::sin(nmod * v12.Phi()))}, qvectors[nmod][cfgQvecEstimator]) / getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange());
924962
if (t1.sign() * t2.sign() < 0) { // ULS
925963
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, sp, weight);
926-
// fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hPrfUQ"), v12.M(), v12.Pt(), pair_dca, sp, weight);
927964
} else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
928965
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, sp, weight);
929-
// fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hPrfUQ"), v12.M(), v12.Pt(), pair_dca, sp, weight);
930966
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
931967
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, sp, weight);
932-
// fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hPrfUQ"), v12.M(), v12.Pt(), pair_dca, sp, weight);
933968
}
934969
} else if constexpr (ev_id == 1) {
935970
if (t1.sign() * t2.sign() < 0) { // ULS
@@ -1316,6 +1351,32 @@ struct Dilepton {
13161351
if (!cut.template IsSelectedPair(t1, t2, d_bz)) {
13171352
return false;
13181353
}
1354+
1355+
if (dielectroncuts.cfg_x_to_go) {
1356+
auto track_par_cov1 = getTrackParCov(t1);
1357+
track_par_cov1.setPID(o2::track::PID::Electron);
1358+
o2::base::Propagator::Instance()->propagateToX(track_par_cov1, dielectroncuts.cfg_x_to_go, d_bz, o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, matCorr);
1359+
auto xyz1 = track_par_cov1.getXYZGlo();
1360+
float eta1 = RecoDecay::eta(std::array{xyz1.X(), xyz1.Y(), xyz1.Z()});
1361+
float phi1 = RecoDecay::phi(std::array{xyz1.X(), xyz1.Y()});
1362+
o2::math_utils::bringTo02Pi(phi1);
1363+
1364+
auto track_par_cov2 = getTrackParCov(t2);
1365+
track_par_cov2.setPID(o2::track::PID::Electron);
1366+
o2::base::Propagator::Instance()->propagateToX(track_par_cov2, dielectroncuts.cfg_x_to_go, d_bz, o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, matCorr);
1367+
auto xyz2 = track_par_cov2.getXYZGlo();
1368+
float eta2 = RecoDecay::eta(std::array{xyz2.X(), xyz2.Y(), xyz2.Z()});
1369+
float phi2 = RecoDecay::phi(std::array{xyz2.X(), xyz2.Y()});
1370+
o2::math_utils::bringTo02Pi(phi2);
1371+
1372+
float deta_geom = t1.sign() * t1.pt() > t2.sign() * t2.pt() ? eta1 - eta2 : eta2 - eta1;
1373+
float dphi_geom = t1.sign() * t1.pt() > t2.sign() * t2.pt() ? phi1 - phi2 : phi2 - phi1;
1374+
o2::math_utils::bringToPMPi(dphi_geom);
1375+
1376+
if (dielectroncuts.cfg_apply_detadphi_geom && std::pow(deta_geom / dielectroncuts.cfg_min_deta_geom, 2) + std::pow(dphi_geom / dielectroncuts.cfg_min_dphi_geom, 2) < 1.f) {
1377+
return false;
1378+
}
1379+
}
13191380
} else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) {
13201381
if (!cut.template IsSelectedPair(t1, t2)) {
13211382
return false;

0 commit comments

Comments
 (0)