Add task for flow with dihadron correlation study

Author: zhaomr13

PWGUD/Tasks/CMakeLists.txt

@@ -244,4 +244,7 @@ o2physics_add_dpl_workflow(flow-cumulants-upc
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::GFWCore
                            COMPONENT_NAME Analysis)
 
-
+o2physics_add_dpl_workflow(flow-correlations-upc
+                           SOURCES flowCorrelationsUpc.cxx
+			   PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGCFCore
+                           COMPONENT_NAME Analysis)

PWGUD/Tasks/flowCorrelationsUpc.cxx

@@ -0,0 +1,264 @@
+// 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 flowCorrelationsUpc.cxx
+/// \brief Provides a sparse with usefull two particle correlation info
+/// \author Mingrui Zhao (mingrui.zhao@cern.ch, mingrui.zhao@mail.labz0.org)
+/// copied from Thor Jensen (thor.kjaersgaard.jensen@cern.ch) and Debojit Sarkar (debojit.sarkar@cern.ch)
+
+#include <vector>
+
+#include "Framework/runDataProcessing.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/ASoAHelpers.h"
+#include "Framework/ASoA.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/RunningWorkflowInfo.h"
+#include "CommonConstants/MathConstants.h"
+#include "CCDB/BasicCCDBManager.h"
+#include "Common/Core/RecoDecay.h"
+
+#include "PWGUD/DataModel/UDTables.h"
+#include "PWGUD/Core/SGSelector.h"
+
+#include "Common/DataModel/EventSelection.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+#include "Common/DataModel/Centrality.h"
+#include "Common/DataModel/Multiplicity.h"
+#include "PWGCF/DataModel/CorrelationsDerived.h"
+#include "PWGCF/Core/CorrelationContainer.h"
+#include "PWGCF/Core/PairCuts.h"
+
+namespace o2::aod
+{
+namespace flowcorrupc
+{
+DECLARE_SOA_COLUMN(Multiplicity, multiplicity, int);
+}
+DECLARE_SOA_TABLE(Multiplicity, "AOD", "MULTIPLICITY",
+                  flowcorrupc::Multiplicity);
+
+} // namespace o2::aod
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+
+// define the filtered collisions and tracks
+#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};
+
+struct CalcNchUpc {
+  O2_DEFINE_CONFIGURABLE(cfgZVtxCut, float, 10.0f, "Accepted z-vertex range")
+  O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.2f, "minimum accepted track pT")
+  O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 10.0f, "maximum accepted track pT")
+  O2_DEFINE_CONFIGURABLE(cfgEtaCut, float, 0.8f, "Eta cut")
+  O2_DEFINE_CONFIGURABLE(cfgMinMixEventNum, int, 5, "Minimum number of events to mix")
+
+  // Filter trackFilter = (nabs(aod::track::eta) < cfgEtaCut) && (aod::track::pt > cfgPtCutMin) && (aod::track::pt < cfgPtCutMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true));
+
+
+  using UdTracks = soa::Join<aod::UDTracks, aod::UDTracksExtra, aod::UDTracksPID>;
+  using UdTracksFull = soa::Join<aod::UDTracks, aod::UDTracksPID, aod::UDTracksExtra, aod::UDTracksFlags, aod::UDTracksDCA>;
+  using UDCollisionsFull = soa::Join<aod::UDCollisions, aod::SGCollisions, aod::UDCollisionsSels, aod::UDZdcsReduced>;
+
+  Produces<aod::Multiplicity> multiplicityNch;
+
+  HistogramRegistry registry{"registry"};
+
+  void init(InitContext&)
+  {
+    AxisSpec axisNch = {100, 0, 100};
+    AxisSpec axisVrtx = {10, -10, 10};
+
+    registry.add("Ncharge", "N_{charge}", {HistType::kTH1D, {axisNch}});
+    registry.add("zVtx_all", "zVtx_all", {HistType::kTH1D, {axisVrtx}});
+  }
+
+  void process(UDCollisionsFull::iterator const& collision, UdTracksFull const& tracks)
+  {
+    multiplicityNch(tracks.size());
+    registry.fill(HIST("Ncharge"), tracks.size());
+    registry.fill(HIST("zVtx_all"), collision.posZ());
+  }
+};
+
+struct FlowCorrelationsUpc {
+  O2_DEFINE_CONFIGURABLE(cfgZVtxCut, float, 10.0f, "Accepted z-vertex range")
+  O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.2f, "minimum accepted track pT")
+  O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 10.0f, "maximum accepted track pT")
+  O2_DEFINE_CONFIGURABLE(cfgEtaCut, float, 0.8f, "Eta cut")
+  O2_DEFINE_CONFIGURABLE(cfgMinMixEventNum, int, 5, "Minimum number of events to mix")
+  O2_DEFINE_CONFIGURABLE(cfgMinMult, int, 0, "Minimum multiplicity for collision")
+  O2_DEFINE_CONFIGURABLE(cfgMaxMult, int, 10, "Maximum multiplicity for collision")
+
+  ConfigurableAxis axisVertex{"axisVertex", {10, -10, 10}, "vertex axis for histograms"};
+  ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"};
+  ConfigurableAxis axisPhi{"axisPhi", {72, 0.0, constants::math::TwoPI}, "phi axis for histograms"};
+  ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt axis for histograms"};
+  ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"};
+  ConfigurableAxis axisDeltaEta{"axisDeltaEta", {40, -2, 2}, "delta eta axis for histograms"};
+  ConfigurableAxis axisPtTrigger{"axisPtTrigger", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt trigger axis for histograms"};
+  ConfigurableAxis axisPtAssoc{"axisPtAssoc", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt associated axis for histograms"};
+  ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100}, "multiplicity / centrality axis for histograms"};
+  ConfigurableAxis vtxMix{"vtxMix", {VARIABLE_WIDTH, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, "vertex axis for mixed event histograms"};
+  ConfigurableAxis multMix{"multMix", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100}, "multiplicity / centrality axis for mixed event histograms"};
+
+  ConfigurableAxis axisVertexEfficiency{"axisVertexEfficiency", {10, -10, 10}, "vertex axis for efficiency histograms"};
+  ConfigurableAxis axisEtaEfficiency{"axisEtaEfficiency", {20, -1.0, 1.0}, "eta axis for efficiency histograms"};
+  ConfigurableAxis axisPtEfficiency{"axisPtEfficiency", {VARIABLE_WIDTH, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0}, "pt axis for efficiency histograms"};
+
+  // make the filters and cuts.
+  // Filter collisionFilter = (nabs(aod::collision::posZ) < cfgZVtxCut) && (aod::flowcorrupc::multiplicity) > cfgMinMult && (aod::flowcorrupc::multiplicity) < cfgMaxMult && (aod::evsel::sel8) == true;
+  // Filter trackFilter = (nabs(aod::track::eta) < cfgEtaCut) && (aod::track::pt > cfgPtCutMin) && (aod::track::pt < cfgPtCutMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true));
+
+  using UdTracks = soa::Join<aod::UDTracks, aod::UDTracksExtra, aod::UDTracksPID>;
+  using UdTracksFull = soa::Join<aod::UDTracks, aod::UDTracksPID, aod::UDTracksExtra, aod::UDTracksFlags, aod::UDTracksDCA>;
+  using UDCollisionsFull = soa::Join<aod::UDCollisions, aod::SGCollisions, aod::UDCollisionsSels, aod::UDZdcsReduced, aod::Multiplicity>;
+
+  // Define the outputs
+  OutputObj<CorrelationContainer> same{Form("sameEvent_%i_%i", static_cast<int>(cfgMinMult), static_cast<int>(cfgMaxMult))};
+  OutputObj<CorrelationContainer> mixed{Form("mixedEvent_%i_%i", static_cast<int>(cfgMinMult), static_cast<int>(cfgMaxMult))};
+
+  HistogramRegistry registry{"registry"};
+
+  void init(InitContext&)
+  {
+    LOGF(info, "Starting init");
+    // Make histograms to check the distributions after cuts
+    registry.add("deltaEta_deltaPhi_same", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}}); // check to see the delta eta and delta phi distribution
+    registry.add("deltaEta_deltaPhi_mixed", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}});
+    registry.add("Phi", "Phi", {HistType::kTH1D, {axisPhi}});
+    registry.add("Eta", "Eta", {HistType::kTH1D, {axisEta}});
+    registry.add("pT", "pT", {HistType::kTH1D, {axisPtTrigger}});
+    registry.add("Nch", "N_{ch}", {HistType::kTH1D, {axisMultiplicity}});
+    registry.add("zVtx", "zVtx", {HistType::kTH1D, {axisVertex}});
+
+    registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisMultiplicity, axisVertex, axisPtTrigger}}});
+
+    registry.add("eventcount", "bin", {HistType::kTH1F, {{3, 0, 3, "bin"}}}); // histogram to see how many events are in the same and mixed event
+
+    std::vector<AxisSpec> corrAxis = {{axisMultiplicity, "Nch"},
+                                      {axisVertex, "z-vtx (cm)"},
+                                      {axisPtTrigger, "p_{T} (GeV/c)"},
+                                      {axisPtAssoc, "p_{T} (GeV/c)"},
+                                      {axisDeltaPhi, "#Delta#varphi (rad)"},
+                                      {axisDeltaEta, "#Delta#eta"}};
+    std::vector<AxisSpec> effAxis = {
+      {axisVertexEfficiency, "z-vtx (cm)"},
+      {axisPtEfficiency, "p_{T} (GeV/c)"},
+      {axisEtaEfficiency, "#eta"},
+    };
+    std::vector<AxisSpec> userAxis;
+
+    same.setObject(new CorrelationContainer(Form("sameEvent_%i_%i", static_cast<int>(cfgMinMult), static_cast<int>(cfgMaxMult)), Form("sameEvent_%i_%i", static_cast<int>(cfgMinMult), static_cast<int>(cfgMaxMult)), corrAxis, effAxis, userAxis));
+    mixed.setObject(new CorrelationContainer(Form("mixedEvent_%i_%i", static_cast<int>(cfgMinMult), static_cast<int>(cfgMaxMult)), Form("mixedEvent_%i_%i", static_cast<int>(cfgMinMult), static_cast<int>(cfgMaxMult)), corrAxis, effAxis, userAxis));
+  }
+  enum EventType {
+    SameEvent = 1,
+    MixedEvent = 3
+  };
+
+  // fill multiple histograms
+  template <typename TCollision, typename TTracks>
+  void fillYield(TCollision collision, TTracks tracks) // function to fill the yield and etaphi histograms.
+  {
+    registry.fill(HIST("Nch"), tracks.size());
+    registry.fill(HIST("zVtx"), collision.posZ());
+
+    for (auto const& track1 : tracks) {
+      auto momentum1 = std::array<double, 3>{track1.px(), track1.py(), track1.pz()};
+      registry.fill(HIST("Phi"), RecoDecay::phi(momentum1));
+      registry.fill(HIST("Eta"), RecoDecay::eta(momentum1));
+      registry.fill(HIST("pT"), track1.pt());
+    }
+  }
+
+  template <CorrelationContainer::CFStep step, typename TTracks>
+  void fillCorrelations(TTracks tracks1, TTracks tracks2, float posZ, int system, float Nch) // function to fill the Output functions (sparse) and the delta eta and delta phi histograms
+  {
+    // loop over all tracks
+    for (auto const& track1 : tracks1) {
+
+      if (system == SameEvent) {
+        registry.fill(HIST("Trig_hist"), Nch, posZ, track1.pt());
+      }
+
+      for (auto const& track2 : tracks2) {
+
+        if (track1.pt() <= track2.pt())
+          continue; // skip if the trigger pt is less than the associate p
+
+	auto momentum1 = std::array<double, 3>{track1.px(), track1.py(), track1.pz()};
+	auto momentum2 = std::array<double, 3>{track2.px(), track2.py(), track2.pz()};
+        double phi1 = RecoDecay::phi(momentum1); 
+	double phi2 = RecoDecay::phi(momentum2);
+        float deltaPhi = RecoDecay::constrainAngle(phi1-phi2, -PIHalf);
+        float deltaEta = RecoDecay::eta(momentum1) - RecoDecay::eta(momentum2);
+
+        // fill the right sparse and histograms
+        if (system == SameEvent) {
+          same->getPairHist()->Fill(step, Nch, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta);
+          registry.fill(HIST("deltaEta_deltaPhi_same"), deltaPhi, deltaEta);
+        } else if (system == MixedEvent) {
+          mixed->getPairHist()->Fill(step, Nch, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta);
+          registry.fill(HIST("deltaEta_deltaPhi_mixed"), deltaPhi, deltaEta);
+        }
+      }
+    }
+  }
+
+  void processSame(UDCollisionsFull::iterator const& collision, UdTracksFull const& tracks)
+  {
+    int gapSide = collision.gapSide();
+    if (gapSide < 0 || gapSide > 2) {
+      return;
+    }
+
+    int trueGapSide = sgSelector.trueGap(collision, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC);
+    gapSide = trueGapSide;
+    if (gapSide == 2) {
+      return;
+    }
+
+    registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin
+    fillYield(collision, tracks);
+    fillCorrelations<CorrelationContainer::kCFStepReconstructed>(tracks, tracks, collision.posZ(), SameEvent, tracks.size()); // fill the SE histogram and Sparse
+  }
+  PROCESS_SWITCH(FlowCorrelationsUpc, processSame, "Process same event", true);
+
+  // event mixing
+
+  SliceCache cache;
+  using MixedBinning = ColumnBinningPolicy<aod::collision::PosZ, aod::flowcorrupc::Multiplicity>;
+
+  // the process for filling the mixed events
+  void processMixed(UDCollisionsFull const& collisions, UdTracksFull const& tracks)
+  {
+    MixedBinning binningOnVtxAndMult{{vtxMix, multMix}, true}; // true is for 'ignore overflows' (true by default)
+    auto tracksTuple = std::make_tuple(tracks);
+    SameKindPair<UDCollisionsFull, UdTracksFull, MixedBinning> pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip
+
+    for (auto const& [collision1, tracks1, collision2, tracks2] : pairs) {
+      registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin
+      fillCorrelations<CorrelationContainer::kCFStepReconstructed>(tracks1, tracks2, collision1.posZ(), MixedEvent, tracks1.size());
+    }
+  }
+  PROCESS_SWITCH(FlowCorrelationsUpc, processMixed, "Process mixed events", true);
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<CalcNchUpc>(cfgc),
+    adaptAnalysisTask<FlowCorrelationsUpc>(cfgc),
+  };
+}