diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 69c29570f2e..f9752b597a5 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -33,6 +33,7 @@ #include #include +#include #include #include #include @@ -75,7 +76,9 @@ int VarManager::fgCalibrationType = 0; // 0 - no calibration, 1 - bool VarManager::fgUseInterpolatedCalibration = true; // use interpolated calibration histograms (default: true) int VarManager::fgEfficiencyType = 0; // type of efficiency to be applied, default is no efficiency TObject* VarManager::fgEfficiencyHist = nullptr; // histogram for efficiency - +TObject* VarManager::fgPosiPhiMap = nullptr; +TObject* VarManager::fgNegaPhiMap = nullptr; +bool VarManager::fgUsePhiCorrection = false; //__________________________________________________________________ VarManager::VarManager() : TObject() { @@ -453,6 +456,76 @@ void VarManager::FillEfficiency(float* values) } } +void VarManager::SetPhiMap(TObject* hposi, TObject* hnega, bool option) +{ + fgPosiPhiMap = hposi; + fgNegaPhiMap = hnega; + fgUsePhiCorrection = option; +} + +double VarManager::SampleRotationPhi(double pt, double eta, int charge) +{ + // each type only alarm once + static bool warnedEmptyPhi = false; + + if (!fgUsePhiCorrection) { + return gRandom->Uniform(0., o2::constants::math::TwoPI); + } else { + + TH3D* hMap = nullptr; + if (charge > 0) { + hMap = dynamic_cast(fgPosiPhiMap); + } else { + hMap = dynamic_cast(fgNegaPhiMap); + } + + if (!hMap) { + LOGF(fatal, "Phi map is not a TH3D"); + } + // TH3 axes: X=pT, Y=phi, Z=eta + int ptBin = hMap->GetXaxis()->FindBin(pt); + int etaBin = hMap->GetZaxis()->FindBin(eta); + + ptBin = std::clamp(ptBin, 1, hMap->GetNbinsX()); + etaBin = std::clamp(etaBin, 1, hMap->GetNbinsZ()); + + TH1D* hPhi = hMap->ProjectionY( + Form("hTRPhi_tmp_charge%d_ptbin%d_etabin%d", + charge, ptBin, etaBin), + ptBin, + ptBin, + etaBin, + etaBin); + + if (!hPhi || hPhi->Integral(1, hPhi->GetNbinsX()) <= 0.) { + if (!warnedEmptyPhi) { + LOGF(warn, + "Empty phi distribution for " + "pt=%f, eta=%f, charge=%d, ptBin=%d, etaBin=%d. " + "Falling back to uniform phi sampling.", + pt, + eta, + charge, + ptBin, + etaBin); + + warnedEmptyPhi = true; + } + + delete hPhi; + + return gRandom->Uniform( + 0., o2::constants::math::TwoPI); + } + + const double phi = + RecoDecay::constrainAngle(hPhi->GetRandom()); + + delete hPhi; + return phi; + } +} + //__________________________________________________________________ std::tuple VarManager::BimodalityCoefficientUnbinned(const std::vector& data) { diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 60d9373e496..12857ca9766 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1530,6 +1530,8 @@ class VarManager : public TObject static void SetEfficiencyObject(int type, TObject* obj); static void FillEfficiency(float* values = nullptr); + static void SetPhiMap(TObject* h1, TObject* h2, bool option); + static double SampleRotationPhi(double pT, double eta, int charge); static TObject* GetCalibrationObject(CalibObjects calib) { auto obj = fgCalibs.find(calib); @@ -1606,14 +1608,18 @@ class VarManager : public TObject static o2::vertexing::FwdDCAFitterN<3> fgFitterThreeProngFwd; static o2::globaltracking::MatchGlobalFwd mMatching; - static std::map fgCalibs; // map of calibration histograms + static std::map fgCalibs; // map of calibration histograms static std::array fgRunTPCPostCalibration; // 0-electron, 1-pion, 2-kaon, 3-proton - static int fgCalibrationType; // 0 - no calibration, 1 - calibration vs (TPCncls,pIN,eta) typically for pp, 2 - calibration vs (eta,nPV,nLong,tLong) typically for PbPb - static bool fgUseInterpolatedCalibration; // use interpolated calibration histograms (default: true) + static int fgCalibrationType; // 0 - no calibration, 1 - calibration vs (TPCncls,pIN,eta) typically for pp, 2 - calibration vs (eta,nPV,nLong,tLong) typically for PbPb + static bool fgUseInterpolatedCalibration; // use interpolated calibration histograms (default: true) static int fgEfficiencyType; // type of efficiency correction to apply static TObject* fgEfficiencyHist; // histogram for efficiency correction + static TObject* fgPosiPhiMap; // phi map to correct track rotation + static TObject* fgNegaPhiMap; + static bool fgUsePhiCorrection; + VarManager& operator=(const VarManager& c); VarManager(const VarManager& c); }; @@ -3956,8 +3962,7 @@ void VarManager::FillPairRotation(T1 const& t1, T2 const& t2, float* values) m2 = o2::constants::physics::MassMuon; } - double dphi = gRandom->Uniform(0., o2::constants::math::TwoPI); - double rotationphi2 = RecoDecay::constrainAngle(t2.phi() + dphi); + double rotationphi2 = SampleRotationPhi(t2.pt(), t2.eta(), t2.sign()); values[kCharge] = t1.sign() + t2.sign(); values[kCharge1] = t1.sign(); diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index aba3686f2e8..7bc5d77d76f 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -1362,6 +1362,7 @@ struct AnalysisSameEventPairing { Configurable GrpLhcIfPath{"grplhcif", "GLO/Config/GRPLHCIF", "Path on the CCDB for the GRPLHCIF object"}; Configurable efficiencyPath{"effHistPath", "Users/z/zhxiong/efficiency", "Path on the CCDB for the efficiency histograms"}; Configurable flowPath{"flowPath", "Users/y/yiping/FlowResolution", "Path to the flow resolution object"}; + Configurable phiPath{"phiPath", "Users/h/hxiaoyu/TrackPhi/LHC23", "Path to load phi distribution for track rotation"}; } fConfigCCDB; struct : ConfigurableGroup { @@ -1381,6 +1382,7 @@ struct AnalysisSameEventPairing { Configurable useEfficiencyWeighting{"cfgUseEfficiencyWeighting", false, "Apply efficiency weighting to the pairs from CCDB"}; Configurable efficiencyType{"cfgEfficiencyType", 0, "Type of efficiency to apply from CCDB: 0 no efficiency, 1 pt-cent-costhetastar"}; Configurable useFlowReso{"cfgUseFlowReso", false, "Use remote flow information from CCDB"}; + Configurable usePhiDistribution{"cfgUsePhiDistribution", false, "Use phi distribution to correct track rotation"}; } fConfigOptions; struct : ConfigurableGroup { Configurable applyBDT{"applyBDT", false, "Flag to apply ML selections"}; @@ -1855,6 +1857,18 @@ struct AnalysisSameEventPairing { LOGF(fatal, "Flow resolution histograms not available in CCDB at timestamp=%llu", timestamp); } } + + if (fConfigOptions.usePhiDistribution) { + TString PathPhi = fConfigCCDB.phiPath.value; + TString ccdbPathPhiPosi = Form("%s/hPtPhiPositive", PathPhi.Data()); + TString ccdbPathPhiNega = Form("%s/hPtPhiNegative", PathPhi.Data()); + auto phiPosi = fCCDB->getForTimeStamp(ccdbPathPhiPosi.Data(), timestamp); + auto phiNega = fCCDB->getForTimeStamp(ccdbPathPhiNega.Data(), timestamp); + if (phiPosi == nullptr || phiNega == nullptr) { + LOGF(fatal, "Phi distribution histograms not available in CCDB at timestamp=%llu", timestamp); + } + VarManager::SetPhiMap(phiPosi, phiNega, fConfigOptions.usePhiDistribution.value); + } } // Template function to run same event pairing (barrel-barrel, muon-muon, barrel-muon)