[PWGDQ] Added a process function in dqFlow to compute FT0C cumulant for flow analysis

PWGDQ/Tasks/dqFlow.cxx

@@ -1,4 +1,4 @@
 // 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.
 //
@@ -142,8 +142,8 @@
   ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1}, "multiplicity / centrality axis for histograms"};
 
   // Define the filter for barrel tracks and forward tracks
-  Filter trackFilter = (requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true);
+  Filter trackFilter = (requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true);
   Filter fwdFilter = (aod::fwdtrack::eta < -2.45f) && (aod::fwdtrack::eta > -3.6f);
 
   // Histograms used for optionnal efficiency and non-uniform acceptance corrections
   struct Config {
@@ -352,6 +352,52 @@
     }
   }
 
+  // Templated function instantianed for all of the process functions
+  template <uint32_t TEventFillMap, typename TEvent, typename TFT0s>
+  void runFillQvectorFromCentralFWWithFT0CCumulants(TEvent const& collision, TFT0s const& /*ft0s*/)
+  {
+    VarManager::ResetValues(0, VarManager::kNVars);
+    VarManager::FillEvent<TEventFillMap>(collision);
+
+    float S11C = collision.sumAmplFT0C();
+    float S12C = 0.f;
+    float S21C = S11C * S11C;
+
+    // Compute sum of squares of amplitudes from FT0C for flow analysis
+    if (collision.has_foundFT0()) {
+      auto ft0 = collision.foundFT0();
+      auto const& ampC = ft0.amplitudeC();
+      for (auto amp : ampC) {
+        if (amp > 0.f) {
+          S12C += amp * amp;
+        }
+      }
+    }
+    VarManager::FillQVectorFromCentralFW(collision);
+    std::complex<double> Q21C(collision.qvecFT0CRe() * S11C, collision.qvecFT0CIm() * S11C);
+
+    // Fill necessary quantities for cumulant calculations with weighted Q-vectors
+    float M11REF = S21C - S12C;
+    float CORR2REF = (norm(Q21C) - S12C) / M11REF;
+    CORR2REF = std::isnan(M11REF) || std::isinf(M11REF) || std::isnan(CORR2REF) || std::isinf(CORR2REF) ? 0 : CORR2REF;
+    M11REF = std::isnan(M11REF) || std::isinf(M11REF) || std::isnan(CORR2REF) || std::isinf(CORR2REF) ? 0 : M11REF;
+
+    if (fConfigQA) {
+      fHistMan->FillHistClass("Event_BeforeCuts_centralFW", VarManager::fgValues);
+      if (fEventCut->IsSelected(VarManager::fgValues)) {
+        fHistMan->FillHistClass("Event_AfterCuts_centralFW", VarManager::fgValues);
+      }
+    }
+
+    if (fEventCut->IsSelected(VarManager::fgValues)) {
+      eventQvectorCentr(collision.qvecFT0ARe(), collision.qvecFT0AIm(), collision.qvecFT0CRe(), collision.qvecFT0CIm(), collision.qvecFT0MRe(), collision.qvecFT0MIm(), collision.qvecFV0ARe(), collision.qvecFV0AIm(), collision.qvecTPCposRe(), collision.qvecTPCposIm(), collision.qvecTPCnegRe(), collision.qvecTPCnegIm(),
+                        collision.sumAmplFT0A(), collision.sumAmplFT0C(), collision.sumAmplFT0M(), collision.sumAmplFV0A(), collision.nTrkTPCpos(), collision.nTrkTPCneg());
+      eventQvectorCentrExtra(collision.qvecTPCallRe(), collision.qvecTPCallIm(), collision.nTrkTPCall());
+      eventRefFlow(M11REF, -9999, -9999, CORR2REF, -9999, -9999, VarManager::fgValues[VarManager::kCentFT0C]);
+      eventQvectorExtra(-9999, -9999, -9999, -9999, S11C, S12C, -9999, -9999);
+    }
+  }
+
   // Templated function instantianed for all of the process functions
   template <uint32_t TEventFillMap, uint32_t TTrackFillMap, typename TEvent, typename TTracks>
   void runFillQvector(TEvent const& collision, MyBcs const&, TTracks const& tracks1, aod::Zdcs const&)
@@ -385,7 +431,7 @@
     float weff = 1.0, wacc = 1.0;
 
     // Fill the GFW object in the track loop
     for (auto& track : tracks1) {
 
       // Selections for barrel tracks
       if (!isTrackSelected<TTrackFillMap>(track)) {
@@ -475,10 +521,10 @@
       S12Full = gfwCumFull.Vec(0, 2).real();
       S13Full = gfwCumFull.Vec(0, 3).real();
       S14Full = gfwCumFull.Vec(0, 4).real();
       S21Full = pow(gfwCumFull.Vec(0, 1).real(), 2);
       S22Full = pow(gfwCumFull.Vec(0, 2).real(), 2);
       S31Full = pow(gfwCumFull.Vec(0, 1).real(), 3);
       S41Full = pow(gfwCumFull.Vec(0, 1).real(), 4);
 
       // Get the Q vector for selected harmonic, power (for minPt=0)
       Q11N = gfwCumN.Vec(1, 1);
@@ -574,6 +620,12 @@
     runFillQvectorFromCentralFW<gkEventFillMapRun3>(collisions);
   }
 
+  // Process to fill Q vector using barrel tracks in a reduced event table for barrel/muon tracks flow related analyses Run 3
+  void processCentralQvectorWithFT0CCumulants(MyEventsWithCentQvectRun3::iterator const& collision, aod::FT0s const& ft0s)
+  {
+    runFillQvectorFromCentralFWWithFT0CCumulants<gkEventFillMapRun3>(collision, ft0s);
+  }
+
   // Process to fill Q vector using forward tracks in a reduced event table for barrel/muon tracks flow related analyses Run 3
   void processForwardQvector(MyEventsWithCentRun3::iterator const& collisions, MyBcs const& bcs, soa::Filtered<aod::MFTTracks> const& tracks, aod::Zdcs const& zdcs)
   {
@@ -581,7 +633,7 @@
   }
 
   // TODO: dummy function for the case when no process function is enabled
   void processDummy(MyEvents&)
   {
     // do nothing
   }
@@ -590,6 +642,7 @@
   PROCESS_SWITCH(DQEventQvector, processBarrelQvector, "Run q-vector task on barrel tracks for Run3", false);
   PROCESS_SWITCH(DQEventQvector, processAllQvector, "Run q-vector task on barrel tracks for Run3 and using central q-vector", false);
   PROCESS_SWITCH(DQEventQvector, processCentralQvector, "Run q-vector task using central q-vector", false);
+  PROCESS_SWITCH(DQEventQvector, processCentralQvectorWithFT0CCumulants, "Run q-vector task using central q-vector with FT0C cumulants", false);
   PROCESS_SWITCH(DQEventQvector, processForwardQvector, "Run q-vector task on forward tracks for Run3", false);
   PROCESS_SWITCH(DQEventQvector, processDummy, "Dummy function", false);
 };