diff --git a/Tutorials/PWGCF/EventPlane/src/CMakeLists.txt b/Tutorials/PWGCF/EventPlane/src/CMakeLists.txt index 86e39d138cf..bf15a1a13c9 100644 --- a/Tutorials/PWGCF/EventPlane/src/CMakeLists.txt +++ b/Tutorials/PWGCF/EventPlane/src/CMakeLists.txt @@ -13,3 +13,8 @@ o2physics_add_dpl_workflow(qvectors-tutorial SOURCES qVectorstutorial.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(spectator-plane-tutorial + SOURCES spectatorPlaneTutorial.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/Tutorials/PWGCF/EventPlane/src/spectatorPlaneTutorial.cxx b/Tutorials/PWGCF/EventPlane/src/spectatorPlaneTutorial.cxx new file mode 100644 index 00000000000..6790dd30ae5 --- /dev/null +++ b/Tutorials/PWGCF/EventPlane/src/spectatorPlaneTutorial.cxx @@ -0,0 +1,237 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file SpectatorPlaneTutorial.cxx +/// \author Noor Koster +/// \since 11/2025 +/// \brief This is a tutorial task to show how to use the ZDC q-vectors and the spectator plane resolution. + +#include "PWGCF/DataModel/SPTableZDC.h" + +#include "Common/Core/EventPlaneHelper.h" +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/Qvectors.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "DataFormatsParameters/GRPLHCIFData.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/runDataProcessing.h" + +#include "TF1.h" +#include "TPDGCode.h" + +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod::rctsel; +// using namespace o2::analysis; + +#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; + +struct SpectatorPlaneTutorial { + RCTFlagsChecker rctChecker; + + struct : ConfigurableGroup { + O2_DEFINE_CONFIGURABLE(cfgEvtUseRCTFlagChecker, bool, false, "Evt sel: use RCT flag checker"); + O2_DEFINE_CONFIGURABLE(cfgEvtRCTFlagCheckerLabel, std::string, "CBT_hadronPID", "Evt sel: RCT flag checker label (CBT, CBT_hadronPID)"); // all Labels can be found in Common/CCDB/RCTSelectionFlags.h + O2_DEFINE_CONFIGURABLE(cfgEvtRCTFlagCheckerZDCCheck, bool, false, "Evt sel: RCT flag checker ZDC check"); + O2_DEFINE_CONFIGURABLE(cfgEvtRCTFlagCheckerLimitAcceptAsBad, bool, false, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"); + } rctFlags; + + // settings for CCDB data + O2_DEFINE_CONFIGURABLE(cfgCCDBdir_QQ, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/meanQQ/Default", "ccdb dir for average QQ values in 1% centrality bins"); + O2_DEFINE_CONFIGURABLE(cfgCCDBdir_SP, std::string, "Users/c/ckoster/ZDC/LHC23_zzh_pass4/SPPlaneRes", "ccdb dir for average event plane resolution in 1% centrality bins"); + // Confogirable axis + ConfigurableAxis axisCentrality{"axisCentrality", {20, 0, 100}, "Centrality bins for vn "}; + + // Configurables containing vector + float vtxZ = 10.0; + float etaMax = 0.8; + float ptMin = 0.2; + float ptMax = 10.0; + float dcaXYMax = 0.2; + float dcaZMax = 2.0; + + Filter collisionFilter = nabs(aod::collision::posZ) < vtxZ; + Filter trackFilter = nabs(aod::track::eta) < etaMax && aod::track::pt > ptMin&& aod::track::pt < ptMax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && nabs(aod::track::dcaXY) < dcaXYMax&& nabs(aod::track::dcaZ) < dcaZMax; + using GeneralCollisions = soa::Join; + using UnfilteredTracks = soa::Join; + + using UsedTracks = soa::Filtered; + using ZDCCollisions = soa::Filtered>; // IMPORTANT: ZDCCollisions must be used to get access to ZDC q-vectors --> PWGCF/DataModel/SPTableZDC.h & PWGCF/Flow/TableProducer/zdcQVectors.cxx + + // Connect to ccdb + Service ccdb; + + HistogramRegistry registry{"registry"}; + + void init(InitContext const&) + { + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + + int64_t now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); + ccdb->setCreatedNotAfter(now); + + AxisSpec axisPhi = {60, 0, constants::math::TwoPI, "#varphi"}; + AxisSpec axisEta = {64, -1.6, 1.6, "#eta"}; + AxisSpec axisEtaVn = {8, -.8, .8, "#eta"}; + AxisSpec axisCent = {90, 0, 90, "Centrality(%)"}; + AxisSpec axisPhiPlane = {100, -constants::math::PI, constants::math::PI, "#Psi"}; + AxisSpec axisQQ = {100, -0.2, 0.2, "#LT Q_{X}^{A}Q_{Y}^{C} #GT"}; + + std::vector ptbinning = {0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.5, 4, 5, 6, 8, 10}; + AxisSpec axisPt = {ptbinning, "#it{p}_{T} GeV/#it{c}"}; + + rctChecker.init(rctFlags.cfgEvtRCTFlagCheckerLabel, rctFlags.cfgEvtRCTFlagCheckerZDCCheck, rctFlags.cfgEvtRCTFlagCheckerLimitAcceptAsBad); + + registry.add("hCentrality", "Centrality; #Events; Centrality (%)", {HistType::kTH1D, {axisCent}}); + registry.add("hSPplaneA", "Spectator Plane Angle A; #Events; #Psi_{A}", {HistType::kTH1D, {axisPhiPlane}}); + registry.add("qAqCXY", "Qx ZDCA x Qy ZDCC vs cent ; #ltQ_{X}^{A}Q_{Y}^{C}#GT; Centrality (%)", {HistType::kTH2D, {axisCent, axisQQ}}); + // Add here the histograms for the spectator plane angle C + // registry.add(.........) + registry.add("hSPplaneAvsC", "Spectator Plane Angle A vs C; #Events; #Psi_{A}; #Psi_{C}", {HistType::kTH2D, {axisPhiPlane, axisPhiPlane}}); + registry.add("hSPplaneFull", "Spectator Plane Angle Full; #Events; #Psi_{Full}", {HistType::kTH1D, {axisPhiPlane}}); + + // Note: we will fill these with data from the CCDB, just to take a look! + registry.add("CalibHistos/hQQx", "QQx; #Events; QQx", {HistType::kTProfile, {axisCent}}); + registry.add("CalibHistos/hQQy", "QQy; #Events; QQy", {HistType::kTProfile, {axisCent}}); + registry.add("CalibHistos/hQQ", "QQ; #Events; QQx + QQy", {HistType::kTProfile, {axisCent}}); + // Add here the histograms for the cross-terms of the Q-vectors from the ZDC + + registry.add("CalibHistos/hEvPlaneRes", "Event Plane Resolution; #Events; Event Plane Resolution", {HistType::kTProfile, {axisCent}}); + + // Flow Histograms + registry.add("flow/v1A", "", {HistType::kTProfile, {axisPt}}); + registry.add("flow/v1C", "", {HistType::kTProfile, {axisPt}}); + + registry.add("flow/vnAxCxUxMH", "", {HistType::kTProfile, {axisCent}}); + registry.add("flow/vnAyCyUxMH", "", {HistType::kTProfile, {axisCent}}); + registry.add("flow/vnAxCyUyMH", "", {HistType::kTProfile, {axisCent}}); + registry.add("flow/vnAyCxUyMH", "", {HistType::kTProfile, {axisCent}}); + } + + void process(ZDCCollisions::iterator const& collision, aod::BCsWithTimestamps const&, UsedTracks const& tracks) + { + + auto bc = collision.bc_as(); + + float centrality = collision.centFT0C(); + registry.fill(HIST("hCentrality"), centrality); + + float centMin = 0; + float centMax = 80; + + if (centrality > centMax || centrality < centMin) + return; + + if (collision.isSelected() == false) + return; + + // Get the ZDC q-vectors + double qxA = collision.qxA(); + double qyA = collision.qyA(); + double qxC = collision.qxC(); + double qyC = collision.qyC(); + + // Calculate the spectator plane angles + double psiA = 1.0 * std::atan2(qyA, qxA); + registry.fill(HIST("hSPplaneA"), psiA); + + // Add the PsiA vs PsiC as a TH2D + + double psiC = 1.0 * std::atan2(qyC, qxC); + double psiFull = 1.0 * std::atan2(qyA + qyC, qxA + qxC); + registry.fill(HIST("hSPplaneFull"), psiFull); + + registry.fill(HIST("hSPplaneAvsC"), psiA, psiC); + + // Fill the q-vector correlations + registry.fill(HIST("qAqCXY"), centrality, qxA * qxC + qyA * qyC); + + double corrQQx = 1; + double corrQQy = 1; + double corrQQ = 1; + double evPlaneRes = 1; + + // Get QQ-correlations from CCDB + if (cfgCCDBdir_QQ.value.empty() == false) { + TList* list = ccdb->getForTimeStamp(cfgCCDBdir_QQ.value, bc.timestamp()); + TProfile* qAqCX = reinterpret_cast(list->FindObject("qAqCX")); + TProfile* qAqCY = reinterpret_cast(list->FindObject("qAqCY")); + TProfile* qAqCXY = reinterpret_cast(list->FindObject("qAqCXY")); + // The sum is qAqCXY + corrQQx = qAqCX->GetBinContent(centrality); + corrQQy = qAqCY->GetBinContent(centrality); + corrQQ = qAqCXY->GetBinContent(centrality); + registry.fill(HIST("CalibHistos/hQQx"), centrality, corrQQx); + registry.fill(HIST("CalibHistos/hQQy"), centrality, corrQQy); + registry.fill(HIST("CalibHistos/hQQ"), centrality, corrQQ); + } + // Get event plane resolution from CCDB + if (cfgCCDBdir_SP.value.empty() == false) { + evPlaneRes = ccdb->getForTimeStamp(cfgCCDBdir_SP.value, bc.timestamp())->GetBinContent(centrality); + registry.fill(HIST("CalibHistos/hEvPlaneRes"), centrality, evPlaneRes); + } + + for (const auto& track : tracks) { + + // constrain angle to 0 -> [0,0+2pi] + auto phi = RecoDecay::constrainAngle(track.phi(), 0); + + double ux = std::cos(phi); + double uy = std::sin(phi); + + double uxMH = std::cos(2 * phi); + double uyMH = std::sin(2 * phi); + + double v1A = (uy * qyA + ux * qxA) / std::sqrt(std::fabs(corrQQ)); + double v1C = (uy * qyC + ux * qxC) / std::sqrt(std::fabs(corrQQ)); + + double v2AxCxUxMH = (uxMH * qxA * qxC) / corrQQx; + double v2AyCyUxMH = (uxMH * qyA * qyC) / corrQQy; + double v2AxCyUyMH = (uyMH * qxA * qyC) / corrQQx; + double v2AyCxUyMH = (uyMH * qyA * qxC) / corrQQy; + + registry.fill(HIST("flow/v1A"), track.eta(), v1A); + registry.fill(HIST("flow/v1C"), track.eta(), v1C); + + registry.fill(HIST("flow/v2AxCxUxMH"), centrality, v2AxCxUxMH); + registry.fill(HIST("flow/v2AyCyUxMH"), centrality, v2AyCyUxMH); + registry.fill(HIST("flow/v2AxCyUyMH"), centrality, v2AxCyUyMH); + registry.fill(HIST("flow/v2AyCxUyMH"), centrality, v2AyCxUyMH); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +}