diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index 92cbaf7b7bc..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 @@ -3001,6 +3029,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 @@ -4141,6 +4180,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); @@ -4597,6 +4645,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); @@ -5080,6 +5134,12 @@ 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; @@ -5298,6 +5358,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..89b9a8401db 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -2484,6 +2484,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..adadf97125b 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -60,6 +60,7 @@ #include #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include +#include #include #include #include @@ -384,6 +385,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 +809,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, + kSel1, + kSel2, kDeltaPhiPair, kOpeningAngle, kQuadDCAabsXY, @@ -820,6 +840,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, @@ -2388,6 +2414,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(); @@ -2983,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) { @@ -5945,6 +6069,116 @@ 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]; + + // checkTrack(t1); + 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; + Q2Y0A = Q2Y0A - qy1; + nNorm = nNorm - 1.; + } + // checkTrack(t2); + 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; + 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 +7024,7 @@ float VarManager::LorentzTransformJpsihadroncosChi(TString Option, T1 const& v1, } return value; } + 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 0bf746ea95d..fe1c961212e 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -49,6 +49,7 @@ #include #include #include +#include #include #include #include @@ -240,7 +241,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; @@ -1302,6 +1303,8 @@ struct AnalysisSameEventPairing { Produces dileptonEventInfoList; o2::base::MatLayerCylSet* fLUT = nullptr; + TH1D* ResoFlowSP = nullptr; + TH1D* ResoFlowEP = nullptr; int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. OutputObj fOutputList{"output"}; @@ -1310,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; @@ -1335,6 +1339,7 @@ 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"}; } fConfigCCDB; struct : ConfigurableGroup { @@ -1353,6 +1358,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 useFlowReso{"cfgUseFlowReso", false, "Use remote flow information from CCDB"}; } fConfigOptions; struct : ConfigurableGroup { Configurable applyBDT{"applyBDT", false, "Flag to apply ML selections"}; @@ -1385,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; @@ -1457,6 +1464,7 @@ struct AnalysisSameEventPairing { if (!muonCutsStr.IsNull()) { objArrayMuonCuts = muonCutsStr.Tokenize(","); } + TString qVectorCutStr = fConfigCuts.qVector.value; if (fConfigML.applyBDT) { // BDT cuts via JSON @@ -1524,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); @@ -1807,6 +1818,17 @@ struct AnalysisSameEventPairing { auto effList = fCCDB->getForTimeStamp(fConfigCCDB.efficiencyPath.value, timestamp); VarManager::SetEfficiencyObject(fConfigOptions.efficiencyType.value, effList->FindObject("efficiency")); } + + if (fConfigOptions.useFlowReso) { + 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); + } + } } // Template function to run same event pairing (barrel-barrel, muon-muon, barrel-muon) @@ -1891,6 +1913,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; @@ -1912,6 +1935,13 @@ struct AnalysisSameEventPairing { continue; } + if (fillFlowReso) { + if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) { + LOGF(fatal, "Flow resolution histograms are not available, cannot fill flow variables!"); + } + VarManager::FillEventFlowResoFactor(ResoFlowSP, ResoFlowEP); + } + bool isFirst = true; for (auto& [a1, a2] : o2::soa::combinations(groupedAssocs, groupedAssocs)) { if constexpr (TPairType == VarManager::kDecayToEE) { @@ -1941,6 +1971,19 @@ struct AnalysisSameEventPairing { } fNPairPerEvent++; + + 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()) { + 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 if (fConfigOptions.propTrack) { @@ -2593,6 +2636,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(); @@ -2607,6 +2656,10 @@ struct AnalysisSameEventPairing { assocs2.bindExternalIndices(&events); fNPairPerEvent = 0; + VarManager::FillTwoMixEvents(event1, event2, assocs1, assocs2); + if (fConfigOptions.useFlowReso) { + 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..cca9d775d3f 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,8 @@ struct AnalysisSameEventPairing { Produces electronmuonList; o2::base::MatLayerCylSet* fLUT = nullptr; + TH1D* ResoFlowSP = nullptr; + TH1D* ResoFlowEP = nullptr; int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. OutputObj fOutputList{"output"}; @@ -1309,6 +1323,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 useFlowReso{"cfgUseFlowReso", false, "Use remote flow information from CCDB"}; } fConfigOptions; struct : ConfigurableGroup { @@ -1316,6 +1331,7 @@ 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"}; } fConfigCCDB; Service fCCDB; @@ -1358,7 +1374,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 @@ -1563,6 +1579,17 @@ 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.useFlowReso) { + 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); + } + } } template @@ -1662,8 +1689,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 +1708,13 @@ 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); + } + for (auto& [a1, a2] : o2::soa::combinations(groupedAssocs, groupedAssocs)) { if constexpr (TPairType == VarManager::kDecayToEE) { @@ -1755,7 +1790,7 @@ struct AnalysisSameEventPairing { } if constexpr (eventHasQvector) { - VarManager::FillPairVn(t1, t2); + VarManager::FillPairVn(t1, t2); } } // Fill normal histograms @@ -1947,10 +1982,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); };