diff --git a/PWGCF/Flow/Tasks/flowSP.cxx b/PWGCF/Flow/Tasks/flowSP.cxx index a8ffd936110..d3358c22496 100644 --- a/PWGCF/Flow/Tasks/flowSP.cxx +++ b/PWGCF/Flow/Tasks/flowSP.cxx @@ -47,7 +47,14 @@ using namespace o2::framework::expressions; #define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; struct FlowSP { - + // QA Plots + O2_DEFINE_CONFIGURABLE(cfgFillEventQA, bool, true, "Fill histograms for event QA"); + // Centrality Estimators -> standard is FT0C + O2_DEFINE_CONFIGURABLE(cfgFT0Cvariant1, bool, false, "Set centrality estimator to cfgFT0Cvariant1"); + O2_DEFINE_CONFIGURABLE(cfgFT0M, bool, false, "Set centrality estimator to cfgFT0M"); + O2_DEFINE_CONFIGURABLE(cfgFV0A, bool, false, "Set centrality estimator to cfgFV0A"); + O2_DEFINE_CONFIGURABLE(cfgNGlobal, bool, false, "Set centrality estimator to cfgNGlobal"); + // Standard selections O2_DEFINE_CONFIGURABLE(cfgDCAxy, float, 0.2, "Cut on DCA in the transverse direction (cm)"); O2_DEFINE_CONFIGURABLE(cfgDCAz, float, 2, "Cut on DCA in the longitudinal direction (cm)"); O2_DEFINE_CONFIGURABLE(cfgNcls, float, 70, "Cut on number of TPC clusters found"); @@ -56,33 +63,37 @@ struct FlowSP { O2_DEFINE_CONFIGURABLE(cfgEta, float, 0.8, "eta cut"); O2_DEFINE_CONFIGURABLE(cfgVtxZ, float, 10, "vertex cut (cm)"); O2_DEFINE_CONFIGURABLE(cfgMagField, float, 99999, "Configurable magnetic field;default CCDB will be queried"); - O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, true, "Bool to enable Additional Event Cut"); - O2_DEFINE_CONFIGURABLE(cfgUseAdditionalTrackCut, bool, true, "Bool to enable Additional Track Cut"); O2_DEFINE_CONFIGURABLE(cfgCentMin, float, 0, "Minimum cenrality for selected events"); O2_DEFINE_CONFIGURABLE(cfgCentMax, float, 90, "Maximum cenrality for selected events"); + // NUA and NUE weights O2_DEFINE_CONFIGURABLE(cfgFillWeights, bool, true, "Fill NUA weights"); O2_DEFINE_CONFIGURABLE(cfgFillWeightsPOS, bool, false, "Fill NUA weights only for positive charges"); O2_DEFINE_CONFIGURABLE(cfgFillWeightsNEG, bool, false, "Fill NUA weights only for negative charges"); O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "ccdb dir for NUA corrections"); O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "ccdb dir for NUE corrections"); + // Additional track Selections + O2_DEFINE_CONFIGURABLE(cfgUseAdditionalTrackCut, bool, true, "Bool to enable Additional Track Cut"); O2_DEFINE_CONFIGURABLE(cfgDoubleTrackFunction, bool, true, "Include track cut at low pt"); O2_DEFINE_CONFIGURABLE(cfgTrackCutSize, float, 0.06, "Spread of track cut"); + // Additional event selections + O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, true, "Bool to enable Additional Event Cut"); O2_DEFINE_CONFIGURABLE(cfgMaxOccupancy, int, 10000, "Maximum occupancy of selected events"); O2_DEFINE_CONFIGURABLE(cfgNoSameBunchPileupCut, bool, true, "kNoSameBunchPileupCut"); O2_DEFINE_CONFIGURABLE(cfgIsGoodZvtxFT0vsPV, bool, true, "kIsGoodZvtxFT0vsPV"); O2_DEFINE_CONFIGURABLE(cfgNoCollInTimeRangeStandard, bool, true, "kNoCollInTimeRangeStandard"); O2_DEFINE_CONFIGURABLE(cfgDoOccupancySel, bool, true, "Bool for event selection on detector occupancy"); - O2_DEFINE_CONFIGURABLE(cfgMultCut, bool, true, "Use additional evenr cut on mult correlations"); O2_DEFINE_CONFIGURABLE(cfgTVXinTRD, bool, false, "Use kTVXinTRD (reject TRD triggered events)"); O2_DEFINE_CONFIGURABLE(cfgIsVertexITSTPC, bool, true, "Selects collisions with at least one ITS-TPC track"); O2_DEFINE_CONFIGURABLE(cfgIsGoodITSLayersAll, bool, true, "Cut time intervals with dead ITS staves"); - O2_DEFINE_CONFIGURABLE(cfgCCDBdir, std::string, "Users/c/ckoster/ZDC/LHC23_zzh_pass4_small/meanQQ", "ccdb dir for average QQ values in 1% centrality bins"); - O2_DEFINE_CONFIGURABLE(cfgLoadAverageQQ, bool, true, "Load average values for QQ (in centrality bins)"); + // harmonics for v coefficients O2_DEFINE_CONFIGURABLE(cfgHarm, int, 1, "Flow harmonic n for ux and uy: (Cos(n*phi), Sin(n*phi))"); O2_DEFINE_CONFIGURABLE(cfgHarmMixed, int, 2, "Flow harmonic n for ux and uy in mixed harmonics (MH): (Cos(n*phi), Sin(n*phi))"); + // settings for CCDB data + O2_DEFINE_CONFIGURABLE(cfgLoadAverageQQ, bool, true, "Load average values for QQ (in centrality bins)"); + O2_DEFINE_CONFIGURABLE(cfgCCDBdir, std::string, "Users/c/ckoster/ZDC/LHC23_zzh_pass4_small/meanQQ", "ccdb dir for average QQ values in 1% centrality bins"); O2_DEFINE_CONFIGURABLE(cfgLoadSPPlaneRes, bool, false, "Load ZDC spectator plane resolution"); O2_DEFINE_CONFIGURABLE(cfgCCDBdir_SP, std::string, "Users/c/ckoster/ZDC/LHC23_zzh_pass4_small/SPPlaneRes", "ccdb dir for average event plane resolution in 1% centrality bins"); - + // axis ConfigurableAxis axisDCAz{"axisDCAz", {200, -.5, .5}, "DCA_{z} (cm)"}; ConfigurableAxis axisDCAxy{"axisDCAxy", {200, -.5, .5}, "DCA_{xy} (cm)"}; ConfigurableAxis axisPhiMod = {"axisPhiMod", {100, 0, constants::math::PI / 9}, "fmod(#varphi,#pi/9)"}; @@ -97,13 +108,14 @@ struct FlowSP { Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ; Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && nabs(aod::track::dcaXY) < cfgDCAxy&& nabs(aod::track::dcaZ) < cfgDCAz; - using UsedCollisions = soa::Filtered>; + using UsedCollisions = soa::Filtered>; using UsedTracks = soa::Filtered>; // Connect to ccdb Service ccdb; // from Generic Framework + // Adapted to hold weights for: inclusive, positive charged, negative charged struct Config { std::vector mEfficiency = {}; std::vector mAcceptance = {}; @@ -111,6 +123,7 @@ struct FlowSP { int lastRunNumber = 0; } cfg; + // define output objects OutputObj fWeights{GFWWeights("weights")}; OutputObj fWeightsPOS{GFWWeights("weights_positive")}; OutputObj fWeightsNEG{GFWWeights("weights_negative")}; @@ -147,6 +160,13 @@ struct FlowSP { kNegative }; + enum FillType { + kBefore, + kAfter + }; + + static constexpr std::string_view Charge[] = {"incl/", "pos/", "neg/"}; + void init(InitContext const&) { ccdb->setURL("http://alice-ccdb.cern.ch"); @@ -176,102 +196,120 @@ struct FlowSP { fWeightsNEG->Init(true, false); } - registry.add("hSPplaneA", "hSPplaneA", kTH1D, {axisPhiPlane}); - registry.add("hSPplaneC", "hSPplaneC", kTH1D, {axisPhiPlane}); - registry.add("hSPplaneFull", "hSPplaneFull", kTH1D, {axisPhiPlane}); - - registry.add("hCosPhiACosPhiC", "hCosPhiACosPhiC; Centrality(%); #LT Cos(#Psi^{A})Cos(#Psi^{C})#GT", kTProfile, {axisCent}); - registry.add("hSinPhiASinPhiC", "hSinPhiASinPhiC; Centrality(%); #LT Sin(#Psi^{A})Sin(#Psi^{C})#GT", kTProfile, {axisCent}); - registry.add("hSinPhiACosPhiC", "hSinPhiACosPhiC; Centrality(%); #LT Sin(#Psi^{A})Cos(#Psi^{C})#GT", kTProfile, {axisCent}); - registry.add("hCosPhiASinsPhiC", "hCosPhiASinsPhiC; Centrality(%); #LT Cos(#Psi^{A})Sin(#Psi^{C})#GT", kTProfile, {axisCent}); - registry.add("hFullEvPlaneRes", "hFullEvPlaneRes; Centrality(%); -#LT Cos(#Psi^{A} - #Psi^{C})#GT ", kTProfile, {axisCent}); - - registry.add("QA/after/hCent", "", {HistType::kTH1D, {axisCent}}); - registry.add("QA/after/pt_phi", "", {HistType::kTH2D, {axisPt, axisPhiMod}}); - registry.add("QA/after/hPt_inclusive", "", {HistType::kTH1D, {axisPt}}); - registry.add("QA/after/globalTracks_centT0C", "", {HistType::kTH2D, {axisCent, nchAxis}}); - registry.add("QA/after/PVTracks_centT0C", "", {HistType::kTH2D, {axisCent, multpvAxis}}); - registry.add("QA/after/globalTracks_PVTracks", "", {HistType::kTH2D, {multpvAxis, nchAxis}}); - registry.add("QA/after/globalTracks_multT0A", "", {HistType::kTH2D, {t0aAxis, nchAxis}}); - registry.add("QA/after/globalTracks_multV0A", "", {HistType::kTH2D, {t0aAxis, nchAxis}}); - registry.add("QA/after/multV0A_multT0A", "", {HistType::kTH2D, {t0aAxis, t0aAxis}}); - registry.add("QA/after/multT0C_centT0C", "", {HistType::kTH2D, {axisCent, t0cAxis}}); - - registry.add("QA/after/PsiA_vs_Cent", "", {HistType::kTH2D, {axisPhiPlane, axisCent}}); - registry.add("QA/after/PsiC_vs_Cent", "", {HistType::kTH2D, {axisPhiPlane, axisCent}}); - registry.add("QA/after/PsiFull_vs_Cent", "", {HistType::kTH2D, {axisPhiPlane, axisCent}}); - - registry.add("QA/after/PsiA_vs_Vx", "", {HistType::kTH2D, {axisPhiPlane, axisVx}}); - registry.add("QA/after/PsiC_vs_Vx", "", {HistType::kTH2D, {axisPhiPlane, axisVx}}); - registry.add("QA/after/PsiFull_vs_Vx", "", {HistType::kTH2D, {axisPhiPlane, axisVx}}); - - registry.add("QA/after/PsiA_vs_Vy", "", {HistType::kTH2D, {axisPhiPlane, axisVy}}); - registry.add("QA/after/PsiC_vs_Vy", "", {HistType::kTH2D, {axisPhiPlane, axisVy}}); - registry.add("QA/after/PsiFull_vs_Vy", "", {HistType::kTH2D, {axisPhiPlane, axisVy}}); - - registry.add("QA/after/PsiA_vs_Vz", "", {HistType::kTH2D, {axisPhiPlane, axisVz}}); - registry.add("QA/after/PsiC_vs_Vz", "", {HistType::kTH2D, {axisPhiPlane, axisVz}}); - registry.add("QA/after/PsiFull_vs_Vz", "", {HistType::kTH2D, {axisPhiPlane, axisVz}}); - - registry.addClone("QA/after/", "QA/before/"); - - // track properties per centrality and per eta, pt bin - registry.add("incl/vnAx_eta", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnAy_eta", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnCx_eta", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnCy_eta", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnC_eta", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnA_eta", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnA_eta_EP", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnC_eta_EP", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnFull_eta_EP", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnAxCxUx_eta_MH", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnAxCyUx_eta_MH", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnAxCyUy_eta_MH", "", kTProfile, {axisEtaVn}); - registry.add("incl/vnAyCxUy_eta_MH", "", kTProfile, {axisEtaVn}); - - registry.add("incl/vnAx_pt", "", kTProfile, {axisPt}); - registry.add("incl/vnAy_pt", "", kTProfile, {axisPt}); - registry.add("incl/vnCx_pt", "", kTProfile, {axisPt}); - registry.add("incl/vnCy_pt", "", kTProfile, {axisPt}); - registry.add("incl/vnC_pt", "", kTProfile, {axisPt}); - registry.add("incl/vnA_pt", "", kTProfile, {axisPt}); - registry.add("incl/vnC_pt_odd", "", kTProfile, {axisPt}); - registry.add("incl/vnA_pt_odd", "", kTProfile, {axisPt}); - registry.add("incl/vnA_pt_EP", "", kTProfile, {axisPt}); - registry.add("incl/vnC_pt_EP", "", kTProfile, {axisPt}); - registry.add("incl/vnFull_pt_EP", "", kTProfile, {axisPt}); - registry.add("incl/vnAxCxUx_pt_MH", "", kTProfile, {axisPt}); - registry.add("incl/vnAxCyUx_pt_MH", "", kTProfile, {axisPt}); - registry.add("incl/vnAxCyUy_pt_MH", "", kTProfile, {axisPt}); - registry.add("incl/vnAyCxUy_pt_MH", "", kTProfile, {axisPt}); - - registry.add("incl/vnC_cent_minEta", "", kTProfile, {axisCent}); - registry.add("incl/vnA_cent_minEta", "", kTProfile, {axisCent}); - registry.add("incl/vnC_cent_plusEta", "", kTProfile, {axisCent}); - registry.add("incl/vnA_cent_plusEta", "", kTProfile, {axisCent}); - registry.add("incl/vnA_cent_EP", "", kTProfile, {axisCent}); - registry.add("incl/vnC_cent_EP", "", kTProfile, {axisCent}); - registry.add("incl/vnFull_cent_EP", "", kTProfile, {axisCent}); - registry.add("incl/vnAxCxUx_cent_MH", "", kTProfile, {axisCent}); - registry.add("incl/vnAxCyUx_cent_MH", "", kTProfile, {axisCent}); - registry.add("incl/vnAxCyUy_cent_MH", "", kTProfile, {axisCent}); - registry.add("incl/vnAyCxUy_cent_MH", "", kTProfile, {axisCent}); - - // track QA for pos, neg, incl - registry.add("incl/QA/hPt", "", kTH1D, {axisPt}); - registry.add("incl/QA/hPhi", "", kTH1D, {axisPhi}); - registry.add("incl/QA/hEta", "", kTH1D, {axisEta}); - registry.add("incl/QA/hPhi_Eta_vz", "", kTH3D, {axisPhi, axisEta, axisVz}); - registry.add("incl/QA/hDCAxy", "", kTH1D, {axisDCAxy}); - registry.add("incl/QA/hDCAz", "", kTH1D, {axisDCAz}); - - registry.addClone("incl/", "pos/"); - registry.addClone("incl/", "neg/"); - - registry.add("qAqCX", "", kTProfile, {axisCent}); - registry.add("qAqCY", "", kTProfile, {axisCent}); - registry.add("qAqCXY", "", kTProfile, {axisCent}); + if ((doprocessData || doprocessMCReco)) { + if (cfgFillEventQA) { + registry.add("QA/after/hCent", "", {HistType::kTH1D, {axisCent}}); + registry.add("QA/after/pt_phi", "", {HistType::kTH2D, {axisPt, axisPhiMod}}); + registry.add("QA/after/hPt_inclusive", "", {HistType::kTH1D, {axisPt}}); + registry.add("QA/after/globalTracks_centT0C", "", {HistType::kTH2D, {axisCent, nchAxis}}); + registry.add("QA/after/PVTracks_centT0C", "", {HistType::kTH2D, {axisCent, multpvAxis}}); + registry.add("QA/after/globalTracks_PVTracks", "", {HistType::kTH2D, {multpvAxis, nchAxis}}); + registry.add("QA/after/globalTracks_multT0A", "", {HistType::kTH2D, {t0aAxis, nchAxis}}); + registry.add("QA/after/globalTracks_multV0A", "", {HistType::kTH2D, {t0aAxis, nchAxis}}); + registry.add("QA/after/multV0A_multT0A", "", {HistType::kTH2D, {t0aAxis, t0aAxis}}); + registry.add("QA/after/multT0C_centT0C", "", {HistType::kTH2D, {axisCent, t0cAxis}}); + } + + if (doprocessData) { + registry.add("hSPplaneA", "hSPplaneA", kTH1D, {axisPhiPlane}); + registry.add("hSPplaneC", "hSPplaneC", kTH1D, {axisPhiPlane}); + registry.add("hSPplaneFull", "hSPplaneFull", kTH1D, {axisPhiPlane}); + + registry.add("hCosPhiACosPhiC", "hCosPhiACosPhiC; Centrality(%); #LT Cos(#Psi^{A})Cos(#Psi^{C})#GT", kTProfile, {axisCent}); + registry.add("hSinPhiASinPhiC", "hSinPhiASinPhiC; Centrality(%); #LT Sin(#Psi^{A})Sin(#Psi^{C})#GT", kTProfile, {axisCent}); + registry.add("hSinPhiACosPhiC", "hSinPhiACosPhiC; Centrality(%); #LT Sin(#Psi^{A})Cos(#Psi^{C})#GT", kTProfile, {axisCent}); + registry.add("hCosPhiASinsPhiC", "hCosPhiASinsPhiC; Centrality(%); #LT Cos(#Psi^{A})Sin(#Psi^{C})#GT", kTProfile, {axisCent}); + registry.add("hFullEvPlaneRes", "hFullEvPlaneRes; Centrality(%); -#LT Cos(#Psi^{A} - #Psi^{C})#GT ", kTProfile, {axisCent}); + + // track properties per centrality and per eta, pt bin + registry.add("incl/vnAx_eta", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnAy_eta", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnCx_eta", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnCy_eta", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnC_eta", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnA_eta", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnA_eta_EP", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnC_eta_EP", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnFull_eta_EP", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnAxCxUx_eta_MH", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnAxCyUx_eta_MH", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnAxCyUy_eta_MH", "", kTProfile, {axisEtaVn}); + registry.add("incl/vnAyCxUy_eta_MH", "", kTProfile, {axisEtaVn}); + + registry.add("incl/vnAx_pt", "", kTProfile, {axisPt}); + registry.add("incl/vnAy_pt", "", kTProfile, {axisPt}); + registry.add("incl/vnCx_pt", "", kTProfile, {axisPt}); + registry.add("incl/vnCy_pt", "", kTProfile, {axisPt}); + registry.add("incl/vnC_pt", "", kTProfile, {axisPt}); + registry.add("incl/vnA_pt", "", kTProfile, {axisPt}); + registry.add("incl/vnC_pt_odd", "", kTProfile, {axisPt}); + registry.add("incl/vnA_pt_odd", "", kTProfile, {axisPt}); + registry.add("incl/vnA_pt_EP", "", kTProfile, {axisPt}); + registry.add("incl/vnC_pt_EP", "", kTProfile, {axisPt}); + registry.add("incl/vnFull_pt_EP", "", kTProfile, {axisPt}); + registry.add("incl/vnAxCxUx_pt_MH", "", kTProfile, {axisPt}); + registry.add("incl/vnAxCyUx_pt_MH", "", kTProfile, {axisPt}); + registry.add("incl/vnAxCyUy_pt_MH", "", kTProfile, {axisPt}); + registry.add("incl/vnAyCxUy_pt_MH", "", kTProfile, {axisPt}); + + registry.add("incl/vnC_cent_minEta", "", kTProfile, {axisCent}); + registry.add("incl/vnA_cent_minEta", "", kTProfile, {axisCent}); + registry.add("incl/vnC_cent_plusEta", "", kTProfile, {axisCent}); + registry.add("incl/vnA_cent_plusEta", "", kTProfile, {axisCent}); + registry.add("incl/vnA_cent_EP", "", kTProfile, {axisCent}); + registry.add("incl/vnC_cent_EP", "", kTProfile, {axisCent}); + registry.add("incl/vnFull_cent_EP", "", kTProfile, {axisCent}); + registry.add("incl/vnAxCxUx_cent_MH", "", kTProfile, {axisCent}); + registry.add("incl/vnAxCyUx_cent_MH", "", kTProfile, {axisCent}); + registry.add("incl/vnAxCyUy_cent_MH", "", kTProfile, {axisCent}); + registry.add("incl/vnAyCxUy_cent_MH", "", kTProfile, {axisCent}); + + registry.add("qAqCX", "", kTProfile, {axisCent}); + registry.add("qAqCY", "", kTProfile, {axisCent}); + registry.add("qAqCXY", "", kTProfile, {axisCent}); + registry.add("qAXqCY", "", kTProfile, {axisCent}); + registry.add("qAYqCX", "", kTProfile, {axisCent}); + + if (cfgFillEventQA) { + registry.add("QA/after/PsiA_vs_Cent", "", {HistType::kTH2D, {axisPhiPlane, axisCent}}); + registry.add("QA/after/PsiC_vs_Cent", "", {HistType::kTH2D, {axisPhiPlane, axisCent}}); + registry.add("QA/after/PsiFull_vs_Cent", "", {HistType::kTH2D, {axisPhiPlane, axisCent}}); + + registry.add("QA/after/PsiA_vs_Vx", "", {HistType::kTH2D, {axisPhiPlane, axisVx}}); + registry.add("QA/after/PsiC_vs_Vx", "", {HistType::kTH2D, {axisPhiPlane, axisVx}}); + registry.add("QA/after/PsiFull_vs_Vx", "", {HistType::kTH2D, {axisPhiPlane, axisVx}}); + + registry.add("QA/after/PsiA_vs_Vy", "", {HistType::kTH2D, {axisPhiPlane, axisVy}}); + registry.add("QA/after/PsiC_vs_Vy", "", {HistType::kTH2D, {axisPhiPlane, axisVy}}); + registry.add("QA/after/PsiFull_vs_Vy", "", {HistType::kTH2D, {axisPhiPlane, axisVy}}); + + registry.add("QA/after/PsiA_vs_Vz", "", {HistType::kTH2D, {axisPhiPlane, axisVz}}); + registry.add("QA/after/PsiC_vs_Vz", "", {HistType::kTH2D, {axisPhiPlane, axisVz}}); + registry.add("QA/after/PsiFull_vs_Vz", "", {HistType::kTH2D, {axisPhiPlane, axisVz}}); + + registry.add("QA/after/CentFT0C_vs_CentFT0Cvariant1", "", {HistType::kTH2D, {axisCent, axisCent}}); + registry.add("QA/after/CentFT0C_vs_CentFT0M", "", {HistType::kTH2D, {axisCent, axisCent}}); + registry.add("QA/after/CentFT0C_vs_CentFV0A", "", {HistType::kTH2D, {axisCent, axisCent}}); + // registry.add("QA/after/CentFT0C_vs_CentNGlobal", "", {HistType::kTH2D, {axisCent, axisCent}}; + } + } + registry.addClone("QA/after/", "QA/before/"); + // track QA for pos, neg, incl + registry.add("incl/QA/hPt", "", kTH1D, {axisPt}); + registry.add("incl/QA/hPhi", "", kTH1D, {axisPhi}); + registry.add("incl/QA/hPhiCorrected", "", kTH1D, {axisPhi}); + registry.add("incl/QA/hEta", "", kTH1D, {axisEta}); + registry.add("incl/QA/hPhi_Eta_vz", "", kTH3D, {axisPhi, axisEta, axisVz}); + registry.add("incl/QA/hDCAxy", "", kTH1D, {axisDCAxy}); + registry.add("incl/QA/hDCAz", "", kTH1D, {axisDCAz}); + + registry.addClone("incl/", "pos/"); + registry.addClone("incl/", "neg/"); + } else if (doprocessMCGen) { + registry.add("trackMCGen/before/pt_gen_incl", "", {HistType::kTH1D, {axisPt}}); + registry.add("trackMCGen/before/phi_eta_vtxZ_gen", "", {HistType::kTH3D, {axisPhi, axisEta, axisVz}}); + registry.addClone("trackMCGen/before/", "trackMCGen/after/"); + } registry.add("hEventCount", "Number of Event; Cut; #Events Passed Cut", {HistType::kTH1D, {{nEventSelections, 0, nEventSelections}}}); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_FilteredEvent + 1, "Filtered event"); @@ -476,8 +514,8 @@ struct FlowSP { return 1; } - template - bool trackSelected(TTrack track, const int& field) + template + bool trackSelected(TrackObject track, const int& field) { if (track.tpcNClsFound() < cfgNcls) @@ -503,60 +541,49 @@ struct FlowSP { return true; } - template - inline void fillEventQA(CollisionObject collision, TracksObject tracks, bool before) + template + inline void fillEventQA(CollisionObject collision, TracksObject tracks) { - if (before) { - registry.fill(HIST("QA/before/hCent"), collision.centFT0C()); - registry.fill(HIST("QA/before/globalTracks_centT0C"), collision.centFT0C(), tracks.size()); - registry.fill(HIST("QA/before/PVTracks_centT0C"), collision.centFT0C(), collision.multNTracksPV()); - registry.fill(HIST("QA/before/globalTracks_PVTracks"), collision.multNTracksPV(), tracks.size()); - registry.fill(HIST("QA/before/globalTracks_multT0A"), collision.multFT0A(), tracks.size()); - registry.fill(HIST("QA/before/globalTracks_multV0A"), collision.multFV0A(), tracks.size()); - registry.fill(HIST("QA/before/multV0A_multT0A"), collision.multFT0A(), collision.multFV0A()); - registry.fill(HIST("QA/before/multT0C_centT0C"), collision.centFT0C(), collision.multFT0C()); - } else { - registry.fill(HIST("QA/after/hCent"), collision.centFT0C()); - registry.fill(HIST("QA/after/globalTracks_centT0C"), collision.centFT0C(), tracks.size()); - registry.fill(HIST("QA/after/PVTracks_centT0C"), collision.centFT0C(), collision.multNTracksPV()); - registry.fill(HIST("QA/after/globalTracks_PVTracks"), collision.multNTracksPV(), tracks.size()); - registry.fill(HIST("QA/after/globalTracks_multT0A"), collision.multFT0A(), tracks.size()); - registry.fill(HIST("QA/after/globalTracks_multV0A"), collision.multFV0A(), tracks.size()); - registry.fill(HIST("QA/after/multV0A_multT0A"), collision.multFT0A(), collision.multFV0A()); - registry.fill(HIST("QA/after/multT0C_centT0C"), collision.centFT0C(), collision.multFT0C()); + static constexpr std::string_view Time[] = {"before", "after"}; + + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/hCent"), collision.centFT0C()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/globalTracks_centT0C"), collision.centFT0C(), tracks.size()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PVTracks_centT0C"), collision.centFT0C(), collision.multNTracksPV()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/globalTracks_PVTracks"), collision.multNTracksPV(), tracks.size()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/globalTracks_multT0A"), collision.multFT0A(), tracks.size()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/globalTracks_multV0A"), collision.multFV0A(), tracks.size()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/multV0A_multT0A"), collision.multFT0A(), collision.multFV0A()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/multT0C_centT0C"), collision.centFT0C(), collision.multFT0C()); + + if constexpr (framework::has_type_v) { + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/CentFT0C_vs_CentFT0Cvariant1"), collision.centFT0C(), collision.centFT0CVariant1()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/CentFT0C_vs_CentFT0M"), collision.centFT0C(), collision.centFT0M()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/CentFT0C_vs_CentFV0A"), collision.centFT0C(), collision.centFV0A()); + // registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/CentFT0C_vs_CentNGlobal"), collision.centFT0C(), collision.centNGlobal()); double psiA = 1.0 * std::atan2(collision.qyA(), collision.qxA()); double psiC = 1.0 * std::atan2(collision.qyC(), collision.qxC()); double psiFull = 1.0 * std::atan2(collision.qyA() + collision.qyC(), collision.qxA() + collision.qxC()); - registry.fill(HIST("QA/after/PsiA_vs_Cent"), psiA, collision.centFT0C()); - registry.fill(HIST("QA/after/PsiC_vs_Cent"), psiC, collision.centFT0C()); - registry.fill(HIST("QA/after/PsiFull_vs_Cent"), psiFull, collision.centFT0C()); - registry.fill(HIST("QA/after/PsiA_vs_Vx"), psiA, collision.vx()); - registry.fill(HIST("QA/after/PsiC_vs_Vx"), psiC, collision.vx()); - registry.fill(HIST("QA/after/PsiFull_vs_Vx"), psiFull, collision.vx()); - registry.fill(HIST("QA/after/PsiA_vs_Vy"), psiA, collision.vy()); - registry.fill(HIST("QA/after/PsiC_vs_Vy"), psiC, collision.vy()); - registry.fill(HIST("QA/after/PsiFull_vs_Vy"), psiFull, collision.vy()); - registry.fill(HIST("QA/after/PsiA_vs_Vz"), psiA, collision.posZ()); - registry.fill(HIST("QA/after/PsiC_vs_Vz"), psiC, collision.posZ()); - registry.fill(HIST("QA/after/PsiFull_vs_Vz"), psiFull, collision.posZ()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiA_vs_Cent"), psiA, collision.centFT0C()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiC_vs_Cent"), psiC, collision.centFT0C()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiFull_vs_Cent"), psiFull, collision.centFT0C()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiA_vs_Vx"), psiA, collision.vx()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiC_vs_Vx"), psiC, collision.vx()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiFull_vs_Vx"), psiFull, collision.vx()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiA_vs_Vy"), psiA, collision.vy()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiC_vs_Vy"), psiC, collision.vy()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiFull_vs_Vy"), psiFull, collision.vy()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiA_vs_Vz"), psiA, collision.posZ()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiC_vs_Vz"), psiC, collision.posZ()); + registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/PsiFull_vs_Vz"), psiFull, collision.posZ()); } return; } template - inline void fillHistograms(TrackObject track, float wacc, float weff, double ux, double uy, double uxMH, double uyMH, double qxA, double qyA, double qxC, double qyC, double corrQQx, double corrQQy, double corrQQ, double vnA, double vnC, double vnFull, double centrality, double vz) + inline void fillHistograms(TrackObject track, float wacc, float weff, double ux, double uy, double uxMH, double uyMH, double qxA, double qyA, double qxC, double qyC, double corrQQx, double corrQQy, double corrQQ, double vnA, double vnC, double vnFull, double centrality) { - static constexpr std::string_view Charge[] = {"incl/", "pos/", "neg/"}; - - registry.fill(HIST(Charge[ct]) + HIST("QA/hPt"), track.pt()); - registry.fill(HIST(Charge[ct]) + HIST("QA/hPhi"), track.phi()); - registry.fill(HIST(Charge[ct]) + HIST("QA/hEta"), track.eta()); - registry.fill(HIST(Charge[ct]) + HIST("QA/hPhi_Eta_vz"), track.phi(), track.eta(), vz); - registry.fill(HIST(Charge[ct]) + HIST("QA/hDCAxy"), track.dcaXY()); - registry.fill(HIST(Charge[ct]) + HIST("QA/hDCAz"), track.dcaZ()); - registry.fill(HIST(Charge[ct]) + HIST("vnAx_eta"), track.eta(), (ux * qxA) / std::sqrt(std::fabs(corrQQx)), wacc * weff); registry.fill(HIST(Charge[ct]) + HIST("vnAy_eta"), track.eta(), (uy * qyA) / std::sqrt(std::fabs(corrQQy)), wacc * weff); registry.fill(HIST(Charge[ct]) + HIST("vnCx_eta"), track.eta(), (ux * qxC) / std::sqrt(std::fabs(corrQQx)), wacc * weff); @@ -615,7 +642,19 @@ struct FlowSP { registry.fill(HIST(Charge[ct]) + HIST("vnFull_cent_EP"), centrality, vnFull, wacc * weff); } - void process(UsedCollisions::iterator const& collision, aod::BCsWithTimestamps const&, UsedTracks const& tracks) + template + inline void fillTrackQA(TrackObject track, double vz, float wacc = 1, float weff = 1) + { + registry.fill(HIST(Charge[ct]) + HIST("QA/hPt"), track.pt()); + registry.fill(HIST(Charge[ct]) + HIST("QA/hPhi"), track.phi()); + registry.fill(HIST(Charge[ct]) + HIST("QA/hPhiCorrected"), track.phi(), wacc * weff); + registry.fill(HIST(Charge[ct]) + HIST("QA/hEta"), track.eta()); + registry.fill(HIST(Charge[ct]) + HIST("QA/hPhi_Eta_vz"), track.phi(), track.eta(), vz); + registry.fill(HIST(Charge[ct]) + HIST("QA/hDCAxy"), track.dcaXY()); + registry.fill(HIST(Charge[ct]) + HIST("QA/hDCAz"), track.dcaZ()); + } + + void processData(UsedCollisions::iterator const& collision, aod::BCsWithTimestamps const&, UsedTracks const& tracks) { registry.fill(HIST("hEventCount"), evSel_FilteredEvent); @@ -627,16 +666,23 @@ struct FlowSP { cfg.correctionsLoaded = false; cfg.lastRunNumber = bc.runNumber(); } + fillEventQA(collision, tracks); loadCorrections(bc.timestamp()); - auto centrality = collision.centFT0C(); + float centrality = collision.centFT0C(); + + if (cfgFT0Cvariant1) + centrality = collision.centFT0CVariant1(); + if (cfgFT0M) + centrality = collision.centFT0M(); + if (cfgFV0A) + centrality = collision.centFV0A(); + // if (cfgNGlobal) centrality = collision.centNGlobal(); if (!eventSelected(collision, tracks.size(), centrality)) return; - fillEventQA(collision, tracks, true); - if (collision.isSelected()) { registry.fill(HIST("hEventCount"), evSel_isSelectedZDC); @@ -658,7 +704,7 @@ struct FlowSP { double psiFull = 1.0 * std::atan2(qyA + qyC, qxA + qxC); registry.fill(HIST("hSPplaneFull"), psiFull, 1); - fillEventQA(collision, tracks, false); + fillEventQA(collision, tracks); registry.fill(HIST("hCosPhiACosPhiC"), centrality, std::cos(psiA) * std::cos(psiC)); registry.fill(HIST("hSinPhiASinPhiC"), centrality, std::sin(psiA) * std::sin(psiC)); @@ -668,13 +714,17 @@ struct FlowSP { registry.fill(HIST("hFullEvPlaneRes"), centrality, -1 * std::cos(psiA - psiC)); registry.fill(HIST("qAqCXY"), centrality, qxA * qxC + qyA * qyC); + + registry.fill(HIST("qAXqCY"), centrality, qxA * qyC); + registry.fill(HIST("qAYqCX"), centrality, qyA * qxC); + registry.fill(HIST("qAqCX"), centrality, qxA * qxC); registry.fill(HIST("qAqCY"), centrality, qyA * qyC); - double corrQQ = 1.; - double corrQQx = 1.; - double corrQQy = 1.; + double corrQQ = 1., corrQQx = 1., corrQQy = 1.; + // Load correlations and SP resolution needed for Scalar Product and event plane methods. + // If not loaded set to 1 if (cfgLoadAverageQQ) { TList* hcorrList = ccdb->getForTimeStamp(cfgCCDBdir.value, bc.timestamp()); TProfile* hcorrQQ = reinterpret_cast(hcorrList->FindObject("qAqCXY")); @@ -702,12 +752,13 @@ struct FlowSP { float weffN = 1, waccN = 1; if (!trackSelected(track, field)) - continue; + return; if (track.sign() == 0.0) - continue; + return; bool pos = (track.sign() > 0) ? true : false; + // Fill NUA weights if (cfgFillWeights) { fWeights->Fill(track.phi(), track.eta(), vtxz, track.pt(), centrality, 0); } else if (cfgFillWeightsPOS) { @@ -718,16 +769,16 @@ struct FlowSP { fWeightsNEG->Fill(track.phi(), track.eta(), vtxz, track.pt(), centrality, 0); } + // Set weff and wacc for inclusice, negative and positive hadrons if (!setCurrentParticleWeights(kInclusive, weff, wacc, track.phi(), track.eta(), track.pt(), vtxz)) - continue; - + return; if (pos && !setCurrentParticleWeights(kPositive, weffP, waccP, track.phi(), track.eta(), track.pt(), vtxz)) - continue; - + return; if (!pos && !setCurrentParticleWeights(kNegative, weffN, waccN, track.phi(), track.eta(), track.pt(), vtxz)) - continue; + return; registry.fill(HIST("QA/after/hPt_inclusive"), track.pt()); + // // constrain angle to 0 -> [0,0+2pi] auto phi = RecoDecay::constrainAngle(track.phi(), 0); @@ -743,15 +794,90 @@ struct FlowSP { double vnC = std::cos(cfgHarm * (phi - psiC)) / evPlaneRes; double vnFull = std::cos(cfgHarm * (phi - psiFull)) / evPlaneRes; - fillHistograms(track, wacc, weff, ux, uy, uxMH, uyMH, qxA, qyA, qxC, qyC, corrQQx, corrQQy, corrQQ, vnA, vnC, vnFull, centrality, vtxz); + fillHistograms(track, wacc, weff, ux, uy, uxMH, uyMH, qxA, qyA, qxC, qyC, corrQQx, corrQQy, corrQQ, vnA, vnC, vnFull, centrality); + fillTrackQA(track, vtxz, wacc, weff); if (pos) { - fillHistograms(track, waccP, weffP, ux, uy, uxMH, uyMH, qxA, qyA, qxC, qyC, corrQQx, corrQQy, corrQQ, vnA, vnC, vnFull, centrality, vtxz); + fillHistograms(track, waccP, weffP, ux, uy, uxMH, uyMH, qxA, qyA, qxC, qyC, corrQQx, corrQQy, corrQQ, vnA, vnC, vnFull, centrality); + fillTrackQA(track, vtxz, wacc, weff); } else { - fillHistograms(track, waccN, weffN, ux, uy, uxMH, uyMH, qxA, qyA, qxC, qyC, corrQQx, corrQQy, corrQQ, vnA, vnC, vnFull, centrality, vtxz); + fillHistograms(track, waccN, weffN, ux, uy, uxMH, uyMH, qxA, qyA, qxC, qyC, corrQQx, corrQQy, corrQQ, vnA, vnC, vnFull, centrality); + fillTrackQA(track, vtxz, wacc, weff); } - } + } // end of track loop + } // end of collision isSelected loop + } + PROCESS_SWITCH(FlowSP, processData, "Process analysis for non-derived data", true); + + void processMCReco(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, soa::Filtered> const& tracks, aod::McParticles const&) + { + auto bc = collision.template bc_as(); + auto field = (cfgMagField == 99999) ? getMagneticField(bc.timestamp()) : cfgMagField; + + double vtxz = collision.posZ(); + float centrality = collision.centFT0C(); + + fillEventQA(collision, tracks); + + if (!eventSelected(collision, tracks.size(), centrality)) + return; + + fillEventQA(collision, tracks); + + for (const auto& track : tracks) { + + auto mcParticle = track.mcParticle(); + if (!mcParticle.isPhysicalPrimary()) + return; + + if (mcParticle.eta() < -cfgEta || mcParticle.eta() > cfgEta || mcParticle.pt() < cfgPtmin || mcParticle.pt() > cfgPtmax || track.tpcNClsFound() < cfgNcls) + return; + + registry.fill(HIST("QA/before/hPt_inclusive"), track.pt()); + + if (!trackSelected(track, field)) + return; + + registry.fill(HIST("QA/after/hPt_inclusive"), track.pt()); + + fillTrackQA(track, vtxz); + + } // end of track loop + } + PROCESS_SWITCH(FlowSP, processMCReco, "Process analysis for MC reconstructed events", false); + + Filter mcCollFilter = nabs(aod::mccollision::posZ) < cfgVtxZ; + void processMCGen(soa::Filtered::iterator const& mcCollision, soa::SmallGroups> const& collisions, aod::McParticles const& particles) + { + if (collisions.size() != 1) + return; + float centrality = -1; + for (const auto& collision : collisions) { + centrality = collision.centFT0C(); + } + + if (particles.size() < 1) + return; + if (centrality < cfgCentMin || centrality > cfgCentMax) + return; + + float vtxz = mcCollision.posZ(); + + for (const auto& track : particles) { + + if (!track.isPhysicalPrimary()) + continue; + + registry.fill(HIST("trackMCGen/before/pt_gen_incl"), track.pt()); + registry.fill(HIST("trackMCGen/before/phi_eta_vtxZ_gen"), track.phi(), track.eta(), vtxz); + + if (track.eta() < -cfgEta || track.eta() > cfgEta || track.pt() < cfgPtmin || track.pt() > cfgPtmax) + return; + + registry.fill(HIST("trackMCGen/after/pt_gen_incl"), track.pt()); + registry.fill(HIST("trackMCGen/after/phi_eta_vtxZ_gen"), track.phi(), track.eta(), vtxz); } } + PROCESS_SWITCH(FlowSP, processMCGen, "Process analysis for MC generated events", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)