From 89cdbf99a006d5d230804abf0e96bd6859adcd22 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Mon, 15 Jun 2026 01:33:50 +0800 Subject: [PATCH 1/8] Coherent Jpsi A2 definition using event plane at midrapidity, same event only --- PWGDQ/Core/CutsLibrary.cxx | 48 +++ PWGDQ/Core/HistogramsLibrary.cxx | 39 +++ PWGDQ/Core/VarManager.cxx | 8 + PWGDQ/Core/VarManager.h | 294 +++++++++++++++++++ PWGDQ/Tasks/tableReader_withAssoc.cxx | 71 ++++- PWGDQ/Tasks/tableReader_withAssoc_direct.cxx | 98 ++++++- 6 files changed, 551 insertions(+), 7 deletions(-) diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index bcc8a1ba788..d5eaa426ef1 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -2973,6 +2973,17 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) return cut; } + if (!nameStr.compare("jpsi_debug_TPCTOF3_rejBadTOF")) { + cut->AddCut(GetAnalysisCut("jpsiStandardKine5")); + cut->AddCut(GetAnalysisCut("electronStandardQualityTPCOnly3")); + cut->AddCut(GetAnalysisCut("SPDfirst")); + cut->AddCut(GetAnalysisCut("dcaCut1_ionut")); + cut->AddCut(GetAnalysisCut("pidJpsi_TPCpion0")); + cut->AddCut(GetAnalysisCut("pidJpsi_beta")); + cut->AddCut(GetAnalysisCut("pidJpsi_noTOF_prot")); + return cut; + } + // ------------------------------------------------------------------------------------------------- // lmee pair cuts @@ -4113,6 +4124,15 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (!nameStr.compare("eventStandardSel8NoPileup")) { + cut->AddCut(VarManager::kVtxZ, -10.0, 10.0); + cut->AddCut(VarManager::kIsSel8, 0.5, 1.5); + cut->AddCut(VarManager::kIsNoSameBunch, 0.5, 1.5); + cut->AddCut(VarManager::kIsGoodZvtxFT0vsPV, 0.5, 1.5); + cut->AddCut(VarManager::kNoCollInTimeRangeStandard, 0.5, 1.5); + return cut; + } + if (!nameStr.compare("eventStandardSel8PbPbQualityCent90")) { cut->AddCut(VarManager::kVtxZ, -10.0, 10.0); cut->AddCut(VarManager::kIsSel8, 0.5, 1.5); @@ -4569,6 +4589,12 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (!nameStr.compare("jpsiStandardKine5")) { + cut->AddCut(VarManager::kP, 1.0, 1000.0); + cut->AddCut(VarManager::kEta, -0.9, 0.9); + return cut; + } + if (!nameStr.compare("jpsiKineSkimmed")) { cut->AddCut(VarManager::kPt, 0.7, 1000.0); cut->AddCut(VarManager::kEta, -0.9, 0.9); @@ -5052,6 +5078,13 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (!nameStr.compare("electronStandardQualityTPCOnly3")) { + cut->AddCut(VarManager::kTPCchi2, 0.0, 4.0); + cut->AddCut(VarManager::kTPCncls, 120, 161.); + return cut; + } + + if (!nameStr.compare("NoelectronStandardQualityTPCOnly")) { cut->AddCut(VarManager::kTPCchi2, 0.0, 4.0, true, VarManager::kTPCncls, 70, 161.); return cut; @@ -5270,6 +5303,21 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (!nameStr.compare("pidJpsi_TPCpion0")) { + cut->AddCut(VarManager::kTPCnSigmaPi, 4.0, 1000.0); + return cut; + } + + if (!nameStr.compare("pidJpsi_noTOF_prot")) { + cut->AddCut(VarManager::kTPCnSigmaPr, 3.5, 1000.0, false, VarManager::kHasTOF, -0.5, 0.5); + return cut; + } + + if (!nameStr.compare("pidJpsi_beta")) { + cut->AddCut(VarManager::kTOFbeta, 0.98, 1.02, false, VarManager::kHasTOF, 0.5, 1.5); + return cut; + } + // Magnus cuts ---------------------------------------------------------- if (!nameStr.compare("pidJpsi_magnus_ele1")) { diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index edfcf3d4985..189f488d5a9 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -1389,6 +1389,45 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "CosThetaStarMC", "", false, 100, -1.0, 1.0, VarManager::kMCCosThetaStar); } } + if (subGroupStr.Contains("flow-jpsi-ep")) { + int bins_A2[5] = {50, 20, 20, 9, 200}; + double minBins_A2[5] = {2.0, 0.0, -1., 0.0, -20.0}; + double maxBins_A2[5] = {4.0, 2.0, 1.0, 90.0, 20.0}; + int bins_DeltaPhi[5] = {50, 20, 20, 9, 10}; + double minBins_DeltaPhi[5] = {2.0, 0.0, -1., 0.0, 0}; + double maxBins_DeltaPhi[5] = {4.0, 2.0, 1.0, 90.0, 3.14}; + TString labels[5] = {"kMass", "kPt", "kRapidity", "kCentFT0C", "kA2EP"}; + if (subGroupStr.Contains("tpc")) { + int varA2_TPC_PP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kA2EP_PP_TPC}; + int varA2_TPC_RP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kA2EP_RP_TPC}; + int varDeltaPhi_TPC_PP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kDeltaPhiPP_TPC}; + int varDeltaPhi_TPC_RP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kDeltaPhiRP_TPC}; + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_A2PP_TPC", "", 5, varA2_TPC_PP, bins_A2, minBins_A2, maxBins_A2, 0, -1, kTRUE); + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_A2RP_TPC", "", 5, varA2_TPC_RP, bins_A2, minBins_A2, maxBins_A2, 0, -1, kTRUE); + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_DeltaPhiPP_TPC", "", 5, varDeltaPhi_TPC_PP, bins_DeltaPhi, minBins_DeltaPhi, maxBins_DeltaPhi, 0, -1, kTRUE); + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_DeltaPhiRP_TPC", "", 5, varDeltaPhi_TPC_RP, bins_DeltaPhi, minBins_DeltaPhi, maxBins_DeltaPhi, 0, -1, kTRUE); + } + if (subGroupStr.Contains("ft0c")) { + int varA2_FT0C_PP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kA2EP_PP_FT0C}; + int varA2_FT0C_RP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kA2EP_RP_FT0C}; + int varDeltaPhi_FT0C_PP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kDeltaPhiPP_FT0C}; + int varDeltaPhi_FT0C_RP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kDeltaPhiRP_FT0C}; + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_A2PP_FT0C", "", 5, varA2_FT0C_PP, bins_A2, minBins_A2, maxBins_A2, 0, -1, kTRUE); + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_A2RP_FT0C", "", 5, varA2_FT0C_RP, bins_A2, minBins_A2, maxBins_A2, 0, -1, kTRUE); + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_DeltaPhiPP_FT0C", "", 5, varDeltaPhi_FT0C_PP, bins_DeltaPhi, minBins_DeltaPhi, maxBins_DeltaPhi, 0, -1, kTRUE); + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_DeltaPhiRP_FT0C", "", 5, varDeltaPhi_FT0C_RP, bins_DeltaPhi, minBins_DeltaPhi, maxBins_DeltaPhi, 0, -1, kTRUE); + } + if (subGroupStr.Contains("ft0a")) { + int varA2_FT0A_PP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kA2EP_PP_FT0A}; + int varA2_FT0A_RP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kA2EP_RP_FT0A}; + int varDeltaPhi_FT0A_PP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kDeltaPhiPP_FT0A}; + int varDeltaPhi_FT0A_RP[5] = {VarManager::kMass, VarManager::kPt, VarManager::kRap, VarManager::kCentFT0C, VarManager::kDeltaPhiRP_FT0A}; + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_A2PP_FT0A", "", 5, varA2_FT0A_PP, bins_A2, minBins_A2, maxBins_A2, 0, -1, kTRUE); + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_A2RP_FT0A", "", 5, varA2_FT0A_RP, bins_A2, minBins_A2, maxBins_A2, 0, -1, kTRUE); + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_DeltaPhiPP_FT0A", "", 5, varDeltaPhi_FT0A_PP, bins_DeltaPhi, minBins_DeltaPhi, maxBins_DeltaPhi, 0, -1, kTRUE); + hm->AddHistogram(histClass, "Mass_Pt_centrFT0C_DeltaPhiRP_FT0A", "", 5, varDeltaPhi_FT0A_RP, bins_DeltaPhi, minBins_DeltaPhi, maxBins_DeltaPhi, 0, -1, kTRUE); + } + } if (subGroupStr.Contains("upsilon")) { hm->AddHistogram(histClass, "MassUpsilon_Pt", "", false, 500, 7.0, 12.0, VarManager::kMass, 400, 0.0, 40.0, VarManager::kPt); } diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 1faf6c8731e..e94e703748c 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -76,6 +76,8 @@ 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 +TH3F* VarManager::fgObjQvec = nullptr; +bool VarManager::fgApplyQVectorCorrection = false; //__________________________________________________________________ VarManager::VarManager() : TObject() @@ -2484,6 +2486,12 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kDCATrackVtxProd"] = kDCATrackVtxProd; fgVarNamesMap["kV2SP"] = kV2SP; fgVarNamesMap["kV2EP"] = kV2EP; + fgVarNamesMap["kA2EP_PP_TPC"] = kA2EP_PP_TPC; + fgVarNamesMap["kA2EP_PP_FT0A"] = kA2EP_PP_FT0A; + fgVarNamesMap["kA2EP_PP_FT0C"] = kA2EP_PP_FT0C; + fgVarNamesMap["kA2EP_RP_TPC"] = kA2EP_RP_TPC; + fgVarNamesMap["kA2EP_RP_FT0A"] = kA2EP_RP_FT0A; + fgVarNamesMap["kA2EP_RP_FT0C"] = kA2EP_RP_FT0C; fgVarNamesMap["kWV2SP"] = kWV2SP; fgVarNamesMap["kWV2EP"] = kWV2EP; fgVarNamesMap["kU2Q2"] = kU2Q2; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index f2b48e3d71f..d044ec90846 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -384,6 +384,9 @@ class VarManager : public TObject kMultA, // Multiplicity of the sub-event A kMultAPOS, // Multiplicity of the sub-event A kMultANEG, // Multiplicity of the sub-event A + kMultA1, + kMultA2, + kNnorm, kMultB, kMultC, kQ3X0A, // q-vector (e.g. from TPC) with x component (harmonic 3 and power 0), sub-event A @@ -805,6 +808,22 @@ class VarManager : public TObject kDeltaPhiPair2, kDeltaEtaPair2, kPsiPair, + kDeltaPhiRP_TPC, + kDeltaPhiRP_FT0A, + kDeltaPhiRP_FT0C, + kCos2DeltaPhiRP_TPC, + kCos2DeltaPhiRP_FT0A, + kCos2DeltaPhiRP_FT0C, + kDeltaPhiPP_TPC, + kDeltaPhiPP_FT0A, + kDeltaPhiPP_FT0C, + kCos2DeltaPhiPP_TPC, + kCos2DeltaPhiPP_FT0A, + kCos2DeltaPhiPP_FT0C, + kNullA2, + kInfA2, + kAmbi1, + kAmbi2, kDeltaPhiPair, kOpeningAngle, kQuadDCAabsXY, @@ -820,6 +839,12 @@ class VarManager : public TObject kDCATrackVtxProd, kV2SP, kV2EP, + kA2EP_RP_TPC, + kA2EP_RP_FT0A, + kA2EP_RP_FT0C, + kA2EP_PP_TPC, + kA2EP_PP_FT0A, + kA2EP_PP_FT0C, kWV2SP, kWV2EP, kU2Q2, @@ -1523,6 +1548,11 @@ class VarManager : public TObject fgEOR = eor; } + static void SetEventQVectorCorrection(TH3F* qvec) { + fgObjQvec = qvec; + fgApplyQVectorCorrection = true; + } + public: VarManager(); ~VarManager() override; @@ -1566,6 +1596,8 @@ class VarManager : public TObject static float calculatePhiV(const T1& t1, const T2& t2); template static float LorentzTransformJpsihadroncosChi(TString Option, const T1& v1, const T2& v2); + template + static bool isSelectedinTPC(const T& t); static o2::vertexing::DCAFitterN<2> fgFitterTwoProngBarrel; static o2::vertexing::DCAFitterN<3> fgFitterThreeProngBarrel; @@ -1581,6 +1613,8 @@ class VarManager : public TObject static int fgEfficiencyType; // type of efficiency correction to apply static TObject* fgEfficiencyHist; // histogram for efficiency correction + static TH3F* fgObjQvec; + static bool fgApplyQVectorCorrection; VarManager& operator=(const VarManager& c); VarManager(const VarManager& c); @@ -2388,6 +2422,101 @@ void VarManager::FillEvent(T const& event, float* values) values[kWR2EP_BC_Im] = std::isnan(R2EP_BC_Im) || std::isinf(R2EP_BC_Im) ? 0. : 1.0; } + if constexpr ((fillMap & CollisionQvectCentr) > 0) { + values[kQ1X0A] = -999; + values[kQ1Y0A] = -999; + values[kQ1X0B] = -999; + values[kQ1Y0B] = -999; + values[kQ1X0C] = -999; + values[kQ1Y0C] = -999; + values[kQ2X0A] = event.qvecTPCallRe() / event.nTrkTPCall(); + values[kQ2Y0A] = event.qvecTPCallIm() / event.nTrkTPCall(); + values[kQ2X0APOS] = event.qvecTPCposRe(); + values[kQ2Y0APOS] = event.qvecTPCposIm(); + values[kQ2X0ANEG] = event.qvecTPCnegRe(); + values[kQ2Y0ANEG] = event.qvecTPCnegIm(); + values[kQ2X0B] = event.qvecFT0ARe(); + values[kQ2Y0B] = event.qvecFT0AIm(); + values[kQ2X0C] = event.qvecFT0CRe(); + values[kQ2Y0C] = event.qvecFT0CIm(); + values[kMultA] = event.nTrkTPCall(); + values[kMultAPOS] = event.nTrkTPCpos(); + values[kMultANEG] = event.nTrkTPCneg(); + values[kMultB] = event.sumAmplFT0A(); + values[kMultC] = event.sumAmplFT0C(); + values[kQ3X0A] = -999; + values[kQ3Y0A] = -999; + values[kQ3X0B] = -999; + values[kQ3Y0B] = -999; + values[kQ3X0C] = -999; + values[kQ3Y0C] = -999; + values[kQ4X0A] = -999; + values[kQ4Y0A] = -999; + values[kQ4X0B] = -999; + values[kQ4Y0B] = -999; + values[kQ4X0C] = -999; + values[kQ4Y0C] = -999; + + EventPlaneHelper epHelper; + float Psi2A = epHelper.GetEventPlane(values[kQ2X0A], values[kQ2Y0A], 2); + float Psi2APOS = epHelper.GetEventPlane(values[kQ2X0APOS], values[kQ2Y0APOS], 2); + float Psi2ANEG = epHelper.GetEventPlane(values[kQ2X0ANEG], values[kQ2Y0ANEG], 2); + float Psi2B = epHelper.GetEventPlane(values[kQ2X0B], values[kQ2Y0B], 2); + float Psi2C = epHelper.GetEventPlane(values[kQ2X0C], values[kQ2Y0C], 2); + + values[kPsi2A] = Psi2A; + values[kPsi2APOS] = Psi2APOS; + values[kPsi2ANEG] = Psi2ANEG; + values[kPsi2B] = Psi2B; + values[kPsi2C] = Psi2C; + + float R2SP_AB = (values[kQ2X0A] * values[kQ2X0B] + values[kQ2Y0A] * values[kQ2Y0B]); + float R2SP_APOSB = (values[kQ2X0APOS] * values[kQ2X0B] + values[kQ2Y0APOS] * values[kQ2Y0B]); + float R2SP_ANEGB = (values[kQ2X0ANEG] * values[kQ2X0B] + values[kQ2Y0ANEG] * values[kQ2Y0B]); + float R2SP_AC = (values[kQ2X0A] * values[kQ2X0C] + values[kQ2Y0A] * values[kQ2Y0C]); + float R2SP_APOSC = (values[kQ2X0APOS] * values[kQ2X0C] + values[kQ2Y0APOS] * values[kQ2Y0C]); + float R2SP_ANEGC = (values[kQ2X0ANEG] * values[kQ2X0C] + values[kQ2Y0ANEG] * values[kQ2Y0C]); + float R2SP_BC = (values[kQ2X0B] * values[kQ2X0C] + values[kQ2Y0B] * values[kQ2Y0C]); + float R2SP_AB_Im = (values[kQ2Y0A] * values[kQ2X0B] - values[kQ2X0A] * values[kQ2Y0B]); + float R2SP_AC_Im = (values[kQ2Y0A] * values[kQ2X0C] - values[kQ2X0A] * values[kQ2Y0C]); + float R2SP_BC_Im = (values[kQ2Y0B] * values[kQ2X0C] - values[kQ2X0B] * values[kQ2Y0C]); + values[kR2SP_AB] = std::isnan(R2SP_AB) || std::isinf(R2SP_AB) ? 0. : R2SP_AB; + values[kR2SP_FT0CTPCPOS] = std::isnan(R2SP_APOSB) || std::isinf(R2SP_APOSB) ? 0. : R2SP_APOSB; + values[kR2SP_FT0CTPCNEG] = std::isnan(R2SP_ANEGB) || std::isinf(R2SP_ANEGB) ? 0. : R2SP_ANEGB; + values[kWR2SP_AB] = std::isnan(R2SP_AB) || std::isinf(R2SP_AB) ? 0. : 1.0; + values[kR2SP_AC] = std::isnan(R2SP_AC) || std::isinf(R2SP_AC) ? 0. : R2SP_AC; + values[kR2SP_FT0ATPCPOS] = std::isnan(R2SP_APOSC) || std::isinf(R2SP_APOSC) ? 0. : R2SP_APOSC; + values[kR2SP_FT0ATPCNEG] = std::isnan(R2SP_ANEGC) || std::isinf(R2SP_ANEGC) ? 0. : R2SP_ANEGC; + values[kWR2SP_AC] = std::isnan(R2SP_AC) || std::isinf(R2SP_AC) ? 0. : 1.0; + values[kR2SP_BC] = std::isnan(R2SP_BC) || std::isinf(R2SP_BC) ? 0. : R2SP_BC; + values[kWR2SP_BC] = std::isnan(R2SP_BC) || std::isinf(R2SP_BC) ? 0. : 1.0; + values[kR2SP_AB_Im] = std::isnan(R2SP_AB_Im) || std::isinf(R2SP_AB_Im) ? 0. : R2SP_AB_Im; + values[kWR2SP_AB_Im] = std::isnan(R2SP_AB_Im) || std::isinf(R2SP_AB_Im) ? 0. : 1.0; + values[kR2SP_AC_Im] = std::isnan(R2SP_AC_Im) || std::isinf(R2SP_AC_Im) ? 0. : R2SP_AC_Im; + values[kWR2SP_AC_Im] = std::isnan(R2SP_AC_Im) || std::isinf(R2SP_AC_Im) ? 0. : 1.0; + values[kR2SP_BC_Im] = std::isnan(R2SP_BC_Im) || std::isinf(R2SP_BC_Im) ? 0. : R2SP_BC_Im; + values[kWR2SP_BC_Im] = std::isnan(R2SP_BC_Im) || std::isinf(R2SP_BC_Im) ? 0. : 1.0; + + float R2EP_AB = std::isnan(Psi2A) || std::isinf(Psi2A) || std::isnan(Psi2B) || std::isinf(Psi2B) ? 0. : TMath::Cos(2 * (Psi2A - Psi2B)); + float R2EP_AC = std::isnan(Psi2A) || std::isinf(Psi2A) || std::isnan(Psi2C) || std::isinf(Psi2C) ? 0. : TMath::Cos(2 * (Psi2A - Psi2C)); + float R2EP_BC = std::isnan(Psi2B) || std::isinf(Psi2B) || std::isnan(Psi2C) || std::isinf(Psi2C) ? 0. : TMath::Cos(2 * (Psi2B - Psi2C)); + float R2EP_AB_Im = std::isnan(Psi2A) || std::isinf(Psi2A) || std::isnan(Psi2B) || std::isinf(Psi2B) ? 0. : TMath::Sin(2 * (Psi2A - Psi2B)); + float R2EP_AC_Im = std::isnan(Psi2A) || std::isinf(Psi2A) || std::isnan(Psi2C) || std::isinf(Psi2C) ? 0. : TMath::Sin(2 * (Psi2A - Psi2C)); + float R2EP_BC_Im = std::isnan(Psi2B) || std::isinf(Psi2B) || std::isnan(Psi2C) || std::isinf(Psi2C) ? 0. : TMath::Sin(2 * (Psi2B - Psi2C)); + values[kR2EP_AB] = std::isnan(R2EP_AB) || std::isinf(R2EP_AB) ? 0. : R2EP_AB; + values[kWR2EP_AB] = std::isnan(R2EP_AB) || std::isinf(R2EP_AB) ? 0. : 1.0; + values[kR2EP_AC] = std::isnan(R2EP_AC) || std::isinf(R2EP_AC) ? 0. : R2EP_AC; + values[kWR2EP_AC] = std::isnan(R2EP_AC) || std::isinf(R2EP_AC) ? 0. : 1.0; + values[kR2EP_BC] = std::isnan(R2EP_BC) || std::isinf(R2EP_BC) ? 0. : R2EP_BC; + values[kWR2EP_BC] = std::isnan(R2EP_BC) || std::isinf(R2EP_BC) ? 0. : 1.0; + values[kR2EP_AB_Im] = std::isnan(R2EP_AB_Im) || std::isinf(R2EP_AB_Im) ? 0. : R2EP_AB_Im; + values[kWR2EP_AB_Im] = std::isnan(R2EP_AB_Im) || std::isinf(R2EP_AB_Im) ? 0. : 1.0; + values[kR2EP_AC_Im] = std::isnan(R2EP_AC_Im) || std::isinf(R2EP_AC_Im) ? 0. : R2EP_AC_Im; + values[kWR2EP_AC_Im] = std::isnan(R2EP_AC_Im) || std::isinf(R2EP_AC_Im) ? 0. : 1.0; + values[kR2EP_BC_Im] = std::isnan(R2EP_BC_Im) || std::isinf(R2EP_BC_Im) ? 0. : R2EP_BC_Im; + values[kWR2EP_BC_Im] = std::isnan(R2EP_BC_Im) || std::isinf(R2EP_BC_Im) ? 0. : 1.0; + } + if constexpr ((fillMap & CollisionMC) > 0) { values[kMCEventGeneratorId] = event.generatorsID(); values[kMCEventSubGeneratorId] = event.getSubGeneratorId(); @@ -5945,6 +6074,129 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) values[kCos2ThetaStarFT0C] = values[kCosThetaStarFT0C] * values[kCosThetaStarFT0C]; } + // Coherent Jpsi A2 + bool useCoherentJpsiA2 = fgUsedVars[kA2EP_RP_TPC] || fgUsedVars[kA2EP_RP_FT0A] || fgUsedVars[kA2EP_RP_FT0C] || fgUsedVars[kA2EP_PP_TPC] || fgUsedVars[kA2EP_PP_FT0A] || fgUsedVars[kA2EP_PP_FT0C]; + if (useCoherentJpsiA2) { + // remove daughter from TPC Q-vector + // TODO: remove based on track cut in qVectorTable + float Q2X0A = values[kQ2X0A]*values[kMultA]; + float Q2Y0A = values[kQ2Y0A]*values[kMultA]; + float nNorm = values[kMultA]; + + int centBin = static_cast(values[kCentFT0C]) + 1; // centrality bin + int detIdx = 6; // TPC all + EventPlaneHelper epHelper; + // checkTrack(t1); + if (isSelectedinTPC(t1) && values[kAmbi1] > 0) { + float qx1 = t1.pt()*TMath::Cos(2. * t1.phi()); + float qy1 = t1.pt()*TMath::Sin(2. * t1.phi()); + if (fgApplyQVectorCorrection) { + epHelper.DoRecenter(qx1, qy1, fgObjQvec->GetBinContent(centBin, 1, detIdx + 1), fgObjQvec->GetBinContent(centBin, 2, detIdx + 1)); + epHelper.DoTwist(qx1, qy1, fgObjQvec->GetBinContent(centBin, 3, detIdx + 1), fgObjQvec->GetBinContent(centBin, 4, detIdx + 1)); + epHelper.DoRescale(qx1, qy1, fgObjQvec->GetBinContent(centBin, 5, detIdx + 1), fgObjQvec->GetBinContent(centBin, 6, detIdx + 1)); + } + Q2X0A = Q2X0A - qx1; + Q2Y0A = Q2Y0A - qy1; + nNorm = nNorm - 1.; + } + // checkTrack(t2); + if (isSelectedinTPC(t2) && values[kAmbi2] > 0) { + float qx2 = t2.pt()*TMath::Cos(2. * t2.phi()); + float qy2 = t2.pt()*TMath::Sin(2. * t2.phi()); + if (fgApplyQVectorCorrection) { + epHelper.DoRecenter(qx2, qy2, fgObjQvec->GetBinContent(centBin, 1, detIdx + 1), fgObjQvec->GetBinContent(centBin, 2, detIdx + 1)); + epHelper.DoTwist(qx2, qy2, fgObjQvec->GetBinContent(centBin, 3, detIdx + 1), fgObjQvec->GetBinContent(centBin, 4, detIdx + 1)); + epHelper.DoRescale(qx2, qy2, fgObjQvec->GetBinContent(centBin, 5, detIdx + 1), fgObjQvec->GetBinContent(centBin, 6, detIdx + 1)); + } + Q2X0A = Q2X0A - qx2; + Q2Y0A = Q2Y0A - qy2; + nNorm = nNorm - 1.; + } + values[kNnorm] = nNorm; + if (nNorm <= 0) { + values[kQ2X0A] = -999.; + values[kQ2Y0A] = -999.; + return; + } + + Q2X0A = nNorm > 0 ? Q2X0A/nNorm : NAN; + Q2Y0A = nNorm > 0 ? Q2Y0A/nNorm : NAN; + + float Psi2A = getEventPlane(2, Q2X0A, Q2Y0A); + values[kPsi2A] = Psi2A; + + // pT ~ 0.2, non-relativistic + // ROOT::Math::PtEtaPhiMVector v_daughter = t1.sign() > 0 ? v1 - v2 : v2 - v1; + // boost to Jpsi rest frame, then calculate the angle with respect to the event plane + ROOT::Math::Boost boostv12{v12.BoostToCM()}; + ROOT::Math::PtEtaPhiMVector v_daughter = boostv12(t1.sign() > 0 ? v1 : v2); + + // production plane + ROOT::Math::XYZVectorF zAxis_RF{(v12.Vect()).Unit()}; + ROOT::Math::XYZVectorF zAxis{fgBeamA.Vect().Unit()}; + ROOT::Math::XYZVectorF yAxis_RF = zAxis_RF.Cross(zAxis).Unit(); + ROOT::Math::XYZVectorF xAxis_RF = yAxis_RF.Cross(zAxis_RF).Unit(); + ROOT::Math::XYZVectorF daughterVec_RF{(v_daughter.Vect()).Unit()}; + ROOT::Math::XYZVectorF b_TPC_RF = ROOT::Math::XYZVectorF(std::cos(Psi2A), std::sin(Psi2A), 0.f); + ROOT::Math::XYZVectorF b_FT0A_RF = ROOT::Math::XYZVectorF(std::cos(Psi2B), std::sin(Psi2B), 0.f); + ROOT::Math::XYZVectorF b_FT0C_RF = ROOT::Math::XYZVectorF(std::cos(Psi2C), std::sin(Psi2C), 0.f); + float cosPhi = yAxis_RF.Dot(zAxis_RF.Cross(daughterVec_RF)); + float sinPhi = -1. * xAxis_RF.Dot(zAxis_RF.Cross(daughterVec_RF)); + float phi_PP = sinPhi > 0 ? TMath::ACos(cosPhi) : -1. * TMath::ACos(cosPhi); + float cosPsi2APP = b_TPC_RF.Dot(xAxis_RF.Cross(daughterVec_RF)); + float sinPsi2APP = b_TPC_RF.Dot(yAxis_RF.Cross(daughterVec_RF)); + float Psi2APP = sinPsi2APP > 0 ? TMath::ACos(cosPsi2APP) : -1. * TMath::ACos(cosPsi2APP); + float cosPsi2BPP = b_FT0A_RF.Dot(xAxis_RF.Cross(daughterVec_RF)); + float sinPsi2BPP = b_FT0A_RF.Dot(yAxis_RF.Cross(daughterVec_RF)); + float Psi2BPP = sinPsi2BPP > 0 ? TMath::ACos(cosPsi2BPP) : -1. * TMath::ACos(cosPsi2BPP); + float cosPsi2CPP = b_FT0C_RF.Dot(xAxis_RF.Cross(daughterVec_RF)); + float sinPsi2CPP = b_FT0C_RF.Dot(yAxis_RF.Cross(daughterVec_RF)); + float Psi2CPP = sinPsi2CPP > 0 ? TMath::ACos(cosPsi2CPP) : -1. * TMath::ACos(cosPsi2CPP); + values[kDeltaPhiPP_TPC] = phi_PP > Psi2APP ? phi_PP - Psi2APP : Psi2APP - phi_PP; + values[kDeltaPhiPP_TPC] = values[kDeltaPhiPP_TPC] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_TPC] : values[kDeltaPhiPP_TPC]; + values[kDeltaPhiPP_FT0A] = phi_PP > Psi2BPP ? phi_PP - Psi2BPP : Psi2BPP - phi_PP; + values[kDeltaPhiPP_FT0A] = values[kDeltaPhiPP_FT0A] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_FT0A] : values[kDeltaPhiPP_FT0A]; + values[kDeltaPhiPP_FT0C] = phi_PP > Psi2CPP ? phi_PP - Psi2CPP : Psi2CPP - phi_PP; + values[kDeltaPhiPP_FT0C] = values[kDeltaPhiPP_FT0C] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_FT0C] : values[kDeltaPhiPP_FT0C]; + // fold delta phi to [0, pi/2] + values[kDeltaPhiPP_TPC] = values[kDeltaPhiPP_TPC] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiPP_TPC] : values[kDeltaPhiPP_TPC]; + values[kDeltaPhiPP_FT0A] = values[kDeltaPhiPP_FT0A] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiPP_FT0A] : values[kDeltaPhiPP_FT0A]; + values[kDeltaPhiPP_FT0C] = values[kDeltaPhiPP_FT0C] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiPP_FT0C] : values[kDeltaPhiPP_FT0C]; + values[kCos2DeltaPhiPP_TPC] = TMath::Cos(2. * (phi_PP - Psi2A)); + values[kCos2DeltaPhiPP_FT0A] = TMath::Cos(2. * (phi_PP - Psi2B)); + values[kCos2DeltaPhiPP_FT0C] = TMath::Cos(2. * (phi_PP - Psi2C)); + + float A2PP_TPC = values[kCos2DeltaPhiPP_TPC] / values[kR2EP]; + float A2PP_FT0A = values[kCos2DeltaPhiPP_FT0A] / values[kR2EP]; + float A2PP_FT0C = values[kCos2DeltaPhiPP_FT0C] / values[kR2EP]; + values[kA2EP_PP_TPC] = std::isnan(A2PP_TPC) || std::isinf(A2PP_TPC) ? -999. : A2PP_TPC; + values[kA2EP_PP_FT0A] = std::isnan(A2PP_FT0A) || std::isinf(A2PP_FT0A) ? -999. : A2PP_FT0A; + values[kA2EP_PP_FT0C] = std::isnan(A2PP_FT0C) || std::isinf(A2PP_FT0C) ? -999. : A2PP_FT0C; + + // reaction plane + float phi = v_daughter.Phi() > TMath::Pi() ? 2. * TMath::Pi() - v_daughter.Phi() : v_daughter.Phi(); + values[kDeltaPhiRP_TPC] = phi > Psi2A ? phi - Psi2A : Psi2A - phi; + values[kDeltaPhiRP_TPC] = values[kDeltaPhiRP_TPC] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_TPC] : values[kDeltaPhiRP_TPC]; + values[kDeltaPhiRP_FT0A] = phi > Psi2B ? phi - Psi2B : Psi2B - phi; + values[kDeltaPhiRP_FT0A] = values[kDeltaPhiRP_FT0A] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_FT0A] : values[kDeltaPhiRP_FT0A]; + values[kDeltaPhiRP_FT0C] = phi > Psi2C ? phi - Psi2C : Psi2C - phi; + values[kDeltaPhiRP_FT0C] = values[kDeltaPhiRP_FT0C] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_FT0C] : values[kDeltaPhiRP_FT0C]; + // fold delta phi to [0, pi/2] + values[kDeltaPhiRP_TPC] = values[kDeltaPhiRP_TPC] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_TPC] : values[kDeltaPhiRP_TPC]; + values[kDeltaPhiRP_FT0A] = values[kDeltaPhiRP_FT0A] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_FT0A] : values[kDeltaPhiRP_FT0A]; + values[kDeltaPhiRP_FT0C] = values[kDeltaPhiRP_FT0C] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_FT0C] : values[kDeltaPhiRP_FT0C]; + values[kCos2DeltaPhiRP_TPC] = TMath::Cos(2. * (phi - Psi2A)); + values[kCos2DeltaPhiRP_FT0A] = TMath::Cos(2. * (phi - Psi2B)); + values[kCos2DeltaPhiRP_FT0C] = TMath::Cos(2. * (phi - Psi2C)); + + float A2RP_TPC = values[kCos2DeltaPhiRP_TPC] / values[kR2EP]; + float A2RP_FT0A = values[kCos2DeltaPhiRP_FT0A] / values[kR2EP]; + float A2RP_FT0C = values[kCos2DeltaPhiRP_FT0C] / values[kR2EP]; + values[kA2EP_RP_TPC] = std::isnan(A2RP_TPC) || std::isinf(A2RP_TPC) ? -999. : A2RP_TPC; + values[kA2EP_RP_FT0A] = std::isnan(A2RP_FT0A) || std::isinf(A2RP_FT0A) ? -999. : A2RP_FT0A; + values[kA2EP_RP_FT0C] = std::isnan(A2RP_FT0C) || std::isinf(A2RP_FT0C) ? -999. : A2RP_FT0C; + } + // kV4, kC4POI, kC4REF etc. if constexpr ((fillMap & ReducedEventQvectorExtra) > 0) { std::complex Q21(values[kQ2X0A] * values[kS11A], values[kQ2Y0A] * values[kS11A]); @@ -6790,6 +7042,48 @@ float VarManager::LorentzTransformJpsihadroncosChi(TString Option, T1 const& v1, } return value; } + +template +bool VarManager::isSelectedinTPC(const T& t) { + bool trackSelected = false; + if constexpr (requires { t.hasTPC(); }) { + trackSelected = true; + if (!t.hasTPC()) { + trackSelected = false; + return trackSelected; + } + if (std::abs(t.eta()) > 0.8) { // eta cut used in q vector table. Todo: read from config + trackSelected = false; + return trackSelected; + } + if (t.pt() < 0.15 || t.pt() > 5.) { // pt cut used in q vector table. Todo: read from config + trackSelected = false; + return trackSelected; + } + if (t.tpcNClsCrossedRows() / t.tpcNClsFound() < 0.8 || t.tpcChi2NCl() > 4.) { + trackSelected = false; + return trackSelected; + } + if (!((t.itsClusterMap() & (1 << uint8_t(0))) > 0 || (t.itsClusterMap() & (1 << uint8_t(1))) > 0 || (t.itsClusterMap() & (1 << uint8_t(2))) > 0)) { + trackSelected = false; + return trackSelected; + } + if (t.itsChi2NCl() > 36.) { + trackSelected = false; + return trackSelected; + } + if (t.dcaZ() > 2) { + trackSelected = false; + return trackSelected; + } + if (t.dcaXY() > 0.0105 + 0.0350 / std::pow(t.pt(), 1.1)) { + trackSelected = false; + return trackSelected; + } + } + return trackSelected; +} + template void VarManager::FillFIT(T1 const& bc, T2 const& bcs, T3 const& ft0s, T4 const& fv0as, T5 const& fdds, float* values) { diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index 9e61737c6dc..ac99679f996 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -240,7 +240,7 @@ using MyEventsVtxCovZdcFitSelected = soa::Join; using MyEventsQvector = soa::Join; using MyEventsHashSelectedQvector = soa::Join; -using MyEventsQvectorCentr = soa::Join; +using MyEventsQvectorCentr = soa::Join; using MyEventsQvectorCentrSelected = soa::Join; using MyEventsHashSelectedQvectorCentr = soa::Join; @@ -1301,6 +1301,9 @@ struct AnalysisSameEventPairing { Produces dileptonEventInfoList; o2::base::MatLayerCylSet* fLUT = nullptr; + TH1D* ResoFlowSP = nullptr; + TH1D* ResoFlowEP = nullptr; + TH3F* qvecObj = nullptr; int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. OutputObj fOutputList{"output"}; @@ -1334,6 +1337,9 @@ struct AnalysisSameEventPairing { Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; 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 flowPathLocal{"flowPathLocal", "/lustre/alice/users/ywang/calib/FlowReso.root", "Path to the flow resolution object in the local cache"}; + Configurable QvecCalibPath{"cfgQvecCalibPath", "Users/j/junlee/Qvector/Pass5/QvecCalib/v2", "Path to the q vector calibration object"}; } fConfigCCDB; struct : ConfigurableGroup { @@ -1352,6 +1358,10 @@ struct AnalysisSameEventPairing { Configurable useRemoteCollisionInfo{"cfgUseRemoteCollisionInfo", false, "Use remote collision information from CCDB"}; 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 useRemoteFlow{"cfgUseRemoteFlow", false, "Use remote flow information from CCDB"}; + Configurable useLocalFlow{"cfgUseLocalFlow", false, "Use flow information from local cache"}; + Configurable useQvecCalib{"cfgUseQvecCalib", false, "Use flow correction factors for Q-vector recalibration when removing the daughter"}; + Configurable useCorrectionForRun{"cfgUseCorrectionForRun", false, "Apply run-by-run correction factors to the flow vectors"}; } fConfigOptions; struct : ConfigurableGroup { Configurable applyBDT{"applyBDT", false, "Flag to apply ML selections"}; @@ -1797,6 +1807,37 @@ struct AnalysisSameEventPairing { auto effList = fCCDB->getForTimeStamp(fConfigCCDB.efficiencyPath.value, timestamp); VarManager::SetEfficiencyObject(fConfigOptions.efficiencyType.value, effList->FindObject("efficiency")); } + + if (fConfigOptions.useRemoteFlow) { + TString PathFlow = fConfigCCDB.flowPath.value; + TString ccdbPathFlowSP = Form("%s/ScalarProduct", PathFlow.Data()); + TString ccdbPathFlowEP = Form("%s/EventPlane", PathFlow.Data()); + ResoFlowSP = fCCDB->getForTimeStamp(ccdbPathFlowSP.Data(), timestamp); + ResoFlowEP = fCCDB->getForTimeStamp(ccdbPathFlowEP.Data(), timestamp); + if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { + LOGF(fatal, "Flow resolution histograms not available in CCDB at timestamp=%llu", timestamp); + } + } else if (fConfigOptions.useLocalFlow) { + TString pathFlow = fConfigCCDB.flowPathLocal.value; + TFile* fileFlow = TFile::Open(pathFlow.Data(), "READ"); + if (fileFlow == nullptr || fileFlow->IsZombie()) { + LOGF(fatal, "Flow resolution file %s cannot be opened", pathFlow.Data()); + } + fileFlow->GetObject("ScalarProduct", ResoFlowSP); + fileFlow->GetObject("EventPlane", ResoFlowEP); + if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { + LOGF(fatal, "Flow resolution histograms not available in file %s", pathFlow.Data()); + } + } + + if (fConfigOptions.useQvecCalib) { + TString pathQvecCalib = fConfigCCDB.QvecCalibPath.value; + if (fConfigOptions.useCorrectionForRun) { + qvecObj = fCCDB->getForRun(pathQvecCalib.Data(), runNumber); + } else { + qvecObj = fCCDB->getForTimeStamp(pathQvecCalib.Data(), timestamp); + } + } } // Template function to run same event pairing (barrel-barrel, muon-muon, barrel-muon) @@ -1900,6 +1941,19 @@ struct AnalysisSameEventPairing { if (groupedAssocs.size() == 0) { continue; } + + if (fillFlowReso) { + if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { + LOGF(fatal, "Flow resolution histograms are not available, cannot fill flow variables!"); + } + VarManager::FillEventFlowResoFactor(ResoFlowSP, ResoFlowEP); + if (fConfigOptions.useQvecCalib) { + if (qvecObj == nullptr) { + LOGF(fatal, "Q-vector calibration object is not available, cannot fill flow variables!"); + } + VarManager::SetEventQVectorCorrection(qvecObj); + } + } bool isFirst = true; for (auto& [a1, a2] : o2::soa::combinations(groupedAssocs, groupedAssocs)) { @@ -1930,6 +1984,13 @@ struct AnalysisSameEventPairing { } fNPairPerEvent++; + + if (t1.reducedeventId() != event.globalIndex()) { + VarManager::fgValues[VarManager::kAmbi1] = 1.; + } + if (t2.reducedeventId() != event.globalIndex()) { + VarManager::fgValues[VarManager::kAmbi2] = 1.; + } VarManager::FillPair(t1, t2); // compute quantities which depend on the associated collision, such as DCA if (fConfigOptions.propTrack) { @@ -2576,6 +2637,10 @@ struct AnalysisSameEventPairing { template void runSameSideMixing(TEvents& events, TAssocs const& assocs, TTracks const& tracks, Preslice& preSlice) { + LOG(info) << "Flow resolution objects not set, flow will not be filled for mixed events"; + if (events.size() > 0) { + initParamsFromCCDB(events.begin().timestamp(), events.begin().runNumber(), false); + } events.bindExternalIndices(&assocs); int mixingDepth = fConfigMixingDepth.value; fAmbiguousPairs.clear(); @@ -2590,6 +2655,10 @@ struct AnalysisSameEventPairing { assocs2.bindExternalIndices(&events); fNPairPerEvent = 0; + VarManager::FillTwoMixEvents(event1, event2, assocs1, assocs2); + if (fConfigOptions.useRemoteFlow || fConfigOptions.useLocalFlow) { + VarManager::FillTwoMixEventsFlowResoFactor(ResoFlowSP, ResoFlowEP); + } runMixedPairing(assocs1, assocs2, tracks, tracks); VarManager::fgValues[VarManager::kNPairsPerEvent] = fNPairPerEvent; if (fEnableBarrelMixingHistos && fConfigQA) { diff --git a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx index 13e6b82c889..104f213195d 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx @@ -223,8 +223,11 @@ DECLARE_SOA_TABLE(BmesonCandidates, "AOD", "DQBMESONS", // Using definitions (data-only) using MyEvents = soa::Join; +using MyEventsQvectorCentr = soa::Join; using MyEventsSelected = soa::Join; +using MyEventsQvectorCentrSelected = soa::Join; using MyEventsHashSelected = soa::Join; +using MyEventsHashSelectedQvectorCentr = soa::Join; using MyBarrelTracksWithCov = soa::Join(events); } + void processDirectWithQvectorCentr(MyEventsQvectorCentr const& events, BCsWithTimestamps const& bcs) + { + runEventSelection(events, bcs); + publishSelections(events); + } + void processDummy(aod::Collisions&) {} PROCESS_SWITCH(AnalysisEventSelection, processDirect, "Run event selection on framework AO2Ds", false); + PROCESS_SWITCH(AnalysisEventSelection, processDirectWithQvectorCentr, "Run event selection on skimmed data with Q-vector and centrality", false); PROCESS_SWITCH(AnalysisEventSelection, processDummy, "Dummy function", true); }; @@ -857,7 +869,7 @@ struct AnalysisTrackSelection { void processWithCov(TrackAssoc const& assocs, BCsWithTimestamps const& bcs, MyEventsSelected const& events, MyBarrelTracksWithCov const& tracks) { - runTrackSelection(assocs, bcs, events, tracks); + runTrackSelection(assocs, bcs, events, tracks); } void processWithCovTOFService(TrackAssoc const& assocs, BCsWithTimestamps const& bcs, MyEventsSelected const& events, MyBarrelTracksWithCovNoTOF const& tracks) @@ -1275,6 +1287,10 @@ struct AnalysisSameEventPairing { Produces electronmuonList; o2::base::MatLayerCylSet* fLUT = nullptr; + TH1D* ResoFlowSP = nullptr; + TH1D* ResoFlowEP = nullptr; + TH3F* qvecObj = nullptr; + TString pathQvecCalib; int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. OutputObj fOutputList{"output"}; @@ -1309,6 +1325,10 @@ struct AnalysisSameEventPairing { Configurable fConfigMiniTree{"cfgMiniTree", false, "Produce a single flat table with minimal information for analysis"}; Configurable fConfigMiniTreeMinMass{"cfgMiniTreeMinMass", 2, "Min. mass cut for minitree"}; Configurable fConfigMiniTreeMaxMass{"cfgMiniTreeMaxMass", 5, "Max. mass cut for minitree"}; + Configurable useRemoteFlow{"cfgUseRemoteFlow", false, "Use remote flow information from CCDB"}; + Configurable useLocalFlow{"cfgUseLocalFlow", false, "Use flow information from local cache"}; + Configurable useQvecCalib{"cfgUseQvecCalib", false, "Use flow correction factors for Q-vector recalibration when removing the daughter"}; + Configurable useCorrectionForRun{"cfgUseCorrectionForRun", false, "Apply run-by-run correction factors to the flow vectors"}; } fConfigOptions; struct : ConfigurableGroup { @@ -1316,6 +1336,9 @@ struct AnalysisSameEventPairing { Configurable grpMagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + Configurable flowPath{"flowPath", "Users/y/yiping/FlowResolution", "Path to the flow resolution object"}; + Configurable flowPathLocal{"flowPathLocal", "/lustre/alice/users/ywang/calib/FlowReso.root", "Path to the flow resolution object in the local cache"}; + Configurable QvecCalibPath{"cfgQvecCalibPath", "Users/j/junlee/Qvector/Pass5/QvecShift/v2", "Path to the q vector calibration object"}; } fConfigCCDB; Service fCCDB; @@ -1358,7 +1381,7 @@ struct AnalysisSameEventPairing { } VarManager::SetDefaultVarNames(); - fEnableBarrelHistos = context.mOptions.get("processBarrelOnly"); + fEnableBarrelHistos = context.mOptions.get("processBarrelOnly") || context.mOptions.get("processBarrelOnlyWithQvectorCentr"); fEnableBarrelMuonHistos = context.mOptions.get("processElectronMuonDirect"); // Keep track of all the histogram class names to avoid composing strings in the pairing loop @@ -1526,9 +1549,15 @@ struct AnalysisSameEventPairing { dqhistograms::AddHistogramsFromJSON(fHistMan, fConfigOptions.fConfigAddJSONHistograms.value.c_str()); // ad-hoc histograms via JSON VarManager::SetUseVars(fHistMan->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill fOutputList.setObject(fHistMan->GetMainHistogramList()); + + if (fConfigOptions.useQvecCalib) { + string tmppathQvecCalib; + getTaskOptionValue(context, "q-vectors-table", "cfgQvecCalibPath", tmppathQvecCalib, false); + pathQvecCalib = tmppathQvecCalib; + } } - void initParamsFromCCDB(uint64_t timestamp, bool withTwoProngFitter = true) + void initParamsFromCCDB(uint64_t timestamp, int runnumber, bool withTwoProngFitter = true) { if (fConfigOptions.useRemoteField.value) { o2::parameters::GRPMagField* grpmag = fCCDB->getForTimeStamp(fConfigCCDB.grpMagPath, timestamp); @@ -1563,6 +1592,41 @@ struct AnalysisSameEventPairing { VarManager::SetupTwoProngDCAFitter(fConfigOptions.magField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, fConfigOptions.useAbsDCA.value); // needed because take in varmanager Bz from fgFitterTwoProngBarrel for PhiV calculations } } + + if (fConfigOptions.useRemoteFlow) { + TString PathFlow = fConfigCCDB.flowPath.value; + TString ccdbPathFlowSP = Form("%s/ScalarProduct", PathFlow.Data()); + TString ccdbPathFlowEP = Form("%s/EventPlane", PathFlow.Data()); + ResoFlowSP = fCCDB->getForTimeStamp(ccdbPathFlowSP.Data(), timestamp); + ResoFlowEP = fCCDB->getForTimeStamp(ccdbPathFlowEP.Data(), timestamp); + if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { + LOGF(fatal, "Flow resolution histograms not available in CCDB at timestamp=%llu", timestamp); + } + } else if (fConfigOptions.useLocalFlow) { + // LOGP(info, "Developing"); + TString pathFlow = fConfigCCDB.flowPathLocal.value; + TFile* fileFlow = TFile::Open(pathFlow.Data(), "READ"); + if (fileFlow == nullptr || fileFlow->IsZombie()) { + LOGF(fatal, "Flow resolution file %s cannot be opened", pathFlow.Data()); + } + fileFlow->GetObject("ScalarProduct", ResoFlowSP); + fileFlow->GetObject("EventPlane", ResoFlowEP); + if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { + LOGF(fatal, "Flow resolution histograms not available in file %s", pathFlow.Data()); + } + } + + if (fConfigOptions.useQvecCalib) { + TString ccdbPathQvecCalib = Form("%s/v2", pathQvecCalib.Data()); + if (fConfigOptions.useCorrectionForRun) { + qvecObj = fCCDB->getForRun(ccdbPathQvecCalib.Data(), runnumber); // get the object with run dependence if correction for run is needed + } else { + qvecObj = fCCDB->getForTimeStamp(ccdbPathQvecCalib.Data(), timestamp); // get the object without time dependence if no correction for run is needed + } + if (qvecObj == nullptr) { + LOGF(fatal, "Q-vector calibration object not available in CCDB at timestamp=%llu", timestamp); + } + } } template @@ -1625,7 +1689,7 @@ struct AnalysisSameEventPairing { return; } if (fCurrentRun != bcs.begin().runNumber()) { - initParamsFromCCDB(bcs.begin().timestamp(), TTwoProngFitter); + initParamsFromCCDB(bcs.begin().timestamp(), bcs.begin().runNumber(), TTwoProngFitter); fCurrentRun = bcs.begin().runNumber(); } @@ -1662,8 +1726,9 @@ struct AnalysisSameEventPairing { dielectronAllList.reserve(reserveSize); } - constexpr bool eventHasQvector = ((TEventFillMap & VarManager::ObjTypes::CollisionQvect) > 0); + constexpr bool eventHasQvector = ((TEventFillMap & VarManager::ObjTypes::CollisionQvectCentr) > 0); constexpr bool trackHasCov = ((TTrackFillMap & VarManager::ObjTypes::TrackCov) > 0); + constexpr bool fillFlowReso = eventHasQvector; for (auto& event : events) { if (!event.isEventSelected_bit(0)) @@ -1680,6 +1745,19 @@ struct AnalysisSameEventPairing { if (groupedAssocs.size() == 0) continue; + if (fillFlowReso) { + if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { + LOGF(fatal, "Flow resolution histograms are not available, cannot fill flow variables!"); + } + VarManager::FillEventFlowResoFactor(ResoFlowSP, ResoFlowEP); + if (fConfigOptions.useQvecCalib) { + if (qvecObj == nullptr) { + LOGF(fatal, "Q-vector calibration object is not available, cannot fill flow variables!"); + } + VarManager::SetEventQVectorCorrection(qvecObj); + } + } + for (auto& [a1, a2] : o2::soa::combinations(groupedAssocs, groupedAssocs)) { if constexpr (TPairType == VarManager::kDecayToEE) { @@ -1755,7 +1833,7 @@ struct AnalysisSameEventPairing { } if constexpr (eventHasQvector) { - VarManager::FillPairVn(t1, t2); + VarManager::FillPairVn(t1, t2); } } // Fill normal histograms @@ -1947,10 +2025,18 @@ struct AnalysisSameEventPairing { muonAssocsPerCollision, muonAssocs, muons); } + void processBarrelOnlyWithQvectorCentr(MyEventsQvectorCentrSelected const& events, BCsWithTimestamps const& bcs, + soa::Join const& barrelAssocs, + MyBarrelTracksWithCovWithAmbiguities const& barrelTracks) + { + runSameEventPairing(events, bcs, trackAssocsPerCollision, barrelAssocs, barrelTracks); + } + void processDummy(MyEvents&) { /* do nothing */ } PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnly, "Run barrel only pairing", false); PROCESS_SWITCH(AnalysisSameEventPairing, processElectronMuonDirect, "Run electron-muon pairing on AO2D tracks/fwd-tracks", false); + PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlyWithQvectorCentr, "Run barrel only pairing with Q-vector and centrality", false); PROCESS_SWITCH(AnalysisSameEventPairing, processDummy, "Dummy function", true); }; From 27f19696f91c2c4000c6991ba8e38cc8b1e9d789 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Mon, 15 Jun 2026 04:05:47 +0800 Subject: [PATCH 2/8] fix include bug --- PWGDQ/Core/VarManager.h | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index d044ec90846..63590cd8a05 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -65,6 +65,7 @@ #include #include #include +#include #include #include From 4c6b0aee766c5e3f3309f80b87e7367978eae36b Mon Sep 17 00:00:00 2001 From: ypwangg Date: Wed, 17 Jun 2026 00:15:22 +0800 Subject: [PATCH 3/8] fix bugs --- PWGDQ/Tasks/tableReader_withAssoc.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index ac99679f996..849b47a6ab6 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -49,6 +49,7 @@ #include #include #include +#include #include #include #include @@ -1921,6 +1922,7 @@ struct AnalysisSameEventPairing { constexpr bool eventHasQvector = ((TEventFillMap & VarManager::ObjTypes::ReducedEventQvector) > 0); constexpr bool eventHasQvectorCentr = ((TEventFillMap & VarManager::ObjTypes::CollisionQvect) > 0); constexpr bool trackHasCov = ((TTrackFillMap & VarManager::ObjTypes::TrackCov) > 0 || (TTrackFillMap & VarManager::ObjTypes::ReducedTrackBarrelCov) > 0); + constexpr bool fillFlowReso = eventHasQvector || eventHasQvectorCentr; bool isSelectedBDT = false; fNPairPerEvent = 0; @@ -1985,6 +1987,8 @@ struct AnalysisSameEventPairing { fNPairPerEvent++; + VarManager::fgValues[VarManager::kAmbi1] = -999.; + VarManager::fgValues[VarManager::kAmbi2] = -999.; if (t1.reducedeventId() != event.globalIndex()) { VarManager::fgValues[VarManager::kAmbi1] = 1.; } From a9596f377209d70bf332270aacef15cf84ac9456 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Wed, 17 Jun 2026 02:18:07 +0800 Subject: [PATCH 4/8] fix bug --- PWGDQ/Tasks/tableReader_withAssoc.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index ee8cdd59a55..191e2a336f1 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -2659,10 +2659,12 @@ struct AnalysisSameEventPairing { template void runSameSideMixing(TEvents& events, TAssocs const& assocs, TTracks const& tracks, Preslice& preSlice) { + if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { LOG(info) << "Flow resolution objects not set, flow will not be filled for mixed events"; if (events.size() > 0) { initParamsFromCCDB(events.begin().timestamp(), events.begin().runNumber(), false); } + } events.bindExternalIndices(&assocs); int mixingDepth = fConfigMixingDepth.value; fAmbiguousPairs.clear(); From b7f9128fe3b9f86d908a5415c7c6b95ae35f2736 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 16 Jun 2026 18:53:13 +0000 Subject: [PATCH 5/8] Please consider the following formatting changes --- PWGDQ/Core/CutsLibrary.cxx | 1 - PWGDQ/Core/VarManager.h | 44 ++++++++++---------- PWGDQ/Tasks/tableReader_withAssoc.cxx | 4 +- PWGDQ/Tasks/tableReader_withAssoc_direct.cxx | 2 +- 4 files changed, 26 insertions(+), 25 deletions(-) diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index d5eaa426ef1..8c29895c976 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -5084,7 +5084,6 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } - if (!nameStr.compare("NoelectronStandardQualityTPCOnly")) { cut->AddCut(VarManager::kTPCchi2, 0.0, 4.0, true, VarManager::kTPCncls, 70, 161.); return cut; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 63590cd8a05..f25d4efcc60 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -60,12 +60,12 @@ #include #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include +#include #include #include #include #include #include -#include #include #include @@ -1549,8 +1549,9 @@ class VarManager : public TObject fgEOR = eor; } - static void SetEventQVectorCorrection(TH3F* qvec) { - fgObjQvec = qvec; + static void SetEventQVectorCorrection(TH3F* qvec) + { + fgObjQvec = qvec; fgApplyQVectorCorrection = true; } @@ -6080,17 +6081,17 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) if (useCoherentJpsiA2) { // remove daughter from TPC Q-vector // TODO: remove based on track cut in qVectorTable - float Q2X0A = values[kQ2X0A]*values[kMultA]; - float Q2Y0A = values[kQ2Y0A]*values[kMultA]; + float Q2X0A = values[kQ2X0A] * values[kMultA]; + float Q2Y0A = values[kQ2Y0A] * values[kMultA]; float nNorm = values[kMultA]; - int centBin = static_cast(values[kCentFT0C]) + 1; // centrality bin - int detIdx = 6; // TPC all + int centBin = static_cast(values[kCentFT0C]) + 1; // centrality bin + int detIdx = 6; // TPC all EventPlaneHelper epHelper; // checkTrack(t1); if (isSelectedinTPC(t1) && values[kAmbi1] > 0) { - float qx1 = t1.pt()*TMath::Cos(2. * t1.phi()); - float qy1 = t1.pt()*TMath::Sin(2. * t1.phi()); + float qx1 = t1.pt() * TMath::Cos(2. * t1.phi()); + float qy1 = t1.pt() * TMath::Sin(2. * t1.phi()); if (fgApplyQVectorCorrection) { epHelper.DoRecenter(qx1, qy1, fgObjQvec->GetBinContent(centBin, 1, detIdx + 1), fgObjQvec->GetBinContent(centBin, 2, detIdx + 1)); epHelper.DoTwist(qx1, qy1, fgObjQvec->GetBinContent(centBin, 3, detIdx + 1), fgObjQvec->GetBinContent(centBin, 4, detIdx + 1)); @@ -6102,8 +6103,8 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) } // checkTrack(t2); if (isSelectedinTPC(t2) && values[kAmbi2] > 0) { - float qx2 = t2.pt()*TMath::Cos(2. * t2.phi()); - float qy2 = t2.pt()*TMath::Sin(2. * t2.phi()); + float qx2 = t2.pt() * TMath::Cos(2. * t2.phi()); + float qy2 = t2.pt() * TMath::Sin(2. * t2.phi()); if (fgApplyQVectorCorrection) { epHelper.DoRecenter(qx2, qy2, fgObjQvec->GetBinContent(centBin, 1, detIdx + 1), fgObjQvec->GetBinContent(centBin, 2, detIdx + 1)); epHelper.DoTwist(qx2, qy2, fgObjQvec->GetBinContent(centBin, 3, detIdx + 1), fgObjQvec->GetBinContent(centBin, 4, detIdx + 1)); @@ -6119,9 +6120,9 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) values[kQ2Y0A] = -999.; return; } - - Q2X0A = nNorm > 0 ? Q2X0A/nNorm : NAN; - Q2Y0A = nNorm > 0 ? Q2Y0A/nNorm : NAN; + + Q2X0A = nNorm > 0 ? Q2X0A / nNorm : NAN; + Q2Y0A = nNorm > 0 ? Q2Y0A / nNorm : NAN; float Psi2A = getEventPlane(2, Q2X0A, Q2Y0A); values[kPsi2A] = Psi2A; @@ -6154,11 +6155,11 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) float sinPsi2CPP = b_FT0C_RF.Dot(yAxis_RF.Cross(daughterVec_RF)); float Psi2CPP = sinPsi2CPP > 0 ? TMath::ACos(cosPsi2CPP) : -1. * TMath::ACos(cosPsi2CPP); values[kDeltaPhiPP_TPC] = phi_PP > Psi2APP ? phi_PP - Psi2APP : Psi2APP - phi_PP; - values[kDeltaPhiPP_TPC] = values[kDeltaPhiPP_TPC] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_TPC] : values[kDeltaPhiPP_TPC]; + values[kDeltaPhiPP_TPC] = values[kDeltaPhiPP_TPC] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_TPC] : values[kDeltaPhiPP_TPC]; values[kDeltaPhiPP_FT0A] = phi_PP > Psi2BPP ? phi_PP - Psi2BPP : Psi2BPP - phi_PP; - values[kDeltaPhiPP_FT0A] = values[kDeltaPhiPP_FT0A] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_FT0A] : values[kDeltaPhiPP_FT0A]; + values[kDeltaPhiPP_FT0A] = values[kDeltaPhiPP_FT0A] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_FT0A] : values[kDeltaPhiPP_FT0A]; values[kDeltaPhiPP_FT0C] = phi_PP > Psi2CPP ? phi_PP - Psi2CPP : Psi2CPP - phi_PP; - values[kDeltaPhiPP_FT0C] = values[kDeltaPhiPP_FT0C] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_FT0C] : values[kDeltaPhiPP_FT0C]; + values[kDeltaPhiPP_FT0C] = values[kDeltaPhiPP_FT0C] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_FT0C] : values[kDeltaPhiPP_FT0C]; // fold delta phi to [0, pi/2] values[kDeltaPhiPP_TPC] = values[kDeltaPhiPP_TPC] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiPP_TPC] : values[kDeltaPhiPP_TPC]; values[kDeltaPhiPP_FT0A] = values[kDeltaPhiPP_FT0A] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiPP_FT0A] : values[kDeltaPhiPP_FT0A]; @@ -6177,11 +6178,11 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) // reaction plane float phi = v_daughter.Phi() > TMath::Pi() ? 2. * TMath::Pi() - v_daughter.Phi() : v_daughter.Phi(); values[kDeltaPhiRP_TPC] = phi > Psi2A ? phi - Psi2A : Psi2A - phi; - values[kDeltaPhiRP_TPC] = values[kDeltaPhiRP_TPC] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_TPC] : values[kDeltaPhiRP_TPC]; + values[kDeltaPhiRP_TPC] = values[kDeltaPhiRP_TPC] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_TPC] : values[kDeltaPhiRP_TPC]; values[kDeltaPhiRP_FT0A] = phi > Psi2B ? phi - Psi2B : Psi2B - phi; - values[kDeltaPhiRP_FT0A] = values[kDeltaPhiRP_FT0A] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_FT0A] : values[kDeltaPhiRP_FT0A]; + values[kDeltaPhiRP_FT0A] = values[kDeltaPhiRP_FT0A] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_FT0A] : values[kDeltaPhiRP_FT0A]; values[kDeltaPhiRP_FT0C] = phi > Psi2C ? phi - Psi2C : Psi2C - phi; - values[kDeltaPhiRP_FT0C] = values[kDeltaPhiRP_FT0C] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_FT0C] : values[kDeltaPhiRP_FT0C]; + values[kDeltaPhiRP_FT0C] = values[kDeltaPhiRP_FT0C] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_FT0C] : values[kDeltaPhiRP_FT0C]; // fold delta phi to [0, pi/2] values[kDeltaPhiRP_TPC] = values[kDeltaPhiRP_TPC] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_TPC] : values[kDeltaPhiRP_TPC]; values[kDeltaPhiRP_FT0A] = values[kDeltaPhiRP_FT0A] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_FT0A] : values[kDeltaPhiRP_FT0A]; @@ -7045,7 +7046,8 @@ float VarManager::LorentzTransformJpsihadroncosChi(TString Option, T1 const& v1, } template -bool VarManager::isSelectedinTPC(const T& t) { +bool VarManager::isSelectedinTPC(const T& t) +{ bool trackSelected = false; if constexpr (requires { t.hasTPC(); }) { trackSelected = true; diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index 191e2a336f1..443b5a5ac12 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -1830,7 +1830,7 @@ struct AnalysisSameEventPairing { LOGF(fatal, "Flow resolution histograms not available in CCDB at timestamp=%llu", timestamp); } } else if (fConfigOptions.useLocalFlow) { - TString pathFlow = fConfigCCDB.flowPathLocal.value; + TString pathFlow = fConfigCCDB.flowPathLocal.value; TFile* fileFlow = TFile::Open(pathFlow.Data(), "READ"); if (fileFlow == nullptr || fileFlow->IsZombie()) { LOGF(fatal, "Flow resolution file %s cannot be opened", pathFlow.Data()); @@ -1955,7 +1955,7 @@ struct AnalysisSameEventPairing { if (groupedAssocs.size() == 0) { continue; } - + if (fillFlowReso) { if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { LOGF(fatal, "Flow resolution histograms are not available, cannot fill flow variables!"); diff --git a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx index 104f213195d..4df52f24343 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx @@ -1604,7 +1604,7 @@ struct AnalysisSameEventPairing { } } else if (fConfigOptions.useLocalFlow) { // LOGP(info, "Developing"); - TString pathFlow = fConfigCCDB.flowPathLocal.value; + TString pathFlow = fConfigCCDB.flowPathLocal.value; TFile* fileFlow = TFile::Open(pathFlow.Data(), "READ"); if (fileFlow == nullptr || fileFlow->IsZombie()) { LOGF(fatal, "Flow resolution file %s cannot be opened", pathFlow.Data()); From 6c76cd4cf3c7408d65fedeee2243e031efdb7a39 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Fri, 19 Jun 2026 08:54:50 +0800 Subject: [PATCH 6/8] remove local flow --- PWGDQ/Tasks/tableReader_withAssoc.cxx | 18 +++--------------- PWGDQ/Tasks/tableReader_withAssoc_direct.cxx | 17 ++--------------- 2 files changed, 5 insertions(+), 30 deletions(-) diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index 443b5a5ac12..a71646ab98f 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -1360,8 +1360,7 @@ struct AnalysisSameEventPairing { Configurable useRemoteCollisionInfo{"cfgUseRemoteCollisionInfo", false, "Use remote collision information from CCDB"}; 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 useRemoteFlow{"cfgUseRemoteFlow", false, "Use remote flow information from CCDB"}; - Configurable useLocalFlow{"cfgUseLocalFlow", false, "Use flow information from local cache"}; + Configurable useFlowReso{"cfgUseFlowReso", false, "Use remote flow information from CCDB"}; Configurable useQvecCalib{"cfgUseQvecCalib", false, "Use flow correction factors for Q-vector recalibration when removing the daughter"}; Configurable useCorrectionForRun{"cfgUseCorrectionForRun", false, "Apply run-by-run correction factors to the flow vectors"}; } fConfigOptions; @@ -1820,7 +1819,7 @@ struct AnalysisSameEventPairing { VarManager::SetEfficiencyObject(fConfigOptions.efficiencyType.value, effList->FindObject("efficiency")); } - if (fConfigOptions.useRemoteFlow) { + if (fConfigOptions.useFlowReso) { TString PathFlow = fConfigCCDB.flowPath.value; TString ccdbPathFlowSP = Form("%s/ScalarProduct", PathFlow.Data()); TString ccdbPathFlowEP = Form("%s/EventPlane", PathFlow.Data()); @@ -1829,17 +1828,6 @@ struct AnalysisSameEventPairing { if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { LOGF(fatal, "Flow resolution histograms not available in CCDB at timestamp=%llu", timestamp); } - } else if (fConfigOptions.useLocalFlow) { - TString pathFlow = fConfigCCDB.flowPathLocal.value; - TFile* fileFlow = TFile::Open(pathFlow.Data(), "READ"); - if (fileFlow == nullptr || fileFlow->IsZombie()) { - LOGF(fatal, "Flow resolution file %s cannot be opened", pathFlow.Data()); - } - fileFlow->GetObject("ScalarProduct", ResoFlowSP); - fileFlow->GetObject("EventPlane", ResoFlowEP); - if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { - LOGF(fatal, "Flow resolution histograms not available in file %s", pathFlow.Data()); - } } if (fConfigOptions.useQvecCalib) { @@ -2680,7 +2668,7 @@ struct AnalysisSameEventPairing { fNPairPerEvent = 0; VarManager::FillTwoMixEvents(event1, event2, assocs1, assocs2); - if (fConfigOptions.useRemoteFlow || fConfigOptions.useLocalFlow) { + if (fConfigOptions.useFlowReso) { VarManager::FillTwoMixEventsFlowResoFactor(ResoFlowSP, ResoFlowEP); } runMixedPairing(assocs1, assocs2, tracks, tracks); diff --git a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx index 4df52f24343..ae5141a2dc3 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx @@ -1325,8 +1325,7 @@ struct AnalysisSameEventPairing { Configurable fConfigMiniTree{"cfgMiniTree", false, "Produce a single flat table with minimal information for analysis"}; Configurable fConfigMiniTreeMinMass{"cfgMiniTreeMinMass", 2, "Min. mass cut for minitree"}; Configurable fConfigMiniTreeMaxMass{"cfgMiniTreeMaxMass", 5, "Max. mass cut for minitree"}; - Configurable useRemoteFlow{"cfgUseRemoteFlow", false, "Use remote flow information from CCDB"}; - Configurable useLocalFlow{"cfgUseLocalFlow", false, "Use flow information from local cache"}; + Configurable useFlowReso{"cfgUseFlowReso", false, "Use remote flow information from CCDB"}; Configurable useQvecCalib{"cfgUseQvecCalib", false, "Use flow correction factors for Q-vector recalibration when removing the daughter"}; Configurable useCorrectionForRun{"cfgUseCorrectionForRun", false, "Apply run-by-run correction factors to the flow vectors"}; } fConfigOptions; @@ -1593,7 +1592,7 @@ struct AnalysisSameEventPairing { } } - if (fConfigOptions.useRemoteFlow) { + if (fConfigOptions.useFlowReso) { TString PathFlow = fConfigCCDB.flowPath.value; TString ccdbPathFlowSP = Form("%s/ScalarProduct", PathFlow.Data()); TString ccdbPathFlowEP = Form("%s/EventPlane", PathFlow.Data()); @@ -1602,18 +1601,6 @@ struct AnalysisSameEventPairing { if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { LOGF(fatal, "Flow resolution histograms not available in CCDB at timestamp=%llu", timestamp); } - } else if (fConfigOptions.useLocalFlow) { - // LOGP(info, "Developing"); - TString pathFlow = fConfigCCDB.flowPathLocal.value; - TFile* fileFlow = TFile::Open(pathFlow.Data(), "READ"); - if (fileFlow == nullptr || fileFlow->IsZombie()) { - LOGF(fatal, "Flow resolution file %s cannot be opened", pathFlow.Data()); - } - fileFlow->GetObject("ScalarProduct", ResoFlowSP); - fileFlow->GetObject("EventPlane", ResoFlowEP); - if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { - LOGF(fatal, "Flow resolution histograms not available in file %s", pathFlow.Data()); - } } if (fConfigOptions.useQvecCalib) { From c6d03aa562ad75f3284fd69476ceea55170ca3a0 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sat, 20 Jun 2026 06:33:32 +0800 Subject: [PATCH 7/8] discard code for Q vector correction in DQ framework --- PWGDQ/Core/VarManager.cxx | 2 -- PWGDQ/Core/VarManager.h | 21 ------------ PWGDQ/Tasks/tableReader_withAssoc.cxx | 20 ------------ PWGDQ/Tasks/tableReader_withAssoc_direct.cxx | 34 ++------------------ 4 files changed, 2 insertions(+), 75 deletions(-) diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index e94e703748c..89b9a8401db 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -76,8 +76,6 @@ 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 -TH3F* VarManager::fgObjQvec = nullptr; -bool VarManager::fgApplyQVectorCorrection = false; //__________________________________________________________________ VarManager::VarManager() : TObject() diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index f25d4efcc60..ea77d1b076c 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1549,12 +1549,6 @@ class VarManager : public TObject fgEOR = eor; } - static void SetEventQVectorCorrection(TH3F* qvec) - { - fgObjQvec = qvec; - fgApplyQVectorCorrection = true; - } - public: VarManager(); ~VarManager() override; @@ -1615,8 +1609,6 @@ class VarManager : public TObject static int fgEfficiencyType; // type of efficiency correction to apply static TObject* fgEfficiencyHist; // histogram for efficiency correction - static TH3F* fgObjQvec; - static bool fgApplyQVectorCorrection; VarManager& operator=(const VarManager& c); VarManager(const VarManager& c); @@ -6085,18 +6077,10 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) float Q2Y0A = values[kQ2Y0A] * values[kMultA]; float nNorm = values[kMultA]; - int centBin = static_cast(values[kCentFT0C]) + 1; // centrality bin - int detIdx = 6; // TPC all - EventPlaneHelper epHelper; // checkTrack(t1); if (isSelectedinTPC(t1) && values[kAmbi1] > 0) { float qx1 = t1.pt() * TMath::Cos(2. * t1.phi()); float qy1 = t1.pt() * TMath::Sin(2. * t1.phi()); - if (fgApplyQVectorCorrection) { - epHelper.DoRecenter(qx1, qy1, fgObjQvec->GetBinContent(centBin, 1, detIdx + 1), fgObjQvec->GetBinContent(centBin, 2, detIdx + 1)); - epHelper.DoTwist(qx1, qy1, fgObjQvec->GetBinContent(centBin, 3, detIdx + 1), fgObjQvec->GetBinContent(centBin, 4, detIdx + 1)); - epHelper.DoRescale(qx1, qy1, fgObjQvec->GetBinContent(centBin, 5, detIdx + 1), fgObjQvec->GetBinContent(centBin, 6, detIdx + 1)); - } Q2X0A = Q2X0A - qx1; Q2Y0A = Q2Y0A - qy1; nNorm = nNorm - 1.; @@ -6105,11 +6089,6 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) if (isSelectedinTPC(t2) && values[kAmbi2] > 0) { float qx2 = t2.pt() * TMath::Cos(2. * t2.phi()); float qy2 = t2.pt() * TMath::Sin(2. * t2.phi()); - if (fgApplyQVectorCorrection) { - epHelper.DoRecenter(qx2, qy2, fgObjQvec->GetBinContent(centBin, 1, detIdx + 1), fgObjQvec->GetBinContent(centBin, 2, detIdx + 1)); - epHelper.DoTwist(qx2, qy2, fgObjQvec->GetBinContent(centBin, 3, detIdx + 1), fgObjQvec->GetBinContent(centBin, 4, detIdx + 1)); - epHelper.DoRescale(qx2, qy2, fgObjQvec->GetBinContent(centBin, 5, detIdx + 1), fgObjQvec->GetBinContent(centBin, 6, detIdx + 1)); - } Q2X0A = Q2X0A - qx2; Q2Y0A = Q2Y0A - qy2; nNorm = nNorm - 1.; diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index a71646ab98f..f282e253d19 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -1305,7 +1305,6 @@ struct AnalysisSameEventPairing { o2::base::MatLayerCylSet* fLUT = nullptr; TH1D* ResoFlowSP = nullptr; TH1D* ResoFlowEP = nullptr; - TH3F* qvecObj = nullptr; int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. OutputObj fOutputList{"output"}; @@ -1340,8 +1339,6 @@ 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 flowPathLocal{"flowPathLocal", "/lustre/alice/users/ywang/calib/FlowReso.root", "Path to the flow resolution object in the local cache"}; - Configurable QvecCalibPath{"cfgQvecCalibPath", "Users/j/junlee/Qvector/Pass5/QvecCalib/v2", "Path to the q vector calibration object"}; } fConfigCCDB; struct : ConfigurableGroup { @@ -1361,8 +1358,6 @@ 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 useQvecCalib{"cfgUseQvecCalib", false, "Use flow correction factors for Q-vector recalibration when removing the daughter"}; - Configurable useCorrectionForRun{"cfgUseCorrectionForRun", false, "Apply run-by-run correction factors to the flow vectors"}; } fConfigOptions; struct : ConfigurableGroup { Configurable applyBDT{"applyBDT", false, "Flag to apply ML selections"}; @@ -1829,15 +1824,6 @@ struct AnalysisSameEventPairing { LOGF(fatal, "Flow resolution histograms not available in CCDB at timestamp=%llu", timestamp); } } - - if (fConfigOptions.useQvecCalib) { - TString pathQvecCalib = fConfigCCDB.QvecCalibPath.value; - if (fConfigOptions.useCorrectionForRun) { - qvecObj = fCCDB->getForRun(pathQvecCalib.Data(), runNumber); - } else { - qvecObj = fCCDB->getForTimeStamp(pathQvecCalib.Data(), timestamp); - } - } } // Template function to run same event pairing (barrel-barrel, muon-muon, barrel-muon) @@ -1949,12 +1935,6 @@ struct AnalysisSameEventPairing { LOGF(fatal, "Flow resolution histograms are not available, cannot fill flow variables!"); } VarManager::FillEventFlowResoFactor(ResoFlowSP, ResoFlowEP); - if (fConfigOptions.useQvecCalib) { - if (qvecObj == nullptr) { - LOGF(fatal, "Q-vector calibration object is not available, cannot fill flow variables!"); - } - VarManager::SetEventQVectorCorrection(qvecObj); - } } bool isFirst = true; diff --git a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx index ae5141a2dc3..cca9d775d3f 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc_direct.cxx @@ -1289,8 +1289,6 @@ struct AnalysisSameEventPairing { o2::base::MatLayerCylSet* fLUT = nullptr; TH1D* ResoFlowSP = nullptr; TH1D* ResoFlowEP = nullptr; - TH3F* qvecObj = nullptr; - TString pathQvecCalib; int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. OutputObj fOutputList{"output"}; @@ -1326,8 +1324,6 @@ struct AnalysisSameEventPairing { Configurable fConfigMiniTreeMinMass{"cfgMiniTreeMinMass", 2, "Min. mass cut for minitree"}; Configurable fConfigMiniTreeMaxMass{"cfgMiniTreeMaxMass", 5, "Max. mass cut for minitree"}; Configurable useFlowReso{"cfgUseFlowReso", false, "Use remote flow information from CCDB"}; - Configurable useQvecCalib{"cfgUseQvecCalib", false, "Use flow correction factors for Q-vector recalibration when removing the daughter"}; - Configurable useCorrectionForRun{"cfgUseCorrectionForRun", false, "Apply run-by-run correction factors to the flow vectors"}; } fConfigOptions; struct : ConfigurableGroup { @@ -1336,8 +1332,6 @@ struct AnalysisSameEventPairing { Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; Configurable flowPath{"flowPath", "Users/y/yiping/FlowResolution", "Path to the flow resolution object"}; - Configurable flowPathLocal{"flowPathLocal", "/lustre/alice/users/ywang/calib/FlowReso.root", "Path to the flow resolution object in the local cache"}; - Configurable QvecCalibPath{"cfgQvecCalibPath", "Users/j/junlee/Qvector/Pass5/QvecShift/v2", "Path to the q vector calibration object"}; } fConfigCCDB; Service fCCDB; @@ -1548,15 +1542,9 @@ struct AnalysisSameEventPairing { dqhistograms::AddHistogramsFromJSON(fHistMan, fConfigOptions.fConfigAddJSONHistograms.value.c_str()); // ad-hoc histograms via JSON VarManager::SetUseVars(fHistMan->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill fOutputList.setObject(fHistMan->GetMainHistogramList()); - - if (fConfigOptions.useQvecCalib) { - string tmppathQvecCalib; - getTaskOptionValue(context, "q-vectors-table", "cfgQvecCalibPath", tmppathQvecCalib, false); - pathQvecCalib = tmppathQvecCalib; - } } - void initParamsFromCCDB(uint64_t timestamp, int runnumber, bool withTwoProngFitter = true) + void initParamsFromCCDB(uint64_t timestamp, bool withTwoProngFitter = true) { if (fConfigOptions.useRemoteField.value) { o2::parameters::GRPMagField* grpmag = fCCDB->getForTimeStamp(fConfigCCDB.grpMagPath, timestamp); @@ -1602,18 +1590,6 @@ struct AnalysisSameEventPairing { LOGF(fatal, "Flow resolution histograms not available in CCDB at timestamp=%llu", timestamp); } } - - if (fConfigOptions.useQvecCalib) { - TString ccdbPathQvecCalib = Form("%s/v2", pathQvecCalib.Data()); - if (fConfigOptions.useCorrectionForRun) { - qvecObj = fCCDB->getForRun(ccdbPathQvecCalib.Data(), runnumber); // get the object with run dependence if correction for run is needed - } else { - qvecObj = fCCDB->getForTimeStamp(ccdbPathQvecCalib.Data(), timestamp); // get the object without time dependence if no correction for run is needed - } - if (qvecObj == nullptr) { - LOGF(fatal, "Q-vector calibration object not available in CCDB at timestamp=%llu", timestamp); - } - } } template @@ -1676,7 +1652,7 @@ struct AnalysisSameEventPairing { return; } if (fCurrentRun != bcs.begin().runNumber()) { - initParamsFromCCDB(bcs.begin().timestamp(), bcs.begin().runNumber(), TTwoProngFitter); + initParamsFromCCDB(bcs.begin().timestamp(), TTwoProngFitter); fCurrentRun = bcs.begin().runNumber(); } @@ -1737,12 +1713,6 @@ struct AnalysisSameEventPairing { LOGF(fatal, "Flow resolution histograms are not available, cannot fill flow variables!"); } VarManager::FillEventFlowResoFactor(ResoFlowSP, ResoFlowEP); - if (fConfigOptions.useQvecCalib) { - if (qvecObj == nullptr) { - LOGF(fatal, "Q-vector calibration object is not available, cannot fill flow variables!"); - } - VarManager::SetEventQVectorCorrection(qvecObj); - } } for (auto& [a1, a2] : o2::soa::combinations(groupedAssocs, groupedAssocs)) { From 36adc3de0475f3d7d6d1320fb5e611904ccad605 Mon Sep 17 00:00:00 2001 From: ypwangg Date: Sat, 20 Jun 2026 09:43:37 +0800 Subject: [PATCH 8/8] removed hard code for q vector track selection --- PWGDQ/Core/CutsLibrary.cxx | 28 ++++++++++++++ PWGDQ/Core/VarManager.h | 55 ++++----------------------- PWGDQ/Tasks/tableReader_withAssoc.cxx | 22 ++++++++--- 3 files changed, 51 insertions(+), 54 deletions(-) diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index bfd9f405798..a285fde2fca 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -1593,6 +1593,34 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) cut->AddCut(GetAnalysisCut("electronPID1shiftDown")); return cut; } + // ------------------------------------------------------------------------------------------------- + // + // Q vector contributor cut + // + if (!nameStr.compare("selTPCCentral")) { + AnalysisCut* kineCut = new AnalysisCut("kineCut", "kine cut"); + kineCut->AddCut(VarManager::kEta, -0.8, 0.8); + kineCut->AddCut(VarManager::kPt, 0.15, 5); + + AnalysisCut* qualityCuts = new AnalysisCut("qualityCuts", "quality cuts"); + qualityCuts->AddCut(VarManager::kTPCchi2, 0., 4.); + qualityCuts->AddCut(VarManager::kTPCnCRoverFindCls, 0.8, 1.); + qualityCuts->AddCut(VarManager::kIsITSibAny, 0.5, 1.5); + qualityCuts->AddCut(VarManager::kITSchi2, 0., 36.); + + AnalysisCut* dcaCuts = new AnalysisCut("dcaCuts", "DCA cuts"); + std::shared_ptr f1dcaxyHigh = std::make_shared("f1dcaxy", "[0]+[1]/pow(x,[2])", 0., 10.); + f1dcaxyHigh->SetParameters(0.0105, 0.035, 1.1); + std::shared_ptr f1dcaxyLow = std::make_shared("f1dcaxy_low", "[0]+[1]/pow(x,[2])", 0., 10.); + f1dcaxyLow->SetParameters(-0.0105, -0.035, 1.1); + dcaCuts->AddCut(VarManager::kTrackDCAxy, f1dcaxyLow, f1dcaxyHigh); + dcaCuts->AddCut(VarManager::kTrackDCAz, -2., 2.); + + cut->AddCut(kineCut); + cut->AddCut(qualityCuts); + cut->AddCut(dcaCuts); + } + // ------------------------------------------------------------------------------------------------- // // LMee cuts diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index ea77d1b076c..adadf97125b 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -823,8 +823,8 @@ class VarManager : public TObject kCos2DeltaPhiPP_FT0C, kNullA2, kInfA2, - kAmbi1, - kAmbi2, + kSel1, + kSel2, kDeltaPhiPair, kOpeningAngle, kQuadDCAabsXY, @@ -1592,8 +1592,6 @@ class VarManager : public TObject static float calculatePhiV(const T1& t1, const T2& t2); template static float LorentzTransformJpsihadroncosChi(TString Option, const T1& v1, const T2& v2); - template - static bool isSelectedinTPC(const T& t); static o2::vertexing::DCAFitterN<2> fgFitterTwoProngBarrel; static o2::vertexing::DCAFitterN<3> fgFitterThreeProngBarrel; @@ -3106,6 +3104,9 @@ void VarManager::FillTrack(T const& track, float* values) values[kITSncls] += ((track.itsClusterMap() & (1 << i)) ? 1 : 0); } } + if (fgUsedVars[kTPCnCRoverFindCls]) { + values[kTPCnCRoverFindCls] = values[kTPCnclsCR] / values[kTPCncls]; + } values[kTrackDCAxy] = track.dcaXY(); values[kTrackDCAz] = track.dcaZ(); if constexpr ((fillMap & ReducedTrackBarrelCov) > 0) { @@ -6078,7 +6079,7 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) float nNorm = values[kMultA]; // checkTrack(t1); - if (isSelectedinTPC(t1) && values[kAmbi1] > 0) { + if (values[kSel1] > 0) { float qx1 = t1.pt() * TMath::Cos(2. * t1.phi()); float qy1 = t1.pt() * TMath::Sin(2. * t1.phi()); Q2X0A = Q2X0A - qx1; @@ -6086,7 +6087,7 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) nNorm = nNorm - 1.; } // checkTrack(t2); - if (isSelectedinTPC(t2) && values[kAmbi2] > 0) { + if (values[kSel2] > 0) { float qx2 = t2.pt() * TMath::Cos(2. * t2.phi()); float qy2 = t2.pt() * TMath::Sin(2. * t2.phi()); Q2X0A = Q2X0A - qx2; @@ -7024,48 +7025,6 @@ float VarManager::LorentzTransformJpsihadroncosChi(TString Option, T1 const& v1, return value; } -template -bool VarManager::isSelectedinTPC(const T& t) -{ - bool trackSelected = false; - if constexpr (requires { t.hasTPC(); }) { - trackSelected = true; - if (!t.hasTPC()) { - trackSelected = false; - return trackSelected; - } - if (std::abs(t.eta()) > 0.8) { // eta cut used in q vector table. Todo: read from config - trackSelected = false; - return trackSelected; - } - if (t.pt() < 0.15 || t.pt() > 5.) { // pt cut used in q vector table. Todo: read from config - trackSelected = false; - return trackSelected; - } - if (t.tpcNClsCrossedRows() / t.tpcNClsFound() < 0.8 || t.tpcChi2NCl() > 4.) { - trackSelected = false; - return trackSelected; - } - if (!((t.itsClusterMap() & (1 << uint8_t(0))) > 0 || (t.itsClusterMap() & (1 << uint8_t(1))) > 0 || (t.itsClusterMap() & (1 << uint8_t(2))) > 0)) { - trackSelected = false; - return trackSelected; - } - if (t.itsChi2NCl() > 36.) { - trackSelected = false; - return trackSelected; - } - if (t.dcaZ() > 2) { - trackSelected = false; - return trackSelected; - } - if (t.dcaXY() > 0.0105 + 0.0350 / std::pow(t.pt(), 1.1)) { - trackSelected = false; - return trackSelected; - } - } - return trackSelected; -} - template void VarManager::FillFIT(T1 const& bc, T2 const& bcs, T3 const& ft0s, T4 const& fv0as, T5 const& fdds, float* values) { diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index ac502bb036c..fe1c961212e 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -1313,6 +1313,7 @@ struct AnalysisSameEventPairing { Configurable track{"cfgTrackCuts", "jpsiO2MCdebugCuts2", "Comma separated list of barrel track cuts"}; Configurable muon{"cfgMuonCuts", "", "Comma separated list of muon cuts"}; Configurable pair{"cfgPairCuts", "", "Comma separated list of pair cuts"}; + Configurable qVector{"cfgQvectorTrackCut", "", "Track cut of Q-vector, enable this if you want to remove the auto-correlation in TPC"}; Configurable event{"cfgRemoveCollSplittingCandidates", false, "If true, remove collision splitting candidates as determined by the event selection task upstream"}; // TODO: Add pair cuts via JSON } fConfigCuts; @@ -1390,6 +1391,7 @@ struct AnalysisSameEventPairing { uint32_t fTrackFilterMask = 0; // mask for the track cuts required in this task to be applied on the barrel cuts produced upstream uint32_t fMuonFilterMask = 0; // mask for the muon cuts required in this task to be applied on the muon cuts produced upstream + uint32_t fQvectorFilterMask = 0; // mask for the track cuts required to be applied on the tracks used for the Q-vector calculation int fNCutsBarrel = 0; int fNCutsMuon = 0; int fNPairCuts = 0; @@ -1462,6 +1464,7 @@ struct AnalysisSameEventPairing { if (!muonCutsStr.IsNull()) { objArrayMuonCuts = muonCutsStr.Tokenize(","); } + TString qVectorCutStr = fConfigCuts.qVector.value; if (fConfigML.applyBDT) { // BDT cuts via JSON @@ -1529,6 +1532,9 @@ struct AnalysisSameEventPairing { fNCutsBarrel = objArray->GetEntries(); for (int icut = 0; icut < objArray->GetEntries(); ++icut) { TString tempStr = objArray->At(icut)->GetName(); + if (tempStr.CompareTo(qVectorCutStr) == 0) { + fQvectorFilterMask |= (static_cast(1) << icut); + } fTrackCuts.push_back(tempStr); if (objArrayTrackCuts->FindObject(tempStr.Data()) != nullptr) { fTrackFilterMask |= (static_cast(1) << icut); @@ -1966,13 +1972,17 @@ struct AnalysisSameEventPairing { fNPairPerEvent++; - VarManager::fgValues[VarManager::kAmbi1] = -999.; - VarManager::fgValues[VarManager::kAmbi2] = -999.; - if (t1.reducedeventId() != event.globalIndex()) { - VarManager::fgValues[VarManager::kAmbi1] = 1.; + VarManager::fgValues[VarManager::kSel1] = -999.; + VarManager::fgValues[VarManager::kSel2] = -999.; + if (t1.reducedeventId() == event.globalIndex()) { + if ((a1.isBarrelSelected_raw() & fQvectorFilterMask) > 0) { + VarManager::fgValues[VarManager::kSel1] = 1.; + } } - if (t2.reducedeventId() != event.globalIndex()) { - VarManager::fgValues[VarManager::kAmbi2] = 1.; + if (t2.reducedeventId() == event.globalIndex()) { + if ((a2.isBarrelSelected_raw() & fQvectorFilterMask) > 0) { + VarManager::fgValues[VarManager::kSel2] = 1.; + } } VarManager::FillPair(t1, t2); // compute quantities which depend on the associated collision, such as DCA