Upgrades and fixes for coding violations from FC
authorhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 14 Mar 2011 10:05:52 +0000 (10:05 +0000)
committerhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 14 Mar 2011 10:05:52 +0000 (10:05 +0000)
PWG2/CMakelibPWG2forward2.pkg
PWG2/FORWARD/analysis2/AddTaskForwardFlow.C
PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.cxx
PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.h
PWG2/FORWARD/analysis2/AliForwardFlowUtil.cxx [new file with mode: 0644]
PWG2/FORWARD/analysis2/AliForwardFlowUtil.h [new file with mode: 0644]
PWG2/FORWARD/analysis2/MakeFlow.C
PWG2/PWG2forward2LinkDef.h

index 5023a72..22d7cf7 100644 (file)
@@ -57,18 +57,18 @@ set ( SRCS
   FORWARD/analysis2/AliCentralCorrAcceptance.cxx 
   FORWARD/analysis2/AliCentraldNdetaTask.cxx
   FORWARD/analysis2/AliBasedNdetaTask.cxx
-  FORWARD/analysis2/AliForwardFlowBase.cxx
+  FORWARD/analysis2/AliForwardFlowUtil.cxx
   FORWARD/analysis2/AliForwardFlowTaskQC.cxx
 )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
 set ( HDRS ${HDRS} FORWARD/analysis2/AliFMDStripIndex.h )
 
-set ( EINCLUDE  ANALYSIS PWG2/FORWARD/analysis2)
+set ( EINCLUDE  ANALYSIS PWG2/FORWARD/analysis2 STEER)
 
 set ( EXPORT FORWARD/analysis2/AliAODForwardMult.h 
              FORWARD/analysis2/AliAODCentralMult.h 
-             FORWARD/analysis2/AliForwardUtil.h )
+             FORWARD/analysis2/AliForwardUtil.h  )
 
 set ( DHDR  PWG2forward2LinkDef.h)
 
index 227e6e0..fd3df1c 100644 (file)
@@ -1,4 +1,9 @@
-void AddTaskForwardFlow(TString type = "", Int_t etabins = 20)
+void AddTaskForwardFlow(TString type = "", 
+                        Int_t etabins = 40,
+                        Int_t zVertex = 2,
+                        TString addFlow = "",
+                        Int_t addFType = 0,
+                        Int_t addFOrder = 0)
 {
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
   if (!mgr) {
@@ -42,8 +47,12 @@ void AddTaskForwardFlow(TString type = "", Int_t etabins = 20)
 
   qc->SetDoHarmonics(v1, v2, v3, v4);
   qc->SetUseNEtaBins(etabins);
+  qc->AddFlow(addFlow);
+  qc->AddFlowType(addFType);
+  qc->AddFlowOrder(addFOrder);
+  qc->SetVertexRange(zVertex);
   
-  mgr->ConnectInput(qc ,0, mgr->GetCommonInputContainer());
+  mgr->ConnectInput(qc, 0, mgr->GetCommonInputContainer());
   mgr->ConnectOutput(qc, 0, mgr->GetCommonOutputContainer());
   mgr->ConnectOutput(qc, 1, qcout);
 
index 770a529..c6f617e 100644 (file)
@@ -7,10 +7,6 @@
 // Outputs:
 //  - AnalysisResults.root
 //
-// TODO!
-// - Add centrality stuff
-// END OF TODO!
-//
 #include <TROOT.h>
 #include <TSystem.h>
 #include <TInterpreter.h>
 #include <TList.h>
 #include <iostream>
 #include <TMath.h>
+#include "TH3D.h"
 #include "AliLog.h"
 #include "AliForwardFlowTaskQC.h"
 #include "AliAnalysisManager.h"
 #include "AliAODHandler.h"
 #include "AliAODInputHandler.h"
 #include "AliAODMCParticle.h"
-#include "AliForwardFlowBase.h"
 #include "AliAODForwardMult.h"
+#include "AliAODEvent.h"
+
+//
+// Enumeration for adding and retrieving stuff from the histogram
+//
+enum { kW2avg2 = 1, kW2, kW2avg2sq, kW2w2pavg2avg2p, kW2w2p, kW4avg4,
+       kW4, kW4avg4sq, kW2w4, kW2w4avg2avg4, kW2w4p, kW2w4pavg2avg4p, 
+       kW4w2p, kW4w2pavg4avg2p, kW4w4p, kW4w4pavg4avg4p, kQnRe, kQnIm,
+       kM, kCosphi1phi2, kSinphi1phi2, kCosphi1phi2phi3m, kSinphi1phi2phi3m,
+       kMm1m2, kCosphi1phi2phi3p, kSinphi1phi2phi3p,
+       kRW2avg2, kRW2, kRW2avg2sq, kRM, kRW4avg4, kRW4, kRW4avg4sq,
+       kRW2w4, kRW2w4avg2avg4, kRCosphi1phi2, kRSinphi1phi2, 
+       kRCosphi1phi2phi3m, kRSinphi1phi2phi3m, kRMm1m2 };
 
 ClassImp(AliForwardFlowTaskQC)
 #if 0
@@ -34,11 +43,16 @@ ClassImp(AliForwardFlowTaskQC)
 #endif
 
 AliForwardFlowTaskQC::AliForwardFlowTaskQC()
-  : fDebug(0),                 // Debug flag
-    fOutputList(0),    // Output list
+  : fOutputList(0),    // Output list
+    fFlowUtil(0),      // AliForwardFlowUtil
     fAOD(0),           // AOD input event
-    fMC(kFALSE),               // MC flag
-    fEtaBins(20)               // # of eta bins in histograms
+    fMC(kFALSE),       // MC flag
+    fEtaBins(20),      // # of etaBin bins in histograms
+    fAddFlow(0),       // Add flow string
+    fAddType(0),       // Add flow type #
+    fAddOrder(0),      // Add flow order
+    fZvertex(0),       // Z vertex range
+    fCent(0)           // Centrality
 {
   // 
   // Default constructor
@@ -47,11 +61,16 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC()
 //_____________________________________________________________________
 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) :
   AliAnalysisTaskSE(name),
-  fDebug(0),           // Debug flag
   fOutputList(0),      // Output list
+  fFlowUtil(0),                // AliForwardFlowUtil
   fAOD(0),             // AOD input event
   fMC(kFALSE),         // MC flag
-  fEtaBins(20)         // # of Eta bins
+  fEtaBins(20),                // # of Eta bins
+  fAddFlow(0),         // Add flow string
+  fAddType(0),         // Add flow type #
+  fAddOrder(0),                // Add flow order
+  fZvertex(0),         // Z vertex range
+  fCent(0)             // Centrality
 {
   // 
   // Constructor
@@ -65,11 +84,16 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) :
 //_____________________________________________________________________
 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) :
   AliAnalysisTaskSE(o),
-  fDebug(o.fDebug),            // Debug flag
   fOutputList(o.fOutputList),  // Output list
+  fFlowUtil(o.fFlowUtil),      // AliForwardFlowUtil
   fAOD(o.fAOD),                        // AOD input event
   fMC(o.fMC),                  // MC flag
-  fEtaBins(o.fEtaBins)         // # of Eta bins
+  fEtaBins(o.fEtaBins),                // # of Eta bins
+  fAddFlow(o.fAddFlow),                // Add flow string
+  fAddType(o.fAddType),                // Add flow type #
+  fAddOrder(o.fAddOrder),      // Add flow order
+  fZvertex(o.fZvertex),                // Z vertex range
+  fCent(0)                     // Centrality
 {
   // 
   // Copy constructor 
@@ -89,114 +113,126 @@ void AliForwardFlowTaskQC::CreateOutputObjects()
   if (!fOutputList)
     fOutputList = new TList();
   fOutputList->SetName("QCumulants");
+  fOutputList->SetOwner();
+  fFlowUtil = new AliForwardFlowUtil(fOutputList);
+  fFlowUtil->SetVertexRange(fZvertex);
 
   if (fEtaBins % 20) fEtaBins = 20;
 
   // Histograms for cumulants analysis
 
   // We loop over flow histograms here to add different orders of harmonics
-  for (Int_t n = 1; n <= 4; n++) {
-    if (!fv[n]) continue;
-    // Only one flow histogram is needed for each type of data;
-    // x-axis is eta-bins with differential flow, integrated is in underflowbin
-    // y-axis bin 1:  (w_<2> * <2>).Re()
-    // y-axis bin 2:  (w_<2> * <2>).Im()
-    // y-axis bin 3:  w_<2> = M(M-1)
-    // y-axis bin 4:  (w_<2> * <2> * <2>).Re()
-    // y-axis bin 5:  (w_<2> * <2> * <2>).Im()
-    // y-axis bin 6:  (w_<2> * w_<2'> * <2> * <2'>).Re()
-    // y-axis bin 7:  (w_<2> * w_<2'> * <2> * <2'>).Im()
-    // y-axis bin 8:  w_<2> * w_<2'>
-    // y-axis bin 9:  (w_<4> * <4>).Re()
-    // y-axis bin 10:  (w_<4> * <4>).Im()
-    // y-axis bin 11:  w_<4>
-    // y-axis bin 12:  (w_<4> * <4> * <4>).Re()
-    // y-axis bin 13:  (w_<4> * <4> * <4>).Im()
-    // y-axis bin 14:  w_<2> * w_<4>
-    // y-axis bin 15:  (w_<2> * w_<4> * <2> * <4>).Re()
-    // y-axis bin 16:  (w_<2> * w_<4> * <2> * <4>).Im()
-    // y-axis bin 17:  w_<2> * w_<4'>
-    // y-axis bin 18:  (w_<2> * w_<4'> * <2> * <4'>).Re()
-    // y-axis bin 19:  (w_<2> * w_<4'> * <2> * <4'>).Im()
-    // y-axis bin 20:  w_<4> * w_<2'>
-    // y-axis bin 21:  (w_<4> * w_<2'> * <4> * <2'>).Re()
-    // y-axis bin 22:  (w_<4> * w_<2'> * <4> * <2'>).Im()
-    // y-axis bin 23:  w_<4> * w_<4'>
-    // y-axis bin 24:  (w_<4>  * w_<4'> * <4> * <4'>).Re()
-    // y-axis bin 25:  (w_<4>  * w_<4'> * <4> * <4'>).Im()
-    // y-axis bin 26: Qn or pn.Re() = <<cos2phi or psi>>
-    // y-axis bin 27: Qn or pn.Im() = <<sin2phi or psi>>
-    // y-axis bin 28: M or mp
-    // y-axis bin 29: (Qn*Qn-Q2n).Re() = <<cos(2(phi1 or psi1+phi2))>>
-    // y-axis bin 30: (Qn*Qn-Q2n).Im() = <<sin(2(phi1 or psi1+phi2))>>
-    // y-axis bin 31: <<cos(2(phi1 or psi1-phi2-phi3))>>
-    // y-axis bin 32: <<sin(2(phi1 or psi1-phi2-phi3))>>
-    // y-axis bin 33: M*(M-1)*(M-2) or similar for diff
-    // y-axis bin 34: <<cos(2(psi1+phi2-phi3>>
-    // y-axis bin 35: <<sin(2(psi1+phi2-phi3>>
-    TH2D* hFlowHist = new TH2D(Form("hQ%dCumuHist", n), Form("hQ%dCumuHist", n), fEtaBins, -4, 6, 35, 0.5, 35.5);
-    hFlowHist->Sumw2();
-    fOutputList->Add(hFlowHist);
-
-    TH2D* hFlowHistMC = new TH2D(Form("hQ%dCumuHistMC", n), Form("hQ%dCumuHistMC", n), fEtaBins, -4, 6, 35, 0.5, 35.5);
-    hFlowHistMC->Sumw2();
-    fOutputList->Add(hFlowHistMC);
-
-    TH2D* hFlowHistTrRef = new TH2D(Form("hQ%dCumuHistTrRef", n), Form("hQ%dCumuHistTrRef", n), fEtaBins, -4, 6, 35, 0.5, 35.5);
-    hFlowHistTrRef->Sumw2();
-    fOutputList->Add(hFlowHistTrRef);
-
-    // Output histograms
-    TH1D* hCumulant2Flow = new TH1D(Form("hQ%dCumulant2Flow", n), Form("hQ%dCumulant2Flow", n), fEtaBins, -4, 6);
-    hCumulant2Flow->Sumw2();
-    fOutputList->Add(hCumulant2Flow);
-    TH1D* hCumulant2FlowMC = new TH1D(Form("hQ%dCumulant2FlowMC", n),Form("hQ%dCumulant2FlowMC", n), fEtaBins, -4, 6);
-    hCumulant2FlowMC->Sumw2();
-    fOutputList->Add(hCumulant2FlowMC);
-  
-    TH1D* hCumulant2FlowTrRef = new TH1D(Form("hQ%dCumulant2FlowTrRef", n), Form("hQ%dCumulant2FlowTrRef", n), fEtaBins, -4, 6);
-    hCumulant2FlowTrRef->Sumw2();
-    fOutputList->Add(hCumulant2FlowTrRef);
-
-
-    TH1D* hCumulant4Flow = new TH1D(Form("hQ%dCumulant4Flow", n), Form("hQ%dCumulant4Flow", n), fEtaBins, -4, 6);
-    hCumulant4Flow->Sumw2();
-    fOutputList->Add(hCumulant4Flow);
-  
-    TH1D* hCumulant4FlowMC = new TH1D(Form("hQ%dCumulant4FlowMC", n), Form("hQ%dCumulant4FlowMC", n), fEtaBins, -4, 6);
-    hCumulant4FlowMC->Sumw2();
-    fOutputList->Add(hCumulant4FlowMC);
-  
-    TH1D* hCumulant4FlowTrRef = new TH1D(Form("hQ%dCumulant4FlowTrRef", n), Form("hQ%dCumulant4FlowTrRef", n), fEtaBins, -4, 6);
-    hCumulant4FlowTrRef->Sumw2();
-    fOutputList->Add(hCumulant4FlowTrRef);
-  }
-
-  // Single Event histograms
-  TH1D* hdNdphiSE = new TH1D("hdNdphiSE","hdNdphiSE", 20, 0, 2*TMath::Pi());
-  hdNdphiSE->Sumw2();
-  fOutputList->Add(hdNdphiSE);
-
-  TH2D* hdNdetadphiSE = new TH2D("hdNdetadphiSE", "hdNdetadphiSE", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi());
-  hdNdetadphiSE->Sumw2();
-  fOutputList->Add(hdNdetadphiSE);
-
-  TH1D* hdNdphiSEMC = new TH1D("hdNdphiSEMC","hdNdphiSEMC", 20, 0, 2*TMath::Pi());
-  hdNdphiSEMC->Sumw2();
-  fOutputList->Add(hdNdphiSEMC);
+  TString type = "";
+  for (Int_t loop = 1; loop <= 4; loop++) {
+    
+    if (loop == 1) type = "";
+    if (loop == 2) type = "SPD";
+    if (loop == 3) type = "MC";
+    if (loop == 4) type = "TrRef";
 
-  TH2D* hdNdetadphiSEMC = new TH2D("hdNdetadphiSEMC", "hdNdetadphiSEMC", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi());
-  hdNdetadphiSEMC->Sumw2();
-  fOutputList->Add(hdNdetadphiSEMC);
+    for (Int_t n = 1; n <= 4; n++) {
+      if (!fv[n]) continue;
+      // Only one flow histogram is needed for each type of data;
+      // x-axis is etaBin-bins with differential flow, integrated is in underflowbin
+      // z-axis bin 1:  (w_<2> * <2>).Re()
+      // z-axis bin 2:  w_<2> = mqM-mp
+      // z-axis bin 3:  (w_<2> * <2> * <2>).Re()
+      // z-axis bin 4:  (w_<2> * w_<2'> * <2> * <2'>).Re()
+      // z-axis bin 5:  w_<2> * w_<2'>
+      // z-axis bin 6:  (w_<4> * <4>).Re()
+      // z-axis bin 7:  w_<4>
+      // z-axis bin 8:  (w_<4> * <4> * <4>).Re()
+      // z-axis bin 9:  w_<2> * w_<4>
+      // z-axis bin 10:  (w_<2> * w_<4> * <2> * <4>).Re()
+      // z-axis bin 11:  w_<2> * w_<4'>
+      // z-axis bin 12:  (w_<2> * w_<4'> * <2> * <4'>).Re()
+      // z-axis bin 13:  w_<4> * w_<2'>
+      // z-axis bin 14:  (w_<4> * w_<2'> * <4> * <2'>).Re()
+      // z-axis bin 15:  w_<4> * w_<4'>
+      // z-axis bin 16:  (w_<4>  * w_<4'> * <4> * <4'>).Re()
+      // z-axis bin 17: Qn or pn.Re() = <<cos2phi or psi>>
+      // z-axis bin 18: Qn or pn.Im() = <<sin2phi or psi>>
+      // z-axis bin 19: M or mp
+      // z-axis bin 20: (Qn*Qn-Q2n).Re() = <<cos(2(phi1 or psi1+phi2))>>
+      // z-axis bin 21: (Qn*Qn-Q2n).Im() = <<sin(2(phi1 or psi1+phi2))>>
+      // z-axis bin 22: <<cos(2(phi1 or psi1-phi2-phi3))>>
+      // z-axis bin 23: <<sin(2(phi1 or psi1-phi2-phi3))>>
+      // z-axis bin 24: M*(M-1)*(M-2) or similar for diff
+      // z-axis bin 25: <<cos(2(psi1+phi2-phi3>>
+      // z-axis bin 26: <<sin(2(psi1+phi2-phi3>>
+      // z-axis bin 27: ref: W_<2> * <2>
+      // z-axis bin 28: ref: W_<2>
+      // z-axis bin 29: ref: W_<2> * <2> * <2>
+      // z-axis bin 30: ref: Mult
+      // z-axis bin 31: ref: W_<4> * <4>
+      // z-axis bin 32: ref: W_<4>
+      // z-axis bin 33: ref: W_<4> * <4> * <4>
+      // z-axis bin 34: ref: W_<2> * W_<4>
+      // z-axis bin 35: ref: W_<2> * W_<4> * <2> * <4>
+      // z-axis bin 36: ref: <cos(n(phi1phi2))>
+      // z-axis bin 37: ref: <sin(n(phi1phi2))>
+      // z-axis bin 38: ref: <cos(n(phi1-phi2-phi3))>
+      // z-axis bin 39: ref: <sin(n(phi1-phi2-phi3))>
+      // z-axis bin 40: ref: M(M-1)(M-2)
+      Double_t x[41] = { 0. };
+      for (Int_t e = 0; e <=fEtaBins; e++) {
+        x[e] = -4. + e*(10./(Double_t)fEtaBins);
+      }
+//      Double_t x[6] = { 0.0, 1.0, 2.0, 3.0, 4.5, 6.0 };
+      Double_t y[11] = { 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 100. };
+      Double_t z[41] = { 0. };
+      for (Int_t k = 0; k <= 40; k++) {
+        z[k] = 0.5 + k*1.;
+      }
+      TH3D* hFlowHist = new TH3D(Form("hQ%dCumuHist%s", n, type.Data()), 
+                   Form("hQ%dCumuHist%s", n, type.Data()), fEtaBins, x, 10, y, 40, z);
+//      hFlowHist->RebinAxis(40/fEtaBins, hFlowHist->GetXaxis());
+      hFlowHist->Sumw2();
+      fOutputList->Add(hFlowHist);
+      TString tag = TString();
+      for (Int_t c = -2; c < 10; c++) {
+        // Output histograms
+        if (c == -2)      tag = "mb";
+        else if (c == -1) tag = "0_40";
+        else              tag = Form("%d_%d", (Int_t)y[c], (Int_t)y[c+1]);
+    
+        TH1D* hCumulant2RefFlow = new TH1D(Form("hQ%dCumulant2RefFlow%s_%s", n, type.Data(), tag.Data()), 
+                     Form("hQ%dCumulant2RefFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
+        hCumulant2RefFlow->Sumw2();
+        fOutputList->Add(hCumulant2RefFlow);
+   
+        TH1D* hCumulant4RefFlow = new TH1D(Form("hQ%dCumulant4RefFlow%s_%s", n, type.Data(), tag.Data()), 
+                     Form("hQ%dCumulant4RefFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
+        hCumulant4RefFlow->Sumw2();
+        fOutputList->Add(hCumulant4RefFlow);
+        TH1D* hCumulant2DiffFlow = new TH1D(Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), 
+                     Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
+        hCumulant2DiffFlow->Sumw2();
+        fOutputList->Add(hCumulant2DiffFlow);
+   
+        TH1D* hCumulant4DiffFlow = new TH1D(Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), 
+                     Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
+        hCumulant4DiffFlow->Sumw2();
+        fOutputList->Add(hCumulant4DiffFlow);
+      } // end of centrality loop
+
+    } // end of v_{n} loop
+    // Single Event histograms
+    TH1D* hdNdphiSE = new TH1D(Form("hdNdphiSE%s", type.Data()),
+                 Form("hdNdphiSE%s", type.Data()), 20, 0, 2*TMath::Pi());
+    hdNdphiSE->Sumw2();
+    fOutputList->Add(hdNdphiSE);
 
-  TH1D* hdNdphiSETrRef = new TH1D("hdNdphiSETrRef","hdNdphiSETrRef", 20, 0, 2*TMath::Pi());
-  hdNdphiSETrRef->Sumw2();
-  fOutputList->Add(hdNdphiSETrRef);
+    TH2D* hdNdetaBindphiSE = new TH2D(Form("hdNdetaBindphiSE%s", type.Data()), 
+                 Form("hdNdetaBindphiSE%s", type.Data()), fEtaBins, -4, 6, 20, 0, 2*TMath::Pi());
+    hdNdetaBindphiSE->Sumw2();
+    fOutputList->Add(hdNdetaBindphiSE);
+  } // end of type loop
 
-  TH2D* hdNdetadphiSETrRef = new TH2D("hdNdetadphiSETrRef", "hdNdetadphiSETrRef", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi());
-  hdNdetadphiSETrRef->Sumw2();
-  fOutputList->Add(hdNdetadphiSETrRef);
+  TH1D* cent = new TH1D("cent", "cent", 100, 0, 100);
+  fOutputList->Add(cent);
 
   PostData(1, fOutputList);
 }
@@ -214,26 +250,16 @@ void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
   if (!fAOD) return;
 
-  // Load histograms and reset from last event
-  TH1D* dNdphi     = (TH1D*)fOutputList->FindObject("hdNdphiSE");
-  TH2D* dNdetadphi  = (TH2D*)fOutputList->FindObject("hdNdetadphiSE");
-
-  dNdphi->Reset();
-  dNdetadphi->Reset();
 
-  // Initiate FlowCommon and fill histograms
-  AliForwardFlowBase* common = new AliForwardFlowBase(fOutputList);
+  // fill histograms
+  if (!fFlowUtil->LoopAODFMD(fAOD)) return;
 
-  if (!common->LoopAODFMD(fAOD)) return;
-  //  else if (!common->LoopAODSPD(fAOD)) return;
-  //  if (!common->LoopAODFMDandSPD(fAOD)) return;
-
-  // Run analysis
+  // Run analysis on FMD
   for (Int_t n = 1; n <= 4; n++) {
     if (fv[n])
       CumulantsMethod("", n);
   }
-
+  
   // Find out if there's any MC data present
   if (!fMC) {
     TClonesArray* mcArray = 0;
@@ -243,6 +269,12 @@ void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
   if (fMC) 
     ProcessPrimary();
 
+  if (!fFlowUtil->LoopAODSPD(fAOD)) return;
+  for (Int_t n = 1; n <= 4; n++) {
+    if (fv[n])
+      CumulantsMethod("SPD", n);
+  }
+   
 }
 //_____________________________________________________________________
 void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", 
@@ -261,189 +293,252 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "",
   Double_t n = harmonic;
   
   // We get histograms depending on if it's real data or MC truth data
-  TH2D* flowHist   = (TH2D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", harmonic, type.Data()));
-  TH1D* dNdphi     = (TH1D*)fOutputList->FindObject(Form("hdNdphiSE%s", type.Data()));
-  TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject(Form("hdNdetadphiSE%s", type.Data()));
+  TH3D* flowHist   = (TH3D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", harmonic, type.Data()));
+  TH2D* dNdetaBindphi = (TH2D*)fOutputList->FindObject(Form("hdNdetaBindphiSE%s", type.Data()));
+  TH1D* dNdphi = 0;
+  if (!type.Contains("SPD")) dNdphi = (TH1D*)fOutputList->FindObject(Form("hdNdphiSE%s", type.Data()));
+  if ( type.Contains("SPD")) dNdphi = (TH1D*)fOutputList->FindObject("hdNdphiSE");
 
   // We create the objects needed for the analysis
-  Double_t Mult = dNdphi->GetBinContent(0);
+  Double_t mult = dNdphi->GetBinContent(0);
+  if (type.Length() <= 1) fCent = dNdphi->GetBinContent(dNdphi->GetNbinsX()+1);
+
+  TH1D* cent = (TH1D*)fOutputList->FindObject("cent");
+  if (type.Length() <= 1) cent->Fill(fCent);
 
-  Double_t QnRe = 0, Q2nRe = 0, QnIm = 0, Q2nIm = 0;
-  Double_t pnRe = 0, p2nRe = 0, qnRe = 0, qnnRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, qnnIm = 0;
+  Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0;
+  Double_t pnRe = 0, p2nRe = 0, qnRe = 0, q2nRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, q2nIm = 0;
   Double_t avg2 = 0, avg4 = 0, avg2p = 0, avg4p = 0;
   Double_t w2avg2sq = 0, w2pavg2psq = 0, w4avg4sq = 0, w4pavg4psq = 0;
   Double_t w2w2pavg2avg2p = 0, w2w4avg2avg4 = 0, w2pw4pavg2pavg4p = 0;
   Double_t w2w4pavg2avg4p = 0, w4w2pavg4avg2p = 0, w4w4pavg4avg4p = 0;
-  Double_t CosPhi1Phi2 = 0, CosPhi1Phi2Phi3m = 0, CosPhi1Phi2Phi3p = 0;
-  Double_t SinPhi1Phi2 = 0, SinPhi1Phi2Phi3m = 0, SinPhi1Phi2Phi3p = 0;
-  Double_t Phii = 0;
+  Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0, cosPhi1Phi2Phi3p = 0;
+  Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0, sinPhi1Phi2Phi3p = 0;
+  Double_t phi = 0;
   Double_t multi = 0, multp = 0, mp = 0, mq = 0;
-  Double_t W2 = 0, W4 = 0, W2p = 0, W4p = 0;
+  Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
 
   // We loop over the data 1 time!
-  for (Int_t eta = 1; eta <= dNdetadphi->GetNbinsX(); eta++) {
-    // The values for each individual eta bins are reset
+  for (Int_t etaBin = 1; etaBin <= dNdetaBindphi->GetNbinsX(); etaBin++) {
+    // The values for each individual etaBin bins are reset
     mp = 0;
     pnRe = 0;
     p2nRe = 0;
     pnIm = 0;
     p2nIm = 0;
 
-    for (Int_t phii = 1; phii <= dNdphi->GetNbinsX()+1; phii++) {
-      Phii = dNdphi->GetXaxis()->GetBinCenter(phii);
-      multi = dNdphi->GetBinContent(phii);
+    for (Int_t phiN = 1; phiN <= dNdphi->GetNbinsX(); phiN++) {
+      phi = dNdphi->GetXaxis()->GetBinCenter(phiN);
+      multi = dNdphi->GetBinContent(phiN);
 
-      // In the phi loop on the first eta loop the integrated flow
+      // In the phi loop on the first etaBin loop the integrated flow
       // is calculated from the dNdphi histogram
-      if(eta == 1) {
-        QnRe += multi*TMath::Cos(n*Phii);
-        QnIm += multi*TMath::Sin(n*Phii);
-        Q2nRe += multi*TMath::Cos(2.*n*Phii);
-        Q2nIm += multi*TMath::Sin(2.*n*Phii);
+      if(etaBin == 1) {
+        dQnRe += multi*TMath::Cos(n*phi);
+        dQnIm += multi*TMath::Sin(n*phi);
+        dQ2nRe += multi*TMath::Cos(2.*n*phi);
+        dQ2nIm += multi*TMath::Sin(2.*n*phi);
       }
       
-      // For each eta bin the necesarry values for differential flow
+      // For each etaBin bin the necesarry values for differential flow
       // is calculated. Here is the loop over the phi's.
-      multp = dNdetadphi->GetBinContent(eta, phii);
+      multp = dNdetaBindphi->GetBinContent(etaBin, phiN);
       mp += multp;
-      pnRe += multp*TMath::Cos(n*Phii);
-      pnIm += multp*TMath::Sin(n*Phii);
-      p2nRe += multp*TMath::Cos(2.*n*Phii);
-      p2nIm += multp*TMath::Sin(2.*n*Phii);    
+      pnRe += multp*TMath::Cos(n*phi);
+      pnIm += multp*TMath::Sin(n*phi);
+      p2nRe += multp*TMath::Cos(2.*n*phi);
+      p2nIm += multp*TMath::Sin(2.*n*phi);    
     }
 
-    // The integrated flow is calculated
-    if (eta == 1) {
-      Double_t Eta = flowHist->GetXaxis()->GetBinCenter(0);
+    // The reference flow is calculated for the differential flow histograms
+    if (etaBin == 1) {
+      Double_t eta = flowHist->GetXaxis()->GetBinCenter(0);
 
       // 2-particle
-      W2 = Mult * (Mult - 1.);
-      avg2 = QnRe*QnRe + QnIm*QnIm - Mult;
-      avg2 /= W2;
-      w2avg2sq = W2 * avg2 * avg2; 
+      w2 = mult * (mult - 1.);
+      avg2 = dQnRe*dQnRe + dQnIm*dQnIm - mult;
+      avg2 /= w2;
+      w2avg2sq = w2 * avg2 * avg2; 
 
-      flowHist->Fill(Eta, 1, W2 * avg2);
-      flowHist->Fill(Eta, 3, W2);
-      flowHist->Fill(Eta, 4, w2avg2sq);
+      flowHist->Fill(eta, fCent, kRW2avg2, w2 * avg2);
+      flowHist->Fill(eta, fCent, kRW2, w2);
+      flowHist->Fill(eta, fCent, kRW2avg2sq, w2avg2sq);
 
-      flowHist->Fill(Eta, 26, QnRe);
-      flowHist->Fill(Eta, 27, QnIm);
-      flowHist->Fill(Eta, 28, Mult);
+      flowHist->Fill(eta, fCent, kQnRe, dQnRe);
+      flowHist->Fill(eta, fCent, kQnIm, dQnIm);
+      flowHist->Fill(eta, fCent, kRM, mult);
 
       // 4-particle
-      W4 = Mult * (Mult - 1.) * (Mult - 2.) * (Mult - 3.);
-      Double_t real = Q2nRe*QnRe*QnRe - Q2nRe*QnIm*QnIm + 2.*Q2nIm*QnRe*QnIm;
+      w4 = mult * (mult - 1.) * (mult - 2.) * (mult - 3.);
+      Double_t real = dQ2nRe*dQnRe*dQnRe - dQ2nRe*dQnIm*dQnIm + 2.*dQ2nIm*dQnRe*dQnIm;
 
-      avg4 = TMath::Power(QnRe*QnRe + QnIm*QnIm, 2); 
-      avg4 += Q2nRe*Q2nRe + Q2nIm*Q2nIm - 2.*real;
-      avg4 -= 4.*(Mult - 2.)*(QnRe*QnRe + QnIm*QnIm) - 2.*Mult*(Mult - 3.);
+      avg4 = TMath::Power(dQnRe*dQnRe + dQnIm*dQnIm, 2); 
+      avg4 += dQ2nRe*dQ2nRe + dQ2nIm*dQ2nIm - 2.*real;
+      avg4 -= 4.*(mult - 2.)*(dQnRe*dQnRe + dQnIm*dQnIm) - 2.*mult*(mult - 3.);
   
-      avg4 /= W4;
-      w4avg4sq = W4 * avg4 * avg4;
-      w2w4avg2avg4 = W2 * W4 * avg2 * avg4;
-
-      flowHist->Fill(Eta, 9, W4 * avg4);
-      flowHist->Fill(Eta, 11, W4);
-      flowHist->Fill(Eta, 12, w4avg4sq);
-      flowHist->Fill(Eta, 14, W2 * W4);
-      flowHist->Fill(Eta, 15, w2w4avg2avg4);
-
-      CosPhi1Phi2 = QnRe*QnRe - QnIm*QnIm - Q2nRe;
-      SinPhi1Phi2 = 2.*QnRe*QnIm - Q2nIm;
+      avg4 /= w4;
+      w4avg4sq = w4 * avg4 * avg4;
+      w2w4avg2avg4 = w2 * w4 * avg2 * avg4;
+
+      flowHist->Fill(eta, fCent, kRW4avg4, w4 * avg4);
+      flowHist->Fill(eta, fCent, kRW4, w4);
+      flowHist->Fill(eta, fCent, kRW4avg4sq, w4avg4sq);
+      flowHist->Fill(eta, fCent, kRW2w4, w2 * w4);
+      flowHist->Fill(eta, fCent, kRW2w4avg2avg4, w2w4avg2avg4);
+
+      cosPhi1Phi2 = dQnRe*dQnRe - dQnIm*dQnIm - dQ2nRe;
+      sinPhi1Phi2 = 2.*dQnRe*dQnIm - dQ2nIm;
       
-      CosPhi1Phi2Phi3m = TMath::Power(QnRe, 3) + QnRe*QnIm*QnIm; 
-      CosPhi1Phi2Phi3m -= QnRe*Q2nRe + QnIm*Q2nIm + 2.*(Mult - 1.)*QnRe;
-      SinPhi1Phi2Phi3m = -TMath::Power(QnIm, 3) - QnRe*QnRe*QnIm; 
-      SinPhi1Phi2Phi3m -= QnIm*Q2nRe - QnRe*Q2nIm + 2.*(Mult - 1.)*QnIm;
+      cosPhi1Phi2Phi3m = TMath::Power(dQnRe, 3) + dQnRe*dQnIm*dQnIm; 
+      cosPhi1Phi2Phi3m -= dQnRe*dQ2nRe + dQnIm*dQ2nIm + 2.*(mult - 1.)*dQnRe;
+
+      sinPhi1Phi2Phi3m = -TMath::Power(dQnIm, 3) - dQnRe*dQnRe*dQnIm; 
+      sinPhi1Phi2Phi3m -= dQnIm*dQ2nRe - dQnRe*dQ2nIm - 2.*(mult - 1.)*dQnIm;
 
-      flowHist->Fill(Eta, 29, CosPhi1Phi2);
-      flowHist->Fill(Eta, 30, SinPhi1Phi2);
-      flowHist->Fill(Eta, 31, CosPhi1Phi2Phi3m);
-      flowHist->Fill(Eta, 32, SinPhi1Phi2Phi3m);
-      flowHist->Fill(Eta, 33, Mult*(Mult-1.)*(Mult-2.));
+      flowHist->Fill(eta, fCent, kRCosphi1phi2, cosPhi1Phi2);
+      flowHist->Fill(eta, fCent, kRSinphi1phi2, sinPhi1Phi2);
+      flowHist->Fill(eta, fCent, kRCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
+      flowHist->Fill(eta, fCent, kRSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
+      flowHist->Fill(eta, fCent, kRMm1m2, mult*(mult-1.)*(mult-2.));
 
       // Count number of events
-      flowHist->Fill(Eta, 0., 1.);
-    } // end of harmonics loop
+      flowHist->Fill(eta, fCent, 0., 1.);
+    } 
 
-    // Differential flow calculations for each eta bin is done:
+    // Differential flow calculations for each etaBin bin is done:
     if (mp == 0) continue;
-    Double_t Eta = dNdetadphi->GetXaxis()->GetBinCenter(eta);
+    Double_t eta = dNdetaBindphi->GetXaxis()->GetBinCenter(etaBin);
+//    eta = TMath::Abs(eta);
 
     mq = mp;
     qnRe = pnRe;
     qnIm = pnIm;
-    qnnRe = p2nRe;
-    qnnIm = p2nIm;
+    q2nRe = p2nRe;
+    q2nIm = p2nIm;
+    if (type.Contains("SPD") || (type.Contains("") && eta >= 4.)) {
+      mq = 0;
+      qnRe = 0;
+      qnIm = 0;
+      q2nRe = 0;
+      q2nIm = 0;
+    }
+
+    // Then the reference flow is calculated for each etaBin bin also,
+    // TODO: Find smart way to implement in above calculations...
+
+    // 2-particle
+    w2p = mp * (mp - 1.);
+    avg2p = pnRe*pnRe + pnIm*pnIm - mp;
+    avg2p /= w2p;
+    w2pavg2psq = w2p * avg2p * avg2p; 
+
+    flowHist->Fill(eta, fCent, kRW2avg2, w2p * avg2p);
+    flowHist->Fill(eta, fCent, kRW2, w2p);
+    flowHist->Fill(eta, fCent, kRW2avg2sq, w2pavg2psq);
+
+    flowHist->Fill(eta, fCent, kRM, mp);
+
+    // 4-particle
+    w4p = mp * (mp - 1.) * (mp - 2.) * (mp - 3.);
+    Double_t real = p2nRe*pnRe*pnRe - p2nRe*pnIm*pnIm + 2.*p2nIm*pnRe*pnIm;
+
+    avg4p = TMath::Power(pnRe*pnRe + pnIm*pnIm, 2); 
+    avg4p += p2nRe*p2nRe + p2nIm*p2nIm - 2.*real;
+    avg4p -= 4.*(mp - 2.)*(pnRe*pnRe + pnIm*pnIm) - 2.*mp*(mp - 3.);
+  
+    avg4p /= w4p;
+    w4pavg4psq = w4p * avg4p * avg4p;
+    w2w4pavg2avg4p = w2p * w4p * avg2p * avg4p;
+
+    flowHist->Fill(eta, fCent, kRW4avg4, w4p * avg4p);
+    flowHist->Fill(eta, fCent, kRW4, w4p);
+    flowHist->Fill(eta, fCent, kRW4avg4sq, w4pavg4psq);
+    flowHist->Fill(eta, fCent, kRW2w4, w2p * w4p);
+    flowHist->Fill(eta, fCent, kRW2w4avg2avg4, w2w4pavg2avg4p);
+
+    cosPhi1Phi2 = pnRe*pnRe - pnIm*pnIm - p2nRe;
+    sinPhi1Phi2 = 2.*pnRe*pnIm - p2nIm;
+      
+    cosPhi1Phi2Phi3m = TMath::Power(pnRe, 3) + pnRe*pnIm*pnIm; 
+    cosPhi1Phi2Phi3m -= pnRe*p2nRe + pnIm*p2nIm + 2.*(mp - 1.)*pnRe;
+
+    sinPhi1Phi2Phi3m = -TMath::Power(pnIm, 3) - pnRe*pnRe*pnIm; 
+    sinPhi1Phi2Phi3m -= pnIm*p2nRe - pnRe*p2nIm - 2.*(mp - 1.)*pnIm;
+
+    flowHist->Fill(eta, fCent, kRCosphi1phi2, cosPhi1Phi2);
+    flowHist->Fill(eta, fCent, kRSinphi1phi2, sinPhi1Phi2);
+    flowHist->Fill(eta, fCent, kRCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
+    flowHist->Fill(eta, fCent, kRSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
+    flowHist->Fill(eta, fCent, kRMm1m2, mp*(mp-1.)*(mp-2.));
 
     // 2-particle differential flow
-    W2p = mp * Mult - mq;
-    avg2p = pnRe*QnRe + pnIm*QnIm - mq;
-    avg2p /= W2p;
-    w2pavg2psq = W2p * avg2p * avg2p;
-    w2w2pavg2avg2p = W2 * W2p * avg2 * avg2p;
+    w2p = mp * mult - mq;
+    avg2p = pnRe*dQnRe + pnIm*dQnIm - mq;
+    avg2p /= w2p;
+    w2pavg2psq = w2p * avg2p * avg2p;
+    w2w2pavg2avg2p = w2 * w2p * avg2 * avg2p;
     
-    flowHist->Fill(Eta, 1, W2p * avg2p);
-    flowHist->Fill(Eta, 3, W2p);
-    flowHist->Fill(Eta, 4, w2pavg2psq);
-    flowHist->Fill(Eta, 6, w2w2pavg2avg2p);
-    flowHist->Fill(Eta, 8, W2 * W2p);
+    flowHist->Fill(eta, fCent, kW2avg2, w2p * avg2p);
+    flowHist->Fill(eta, fCent, kW2, w2p);
+    flowHist->Fill(eta, fCent, kW2avg2sq, w2pavg2psq);
+    flowHist->Fill(eta, fCent, kW2w2pavg2avg2p, w2w2pavg2avg2p);
+    flowHist->Fill(eta, fCent, kW2w2p, w2 * w2p);
 
-    flowHist->Fill(Eta, 26, pnRe);
-    flowHist->Fill(Eta, 27, pnIm);
-    flowHist->Fill(Eta, 28, mp);
+    flowHist->Fill(eta, fCent, kQnRe, pnRe);
+    flowHist->Fill(eta, fCent, kQnIm, pnIm);
+    flowHist->Fill(eta, fCent, kM, mp);
 
     // 4-particle differential flow
-    W4p = (mp * Mult - 3.*mq)*(Mult - 1.)*(Mult - 2.);
-
-    avg4p =  pnRe*QnRe*(QnRe*QnRe + QnIm*QnIm) + pnIm*QnIm*(QnRe*QnRe + QnIm*QnIm);
-    avg4p -= qnnRe*QnRe*QnRe - qnnRe*QnIm*QnIm + 2.*qnnIm*QnRe*QnIm;
-    avg4p -= pnRe*QnRe*Q2nRe - pnRe*QnIm*Q2nIm + pnIm*QnRe*Q2nIm + pnIm*QnIm*Q2nRe;
-    avg4p -= 2.*Mult*(pnRe*QnRe + pnIm*QnIm);
-
-    avg4p += - 2.*mq*(QnRe*QnRe + QnIm*QnIm) + 7.*(qnRe*QnRe + qnIm*QnIm); 
-    avg4p += - (QnRe*qnRe + QnIm*qnIm) + (qnnRe*Q2nRe + qnnIm*Q2nIm);
-    avg4p += 2.*(pnRe*QnRe + pnIm*QnIm) + 2.*mq*Mult - 6.*mq;
-    avg4p /= W4p;
-
-    w4pavg4psq = W4p * avg4p * avg4p;
-    w2w4pavg2avg4p = W2 * W4p * avg2 * avg4p;
-    w4w2pavg4avg2p = W4 * W2p * avg4 * avg2p;
-    w4w4pavg4avg4p = W4 * W4p * avg4 * avg4p;
-    w2pw4pavg2pavg4p = W2p * W4p * avg2p * avg4p;
-
-    flowHist->Fill(Eta, 9, W4p * avg4p);
-    flowHist->Fill(Eta, 11, W4p);
-    flowHist->Fill(Eta, 12, w4pavg4psq);
-    flowHist->Fill(Eta, 14, W2p * W4p);
-    flowHist->Fill(Eta, 15, w2pw4pavg2pavg4p);
-    flowHist->Fill(Eta, 17, W2 * W4p);
-    flowHist->Fill(Eta, 18, w2w4pavg2avg4p);
-    flowHist->Fill(Eta, 20, W4 * W2p);
-    flowHist->Fill(Eta, 21, w4w2pavg4avg2p);
-    flowHist->Fill(Eta, 23, W4 * W4p);
-    flowHist->Fill(Eta, 24, w4w4pavg4avg4p);
-
-    CosPhi1Phi2 = pnRe*QnRe - pnIm*QnIm - qnnRe;
-    SinPhi1Phi2 = pnRe*QnIm + pnIm*QnRe - qnnIm;
-
-    CosPhi1Phi2Phi3p =  pnRe*(QnRe*QnRe + QnIm*QnIm - Mult);
-    CosPhi1Phi2Phi3p -= qnnRe*QnRe - qnnIm*QnIm + mq*QnRe - 2.*qnRe;
-    SinPhi1Phi2Phi3p =  pnIm*(QnRe*QnRe + QnIm*QnIm - Mult);
-    SinPhi1Phi2Phi3p -= qnnIm*QnRe - qnnRe*QnIm + mq*QnIm - 2.*qnIm;
-
-    CosPhi1Phi2Phi3m =  pnRe*(QnRe*QnRe - QnIm*QnIm) + 2.*pnIm*QnRe*QnIm;
-    CosPhi1Phi2Phi3m -= pnRe*Q2nRe + pnIm*Q2nIm + 2.*mq*QnRe - 2.*qnRe;
-    SinPhi1Phi2Phi3m =  pnIm*(QnRe*QnRe - QnIm*QnIm) - 2.*pnRe*QnRe*QnIm;
-    SinPhi1Phi2Phi3m += - pnIm*Q2nRe + pnRe*Q2nIm + 2.*mq*QnIm - 2.*qnIm;
-
-    flowHist->Fill(Eta, 29, CosPhi1Phi2);
-    flowHist->Fill(Eta, 30, SinPhi1Phi2);
-    flowHist->Fill(Eta, 31, CosPhi1Phi2Phi3m);
-    flowHist->Fill(Eta, 32, SinPhi1Phi2Phi3m);
-    flowHist->Fill(Eta, 33, (mp*Mult-2.*mq)*(Mult-1.));
-    flowHist->Fill(Eta, 34, CosPhi1Phi2Phi3p);
-    flowHist->Fill(Eta, 35, SinPhi1Phi2Phi3p); 
+    w4p = (mp * mult - 3.*mq)*(mult - 1.)*(mult - 2.);
+
+    avg4p =  pnRe*dQnRe*(dQnRe*dQnRe + dQnIm*dQnIm) + pnIm*dQnIm*(dQnRe*dQnRe + dQnIm*dQnIm);
+    avg4p -= q2nRe*dQnRe*dQnRe - q2nRe*dQnIm*dQnIm + 2.*q2nIm*dQnRe*dQnIm;
+    avg4p -= pnRe*dQnRe*dQ2nRe - pnRe*dQnIm*dQ2nIm + pnIm*dQnRe*dQ2nIm + pnIm*dQnIm*dQ2nRe;
+    avg4p -= 2.*mult*(pnRe*dQnRe + pnIm*dQnIm);
+
+    avg4p += - 2.*mq*(dQnRe*dQnRe + dQnIm*dQnIm) + 7.*(qnRe*dQnRe + qnIm*dQnIm); 
+    avg4p += - (dQnRe*qnRe + dQnIm*qnIm) + (q2nRe*dQ2nRe + q2nIm*dQ2nIm);
+    avg4p += 2.*(pnRe*dQnRe + pnIm*dQnIm) + 2.*mq*mult - 6.*mq;
+    avg4p /= w4p;
+
+    w4pavg4psq = w4p * avg4p * avg4p;
+    w2w4pavg2avg4p = w2 * w4p * avg2 * avg4p;
+    w4w2pavg4avg2p = w4 * w2p * avg4 * avg2p;
+    w4w4pavg4avg4p = w4 * w4p * avg4 * avg4p;
+    w2pw4pavg2pavg4p = w2p * w4p * avg2p * avg4p;
+
+    flowHist->Fill(eta, fCent, kW4avg4, w4p * avg4p);
+    flowHist->Fill(eta, fCent, kW4, w4p);
+    flowHist->Fill(eta, fCent, kW4avg4sq, w4pavg4psq);
+    flowHist->Fill(eta, fCent, kW2w4, w2p * w4p);
+    flowHist->Fill(eta, fCent, kW2w4avg2avg4, w2pw4pavg2pavg4p);
+    flowHist->Fill(eta, fCent, kW2w4p, w2 * w4p);
+    flowHist->Fill(eta, fCent, kW2w4pavg2avg4p, w2w4pavg2avg4p);
+    flowHist->Fill(eta, fCent, kW4w2p, w4 * w2p);
+    flowHist->Fill(eta, fCent, kW4w2pavg4avg2p, w4w2pavg4avg2p);
+    flowHist->Fill(eta, fCent, kW4w4p, w4 * w4p);
+    flowHist->Fill(eta, fCent, kW4w4pavg4avg4p, w4w4pavg4avg4p);
+
+    cosPhi1Phi2 = pnRe*dQnRe - pnIm*dQnIm - q2nRe;
+    sinPhi1Phi2 = pnRe*dQnIm + pnIm*dQnRe - q2nIm;
+
+    cosPhi1Phi2Phi3p =  pnRe*(dQnRe*dQnRe + dQnIm*dQnIm - mult);
+    cosPhi1Phi2Phi3p -= q2nRe*dQnRe - q2nIm*dQnIm + mq*dQnRe - 2.*qnRe;
+    sinPhi1Phi2Phi3p =  pnIm*(dQnRe*dQnRe + dQnIm*dQnIm - mult);
+    sinPhi1Phi2Phi3p -= q2nIm*dQnRe - q2nRe*dQnIm + mq*dQnIm - 2.*qnIm;
+
+    cosPhi1Phi2Phi3m =  pnRe*(dQnRe*dQnRe - dQnIm*dQnIm) + 2.*pnIm*dQnRe*dQnIm;
+    cosPhi1Phi2Phi3m -= pnRe*dQ2nRe + pnIm*dQ2nIm + 2.*mq*dQnRe - 2.*qnRe;
+    sinPhi1Phi2Phi3m =  pnIm*(dQnRe*dQnRe - dQnIm*dQnIm) - 2.*pnRe*dQnRe*dQnIm;
+    sinPhi1Phi2Phi3m += - pnIm*dQ2nRe + pnRe*dQ2nIm + 2.*mq*dQnIm - 2.*qnIm;
+
+    flowHist->Fill(eta, fCent, kCosphi1phi2, cosPhi1Phi2);
+    flowHist->Fill(eta, fCent, kSinphi1phi2, sinPhi1Phi2);
+    flowHist->Fill(eta, fCent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
+    flowHist->Fill(eta, fCent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
+    flowHist->Fill(eta, fCent, kMm1m2, (mp*mult-2.*mq)*(mult-1.));
+    flowHist->Fill(eta, fCent, kCosphi1phi2phi3p, cosPhi1Phi2Phi3p);
+    flowHist->Fill(eta, fCent, kSinphi1phi2phi3p, sinPhi1Phi2Phi3p); 
 
   }
 
@@ -458,243 +553,291 @@ void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
   //  Parameters:
   //    option Not used 
   //
+//  Double_t x[6] = { 0., 1., 2., 3., 4.5, 6. };
+
+  TH3D* hcumulantsHist = 0;
+  TH2D* cumulantsHist = new TH2D("tmphist", "tmphist",  fEtaBins, -4, 6, 40, 0.5, 40.5);
+  TH1D* cumulant2refHist = 0;
+  TH1D* cumulant4refHist = 0;
+  TH1D* cumulant2diffHist = 0;
+  TH1D* cumulant4diffHist = 0;
+
+  // For flow calculations
+  Double_t two = 0, qc2 = 0, vTwo2 = 0, four = 0, qc4 = 0, vTwo4 = 0; 
+  Double_t twoPrime = 0, qc2Prime = 0, vTwo2diff = 0, fourPrime = 0, qc4Prime = 0, vTwo4diff = 0;
+  Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0, sqrtW2sq = 0, sqrtW2psq = 0, w2W2p = 0;
+  Double_t w2W4 = 0, w2W4p = 0, w4W2p = 0, w4W4p = 0, w2pW4p = 0;
+  Double_t sqrtW4sq = 0, sqrtW4psq = 0;
+  Double_t w2avg2 = 0, w2pavg2p = 0, w4avg4 = 0, w4pavg4p = 0;
+  Double_t cosP1nPhi = 0, sinP1nPhi = 0, mult = 0, cosP1nPhi1P1nPhi2 = 0, sinP1nPhi1P1nPhi2 = 0;
+  Double_t cosP1nPhi1M1nPhi2M1nPhi3 = 0, sinP1nPhi1M1nPhi2M1nPhi3 = 0, multm1m2 = 0;
+  Double_t cosP1nPsi = 0, sinP1nPsi = 0, mp = 0, cosP1nPsi1P1nPhi2 = 0, sinP1nPsi1P1nPhi2 = 0;
+  Double_t cosP1nPsi1M1nPhi2M1nPhi3 = 0, sinP1nPsi1M1nPhi2M1nPhi3 = 0, mpqMult = 0;
+  Double_t cosP1nPsi1P1nPhi2M1nPhi3 = 0, sinP1nPsi1P1nPhi2M1nPhi3 = 0;
+
+  // For error calculations
+  Double_t w2avg2sq = 0, w2W2pavg2avg2p = 0, w2pavg2psq = 0;
+  Double_t w4avg4sq = 0, w2W4avg2avg4 = 0, w4W4pavg4avg4p = 0, w4pavg4psq = 0;
+  Double_t w2W4pavg2avg4p = 0, w4W2pavg4avg2p = 0;
+  Double_t stwosq = 0, stwoPrimesq = 0, sfoursq = 0, sfourPrimesq = 0;
+  Double_t vTwo2err = 0, vTwo2diffErr = 0, vTwo4err = 0, vTwo4diffErr = 0;
+  Double_t cov22p = 0, cov24 = 0, cov24p = 0, cov42p = 0, cov44p = 0, cov2p2np = 0;
+  Double_t w2pW4pavg2pavg4p = 0;
+  
 
-  TH2D* cumulantsHist;
-  TH1D* cumulant2Hist; 
-  TH1D* cumulant4Hist; 
-
-  Int_t nLoops = (fMC ? 3 : 1);
+  Int_t nLoops = (fMC ? 4 : 2);
+  TString type = "";
 
+  Int_t y[11] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 100 };
+  
   // Do a loop over the difference analysis types, calculating flow
-  // 1 loop for real data, 3 for MC data
+  // 2 loops for real data, 3 for MC data
   // inside each is a nested loop over each harmonic (1, 2, 3 and 4 at the moment)
   for (Int_t loop = 1; loop <= nLoops; loop++) {
 
-    TString type;
     if (loop == 1) type = "";
-    if (loop == 2) type = "MC";
-    if (loop == 3) type = "TrRef";
+    if (loop == 2) type = "SPD";
+    if (loop == 3) type = "MC";
+    if (loop == 4) type = "TrRef";
     
     for (Int_t n = 1; n <= 4; n++) {
       if (!fv[n]) continue;
 
-      cumulantsHist = (TH2D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", n, type.Data()));
-      cumulant2Hist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2Flow%s", n, type.Data()));
-      cumulant4Hist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4Flow%s", n, type.Data()));
-  
-      // For flow calculations
-      Double_t Avg2 = 0, c2 = 0, vTwo2 = 0, Avg4 = 0, c4 = 0, vTwo4 = 0; 
-      Double_t Avg2p = 0, d2 = 0, vTwo2diff = 0, Avg4p = 0, d4 = 0, vTwo4diff = 0;
-      Double_t W2 = 0, W4 = 0, W2p = 0, W4p = 0, sqrtW2sq = 0, sqrtW2psq = 0, W2W2p = 0;
-      Double_t W2W4 = 0, W2W4p = 0, W4W2p = 0, W4W4p = 0, W2pW4p = 0;
-      Double_t sqrtW4sq = 0, sqrtW4psq = 0;
-      Double_t W2avg2 = 0, W2pavg2p = 0, W4avg4 = 0, W4pavg4p = 0;
-      Double_t AvgCos2Phi = 0, AvgSin2Phi = 0, Mult = 0, AvgCos2Phi1Phi2 = 0, AvgSin2Phi1Phi2 = 0;
-      Double_t AvgCos2Phi1Phi2Phi3m = 0, AvgSin2Phi1Phi2Phi3m = 0, Multm1m2 = 0;
-      Double_t AvgCos2Psi = 0, AvgSin2Psi = 0, mp = 0, AvgCos2Psi1Phi2 = 0, AvgSin2Psi1Phi2 = 0;
-      Double_t AvgCos2Psi1Phi2Phi3m = 0, AvgSin2Psi1Phi2Phi3m = 0, mpqMult = 0;
-      Double_t AvgCos2Psi1Phi2Phi3p = 0, AvgSin2Psi1Phi2Phi3p = 0;
-
-      // For error calculations
-      Double_t W2avg2sq = 0, W2W2pavg2avg2p = 0, W2pavg2psq = 0;
-      Double_t W4avg4sq = 0, W2W4avg2avg4 = 0, W4W4pavg4avg4p = 0, W4pavg4psq = 0;
-      Double_t W2W4pavg2avg4p = 0, W4W2pavg4avg2p = 0;
-      Double_t sAvg2sq = 0, sAvg2psq = 0, sAvg4sq = 0, sAvg4psq = 0;
-      Double_t vTwo2err = 0, vTwo2diffErr = 0, vTwo4err = 0, vTwo4diffErr = 0;
-      Double_t Cov22p = 0, Cov24 = 0, Cov24p = 0, Cov42p = 0, Cov44p = 0, Cov2p2np = 0;
-      Double_t W2pW4pavg2pavg4p = 0;
-
-      for (Int_t eta = 0; eta <= cumulantsHist->GetNbinsX(); eta++) {
-        if (eta == 0) {
-          // 2-particle reference flow
-          W2avg2 = cumulantsHist->GetBinContent(eta, 1);
-          if (!W2avg2) continue;
-          W2 = cumulantsHist->GetBinContent(eta, 3);
-          AvgCos2Phi = cumulantsHist->GetBinContent(eta, 26);
-          AvgSin2Phi = cumulantsHist->GetBinContent(eta, 27);
-          Mult = cumulantsHist->GetBinContent(eta, 28);
-          AvgCos2Phi /= Mult;
-          AvgSin2Phi /= Mult;
-          Avg2 = W2avg2 / W2;
-          c2 = Avg2 - TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2); 
-          vTwo2 = TMath::Sqrt(c2);
-          cumulant2Hist->SetBinContent(eta, vTwo2);
-
-          // 4-particle reference flow
-          W4avg4 = cumulantsHist->GetBinContent(eta, 9);
-          W4 = cumulantsHist->GetBinContent(eta, 11);
-          AvgCos2Phi1Phi2 = cumulantsHist->GetBinContent(eta, 29);
-          AvgSin2Phi1Phi2 = cumulantsHist->GetBinContent(eta, 30);
-          AvgCos2Phi1Phi2 /= W2;
-          AvgSin2Phi1Phi2 /= W2;
-          AvgCos2Phi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 31);
-          AvgSin2Phi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 32);
-          Multm1m2 = cumulantsHist->GetBinContent(eta, 33);
-          AvgCos2Phi1Phi2Phi3m /= Multm1m2;
-          AvgSin2Phi1Phi2Phi3m /= Multm1m2;
-          Avg4 = W4avg4 / W4;
-          c4 = Avg4 - 2. * Avg2 * Avg2;
-          c4 -= 4.*AvgCos2Phi*AvgCos2Phi1Phi2Phi3m;
-          c4 += 4.*AvgSin2Phi*AvgSin2Phi1Phi2Phi3m; 
-          c4 -= TMath::Power(AvgCos2Phi1Phi2, 2) + TMath::Power(AvgSin2Phi1Phi2 , 2);
-          c4 += 4.*AvgCos2Phi1Phi2*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2));
-          c4 += 8.*AvgSin2Phi1Phi2*AvgSin2Phi*AvgCos2Phi;  
-          c4 += 8.*Avg2*(TMath::Power(AvgCos2Phi, 2) + TMath::Power(AvgSin2Phi, 2));
-          c4 -= 6.*TMath::Power(TMath::Power(AvgCos2Phi, 2)+TMath::Power(AvgSin2Phi, 2), 2);
-
-          vTwo4 = TMath::Power(-c4, 0.25);
-          cumulant4Hist->SetBinContent(eta, vTwo4);
+      hcumulantsHist = (TH3D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", n, type.Data()));
  
-          // 2-particle reference flow error calculations
-          W2avg2sq = cumulantsHist->GetBinContent(eta, 4);
-          sqrtW2sq = cumulantsHist->GetBinError(eta, 3);
-  
-          sAvg2sq = VarSQ(W2avg2sq, Avg2, W2, W2avg2, sqrtW2sq);
-          vTwo2err = sqrtW2sq * TMath::Sqrt(sAvg2sq) / (2. * TMath::Sqrt(Avg2) * W2);
-          cumulant2Hist->SetBinError(eta, vTwo2err);
-  
-          // 4-particle reference flow error calculations
-          W4avg4sq = cumulantsHist->GetBinContent(eta, 12);
-          sqrtW4sq = cumulantsHist->GetBinError(eta, 11);
-          W2W4 = cumulantsHist->GetBinContent(eta, 14);
-          W2W4avg2avg4 = cumulantsHist->GetBinContent(eta, 15);
-  
-          sAvg4sq = VarSQ(W4avg4sq, Avg4, W4, W4avg4, sqrtW4sq);
-          Cov24 = CovXY(W2W4avg2avg4, W2W4, Avg2*Avg4, W2, W4);
-  
-          vTwo4err = Avg2*Avg2 * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2);
-          vTwo4err += TMath::Power(sqrtW4sq, 2) * sAvg4sq / (16. * W4*W4);
-          vTwo4err -= Avg2 * W2W4 * Cov24 / (2. * W2 * W4);
-          vTwo4err /= TMath::Power(2. * Avg2*Avg2 - Avg4, 1.5);
-          vTwo4err = TMath::Sqrt(vTwo4err);
-          cumulant4Hist->SetBinError(eta, vTwo4err);
-          continue;
+      // Centrality loop
+      for (Int_t c = -2; c < hcumulantsHist->GetNbinsY(); c++) {
+        if (c == -2) {
+          hcumulantsHist->GetYaxis()->SetRange(0, 9);
+          cumulantsHist = (TH2D*)hcumulantsHist->Project3D("zx oe");
+          
+          cumulant2refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2RefFlow%s_mb", n, type.Data()));
+          cumulant4refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4RefFlow%s_mb", n, type.Data()));
+          cumulant2diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2DiffFlow%s_mb", n, type.Data()));
+          cumulant4diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4DiffFlow%s_mb", n, type.Data()));
+        }
+        else if (c == -1) {
+          hcumulantsHist->GetYaxis()->SetRange(1, 5);
+          cumulantsHist = (TH2D*)hcumulantsHist->Project3D("zx oe");
+          
+          cumulant2refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2RefFlow%s_0_40", n, type.Data()));
+          cumulant4refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4RefFlow%s_0_40", n, type.Data()));
+          cumulant2diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2DiffFlow%s_0_40", n, type.Data()));
+          cumulant4diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4DiffFlow%s_0_40", n, type.Data()));
+        }
+        else {
+          hcumulantsHist->GetYaxis()->SetRange(c+1, c+1);
+          cumulantsHist = (TH2D*)hcumulantsHist->Project3D("zx e");
+          cumulant2refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2RefFlow%s_%d_%d", n, type.Data(), y[c], y[c+1]));
+          cumulant4refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4RefFlow%s_%d_%d", n, type.Data(), y[c], y[c+1]));
+          cumulant2diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2DiffFlow%s_%d_%d", n, type.Data(), y[c], y[c+1]));
+          cumulant4diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4DiffFlow%s_%d_%d", n, type.Data(), y[c], y[c+1]));
         }
 
-        // 2-particle differential flow
-        W2pavg2p = cumulantsHist->GetBinContent(eta, 1);
-        if (!W2pavg2p) continue;
-        W2p = cumulantsHist->GetBinContent(eta, 3);
-        AvgCos2Psi = cumulantsHist->GetBinContent(eta, 26);
-        AvgSin2Psi = cumulantsHist->GetBinContent(eta, 27);
-        mp = cumulantsHist->GetBinContent(eta, 28);
-        AvgCos2Psi /= mp;
-        AvgSin2Psi /= mp;
-        Avg2p = W2pavg2p / W2p;
-        d2 = Avg2p - AvgCos2Psi*AvgCos2Phi - AvgSin2Psi*AvgSin2Phi; 
-        vTwo2diff = d2 / TMath::Sqrt(c2);
-        cumulant2Hist->SetBinContent(eta, vTwo2diff);
-        // 4-particle differential flow
-        W4pavg4p = cumulantsHist->GetBinContent(eta, 9);
-        W4p = cumulantsHist->GetBinContent(eta, 11);
-        AvgCos2Psi1Phi2 = cumulantsHist->GetBinContent(eta, 29);
-        AvgSin2Psi1Phi2 = cumulantsHist->GetBinContent(eta, 30);
-        AvgCos2Psi1Phi2 /= W2p;
-        AvgSin2Psi1Phi2 /= W2p;
-        AvgCos2Psi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 31);
-        AvgSin2Psi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 32);
-        mpqMult = cumulantsHist->GetBinContent(eta, 33);
-        AvgCos2Psi1Phi2Phi3m /= mpqMult;
-        AvgSin2Psi1Phi2Phi3m /= mpqMult;
-        AvgCos2Psi1Phi2Phi3p = cumulantsHist->GetBinContent(eta, 34);
-        AvgSin2Psi1Phi2Phi3p = cumulantsHist->GetBinContent(eta, 35);
-        AvgCos2Psi1Phi2Phi3p /= mpqMult;
-        AvgSin2Psi1Phi2Phi3p /= mpqMult;
-
-        Avg4p = W4pavg4p / W4p;
-        d4 = Avg4p - 2. * Avg2p * Avg2;
-        d4 -= AvgCos2Psi*AvgCos2Phi1Phi2Phi3m; 
-        d4 += AvgSin2Psi*AvgSin2Phi1Phi2Phi3m; 
-        d4 -= AvgCos2Phi*AvgCos2Psi1Phi2Phi3m; 
-        d4 += AvgSin2Phi*AvgSin2Psi1Phi2Phi3m; 
-        d4 -= 2.*AvgCos2Phi*AvgCos2Psi1Phi2Phi3p;
-        d4 -= 2.*AvgSin2Phi*AvgSin2Psi1Phi2Phi3p; 
-        d4 -= AvgCos2Psi1Phi2*AvgCos2Phi1Phi2; 
-        d4 -= AvgSin2Psi1Phi2*AvgSin2Phi1Phi2; 
-        d4 += 2.*AvgCos2Phi1Phi2*(AvgCos2Psi*AvgCos2Phi - AvgSin2Psi*AvgSin2Phi);  
-        d4 += 2.*AvgSin2Phi1Phi2*(AvgCos2Psi*AvgSin2Phi + AvgSin2Psi*AvgCos2Phi); 
-        d4 += 4.*Avg2*(AvgCos2Psi*AvgCos2Phi + AvgSin2Psi*AvgSin2Phi);
-        d4 += 2.*AvgCos2Psi1Phi2*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2)); 
-        d4 += 4.*AvgSin2Psi1Phi2*AvgCos2Phi*AvgSin2Phi;
-        d4 += 4.*Avg2p*(TMath::Power(AvgCos2Phi, 2) + TMath::Power(AvgSin2Phi, 2)); 
-        d4 -= 6.*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2))
-         *(AvgCos2Psi*AvgCos2Phi-AvgSin2Psi*AvgSin2Phi); 
-        d4 -= 12.*AvgCos2Phi*AvgSin2Phi*(AvgSin2Psi*AvgCos2Phi+AvgCos2Psi*AvgSin2Phi); 
-        vTwo4diff = - d4 / TMath::Power(-c4, 0.75);
-        cumulant4Hist->SetBinContent(eta, vTwo4diff);    
-      
-        // 2-particle differential flow error calculations
-        W2pavg2psq = cumulantsHist->GetBinContent(eta, 4);
-        sqrtW2psq = cumulantsHist->GetBinError(eta, 3);
-        W2W2pavg2avg2p = cumulantsHist->GetBinContent(eta, 6);
-        W2W2p = cumulantsHist->GetBinContent(eta, 8);
+        Bool_t refDone = kFALSE;
+
+        for (Int_t etaBin = 1; etaBin <= cumulantsHist->GetNbinsX(); etaBin++) {
+          if (refDone == kFALSE || etaBin == 0) {
+            // 2-particle reference flow
+            w2avg2 = cumulantsHist->GetBinContent(etaBin, kRW2avg2);
+            if (etaBin == cumulantsHist->GetNbinsX()) {
+              refDone = kTRUE;
+              etaBin = -1;
+              continue;
+            }
+            if (!w2avg2) continue;
+            w2 = cumulantsHist->GetBinContent(etaBin, kRW2);
+            cosP1nPhi = cumulantsHist->GetBinContent(etaBin, kQnRe);
+            sinP1nPhi = cumulantsHist->GetBinContent(etaBin, kQnIm);
+            mult = cumulantsHist->GetBinContent(etaBin, kRM);
+            cosP1nPhi /= mult;
+            sinP1nPhi /= mult;
+            two = w2avg2 / w2;
+            qc2 = two  - TMath::Power(cosP1nPhi, 2) - TMath::Power(sinP1nPhi, 2); 
+            vTwo2 = TMath::Sqrt(qc2);
+            if (etaBin == 0) cumulant2diffHist->SetBinContent(etaBin, vTwo2);
+            if (etaBin > 0)  cumulant2refHist->SetBinContent(etaBin, vTwo2);
+
+            // 4-particle reference flow
+            w4avg4 = cumulantsHist->GetBinContent(etaBin, kRW4avg4);
+            w4 = cumulantsHist->GetBinContent(etaBin, kRW4);
+            cosP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, kRCosphi1phi2);
+            sinP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, kRSinphi1phi2);
+            cosP1nPhi1P1nPhi2 /= w2;
+            sinP1nPhi1P1nPhi2 /= w2;
+            cosP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kRCosphi1phi2phi3m);
+            sinP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kRSinphi1phi2phi3m);
+            multm1m2 = cumulantsHist->GetBinContent(etaBin, kRMm1m2);
+            cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
+            sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
+            four = w4avg4 / w4;
+            qc4 = four-2.*pow(two,2.)
+               - 4.*cosP1nPhi*cosP1nPhi1M1nPhi2M1nPhi3
+               + 4.*sinP1nPhi*sinP1nPhi1M1nPhi2M1nPhi3-pow(cosP1nPhi1P1nPhi2,2.)-pow(sinP1nPhi1P1nPhi2,2.)
+               + 4.*cosP1nPhi1P1nPhi2*(pow(cosP1nPhi,2.)-pow(sinP1nPhi,2.))+8.*sinP1nPhi1P1nPhi2*sinP1nPhi*cosP1nPhi
+               + 8.*two*(pow(cosP1nPhi,2.)+pow(sinP1nPhi,2.))-6.*pow((pow(cosP1nPhi,2.)+pow(sinP1nPhi,2.)),2.);
+
+
+            vTwo4 = TMath::Power(-qc4, 0.25);
+            if (etaBin == 0) cumulant4diffHist->SetBinContent(etaBin, vTwo4);
+            if (etaBin > 0) cumulant4refHist->SetBinContent(etaBin, vTwo4);
+   
+            // 2-particle reference flow error calculations
+            w2avg2sq = cumulantsHist->GetBinContent(etaBin, kRW2avg2sq);
+            sqrtW2sq = cumulantsHist->GetBinError(etaBin, kRW2);
+   
+            stwosq = VarSQ(w2avg2sq, two, w2, w2avg2, sqrtW2sq);
+            vTwo2err = sqrtW2sq * TMath::Sqrt(stwosq) / (2. * TMath::Sqrt(two) * w2);
+            if (etaBin == 0) cumulant2diffHist->SetBinError(etaBin, vTwo2err);
+            if (etaBin > 0)  cumulant2refHist->SetBinError(etaBin, vTwo2err);
+    
+            // 4-particle reference flow error calculations
+            w4avg4sq = cumulantsHist->GetBinContent(etaBin, kRW4avg4sq);
+            sqrtW4sq = cumulantsHist->GetBinError(etaBin, kRW4);
+            w2W4 = cumulantsHist->GetBinContent(etaBin, kRW2w4);
+            w2W4avg2avg4 = cumulantsHist->GetBinContent(etaBin, kRW2w4avg2avg4);
+   
+            sfoursq = VarSQ(w4avg4sq, four, w4, w4avg4, sqrtW4sq);
+            cov24 = CovXY(w2W4avg2avg4, w2W4, two*four, w2, w4);
   
-        Cov22p = CovXY(W2W2pavg2avg2p, W2W2p, Avg2*Avg2p, W2, W2p);
-        sAvg2psq = VarSQ(W2pavg2psq, Avg2p, W2p, W2pavg2p, sqrtW2psq);
+            vTwo4err = two*two * TMath::Power(sqrtW2sq, 2) * stwosq / (w2*w2);
+            vTwo4err += TMath::Power(sqrtW4sq, 2) * sfoursq / (16. * w4*w4);
+            vTwo4err -= two * w2W4 * cov24 / (2. * w2 * w4);
+            vTwo4err /= TMath::Power(2. * two*two - four, 1.5);
+            vTwo4err = TMath::Sqrt(vTwo4err);
+            if (etaBin == 0) cumulant4diffHist->SetBinError(etaBin, vTwo4err);
+            if (etaBin > 0)  cumulant4refHist->SetBinError(etaBin, vTwo4err);
+
+            continue;
+          }
+
+          // 2-particle differential flow
+          w2pavg2p = cumulantsHist->GetBinContent(etaBin, kW2avg2);
+          if (!w2pavg2p) continue;
+          w2p = cumulantsHist->GetBinContent(etaBin, kW2);
+          cosP1nPsi = cumulantsHist->GetBinContent(etaBin, kQnRe);
+          sinP1nPsi = cumulantsHist->GetBinContent(etaBin, kQnIm);
+          mp = cumulantsHist->GetBinContent(etaBin, kM);
+          cosP1nPsi /= mp;
+          sinP1nPsi /= mp;
+          twoPrime = w2pavg2p / w2p;
+          qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi;
+          vTwo2diff = qc2Prime / TMath::Sqrt(qc2);
+          cumulant2diffHist->SetBinContent(etaBin, vTwo2diff);
  
-        vTwo2diffErr = Avg2p*Avg2p*TMath::Power(sqrtW2psq, 2)*sAvg2sq/(W2*W2);
-        vTwo2diffErr += 4.*Avg2*Avg2*TMath::Power(sqrtW2psq, 2)*sAvg2psq/(W2p*W2p);
-        vTwo2diffErr -= 4.*Avg2*Avg2p*W2W2p*Cov22p/(W2*W2p);
-        vTwo2diffErr /= (4. * TMath::Power(Avg2, 3));
-        vTwo2diffErr = TMath::Sqrt(vTwo2diffErr);
-        cumulant2Hist->SetBinError(eta, vTwo2diffErr);
-
-        // 4-particle differential flow error calculations
-        sqrtW4psq = cumulantsHist->GetBinError(eta, 11);
-        W4pavg4psq = cumulantsHist->GetBinContent(eta, 12);
-        W2pW4p = cumulantsHist->GetBinContent(eta, 14);
-        W2pW4pavg2pavg4p = cumulantsHist->GetBinContent(eta, 15);
-        W2W4p = cumulantsHist->GetBinContent(eta, 17);
-        W2W4pavg2avg4p = cumulantsHist->GetBinContent(eta, 18);
-        W4W2p = cumulantsHist->GetBinContent(eta, 20);
-        W4W2pavg4avg2p = cumulantsHist->GetBinContent(eta, 21);
-        W4W4p = cumulantsHist->GetBinContent(eta, 23);
-        W4W4pavg4avg4p = cumulantsHist->GetBinContent(eta, 24);
-      
-        sAvg4psq = VarSQ(W4pavg4psq, Avg4p, W4p, W4pavg4p, sqrtW4psq);
-        Cov24p = CovXY(W2W4pavg2avg4p, W2W4p, Avg2*Avg4p, W2, W4p);
-        Cov42p = CovXY(W4W2pavg4avg2p, W4W2p, Avg4*Avg2p, W4, W2p);
-        Cov44p = CovXY(W4W4pavg4avg4p, W4W4p, Avg4*Avg4p, W4, W4p);
-        Cov2p2np = CovXY(W2pW4pavg2pavg4p, W2pW4p, Avg2p*Avg4p, W2p, W4p);
+          // 4-particle differential flow
+          w4pavg4p = cumulantsHist->GetBinContent(etaBin, kW4avg4);
+          w4p = cumulantsHist->GetBinContent(etaBin, kW4);
+          cosP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, kCosphi1phi2);
+          sinP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, kSinphi1phi2);
+          cosP1nPsi1P1nPhi2 /= w2p;
+          sinP1nPsi1P1nPhi2 /= w2p;
+          cosP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kCosphi1phi2phi3m);
+          sinP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kSinphi1phi2phi3m);
+          mpqMult = cumulantsHist->GetBinContent(etaBin, kMm1m2);
+          cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
+          sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
+          cosP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kCosphi1phi2phi3p);
+          sinP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kSinphi1phi2phi3p);
+          cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
+          sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
+  
+          fourPrime = w4pavg4p / w4p;
+          qc4Prime = fourPrime-2.*twoPrime*two
+                    - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3
+                    + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3
+                    - cosP1nPhi*cosP1nPsi1M1nPhi2M1nPhi3
+                    + sinP1nPhi*sinP1nPsi1M1nPhi2M1nPhi3
+                    - 2.*cosP1nPhi*cosP1nPsi1P1nPhi2M1nPhi3
+                    - 2.*sinP1nPhi*sinP1nPsi1P1nPhi2M1nPhi3
+                    - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2
+                    - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2
+                    + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi)
+                    + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhi+sinP1nPsi*cosP1nPhi)
+                    + 4.*two*(cosP1nPsi*cosP1nPhi+sinP1nPsi*sinP1nPhi)
+                    + 2.*cosP1nPsi1P1nPhi2*(pow(cosP1nPhi,2.)-pow(sinP1nPhi,2.))
+                    + 4.*sinP1nPsi1P1nPhi2*cosP1nPhi*sinP1nPhi
+                    + 4.*twoPrime*(pow(cosP1nPhi,2.)+pow(sinP1nPhi,2.))
+                    - 6.*(pow(cosP1nPhi,2.)-pow(sinP1nPhi,2.)) 
+                    * (cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi)
+                    - 12.*cosP1nPhi*sinP1nPhi
+                    * (sinP1nPsi*cosP1nPhi+cosP1nPsi*sinP1nPhi);
  
-        // Numbers on the side reference term number in paper (cite needed) loosely
-       /*1*/ vTwo4diffErr =  TMath::Power(2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p, 2) 
-         * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2);
-       /*2*/ vTwo4diffErr += 9. * TMath::Power(2.*Avg2*Avg2p - Avg4p, 2) * TMath::Power(sqrtW4sq, 2)
-         * sAvg4sq / (16. * W4*W4);
-       /*3*/ vTwo4diffErr += 4. * Avg2*Avg2 * TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW2psq, 2)
-         * sAvg2psq / (W2p*W2p);
-       /*4*/ vTwo4diffErr += TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW4psq, 2) * sAvg4psq
-         / (W4p*W4p);
-       /*5*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2p - Avg4p) * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
-         * W2W4 * Cov24 / (W2*W4);
-       /*6*/ vTwo4diffErr -= 4. * Avg2 * (2.*Avg2*Avg2 - Avg4) 
-         * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
-         * W2W2p * Cov22p / (W2 * W2p);
-       /*7*/ vTwo4diffErr += 2. * (2.*Avg2*Avg2 - Avg4)
-         * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
-         * W2W4p * Cov24p / (W2 * W4p);
-       /*8*/ vTwo4diffErr += 3.*Avg2*(2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p)
-         * W4W2p * Cov42p / (W4*W2p);
-       /*9*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p)
-         * W4W4p * Cov44p / (W4 * W4p);
-       /*10*/vTwo4diffErr -= 4.*Avg2*TMath::Power(2.*Avg2*Avg2 - Avg4, 2)
-         * W2pW4p * Cov2p2np / (W2p * W4p);
-       /*11*/vTwo4diffErr /= TMath::Power(2.*Avg2*Avg2 - Avg4, 3.5);
-       vTwo4diffErr = TMath::Sqrt(vTwo4diffErr);
-
-        cumulant4Hist->SetBinError(eta, vTwo4diffErr);
-      } // End of eta loop
-
-      // Number of events:
-      Int_t nEv = (Int_t)cumulantsHist->GetBinContent(0,0);
-      cumulant2Hist->SetBinContent(cumulant2Hist->GetNbinsX() + 1, nEv);
-      cumulant4Hist->SetBinContent(cumulant4Hist->GetNbinsX() + 1, nEv);
+          vTwo4diff = - qc4Prime / TMath::Power(-qc4, 0.75);
+          cumulant4diffHist->SetBinContent(etaBin, vTwo4diff);    
+       
+          // 2-particle differential flow error calculations
+          w2pavg2psq = cumulantsHist->GetBinContent(etaBin, kW2avg2sq);
+          sqrtW2psq = cumulantsHist->GetBinError(etaBin, kW2);
+          w2W2pavg2avg2p = cumulantsHist->GetBinContent(etaBin, kW2w2pavg2avg2p);
+          w2W2p = cumulantsHist->GetBinContent(etaBin, kW2w2p);
+    
+          cov22p = CovXY(w2W2pavg2avg2p, w2W2p, two*twoPrime, w2, w2p);
+          stwoPrimesq = VarSQ(w2pavg2psq, twoPrime, w2p, w2pavg2p, sqrtW2psq);
+   
+          vTwo2diffErr = twoPrime*twoPrime*TMath::Power(sqrtW2sq, 2)*stwosq/(w2*w2);
+          vTwo2diffErr += 4.*two*two*TMath::Power(sqrtW2psq, 2)*stwoPrimesq/(w2p*w2p);
+          vTwo2diffErr -= 4.*two*twoPrime*w2W2p*cov22p/(w2*w2p);
+          vTwo2diffErr /= (4. * TMath::Power(two, 3));
+          vTwo2diffErr = TMath::Sqrt(vTwo2diffErr);
+          cumulant2diffHist->SetBinError(etaBin, vTwo2diffErr);
+  
+          // 4-particle differential flow error calculations
+          sqrtW4psq = cumulantsHist->GetBinError(etaBin, kW4);
+          w4pavg4psq = cumulantsHist->GetBinContent(etaBin, kW4avg4sq);
+          w2pW4p = cumulantsHist->GetBinContent(etaBin, kW2w4);
+          w2pW4pavg2pavg4p = cumulantsHist->GetBinContent(etaBin, kW2w4avg2avg4);
+          w2W4p = cumulantsHist->GetBinContent(etaBin, kW2w4p);
+          w2W4pavg2avg4p = cumulantsHist->GetBinContent(etaBin, kW2w4pavg2avg4p);
+          w4W2p = cumulantsHist->GetBinContent(etaBin, kW4w2p);
+          w4W2pavg4avg2p = cumulantsHist->GetBinContent(etaBin, kW4w2pavg4avg2p);
+          w4W4p = cumulantsHist->GetBinContent(etaBin, kW4w4p);
+          w4W4pavg4avg4p = cumulantsHist->GetBinContent(etaBin, kW4w4pavg4avg4p);
+        
+          sfourPrimesq = VarSQ(w4pavg4psq, fourPrime, w4p, w4pavg4p, sqrtW4psq);
+          cov24p = CovXY(w2W4pavg2avg4p, w2W4p, two*fourPrime, w2, w4p);
+          cov42p = CovXY(w4W2pavg4avg2p, w4W2p, four*twoPrime, w4, w2p);
+          cov44p = CovXY(w4W4pavg4avg4p, w4W4p, four*fourPrime, w4, w4p);
+          cov2p2np = CovXY(w2pW4pavg2pavg4p, w2pW4p, twoPrime*fourPrime, w2p, w4p);
+   
+          // Numbers on the side reference term number in paper (cite needed) loosely
+    /*1*/ vTwo4diffErr =  TMath::Power(2.*two*two*twoPrime - 3.*two*fourPrime + 2.*four*twoPrime, 2) 
+                          * TMath::Power(sqrtW2sq, 2) * stwosq / (w2*w2);
+    /*2*/ vTwo4diffErr += 9. * TMath::Power(2.*two*twoPrime - fourPrime, 2) * TMath::Power(sqrtW4sq, 2)
+                          * sfoursq / (16. * w4*w4);
+    /*3*/ vTwo4diffErr += 4. * two*two * TMath::Power(2.*two*two - four, 2) * TMath::Power(sqrtW2psq, 2)
+                          * stwoPrimesq / (w2p*w2p);
+    /*4*/ vTwo4diffErr += TMath::Power(2.*two*two - four, 2) * TMath::Power(sqrtW4psq, 2) * sfourPrimesq
+                          / (w4p*w4p);
+    /*5*/ vTwo4diffErr -= 1.5 * (2.*two*twoPrime - fourPrime) * (2.*two*two*twoPrime - 3.*two*fourPrime + 2.*four*twoPrime)
+                          * w2W4 * cov24 / (w2*w4);
+    /*6*/ vTwo4diffErr -= 4. * two * (2.*two*two - four) 
+                          * (2.*two*two*twoPrime - 3.*two*fourPrime + 2.*four*twoPrime)
+                          * w2W2p * cov22p / (w2 * w2p);
+    /*7*/ vTwo4diffErr += 2. * (2.*two*two - four)
+                          * (2.*two*two*twoPrime - 3.*two*fourPrime + 2.*four*twoPrime)
+                          * w2W4p * cov24p / (w2 * w4p);
+    /*8*/ vTwo4diffErr += 3.*two*(2.*two*two - four)*(2.*two*twoPrime - fourPrime)
+                          * w4W2p * cov42p / (w4*w2p);
+    /*9*/ vTwo4diffErr -= 1.5 * (2.*two*two - four)*(2.*two*twoPrime - fourPrime)
+                          * w4W4p * cov44p / (w4 * w4p);
+    /*10*/vTwo4diffErr -= 4.*two*TMath::Power(2.*two*two - four, 2)
+                          * w2pW4p * cov2p2np / (w2p * w4p);
+    /*11*/vTwo4diffErr /= TMath::Power(2.*two*two - four, 3.5);
+          vTwo4diffErr = TMath::Sqrt(vTwo4diffErr);
+  
+          cumulant4diffHist->SetBinError(etaBin, vTwo4diffErr);
+        } // End of etaBin loop
+        // Number of events:
+        Int_t nEv = (Int_t)cumulantsHist->GetBinContent(0,0);
+        cumulant2diffHist->SetBinContent(cumulant2diffHist->GetNbinsX() + 1, nEv);
+        cumulant4diffHist->SetBinContent(cumulant4diffHist->GetNbinsX() + 1, nEv);
+      } // End of centrality loop
     } // End of harmonics loop
   }  // End of type loop
+
+  delete cumulantsHist;
+
 }
 // _____________________________________________________________________
 void AliForwardFlowTaskQC::ProcessPrimary() 
@@ -705,45 +848,29 @@ void AliForwardFlowTaskQC::ProcessPrimary()
   // can be run on it.
   //
 
-  // Histograms are loaded and reset
-  TH1D* dNdphi          = (TH1D*)fOutputList->FindObject("hdNdphiSEMC");
-  TH2D* dNdetadphi      = (TH2D*)fOutputList->FindObject("hdNdetadphiSEMC");
-  TH1D* dNdphiTrRef     = (TH1D*)fOutputList->FindObject("hdNdphiSETrRef");
-  TH2D* dNdetadphiTrRef = (TH2D*)fOutputList->FindObject("hdNdetadphiSETrRef");
-
-  dNdphi->Reset();
-  dNdetadphi->Reset();
-  dNdphiTrRef->Reset();
-  dNdetadphiTrRef->Reset();
-
-  // Loads AliFMDFlowCommon and fills histograms and runs analysis.
-  // AOD events also get a TrackRef histogram
-  AliForwardFlowBase* common = new AliForwardFlowBase(fOutputList);
-
-  if (fAOD) {
-    if (!common->LoopAODmc(fAOD)) return;
-    if (!common->LoopAODtrrefHits(fAOD)) return;
-    //    if (!common->LoopMCaddptFlow(fAOD)) return;
-    //    if (!common->LoopMCaddpdgFlow(fAOD)) return;
-    //    if (!common->LoopMCaddetaFlow(fAOD)) return;
+  if (fFlowUtil->LoopAODtrrefHits(fAOD)) {
+   // Run analysis on TrackRefs
+    for (Int_t n = 1; n <= 4; n++) {
+      if (fv[n])
+        CumulantsMethod("TrRef", n);
+    }
   }
 
-  // Run analysis on MC truth
-  for (Int_t n = 1; n <= 4; n++) {
-    if (fv[n])
-      CumulantsMethod("MC", n);
-  }
-  
-  // Run analysis on TrackRefs
-  for (Int_t n = 1; n <= 4; n++) {
-    if (fv[n])
-      CumulantsMethod("TrRef", n);
+  if (fFlowUtil->LoopAODmc(fAOD, fAddFlow, fAddType, fAddOrder)) {
+    // Run analysis on MC truth
+    for (Int_t n = 1; n <= 4; n++) {
+      if (fv[n])
+        CumulantsMethod("MC", n);
+    }
   }
 
-}
+ }
 //_____________________________________________________________________
-Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, Double_t x, Double_t wx, 
-                                    Double_t wxx, Double_t sqrtwx2)
+Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, 
+                                    Double_t x, 
+                                    Double_t wx, 
+                                    Double_t wxx, 
+                                    Double_t sqrtwx2) const
 {
   //
   // Small function to compute the variance squared - used by Terminte()
@@ -757,20 +884,23 @@ Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, Double_t x, Double_t wx,
   return sx;
 }
 //_____________________________________________________________________
-Double_t AliForwardFlowTaskQC::CovXY(Double_t wxwyxy, Double_t wxwy, 
-                                    Double_t xy, Double_t wx, Double_t wy)
+Double_t AliForwardFlowTaskQC::CovXY(Double_t wxwyxy, 
+                                    Double_t wxwy, 
+                                    Double_t xy, 
+                                    Double_t wx, 
+                                    Double_t wy) const
 {
   //
   // Small function to compute the covariance between two numbers
   // - used by Terminate()
   //
-  Double_t Cov, denominator, numerator;
+  Double_t cov, denominator, numerator;
 
   denominator = (wxwyxy / wxwy) - xy;
   numerator = 1 - (wxwy / (wx * wy));
   
-  Cov = denominator / numerator;
-  return Cov;
+  cov = denominator / numerator;
+  return cov;
 }
 //_____________________________________________________________________
 //
index 70532ac..dcee800 100644 (file)
@@ -4,7 +4,8 @@
 #ifndef ALIFORWARDFLOWTASKQC_H
 #define ALIFORWARDFLOWTASKQC_H
 #include "AliAnalysisTaskSE.h"
-#include "AliAODEvent.h"
+#include "AliForwardFlowUtil.h"
+class AliAODEvent;
 
 /**
  * Calculate the flow in the forward regions using the Q cumulants method
@@ -95,10 +96,25 @@ public:
   void SetDoHarmonics(Bool_t v1 = kTRUE, Bool_t v2 = kTRUE, 
                      Bool_t v3 = kTRUE, Bool_t v4 = kTRUE) { 
     fv[1] = v1; fv[2] = v2; fv[3] = v3; fv[4] = v4; }
+  /*
+   * Set string to add flow to MC truth particles
+   */
+  void AddFlow(TString type = "") { fAddFlow = type; }
+  /*
+   * Set which function fAddFlow should use
+   */
+  void AddFlowType(Int_t number = 0) { fAddType = number; }
+  /*
+   * Set which order of flow to add
+   */
+  void AddFlowOrder(Int_t order = 2) { fAddOrder = order; }
+ /**
+   * Set Z vertex range - Used by flow task
+   */
+  void SetVertexRange(Int_t vertex = 2) { fZvertex = vertex; }
   /** 
    * @} 
    */
-  virtual void SetDebugLevel(Int_t level) { fDebug = level; }
   
 protected:
   /*
@@ -123,22 +139,28 @@ protected:
    *
    */
   Double_t VarSQ(Double_t wxx2, Double_t x, Double_t wx, 
-                Double_t wxx, Double_t sqrtwx2);
+                Double_t wxx, Double_t sqrtwx2) const ;
   /**
    * Caclulate the covariance between x and y - used to finalize
    * calculations in Terminate()
    *
    */
   Double_t CovXY(Double_t wxwyxy, Double_t wxwy, Double_t XY, 
-                Double_t wx, Double_t wy);
+                Double_t wx, Double_t wy) const;
 
-  Int_t          fDebug;        //  Debug flag
   TList*         fOutputList;   //  Output list
+  AliForwardFlowUtil* fFlowUtil;//  AliForwardFlowUtil
   AliAODEvent*   fAOD;          //  AOD event
   Bool_t         fMC;           //  Is MC flags
   Int_t          fEtaBins;      //  Number of eta bins in histograms
   Bool_t         fv[5];         //  Calculate v_{n} flag
-   
+  TString        fAddFlow;     //  Add flow string
+  Int_t          fAddType;     //  Add flow type #
+  Int_t          fAddOrder;    //  Add flow order
+  Int_t                 fZvertex;      //  Z vertex range
+  Double_t       fCent;         //  Centrality
+  
+
   ClassDef(AliForwardFlowTaskQC, 1); // Analysis task for FMD analysis
 };
  
diff --git a/PWG2/FORWARD/analysis2/AliForwardFlowUtil.cxx b/PWG2/FORWARD/analysis2/AliForwardFlowUtil.cxx
new file mode 100644 (file)
index 0000000..1425dd0
--- /dev/null
@@ -0,0 +1,363 @@
+//
+// Class used to handle the input from AODs and put it into histograms
+// the Forward Flow tasks can run on.
+// It can also add flow to AliAODMCParticles. 
+//
+#include <iostream>
+#include "AliForwardFlowUtil.h"
+#include "AliAODCentralMult.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "AliLog.h"
+#include "AliAODForwardMult.h"
+#include "AliAODEvent.h"
+
+ClassImp(AliForwardFlowUtil)
+
+//_____________________________________________________________________
+AliForwardFlowUtil::AliForwardFlowUtil() : 
+  fList(0),
+  fZvertex(0) {} 
+  //
+  // Default Constructor
+  //
+//_____________________________________________________________________
+AliForwardFlowUtil::AliForwardFlowUtil(TList* outputList) :
+  fList(0),
+  fZvertex(0)
+{ 
+  // 
+  // Constructor
+  //
+  // Parameters:
+  // TList: list containing histograms for flow analysis
+  //
+  fList = outputList;
+}
+//_____________________________________________________________________
+Bool_t AliForwardFlowUtil::AODCheck(const AliAODForwardMult* aodfm) const
+{
+  // 
+  // Function to check that and AOD event meets the cuts
+  //
+  // Parameters: 
+  //  AliAODForwardMult: forward mult object with trigger and vertex info
+  //
+  if (!aodfm->IsTriggerBits(AliAODForwardMult::kInel)) return kFALSE;
+  if (!aodfm->HasIpZ()) return kFALSE;
+  if (!aodfm->InRange(-fZvertex,fZvertex)) return kFALSE;
+
+  return kTRUE;
+}
+//_____________________________________________________________________
+Bool_t AliForwardFlowUtil::LoopAODFMD(const AliAODEvent* aodevent) const
+{
+  //
+  // Loop over AliAODFowardMult object and fill histograms provided
+  // by flow task.
+  //
+  AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(aodevent->FindListObject("Forward"));
+  if (!aodfmult) return kFALSE;
+  if (!AODCheck(aodfmult)) return kFALSE;
+
+  TH2D hdNdetadphi = aodfmult->GetHistogram();
+  TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSE");
+  TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSE");
+  dNdphi->Reset();
+  dNdetadphi->Reset();
+
+  Double_t mult = 0;
+  Double_t eta = 0;
+  Double_t phi = 0;
+  Double_t weight = 0;
+  for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
+    eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
+    for (Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
+      phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
+      weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
+      dNdetadphi->Fill(eta, phi, weight);
+      if (eta < 4.) dNdphi->Fill(phi, weight);
+      mult += weight;
+    }
+  }
+  dNdphi->SetBinContent(0, mult);
+//  if (aodfmult->HasCentrality()) dNdphi->SetBinContent(dNdphi->GetNbinsX()+1, aodfmult->GetCentrality());
+//  else dNdphi->SetBinContent(dNdphi->GetNbinsX()+1, -1);
+
+  dNdphi->SetBinContent(dNdphi->GetNbinsX()+1, GetCentFromMC(aodevent));
+  return kTRUE;
+}
+//_____________________________________________________________________
+Bool_t AliForwardFlowUtil::LoopAODSPD(const AliAODEvent* aodevent) const
+{
+  // 
+  // Loop over AliAODCentralMult object and fill histograms
+  // provided by flow task
+  //
+  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(aodevent->FindListObject("CentralClusters"));
+  if (!aodcmult) return kFALSE;
+
+  TH2D hdNdetadphi = aodcmult->GetHistogram();
+  TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSESPD");
+  dNdetadphi->Reset();
+
+  Double_t eta = 0;
+  Double_t phi = 0;
+  Double_t weight = 0;
+  for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
+    eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
+    for (Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
+      phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
+      weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
+      dNdetadphi->Fill(eta, phi, weight);
+    }
+  }
+  return kTRUE;
+}
+//_____________________________________________________________________
+Bool_t AliForwardFlowUtil::LoopAODtrrefHits(const AliAODEvent* aodevent) const 
+{
+  //
+  // Loop over AliAODForwardMult object, get MC track ref information
+  // and fill flow histograms
+  //
+
+  AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(aodevent->FindListObject("ForwardMC"));
+  if (!aodfmult) return kFALSE;
+ // if (!AODCheck(aodfmult)) return kFALSE;
+
+  TH2D hdNdetadphi = aodfmult->GetHistogram();
+  TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSETrRef");
+  TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSETrRef");
+  dNdphi->Reset();
+  dNdetadphi->Reset();
+
+  Double_t mult = 0;
+  Double_t eta = 0;
+  Double_t phi = 0;
+  Double_t weight = 0;
+  for(Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
+    eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
+    for(Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
+      phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
+      weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
+      dNdetadphi->Fill(eta, phi, weight);
+      dNdphi->Fill(phi, weight);
+      mult += weight;
+    }
+  }
+  dNdphi->SetBinContent(0, mult);
+
+  return kTRUE;
+}
+//_____________________________________________________________________
+Bool_t AliForwardFlowUtil::LoopAODmc(const AliAODEvent* aodevent, 
+                                     TString addFlow = "", 
+                                     Int_t type = 0, 
+                                     Int_t order = 2) const 
+{
+  // 
+  // Loop over AliAODParticle branch and fill flow histograms
+  //
+
+  //retreive MC particles from event
+  TClonesArray* mcArray = (TClonesArray*)aodevent->FindListObject(AliAODMCParticle::StdBranchName());
+  if(!mcArray){
+    AliWarning("No MC array found in AOD. Try making it again.");
+    return kFALSE;
+  }
+
+  Double_t rp = 0;
+  if (addFlow.Length() > 1) {
+    AliAODMCHeader* header = (AliAODMCHeader*)aodevent->FindListObject(AliAODMCHeader::StdBranchName());
+    rp = header->GetReactionPlaneAngle();
+  }
+
+  TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEMC");
+  TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSEMC");
+  dNdphi->Reset();
+  dNdetadphi->Reset();
+  
+  Int_t ntracks = mcArray->GetEntriesFast();
+
+  // Track loop: chek how many particles will be accepted
+  Double_t mult = 0;
+  Double_t weight = 0;
+  for (Int_t it = 0; it < ntracks; it++) {
+    weight = 0;
+    AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
+    if (!particle) {
+      AliError(Form("Could not receive track %d", it));
+      continue;
+    }
+    if (!particle->IsPhysicalPrimary()) continue;
+    if (particle->Charge() == 0) continue;
+    if (particle->Eta() > -3.4 && particle->Eta() < 5) {
+      // Add flow if it is in the argument
+      if (addFlow.Length() > 1) {
+        if (addFlow.Contains("pt"))
+          weight += AddptFlow(particle->Pt(), type);
+        if (addFlow.Contains("pid"))
+          weight += AddpidFlow(particle->PdgCode(), type);
+        if (addFlow.Contains("eta"))
+          weight += AddetaFlow(particle->Eta(), type);
+        if (addFlow.Contains("flat"))
+          weight = 0.1*type;
+        weight *= 2.*TMath::Cos((Double_t)order*(particle->Phi()-rp)); 
+      }
+      weight += 1;
+      
+      dNdphi->Fill(particle->Phi(), weight);
+      dNdetadphi->Fill(particle->Eta(), particle->Phi(), weight);
+      mult += weight;
+    }
+  }
+
+  dNdphi->SetBinContent(0, mult);
+
+  return kTRUE;
+}
+//_____________________________________________________________________
+Double_t AliForwardFlowUtil::AddptFlow(Double_t Pt = 0, Int_t type = 0) const 
+{
+  //
+  // Add pt dependent flow factor
+  //
+  Double_t weight = 0;
+
+  if (type == 1) weight = 0.125*Pt;
+      
+  if (type == 2) {
+    weight = 0.2*(Pt/2.0);
+    if (Pt > 2.0) weight = 0.2;
+  }
+
+  if (type == 3) {
+      weight = 0.05*(Pt/2.0);
+      if (Pt < 2.0) weight = 0.05;
+  } 
+
+  if (type == 4) {
+      weight = 0.2*(Pt/1.0);
+      if (Pt > 1.0) weight = 0.2;
+  }
+
+  if (type == 5) { 
+      weight = 0.05*(Pt/1.0);
+      if (Pt < 1.0) weight = 0.05;
+  }                                                      
+
+  if (type == 6) {
+      weight = 0.2*(Pt/3.0);
+      if (Pt > 3.0) weight = 0.2;
+  }
+
+  return weight;
+}
+//_____________________________________________________________________
+Double_t AliForwardFlowUtil::AddpidFlow(Int_t ID = 0, Int_t type = 0) const 
+{
+  //
+  // Add pid dependent flow factor 
+  //
+  Double_t weight = 0;
+
+  if (type == 1) {
+    weight = 0.07;
+    if (TMath::Abs(ID) ==  211) // pion flow
+      weight = 0.1;
+    if (TMath::Abs(ID) ==  2212) // proton flow
+      weight = 0.05;
+  }
+
+  if (type == 2) {
+    weight = 0.06;
+    if (TMath::Abs(ID) ==  211) // pion flow
+      weight = 0.1;
+    if (TMath::Abs(ID) ==  2212) // proton flow
+      weight = 0.08;
+  }
+
+  if (type == 3) {
+    weight = 0.05;
+    if (TMath::Abs(ID) ==  211) // pion flow
+      weight = 0.1;
+    if (TMath::Abs(ID) ==  2212) // proton flow
+      weight = 0.07;
+  }
+
+  if (type == 4) {
+    weight = 0.07;
+    if (TMath::Abs(ID) ==  211) // pion flow
+      weight = 0.1;
+    if (TMath::Abs(ID) ==  2212) // proton flow
+      weight = 0.085;
+  }
+
+  return weight;
+}
+//_____________________________________________________________________
+Double_t AliForwardFlowUtil::AddetaFlow(Double_t eta = 0, Int_t type = 0) const 
+{
+  //
+  // Add eta dependent flow factor 
+  //
+  Double_t weight = 0;
+  
+  if (type == 1) weight = 0.03 + 0.07 * (1 - TMath::Abs(eta) / 6.);
+
+  if (type == 2) weight = 0.07 * (1 - TMath::Abs(eta) / 6.);
+  
+  if (type == 3) weight = 0.07 * (1 - TMath::Abs(eta) / 8.);
+  
+  if (type == 4) weight = 0.07 * (1 - TMath::Abs(eta) / 10.);
+  
+  if (type == 5) weight = 0.07 * (1 - TMath::Abs(eta) / 12.); 
+
+  if (type == 6) weight = 0.07 * (1 - TMath::Abs(eta) / 14.); 
+
+  if (type == 7) weight = 0.07 * (1 - TMath::Abs(eta) / 16.); 
+
+  return weight;
+}
+//_____________________________________________________________________
+Double_t AliForwardFlowUtil::GetCentFromMC(const AliAODEvent* aodevent) const
+{
+  //
+  // Get centrality from MC impact parameter.
+  // Values taken from: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
+  //
+  Double_t cent = -1.;
+  Double_t b = -1.;
+  AliAODMCHeader* header = (AliAODMCHeader*)aodevent->FindListObject(AliAODMCHeader::StdBranchName());
+  if (!header) return cent;
+  b = header->GetImpactParameter();
+  
+  if ( 0.00 <= b && b <  3.50)
+    cent = 2.5;
+  if ( 3.50 <= b && b <  4.95)
+    cent = 7.5;
+  if ( 4.95 <= b && b <  6.98)
+    cent = 15.;
+  if ( 6.98 <= b && b <  8.55)
+    cent = 25.;
+  if ( 8.55 <= b && b <  9.88)
+    cent = 35.;
+  if ( 9.88 <= b && b < 11.04)
+    cent = 45.;
+  if (11.04 <= b && b < 12.09)
+    cent = 55.;
+  if (12.09 <= b && b < 13.06)
+    cent = 65.;
+  if (13.06 <= b && b < 13.97)
+    cent = 75.;
+  if (13.97 <= b && b < 19.61)
+    cent = 90.; 
+
+  return cent;
+}
+//_____________________________________________________________________
+//
+//
+// EOF
+
diff --git a/PWG2/FORWARD/analysis2/AliForwardFlowUtil.h b/PWG2/FORWARD/analysis2/AliForwardFlowUtil.h
new file mode 100644 (file)
index 0000000..fecd7f5
--- /dev/null
@@ -0,0 +1,102 @@
+//
+// Class used to handle the input from AODs and put it into histograms
+// the Forward Flow tasks can run on
+//
+#ifndef ALIFORWARDFLOWUTIL_H
+#define ALIFORWARDFLOWUTIL_H
+#include "TNamed.h"
+class AliAODForwardMult;
+class AliAODEvent;
+class TList;
+
+/**
+ * 
+ * Class used to handle the input from AODs and put it into histograms
+ * the Forward Flow tasks can run on.
+ *
+ * @ingroup pwg2_forward_tasks
+ */
+class AliForwardFlowUtil : public TNamed
+{
+public:
+  /**
+   * Constructor
+   */
+  AliForwardFlowUtil();
+  /*
+   * Constructor
+   *
+   * @param fList list of histograms for flow analysis
+   */
+  AliForwardFlowUtil(TList* fList);
+  /*
+   * Copy constructor
+   *
+   * @param o Object to copy from
+   */
+  AliForwardFlowUtil(const AliForwardFlowUtil& o) : TNamed(),
+                                                   fList(o.fList),
+                                                   fZvertex(o.fZvertex) {}
+  /** 
+   * Assignment operator 
+   * 
+   * @param o Object to assign from 
+   * 
+   * @return Reference to this object 
+   */
+  AliForwardFlowUtil& operator=(const AliForwardFlowUtil&) { return *this; }
+  /**
+   * Check that AOD event meet trigger requirements
+   */
+  Bool_t AODCheck(const AliAODForwardMult* aodfm) const;
+  /**
+   * Loop over AliAODForwardMult object and fill flow histograms
+   */
+  Bool_t LoopAODFMD(const AliAODEvent* AODevent) const;
+  /*
+   * Loop over AliAODForwardCentral object and fill flow histograms
+   */
+  Bool_t LoopAODSPD(const AliAODEvent* AODevent) const;
+  /**
+   * Loop over AliAODForwardMult object and fill flow histograms from
+   * track refs
+   */
+  Bool_t LoopAODtrrefHits(const AliAODEvent* AODevent) const;
+ /**
+   * Loop over AliAODMCParticle branch object and fill flow histograms
+   * add flow if arguments are set
+   */
+  Bool_t LoopAODmc(const AliAODEvent* AODevent, TString addFlow, Int_t type, Int_t order) const;
+ /**
+   * Set Z vertex range - Used by flow task
+   */
+  void SetVertexRange(Int_t vertex = 2) { fZvertex = vertex; }
+   
+protected:
+  /**
+   * Add pt dependent flow factor
+   */
+  Double_t AddptFlow(Double_t Pt, Int_t type) const;
+  /**
+   * Add pid dependent flow factor
+   */
+  Double_t AddpidFlow(Int_t ID, Int_t type) const;
+  /**
+   * Add eta dependent flow factor
+   */
+  Double_t AddetaFlow(Double_t Eta, Int_t type) const;
+  /**
+   * Get centrality form MC impact parameter
+   */
+  Double_t GetCentFromMC(const AliAODEvent* AODevent) const;
+
+  TList*       fList;  // List of flow histograms
+  Int_t                fZvertex; // Z vertex range
+
+  ClassDef(AliForwardFlowUtil, 1); 
+};
+#endif
+// Local Variables:
+//   mode: C++ 
+// End:
index 1f401fb..a85d413 100644 (file)
@@ -1,43 +1,62 @@
 /**
- * Script to analyse AOD input for flow 
+ * Script to analyse AOD input for flow
+ * 
+ * Takes either a single (AOD) .root file as input or a .txt
+ * The .txt file is expected to contain the path to the files 
+ * from the current directory or the absolute path.
  * 
  * @par Inputs: 
- * - 
+ *  
  * 
  * @par Outputs: 
  * - 
  * 
  */
-void MakeFlow(TString path      = "", 
-             Bool_t  recursive = true,
-             Int_t   nevents   = 100, 
+void MakeFlow(TString data      = "", 
+             Int_t   nevents   = 0, 
              TString type      = "", 
-             Int_t   etabins)
+             Int_t   etabins   = 40,
+             Int_t   zVertex   = 2,
+             TString addFlow   = "",
+              Int_t   addFType  = 0,
+              Int_t   addFOrder = 0)
 {
-  Bool_t proof = false;
+  Bool_t proof = kFALSE;
 
-  // --- Load libs -------------------------------------------------
+  // --- Load libs ---------------------------------------------------
   gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
 
   // --- Check for proof mode, and possibly upload pars --------------
   if (proof> 0) { 
     gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C");
     if (!LoadPars(proof)) { 
-      Error("MakeAOD", "Failed to load PARs");
+      AliError("MakeFlow", "Failed to load PARs");
       return;
     }
   }
 
-  // --- Add to chain either ESD or AOD ----------------------------
-  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
-  TChain* chain = MakeChain("AOD", path, recursive);
-  // If 0 or less events is select, choose all 
-  if (nevents <= 0) nevents = chain->GetEntries();
+  // --- Add to chain either AOD ------------------------------------
+  if (data.Length() <= 1) {
+    AliError("You didn't add a data file");
+    return;
+  }
+  TChain* chain = new TChain("aodTree");
+
+  if (data.Contains(".txt"))
+    MakeChain(data, chain);
 
-  // --- Initiate the event handlers -------------------------------
+  if (data.Contains(".root")) {
+    if (!TFile::Open(data.Data())) {
+      AliError(Form("AOD file %s not found", data.Data()));
+      return;
+    }
+    chain->Add(data.Data());
+  }
+
+  // --- Initiate the event handlers --------------------------------
   AliAnalysisManager *mgr  = new AliAnalysisManager("Forward Flow", 
-                                                   "Flow in forward region");
-  mgr->SetUseProgressBar(kTRUE);
+                                                   "Flow in the forward region");
+  mgr->SetUseProgressBar(kTRUE, 10);
 
   // --- AOD input handler -------------------------------------------
   AliAODInputHandler *aodInputHandler = new AliAODInputHandler();
@@ -50,7 +69,7 @@ void MakeFlow(TString path      = "",
 
   // --- Add the tasks ---------------------------------------------
   gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskForwardFlow.C");
-  AddTaskForwardFlow(type.Data(), etabins);
+  AddTaskForwardFlow(type, etabins, zVertex, addFlow, addFType, addFOrder);
 
   // --- Run the analysis --------------------------------------------
   TStopwatch t;
@@ -66,6 +85,32 @@ void MakeFlow(TString path      = "",
   t.Stop();
   t.Print();
 }
+//----------------------------------------------------------------
+void MakeChain(TString data = "", TChain* chain = 0)
+{
+  // creates chain of files in a given directory or file containing a list.
+  
+  // Open the input stream
+  ifstream in;
+  in.open(data.Data());
+
+  // Read the input list of files and add them to the chain
+  TString line;
+  while(in.good()) 
+  {
+    in >> line;
+      
+    if (line.Length() == 0)
+      continue;      
+    
+    if (TFile::Open(line))
+      chain->Add(line);
+  }
+
+  in.close();
+
+  return;
+}
 //
 // EOF
 //
index 36f5d47..f9ed7e9 100644 (file)
@@ -71,7 +71,7 @@
 #pragma link C++ class AliCentralCorrSecondaryMap+;
 #pragma link C++ class AliCentralCorrAcceptance+;
 #pragma link C++ class AliCentraldNdetaTask+;
-#pragma link C++ class AliForwardFlowBase+;
+#pragma link C++ class AliForwardFlowUtil+;
 #pragma link C++ class AliForwardFlowTaskQC+;
 
 #else