Mega-commit by Alexander - added Event plane to analysis task, and other flow updates...
authorhansena <hansena@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Feb 2012 12:03:58 +0000 (12:03 +0000)
committerhansena <hansena@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Feb 2012 12:03:58 +0000 (12:03 +0000)
29 files changed:
PWGLF/CMakelibPWGLFforward2.pkg
PWGLF/FORWARD/analysis2/AddTaskForwardFlow.C
PWGLF/FORWARD/analysis2/AddTaskMCParticleFilter.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliAODForwardEP.cxx [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliAODForwardEP.h [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx
PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.cxx [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.h [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.cxx
PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.h
PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx
PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.h
PWGLF/FORWARD/analysis2/AliForwardFlowUtil.cxx [deleted file]
PWGLF/FORWARD/analysis2/AliForwardFlowUtil.h [deleted file]
PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.h [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx
PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.h
PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx
PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h
PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx
PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.h
PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.cxx
PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.h
PWGLF/FORWARD/analysis2/ForwardAODConfig.C
PWGLF/FORWARD/analysis2/MakeFlow.C
PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C
PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C [new file with mode: 0644]
PWGLF/PWGLFforward2LinkDef.h

index 398357d..7589ea8 100644 (file)
@@ -39,6 +39,7 @@ set ( SRCS
   FORWARD/analysis2/AliFMDEnergyFitter.cxx
   FORWARD/analysis2/AliFMDEnergyFitterTask.cxx 
   FORWARD/analysis2/AliFMDEventInspector.cxx
+  FORWARD/analysis2/AliFMDEventPlaneFinder.cxx
   FORWARD/analysis2/AliFMDHistCollector.cxx
   FORWARD/analysis2/AliFMDSharingFilter.cxx
   FORWARD/analysis2/AliFMDMCEventInspector.cxx
@@ -57,13 +58,14 @@ set ( SRCS
   FORWARD/analysis2/AliCentralMCMultiplicityTask.cxx
   FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx
   FORWARD/analysis2/AliAODCentralMult.cxx
+  FORWARD/analysis2/AliAODForwardEP.cxx
   FORWARD/analysis2/AliCentralCorrSecondaryMap.cxx
   FORWARD/analysis2/AliCentralCorrAcceptance.cxx 
   FORWARD/analysis2/AliCentraldNdetaTask.cxx
   FORWARD/analysis2/AliBasedNdetaTask.cxx
   FORWARD/analysis2/AliMCTruthdNdetaTask.cxx
-  FORWARD/analysis2/AliForwardFlowUtil.cxx
   FORWARD/analysis2/AliForwardFlowTaskQC.cxx
+  FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx
   FORWARD/analysis2/AliSPDMCTrackDensity.cxx
   FORWARD/analysis2/AliFMDMultCuts.cxx
   FORWARD/analysis2/AliPoissonCalculator.cxx
@@ -159,6 +161,7 @@ install ( FILES FORWARD/analysis2/AddTaskCentraldNdeta.C
                 FORWARD/analysis2/AddTaskForwardQA.C
                FORWARD/analysis2/AddTaskForwarddNdeta.C
                FORWARD/analysis2/AddTaskMCTruthdNdeta.C
+               FORWARD/analysis2/AddTaskMCParticleFilter.C
                FORWARD/analysis2/DrawdNdeta.C
                FORWARD/analysis2/ForwardAODConfig.C
                FORWARD/analysis2/CentralAODConfig.C
index 5fc568f..9767dbb 100644 (file)
@@ -1,7 +1,7 @@
 /**
  * @file   AddTaskForwardFlow.C
- * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
- * @date   Wed Mar 23 12:14:17 2011
+ * @author Alexander Hansen alexander.hansen@cern.ch 
+ * @date   Wed Sep 07 12:14:17 2011
  * 
  * @brief  
  * 
@@ -25,7 +25,6 @@
  * @ingroup pwglf_forward_flow
  */
 void AddTaskForwardFlow(TString type = "", 
-                        Int_t etabins = 48,
                         Bool_t mc = kFALSE,
                         TString addFlow = "",
                         Int_t addFType = 0,
@@ -64,27 +63,36 @@ void AddTaskForwardFlow(TString type = "",
     if (!type.Contains("6")) v6 = kFALSE;
   }
 
-  // --- Create output containers and find input from fmd task --- //
-
-  TString outputFile = AliAnalysisManager::GetCommonFileName();
-  outputFile += ":FlowResults";
-
-  AliAnalysisDataContainer* qcout = mgr->CreateContainer("QCumulants", TList::Class(), AliAnalysisManager::kOutputContainer, outputFile);
+  // --- Create containers for output --- //
+  AliAnalysisDataContainer* sums = 
+    mgr->CreateContainer("FlowQCSums", TList::Class(), 
+                        AliAnalysisManager::kOutputContainer, 
+                        AliAnalysisManager::GetCommonFileName());
+  AliAnalysisDataContainer* output = 
+    mgr->CreateContainer("FlowQCResults", TList::Class(), 
+                        AliAnalysisManager::kParamContainer, 
+                        AliAnalysisManager::GetCommonFileName());
 
   // --- For the selected flow tasks the input and output is set --- //
   
-  AliForwardFlowTaskQC* qc = new AliForwardFlowTaskQC("QCumulants");
+  AliForwardFlowTaskQC* task = 0;
+  if (mc) 
+    task = new AliForwardMCFlowTaskQC("QCumulants");
+  else
+    task = new AliForwardFlowTaskQC("QCumulants");
+  mgr->AddTask(task); 
 
-  qc->SetDoHarmonics(v1, v2, v3, v4, v5, v6);
-  qc->SetUseNEtaBins(etabins);
-  qc->SetMCinput(mc);
-  qc->AddFlow(addFlow);
-  qc->AddFlowType(addFType);
-  qc->AddFlowOrder(addFOrder);
-  qc->SetRefEtaBins(4);
-  
-  mgr->ConnectInput(qc, 0, mgr->GetCommonInputContainer());
-  mgr->ConnectOutput(qc, 1, qcout);
+  task->SetDoHarmonics(v1, v2, v3, v4, v5, v6);
+  if (mc) {
+    AliForwardMCFlowTaskQC* mcTask = 
+      static_cast<AliForwardMCFlowTaskQC*>task;
+    mcTask->AddFlow(addFlow);
+    mcTask->AddFlowType(addFType);
+    mcTask->AddFlowOrder(addFOrder);
+  }
+  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task, 1, sums);
+  mgr->ConnectOutput(task, 2, output);
 
   return;
 }
diff --git a/PWGLF/FORWARD/analysis2/AddTaskMCParticleFilter.C b/PWGLF/FORWARD/analysis2/AddTaskMCParticleFilter.C
new file mode 100644 (file)
index 0000000..104f72d
--- /dev/null
@@ -0,0 +1,20 @@
+AliAnalysisTask*
+AddTaskMCParticleFilter() 
+{
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("AddTaskFMD", "No analysis manager to connect to.");
+    return NULL;
+  }   
+
+  AliAnalysisTaskMCParticleFilter* mctask = new AliAnalysisTaskMCParticleFilter("mcfilter");
+  mgr->AddTask(mctask);
+  AliAnalysisDataContainer* histOut = 
+    mgr->CreateContainer("mcfilter", TList::Class(), 
+                        AliAnalysisManager::kOutputContainer, "mcParticles.root");
+  mgr->ConnectInput(mctask, 0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(mctask, 0, mgr->GetCommonOutputContainer());
+  mgr->ConnectOutput(mctask, 1, histOut);
+
+  return mctask;
+}
diff --git a/PWGLF/FORWARD/analysis2/AliAODForwardEP.cxx b/PWGLF/FORWARD/analysis2/AliAODForwardEP.cxx
new file mode 100644 (file)
index 0000000..ab50c83
--- /dev/null
@@ -0,0 +1,124 @@
+// 
+// Class that contains results from FMD eventplane calculations
+//
+#include "AliAODForwardEP.h"
+#include <TBrowser.h>
+#include <iostream>
+#include <TMath.h>
+#include <TObjString.h>
+#include <TObjArray.h>
+#include "AliLog.h"
+ClassImp(AliAODForwardEP)
+#ifdef DOXY_INPUT
+; // For Emacs 
+#endif
+
+//____________________________________________________________________
+AliAODForwardEP::AliAODForwardEP()
+  : fIsMC(false),
+    fEpT(-1),
+    fEpA(-1),
+    fEpC(-1),
+    fEp1(-1),
+    fEp2(-1),
+    fHist()
+{
+  // 
+  // Constructor 
+  // 
+}
+
+//____________________________________________________________________
+AliAODForwardEP::AliAODForwardEP(Bool_t isMC) 
+  : fIsMC(isMC),
+    fEpT(-1),
+    fEpA(-1),
+    fEpC(-1),
+    fEp1(-1),
+    fEp2(-1),
+    fHist()
+{
+  // 
+  // Constructor 
+  // 
+  // Parameters: 
+  //  isMC   If set to true this is for MC data (effects branch name)
+  // 
+  fHist.SetXTitle("#eta");
+  fHist.SetDirectory(0);
+  fHist.Sumw2();
+}
+
+//____________________________________________________________________
+void
+AliAODForwardEP::Init(const TAxis& etaAxis)
+{
+  // Initialize the histogram with an eta axis 
+  // 
+  // Parameters: 
+  //   etaAxis       Eta axis to use 
+  // 
+  fHist.SetBins(etaAxis.GetNbins()/10, etaAxis.GetXmin(), etaAxis.GetXmax()); 
+}
+
+//____________________________________________________________________
+void
+AliAODForwardEP::Clear(Option_t* option)
+{
+  // Clear (or reset) internal values 
+  // 
+  // Parameters: 
+  //  option   Passed to TH1::Reset 
+  // 
+  fHist.Reset(option);
+  fEpT = -1;
+  fEpA = -1;
+  fEpC = -1;
+  fEp1 = -1;
+  fEp2 = -1;
+}
+//____________________________________________________________________
+void
+AliAODForwardEP::Browse(TBrowser* /*b*/)
+{
+  // Browse this object 
+  // 
+  // Parameters: 
+  //   b   Browser to use 
+  // TODO: Make nice
+/*  static TObjString ipz;
+  static TObjString trg;
+  static TObjString cnt;
+  static TObjString ncl;
+  ipz = Form("ip_z=%fcm", fIpZ);
+  trg = GetTriggerString(fTriggers);
+  cnt = Form("%+6.1f%%", fCentrality);
+  ncl = Form("%d clusters", fNClusters);
+  b->Add(&fHist);
+  b->Add(&ipz);
+  b->Add(&trg);
+  b->Add(&cnt);
+  b->Add(&ncl);
+*/
+}
+
+//____________________________________________________________________
+void
+AliAODForwardEP::Print(Option_t* option) const
+{
+  // Print this object 
+  // 
+  // Parameters: 
+  //  option   Not used 
+  fHist.Print(option);
+  std::cout << "Total FMD EP: \t" << fEpT <<std::endl; 
+  std::cout << "FMD1+2 EP: \t" << fEpA <<std::endl; 
+  std::cout << "FMD3   EP: \t" << fEpC <<std::endl; 
+  std::cout << "Random Subep 1: \t" << fEp1 <<std::endl; 
+  std::cout << "Random Subep 2: \t" << fEp2 <<std::endl; 
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/PWGLF/FORWARD/analysis2/AliAODForwardEP.h b/PWGLF/FORWARD/analysis2/AliAODForwardEP.h
new file mode 100644 (file)
index 0000000..275aba6
--- /dev/null
@@ -0,0 +1,165 @@
+
+// 
+// Class that contains results from FMD eventplane calculations
+//
+#ifndef ALIAODFORWARDEP_H
+#define ALIAODFORWARDEP_H
+/**
+ * @file AliAODForwardMult.h
+ * @author Alexander Hansen
+ * @date   Tue Feb 14 2012
+ * 
+ * @brief
+ * 
+ * 
+ * @ingroup pwglf_forward
+ */
+#include <TObject.h>
+#include <TH1D.h>
+class TBrowser;
+class TH1D;
+
+class AliAODForwardEP : public TObject
+{
+public:
+  /** 
+   * Default constructor 
+   * 
+   * Used by ROOT I/O sub-system - do not use
+   */
+  AliAODForwardEP();
+  /** 
+   * Constructor 
+   * 
+   * @param isMC Whether this was from MC or not 
+   */
+  AliAODForwardEP(Bool_t isMC);
+  /** 
+   * Destructor 
+   */
+  virtual ~AliAODForwardEP() {} // Destructor 
+  /** 
+   * Initialize 
+   * 
+   * @param etaAxis  Pseudo-rapidity axis
+   */
+  void Init(const TAxis& etaAxis);
+  /** 
+   * Clear all data 
+   * 
+   * @param option  Passed on to TH2::Reset verbatim
+   */
+  void Clear(Option_t* option="");
+  /** 
+   * browse this object 
+   * 
+   * @param b Browser 
+   */
+  void Browse(TBrowser* b);
+  /** 
+   * This is a folder 
+   * 
+   * @return Always true
+   */
+  Bool_t IsFolder() const { return kTRUE; } // Always true 
+  /** 
+   * Print content 
+   * 
+   * @param option Passed verbatim to TH2::Print 
+   */
+  void Print(Option_t* option="") const;
+ /*
+  * Get the name of the AOD object
+  *
+  * @return for now ForwardEP
+  */
+  const Char_t* GetName() const { return (fIsMC ? "ForwardEP" : "ForwardEP"); }
+ /*
+  * Set the overall FMD eventplane
+  * 
+  * @param ep FMD EP
+  */
+  void SetEventplane(Double_t ep) { fEpT = ep; }
+ /* 
+  * Get the overall FMD eventplane
+  *
+  * @return FMD EP
+  */
+  Double_t GetEventplane() { return fEpT; }
+ /*
+  * Set FMD eventplane from A side
+  *
+  * @param epA FMD A side EP
+  */
+  void SetEventplaneA(Double_t epA) { fEpA = epA; }
+ /*
+  * Get FMD eventplane from A side
+  *
+  * @return FMD A side EP
+  */
+  Double_t GetEventplaneA() { return fEpA; }
+ /*
+  * Set FMD eventplane from C side
+  *
+  * @param epC FMD C side EP
+  */
+  void SetEventplaneC(Double_t epC) { fEpC = epC; }
+ /*
+  * Get FMD eventplane from C side
+  *
+  * @return FMD C side EP
+  */
+  Double_t GetEventplaneC() { return fEpC; }
+ /*
+  * Set FMD eventplane 1 using random particles
+  *
+  * @param ep1 FMD EP 1 from random selection
+  */
+  void SetEventplane1(Double_t ep1) { fEp1 = ep1; }
+ /* 
+  * Get FMD eventplane 1 from random particles
+  *
+  * @return FMD EP 1 from random selection
+  */
+  Double_t GetEventplane1() { return fEp1; }
+ /*
+  * Set FMD eventplane 2 using random particles
+  *
+  * @param ep2 FMD EP 2 from random selection
+  */
+  void SetEventplane2(Double_t ep2) { fEp2 = ep2; }
+ /* 
+  * Get FMD eventplane 2 from random particles
+  *
+  * @return FMD EP 2 from random selection
+  */
+  Double_t GetEventplane2() { return fEp2; }
+ /*
+  * Get eta histogram from eta vs. Psi_2 calculatins
+  *
+  * @return Returns eta vs. Psi_2 histogram
+  */
+  const TH1D& GetHistogram() const { return fHist; } // Get histogram 
+ /*
+  * Get eta histogram from eta vs. Psi_2 calculatins
+  *
+  * @return Returns eta vs. Psi_2 histogram
+  */
+  TH1D& GetHistogram() { return fHist; } // Get histogram 
+
+protected: 
+  Bool_t   fIsMC;       // Whether this is from MC 
+  Double_t fEpT;        // Total FMD event plane
+  Double_t fEpA;        // FMD1+2 subevent plane
+  Double_t fEpC;        // FMD3 subevent plane
+  Double_t fEp1;        // Random FMD subevent plane 1
+  Double_t fEp2;        // Random FMD subevent plane 2
+  TH1D     fHist;       // Histogram of subevent planes in eta for this event
+
+  ClassDef(AliAODForwardEP,1); // AOD forward multiplicity 
+};
+
+#endif
+// Local Variables:
+//  mode: C++
+// End:
index ba4d205..38613aa 100644 (file)
@@ -188,7 +188,7 @@ AliCentralMultiplicityTask::GetESDEvent()
 
   fNClusterTracklet = new TH2D("nClusterVsnTracklet", 
                               "Total number of cluster vs number of tracklets",
-                              100, 0, 100, 100, 0, 100);
+                              100, 0, 10000, 100, 0, 10000);
   fNClusterTracklet->SetDirectory(0);
   fNClusterTracklet->SetXTitle("# of free clusters");
   fNClusterTracklet->SetYTitle("# of tracklets");
diff --git a/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.cxx b/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.cxx
new file mode 100644 (file)
index 0000000..31ce9cb
--- /dev/null
@@ -0,0 +1,555 @@
+// This task finds the eventplane
+// using the FMD
+// 
+#include "AliFMDEventPlaneFinder.h"
+#include <TList.h>
+#include <TMath.h>
+#include "AliLog.h"
+#include <TH2D.h>
+#include <TH3D.h>
+#include "TROOT.h"
+#include <iostream>
+#include <iomanip>
+#include <TString.h>
+#include <TFile.h>
+#include "AliOADBContainer.h"
+#include "AliAnalysisManager.h"
+#include "AliVEvent.h"
+#include "AliEventplane.h"
+#include "AliCentrality.h"
+#include "AliAODForwardEP.h"
+#include "AliAODForwardMult.h"
+#include "AliAODEvent.h"
+
+ClassImp(AliFMDEventPlaneFinder)
+#if 0
+; // For Emacs
+#endif 
+
+//____________________________________________________________________
+AliFMDEventPlaneFinder::AliFMDEventPlaneFinder()
+  : TNamed(),
+    fList(0),
+    fEvent(0),
+    fQt(),
+    fQa(),
+    fQc(),
+    fQ1(),
+    fQ2(),
+    fQeta(),
+    fHistFMDEventplane(0),
+    fHistFMDEventplaneA(0),
+    fHistFMDEventplaneC(0),
+    fHistFMDEventplane1(0),
+    fHistFMDEventplane2(0),
+    fHistPhiDist(0),
+    fHistFMDmTPCep(0),
+    fHistFMDvsTPCep(0),
+    fHistFMDmVZEROep(0),
+    fHistFMDvsVZEROep(0),
+    fDebug(0),
+    fOADBFileName(0),
+    fOADBContainer(0),
+    fPhiDist(0),
+    fRunNumber(0),
+    fUsePhiWeights(0)
+{
+  // 
+  // Constructor 
+  //
+   
+}
+
+//____________________________________________________________________
+AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const char* title)
+  : TNamed("fmdEventPlaneFinder", title), 
+    fList(0),
+    fEvent(0),
+    fQt(),
+    fQa(),
+    fQc(),
+    fQ1(),
+    fQ2(),
+    fQeta(),
+    fHistFMDEventplane(0),
+    fHistFMDEventplaneA(0),
+    fHistFMDEventplaneC(0),
+    fHistFMDEventplane1(0),
+    fHistFMDEventplane2(0),
+    fHistPhiDist(0),
+    fHistFMDmTPCep(0),
+    fHistFMDvsTPCep(0),
+    fHistFMDmVZEROep(0),
+    fHistFMDvsVZEROep(0),
+    fDebug(0),
+    fOADBFileName(0),
+    fOADBContainer(0),
+    fPhiDist(0),
+    fRunNumber(0),
+    fUsePhiWeights(1)
+{
+  // 
+  // Constructor 
+  // 
+  // Parameters:
+  //    name Name of object
+  //
+}
+
+//____________________________________________________________________
+AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const 
+                                                AliFMDEventPlaneFinder& o)
+  : TNamed(o),
+    fList(o.fList),
+    fEvent(o.fEvent),
+    fQt(o.fQt),
+    fQa(o.fQa),
+    fQc(o.fQc),
+    fQ1(o.fQ1),
+    fQ2(o.fQ2),
+    fQeta(o.fQeta),
+    fHistFMDEventplane(o.fHistFMDEventplane),
+    fHistFMDEventplaneA(o.fHistFMDEventplaneA),
+    fHistFMDEventplaneC(o.fHistFMDEventplaneC),
+    fHistFMDEventplane1(o.fHistFMDEventplane1),
+    fHistFMDEventplane2(o.fHistFMDEventplane2),
+    fHistPhiDist(o.fHistPhiDist),
+    fHistFMDmTPCep(o.fHistFMDmTPCep),
+    fHistFMDvsTPCep(o.fHistFMDvsTPCep),
+    fHistFMDmVZEROep(o.fHistFMDmVZEROep),
+    fHistFMDvsVZEROep(o.fHistFMDvsVZEROep),
+    fDebug(o.fDebug),
+    fOADBFileName(o.fOADBFileName),
+    fOADBContainer(o.fOADBContainer),
+    fPhiDist(o.fPhiDist),
+    fRunNumber(o.fRunNumber),
+    fUsePhiWeights(o.fUsePhiWeights)
+{
+  // 
+  // Copy constructor 
+  // 
+  // Parameters:
+  //    o Object to copy from 
+  //
+}
+
+//____________________________________________________________________
+AliFMDEventPlaneFinder::~AliFMDEventPlaneFinder()
+{
+  // 
+  // Destructor 
+  //
+}
+
+//____________________________________________________________________
+AliFMDEventPlaneFinder&
+AliFMDEventPlaneFinder::operator=(const AliFMDEventPlaneFinder& o)
+{
+  // 
+  // Assignement operator
+  // 
+  // Parameters:
+  //    o Object to assign from 
+  // 
+  // Return:
+  //    Reference to this object
+  //
+  if (&o == this) return *this; 
+  TNamed::operator=(o);
+
+  fList               = o.fList;
+  fEvent              = o.fEvent;
+  fQt                 = o.fQt;
+  fQa                 = o.fQa;
+  fQc                 = o.fQc;
+  fQ1                 = o.fQ1;
+  fQ2                 = o.fQ2;
+  fQeta               = o.fQeta;
+  fHistFMDEventplane  = o.fHistFMDEventplane;
+  fHistFMDEventplaneA = o.fHistFMDEventplaneA;
+  fHistFMDEventplaneC = o.fHistFMDEventplaneC;
+  fHistFMDEventplane1 = o.fHistFMDEventplane1;
+  fHistFMDEventplane2 = o.fHistFMDEventplane2;
+  fHistPhiDist        = o.fHistPhiDist;
+  fHistFMDmTPCep      = o.fHistFMDmTPCep;
+  fHistFMDvsTPCep     = o.fHistFMDvsTPCep;
+  fHistFMDmVZEROep    = o.fHistFMDmVZEROep;
+  fHistFMDvsVZEROep   = o.fHistFMDvsVZEROep;
+  fDebug              = o.fDebug;
+  fOADBFileName       = o.fOADBFileName;
+  fOADBContainer      = o.fOADBContainer;
+  fPhiDist            = o.fPhiDist;
+  fRunNumber          = o.fRunNumber;
+  fUsePhiWeights      = o.fUsePhiWeights;
+
+  return *this;
+}
+
+//____________________________________________________________________
+void
+AliFMDEventPlaneFinder::Init(const TAxis& etaAxis)
+{
+  // Intialize this sub-algorithm 
+  //
+  // Parameters:
+  //   etaAxis   fmd eta axis binning
+  //
+  fHistFMDEventplane = new TH1D("hFMDEventplane", "FMD eventplane", 
+                                100, 0., TMath::Pi());
+  fHistFMDEventplane->Sumw2();
+  fHistFMDEventplane->SetDirectory(0);
+  fList->Add(fHistFMDEventplane);
+
+  //
+  fHistFMDEventplaneA = new TH1D("hFMDEventplaneA", "FMD eventplaneA", 
+                                100, 0., TMath::Pi());
+  fHistFMDEventplaneA->Sumw2();
+  fHistFMDEventplaneA->SetDirectory(0);
+  fList->Add(fHistFMDEventplaneA);
+
+  //
+  fHistFMDEventplaneC = new TH1D("hFMDEventplaneC", "FMD eventplaneC", 
+                                100, 0., TMath::Pi());
+  fHistFMDEventplaneC->Sumw2();
+  fHistFMDEventplaneC->SetDirectory(0);
+  fList->Add(fHistFMDEventplaneC);
+
+  //
+  fHistFMDEventplane1 = new TH1D("hFMDEventplane1", "FMD eventplane1", 
+                                100, 0., TMath::Pi());
+  fHistFMDEventplane1->Sumw2();
+  fHistFMDEventplane1->SetDirectory(0);
+  fList->Add(fHistFMDEventplane1);
+
+  //
+  fHistFMDEventplane2 = new TH1D("hFMDEventplane2", "FMD eventplane2", 
+                                100, 0., TMath::Pi());
+  fHistFMDEventplane2->Sumw2();
+  fHistFMDEventplane2->SetDirectory(0);
+  fList->Add(fHistFMDEventplane2);
+
+  //
+  fHistPhiDist = new TH2D("hPhiDist", "Phi distribution in FMD",
+       etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(), 
+       20, 0., TMath::TwoPi());
+  fHistPhiDist->Sumw2();
+  fHistPhiDist->SetDirectory(0);
+  fList->Add(fHistPhiDist);
+
+  //
+  fHistFMDmTPCep = new TH1D("hFMDmTPCep", 
+                             "Eventplane from FMD - TPC",
+                             100, -TMath::Pi()/2., TMath::Pi()/2.);
+  fHistFMDmTPCep->Sumw2();
+  fHistFMDmTPCep->SetDirectory(0);
+  fList->Add(fHistFMDmTPCep);
+
+  //
+  fHistFMDvsTPCep = new TH2F("hFMDvsTPCep", 
+                             "Eventplane from FMD vs. TPC",
+                             100, 0., TMath::Pi(),
+                             100, 0., TMath::Pi());
+  fHistFMDvsTPCep->Sumw2();
+  fHistFMDvsTPCep->SetDirectory(0);
+  fList->Add(fHistFMDvsTPCep);
+
+  //
+  fHistFMDmVZEROep = new TH1D("hFMDmVZEROep", 
+                             "Eventplane from FMD - VZERO",
+                             100, -TMath::Pi()/2., TMath::Pi()/2.);
+  fHistFMDmVZEROep->Sumw2();
+  fHistFMDmVZEROep->SetDirectory(0);
+  fList->Add(fHistFMDmVZEROep);
+
+  //
+  fHistFMDvsVZEROep = new TH2F("hFMDvsVZEROep", 
+                             "Eventplane from FMD vs. VZERO",
+                             100, 0., TMath::Pi(),
+                             100, 0., TMath::Pi());
+  fHistFMDvsVZEROep->Sumw2();
+  fHistFMDvsVZEROep->SetDirectory(0);
+  fList->Add(fHistFMDvsVZEROep);
+
+  if (!fUsePhiWeights) return;
+  
+  if (fOADBFileName.Length()==0)
+    fOADBFileName = Form("%s/PWGLF/FORWARD/FMDEVENTPLANE/data/fmdEPoadb.root", 
+                       AliAnalysisManager::GetOADBPath());
+
+  TFile foadb(fOADBFileName);
+  if(!foadb.IsOpen()) {
+    AliError(Form("Cannot open OADB file %s", fOADBFileName.Data()));
+    return;
+  }
+  fOADBContainer = (AliOADBContainer*)foadb.Get("FMDphidist");
+  if (!fOADBContainer) AliError(Form("No OADB container found in %s", fOADBFileName.Data()));
+  foadb.Close();
+
+  return;
+}
+
+//_____________________________________________________________________
+void
+AliFMDEventPlaneFinder::DefineOutput(TList* dir)
+{
+  // 
+  // Output diagnostic histograms to directory 
+  // 
+  // Parameters:
+  //    dir List to write in
+  //  
+  fList = new TList();
+  fList->SetOwner();
+  fList->SetName(GetName());
+  dir->Add(fList);
+
+  return;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDEventPlaneFinder::FindEventplane(AliVEvent* ev,
+                                      AliAODForwardEP& aodep,
+                                       TH2D* h,
+                                       AliForwardUtil::Histos* hists)
+{
+  // 
+  // Do the calculations 
+  // 
+  // Parameters:
+  //    hists    Histogram cache
+  //    ep       calculated eventplane 
+  // 
+  // Return:
+  //    true on successs 
+
+  fQt.Set(0., 0.);
+  fQa.Set(0., 0.);
+  fQc.Set(0., 0.);
+  fQ1.Set(0., 0.);
+  fQ2.Set(0., 0.);
+
+  TH1D epEtaHist = aodep.GetHistogram();
+
+  fEvent = ev;
+  if (hists) {
+    for (UShort_t d=1; d<=3; d++) { 
+      UShort_t nr = (d == 1 ? 1 : 2);
+      for (UShort_t q=0; q<nr; q++) { 
+        Char_t r = (q == 0 ? 'I' : 'O');
+        CalcQVectors(hists->Get(d,r), &epEtaHist);
+      } // for q
+    } // for d
+  }
+  else if (h) {
+    CalcQVectors(h, &epEtaHist);
+  }
+
+  aodep.SetEventplane(CalcEventplane(fQt));
+  aodep.SetEventplaneA(CalcEventplane(fQa));
+  aodep.SetEventplaneC(CalcEventplane(fQc));
+  // aodep.SetEventplane1(CalcEventplane(fQ1));
+  // aodep.SetEventplane2(CalcEventplane(fQ2));
+  
+  FillHists(&aodep);
+
+  return kTRUE;
+}
+
+//_____________________________________________________________________
+void
+AliFMDEventPlaneFinder::CalcQVectors(TH2D* h, TH1D* eHist)
+{
+  //
+  // Calculate the Q vectors
+  //
+  Double_t phi = 0, eta = 0, weight = 0;
+  for (Int_t e = 1; e <= h->GetNbinsX(); e++) {
+    Double_t qx = 0, qy = 0;
+    eta = h->GetXaxis()->GetBinCenter(e);
+    for (Int_t p = 1; p <= h->GetNbinsY(); p++) {
+      phi = h->GetYaxis()->GetBinCenter(p);
+      weight = h->GetBinContent(e, p);
+      if (fUsePhiWeights) weight *= GetPhiWeight(e, p);
+      // fix for FMD1 hole
+      if (e > 168 && p == 14) {
+       h->GetBinContent(e, 4);
+        if (fUsePhiWeights) weight *= GetPhiWeight(e, 4);
+      }
+      fHistPhiDist->Fill(eta, phi, weight);
+      // for underflowbin total Nch/eta
+      fHistPhiDist->Fill(eta, -1., weight);
+      
+      // increment Q vectors
+      qx += weight*TMath::Cos(2.*phi);
+      qy += weight*TMath::Sin(2.*phi);
+    }
+    TVector2 qVec = TVector2(qx, qy);
+    fQt += qVec;
+    if (eta < 0) fQa += qVec;
+    if (eta > 0) fQc += qVec;
+    // TODO: Add random ep increments
+    fQeta += qVec;
+    if (e % 10 == 0) {
+      eHist->Fill(eta, CalcEventplane(fQeta));
+      fQeta.Set(0., 0.);
+    }
+  }
+
+  return;
+}
+
+//_____________________________________________________________________
+Double_t
+AliFMDEventPlaneFinder::CalcEventplane(TVector2 v) const
+{
+  //
+  // Calculate the eventplane
+  //
+  Double_t ep = -1;
+  if (v.Mod() == 0.) return ep;
+
+  ep = v.Phi();
+  ep /= 2.;
+
+  if (fDebug > 0)
+    AliInfo(Form("Eventplane found to be: %f", ep));
+  return ep;
+}
+
+//_____________________________________________________________________
+void 
+AliFMDEventPlaneFinder::FillHists(AliAODForwardEP* fmdEP)
+{
+  //
+  // Fill diagnostics histograms
+  //
+  Double_t fmdEPt = fmdEP->GetEventplane();
+  Double_t fmdEPa = fmdEP->GetEventplaneA();
+  Double_t fmdEPc = fmdEP->GetEventplaneC();
+  // Double_t fmdEP1 = fmdEP->GetEventplane1();
+  // Double_t fmdEP2 = fmdEP->GetEventplane2();
+
+  // FMD hists
+  fHistFMDEventplane->Fill(fmdEPt);
+  fHistFMDEventplaneA->Fill(fmdEPa);
+  fHistFMDEventplaneC->Fill(fmdEPc);
+//  fHistFMDEventplane1->Fill(fmdEP1);
+//  fHistFMDEventplane2->Fill(fmdEP2);
+
+  // TPC comparison
+  AliEventplane* ep = fEvent->GetEventplane();
+  Double_t tpcEP = -1;
+  if (ep) ep->GetEventplane("Q");
+
+  Double_t diffTPC = fmdEPt - tpcEP;
+
+  if (diffTPC <  TMath::Pi()/2.) diffTPC = TMath::Pi() - diffTPC;
+  if (diffTPC >= TMath::Pi()/2.) diffTPC = TMath::Pi() - diffTPC;
+
+  fHistFMDmTPCep->Fill(diffTPC);
+
+//  if (fmdEPt >= TMath::Pi()/2. && fmdEPt - tpcEP >= TMath::Pi()/2.) tpcEP = TMath::Pi()-tpcEP;
+//  if (fmdEPt <  TMath::Pi()/2. && tpcEP - fmdEPt >= TMath::Pi()/2.) tpcEP = TMath::Pi()-tpcEP;
+
+  fHistFMDvsTPCep->Fill(fmdEPt, tpcEP);
+
+  // VZERO comparison
+  Double_t vzeroEP = -1;
+  if (ep) ep->GetEventplane("V0");
+
+  Double_t diffVZERO = fmdEPt - vzeroEP;
+
+  if (diffVZERO <  TMath::Pi()/2.) diffVZERO = TMath::Pi() - diffVZERO;
+  if (diffVZERO >= TMath::Pi()/2.) diffVZERO = TMath::Pi() - diffVZERO;
+
+  fHistFMDmVZEROep->Fill(diffVZERO);
+
+//  if (fmdEPt >= TMath::Pi()/2. && fmdEPt - vzeroEP >= TMath::Pi()/2.) vzeroEP = TMath::Pi()-vzeroEP;
+//  if (fmdEPt <  TMath::Pi()/2. && vzeroEP - fmdEPt >= TMath::Pi()/2.) vzeroEP = TMath::Pi()-vzeroEP;
+
+  fHistFMDvsVZEROep->Fill(fmdEPt, vzeroEP);
+  return;
+}
+
+//_____________________________________________________________________
+Double_t 
+AliFMDEventPlaneFinder::GetPhiWeight(Int_t etaBin, Int_t phiBin) const
+{
+  //
+  // Get phi weight for flattening
+  //
+  Double_t phiWeight = 1;
+  
+  if (fPhiDist) {
+    Double_t nParticles = fPhiDist->GetBinContent(etaBin, 0);
+    Double_t nPhiBins = fPhiDist->GetNbinsY();
+    Double_t phiDistValue = fPhiDist->GetBinContent(etaBin, phiBin);
+    
+    if (phiDistValue > 0) phiWeight = nParticles/nPhiBins/phiDistValue;
+  }
+
+  return phiWeight;
+}
+
+//____________________________________________________________________
+void
+AliFMDEventPlaneFinder::SetRunNumber(Int_t run)
+{
+  // 
+  // Set run number
+  //
+  if (fRunNumber == run) return;
+  
+  fRunNumber = run;
+  if (fUsePhiWeights) GetPhiDist();
+
+  return;
+}
+
+//____________________________________________________________________
+void
+AliFMDEventPlaneFinder::GetPhiDist()
+{
+  //
+  // Get phi dist from OADB
+  //
+  if (!fOADBContainer) return;
+
+  fPhiDist = (TH2D*)fOADBContainer->GetObject(fRunNumber, "Default")->Clone(Form("fPhiDist_%d", fRunNumber));
+  fList->Add(fPhiDist);
+
+  return;
+}
+
+//____________________________________________________________________
+void
+AliFMDEventPlaneFinder::Print(Option_t* /*option*/) const
+{
+  // 
+  // Print information 
+  // 
+  // Parameters:
+  //    option Not used
+  //
+  char ind[gROOT->GetDirLevel()+3];
+  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+  ind[gROOT->GetDirLevel()] = '\0';
+  std::cout << ind << ClassName() << ": " << GetName() << '\n'
+           << std::boolalpha 
+           << ind << "Eventplane finder active!" << '\n'
+           << ind << "Loading OADB object for run number: " << fRunNumber << '\n'
+           << ind << "Using phi weights: " << (fUsePhiWeights ? "true" : "false") << '\n'
+           << std::noboolalpha
+           << std::flush;
+  std::cout << std::endl;
+}
+//____________________________________________________________________
+//
+// EOF
+//
+         
+
+
diff --git a/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.h b/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.h
new file mode 100644 (file)
index 0000000..226f12f
--- /dev/null
@@ -0,0 +1,186 @@
+// This task finds the eventplane
+// using the FMD
+// 
+#ifndef ALIFMDEVENTPLANEFINDER_H 
+#define ALIFMDEVENTPLANEFINDER_H 
+/**
+ * @file AliFMDEventPlaneFinder.h
+ * @author Alexander Hansen
+ * @date   Tue Feb 14 2012
+ * 
+ * @brief
+ * 
+ * 
+ * @ingroup pwglf_forward
+ */
+#include <TNamed.h>
+#include <TVector2.h>
+#include "AliForwardUtil.h"
+class AliVEvent;
+class TH1D;
+class TH2F;
+class TH2D;
+class TString;
+class AliOADBContainer;
+class AliAODForwardEP;
+
+class AliFMDEventPlaneFinder : public TNamed
+{
+public:
+  /** 
+   * Constructor 
+   */
+  AliFMDEventPlaneFinder();
+  /** 
+   * Constructor 
+   * 
+   * @param name Name of object
+   */
+  AliFMDEventPlaneFinder(const char* name);
+  /** 
+   * Copy constructor 
+   * 
+   * @param o Object to copy from 
+   */
+  AliFMDEventPlaneFinder(const AliFMDEventPlaneFinder& o);
+  /** 
+   * Destructor 
+   */
+  virtual ~AliFMDEventPlaneFinder();
+  /** 
+   * Assignement operator
+   * 
+   * @param o Object to assign from 
+   * 
+   * @return Reference to this object
+   */
+  AliFMDEventPlaneFinder& operator=(const AliFMDEventPlaneFinder& o);
+ /** 
+   * Initialize this sub-algorithm
+   * 
+   */
+  virtual void Init(const TAxis& etaAxis);
+  /** 
+   * Do the calculations 
+   * 
+   * @param hists    Histogram cache
+   * 
+   * @return true on successs 
+   */
+  Bool_t FindEventplane(AliVEvent* esd,
+                        AliAODForwardEP& aodEp,
+                        TH2D* h,
+                        AliForwardUtil::Histos* hists);
+  /** 
+   * Output diagnostic histograms to directory 
+   * 
+   * @param dir List to write in
+   */  
+  virtual void DefineOutput(TList* dir);
+  /** 
+   * Print information 
+   * 
+   * @param option Print options 
+   *   - max  Print max weights 
+   */
+  void Print(Option_t* option="") const;
+ /** 
+   * Set the debug level.  The higher the value the more output 
+   * 
+   * @param dbg Debug level 
+   */
+  void SetDebug(Int_t dbg=1) { fDebug = dbg; }
+ /*
+  * Calculate Q vectors
+  *
+  * @param h dN/detadphi histogram
+  * @param eHist histogram for ep vs. eta
+  */
+  void CalcQVectors(TH2D* h, TH1D* eHist);
+ /*
+  * Calculate the eventplane from a vector
+  *
+  * @param v TVector2 of Q-vectors
+  *
+  * @return the eventplane as a double
+  */
+  Double_t CalcEventplane(TVector2 v) const;
+ /*
+  * Set the run number, used for OADB object
+  *
+  * @param run Run number
+  */
+  void SetRunNumber(Int_t run);
+ /*
+  * Get the run number
+  *
+  * @return returns the run number
+  */
+  Int_t GetRunNumber() { return fRunNumber; }
+ /*
+  * Get the OADB phi distribution for flattening
+  */
+  void GetPhiDist();
+ /*
+  * Flag for setting the use of phi weights for flattening
+  * 
+  * @param use true or false
+  */
+  void SetUsePhiWeights(Bool_t use = kTRUE) { fUsePhiWeights = use; }
+ /*
+  * Fill diagnostics hists
+  *
+  * @param fmdEP Object containing results of FMD EP calculations
+  */
+  void FillHists(AliAODForwardEP* fmdEP);
+ /*
+  * Set the OADB path, for using a custom OADB path and file
+  *
+  * @param fname Name of the custom OADB file, including path
+  */
+  void SetOADBPath(Char_t* fname) { fOADBFileName = fname; }
+
+protected:
+ /*
+  * Get the phi weight from OADB histogram for the ep flattening
+  *
+  * @param etaBin which eta bin
+  * @param phiBin which phi bin
+  *
+  * @return phi weight for etaBin, phiBin as double
+  */
+  Double_t GetPhiWeight(Int_t etaBin, Int_t phiBin) const;
+  TList*            fList;              // List for diag. hists.
+  AliVEvent*        fEvent;             // Current event
+  TVector2          fQt;                // Q vector for total ep
+  TVector2          fQa;                // Q vector for sub-ep A
+  TVector2          fQc;                // Q vector for sub-ep C
+  TVector2          fQ1;                // Q vector for sub-ep 1
+  TVector2          fQ2;                // Q vector for sub-ep 2
+  TVector2          fQeta;              // Q vector for psi eta-dependence
+  TH1D*             fHistFMDEventplane; // Diagnostics histogram
+  TH1D*             fHistFMDEventplaneA;// Diagnostics histogram
+  TH1D*             fHistFMDEventplaneC;// Diagnostics histogram
+  TH1D*             fHistFMDEventplane1;// Diagnostics histogram
+  TH1D*             fHistFMDEventplane2;// Diagnostics histogram
+  TH2D*             fHistPhiDist;       // Diagnostics histogram
+  TH1D*             fHistFMDmTPCep;     // Diagnostics histogram
+  TH2F*             fHistFMDvsTPCep;    // Diagnostics histogram
+  TH1D*             fHistFMDmVZEROep;   // Diagnostics histogram
+  TH2F*             fHistFMDvsVZEROep;  // Diagnostics histogram
+  Int_t             fDebug;             // Debug flag
+  TString           fOADBFileName;      // Path to OADB container
+  AliOADBContainer* fOADBContainer;     // OADBContainer object
+  TH2D*             fPhiDist;           // Phi dist. for phi weights
+  Int_t             fRunNumber;         // Run number supplied
+  Bool_t            fUsePhiWeights;     // Flag for phi weights
+
+  ClassDef(AliFMDEventPlaneFinder,1); //  
+};
+
+#endif
+// Local Variables:
+//   mode: C++
+// End:
+
index 63295dc..8972f27 100644 (file)
@@ -12,6 +12,9 @@
 #include <TList.h>
 #include <TROOT.h>
 #include <iostream>
+#include "AliGenHijingEventHeader.h"
+#include <TF1.h>
+#include <TGraph.h>
 
 //____________________________________________________________________
 AliFMDMCTrackDensity::AliFMDMCTrackDensity()
@@ -121,7 +124,8 @@ AliFMDMCTrackDensity::StoreParticle(AliMCParticle*       particle,
                                    Double_t             vz,
                                    UShort_t             nC, 
                                    UShort_t             nT,
-                                   AliESDFMD&           output) const
+                                   AliESDFMD&           output,
+                                    Double_t             w) const
 {
   // Store a particle. 
   if (longest < 0) return;
@@ -141,7 +145,7 @@ AliFMDMCTrackDensity::StoreParticle(AliMCParticle*       particle,
   if (old == AliESDFMD::kInvalidMult) old = 0;
 
   // Increment count 
-  output.SetMultiplicity(d,r,s,t,old+1);
+  output.SetMultiplicity(d,r,s,t,old+w);
 
   // Get track-reference stuff 
   Double_t x      = ref->X();
@@ -211,6 +215,7 @@ AliFMDMCTrackDensity::GetMother(Int_t     iTr,
   return 0;
 }  
 
+// #define USE_FLOW_WEIGHTS 1
 //____________________________________________________________________
 Bool_t
 AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input, 
@@ -247,6 +252,14 @@ AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input,
          output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
     }
   }
+
+#ifdef USE_FLOW_WEIGHTS
+  AliGenHijingEventHeader* hd = dynamic_cast<AliGenHijingEventHeader*>
+                (event.GenEventHeader());
+  Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.);
+  Double_t b = (hd ? hd->ImpactParameter() : -1 );
+#endif
+
   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
   Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
   Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
@@ -332,7 +345,17 @@ AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input,
               fMaxConsequtiveStrips);
        Int_t nnT   = TMath::Abs(oT - ooT) + 1;
        const AliMCParticle* mother = GetMother(iTr, event);
-       StoreParticle(particle, mother, longest, vz, nC, nnT, output);
+       
+#ifdef USE_FLOW_WEIGHTS
+        phi = (mother ? mother->Phi() : particle->Phi());
+        eta = (mother ? mother->Eta() : particle->Eta());
+        Double_t pt = (mother ? mother->Pt() : particle->Pt());
+        Int_t   id = (mother ? mother->PdgCode() : 2212);
+       Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
+#else
+        Double_t weight = 1; // // 
+#endif
+        StoreParticle(particle, mother, longest, vz, nC, nnT, output, weight);
        longest = -1;
        angle   = 0;
        nC  = 1;    // Reset track-ref counter - we have this->1
@@ -374,10 +397,124 @@ AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input,
       Info("Process", "I=%3d L=%3d nT=%3d (out of %3d)", 
           iTr, longest, nT, fMaxConsequtiveStrips);
     const AliMCParticle* mother = GetMother(iTr, event);
-    StoreParticle(particle, mother, longest, vz, nC, nT, output);
+
+#ifdef USE_FLOW_WEIGHTS
+    phi = (mother ? mother->Phi() : particle->Phi());
+    eta = (mother ? mother->Eta() : particle->Eta());
+    Double_t pt = (mother ? mother->Pt() : particle->Pt());
+    Int_t    id = (mother ? mother->PdgCode() : 2212);
+       Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
+#else
+        Double_t weight = 1; // // 
+#endif
+    StoreParticle(particle, mother, longest, vz, nC, nT, output, weight);
   } // Loop over tracks
   return kTRUE;
 }
+
+//____________________________________________________________________
+Double_t
+AliFMDMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, Double_t b, 
+                                     Double_t phi, Double_t rp, Int_t id) const
+{
+  static TF1 gaus = TF1("gaus", "gaus", -6, 6);
+  gaus.SetParameters(0.1, 0., 9);
+  //  gaus.SetParameters(0.1, 0., 3);
+  //  gaus.SetParameters(0.1, 0., 15);
+  
+  const Double_t xCumulant2nd4050ALICE[] = {0.00, 0.25, 0.350,
+                                           0.45, 0.55, 0.650, 
+                                           0.75, 0.85, 0.950,
+                                           1.10, 1.30, 1.500,
+                                           1.70, 1.90, 2.250,
+                                           2.75, 3.25, 3.750,
+                                           4.50};
+  const Double_t yCumulant2nd4050ALICE[] = {0.00000, 0.043400,
+                                           0.059911,0.073516,
+                                           0.089756,0.105486,
+                                           0.117391,0.128199,
+                                           0.138013,0.158271,
+                                           0.177726,0.196383,
+                                           0.208277,0.216648,
+                                           0.242954,0.249961,
+                                           0.240131,0.269006,
+                                           0.207796};
+  const Int_t nPointsCumulant2nd4050ALICE = 
+    sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);                                      
+  static TGraph alicePointsPt2(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
+#if 0
+  const Double_t xCumulant4th3040ALICE[] = {0.00,0.250,0.35,
+                                           0.45,0.550,0.65,
+                                           0.75,0.850,0.95,
+                                           1.10,1.300,1.50,
+                                           1.70,1.900,2.25,
+                                           2.75,3.250,3.75,
+                                           4.50,5.500,7.00,
+                                           9.000000};
+  const Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,
+                                           0.048566,0.061083,
+                                           0.070910,0.078831,
+                                           0.091396,0.102026,
+                                           0.109691,0.124449,
+                                           0.139819,0.155561,
+                                           0.165701,0.173678,
+                                           0.191149,0.202015,
+                                           0.204540,0.212560,
+                                           0.195885,0.000000,
+                                           0.000000,0.000000};
+#endif
+  const Double_t xCumulant4th4050ALICE[] = {0.00,0.25,0.350,
+                                           0.45,0.55,0.650,
+                                           0.75,0.85,0.950,
+                                           1.10,1.30,1.500,
+                                           1.70,1.90,2.250,
+                                           2.75,3.25,3.750,
+                                           4.50};
+  const Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,
+                                           0.049824,0.066662,
+                                           0.075856,0.081583,
+                                           0.099778,0.104674,
+                                           0.118545,0.131874,
+                                           0.152959,0.155348,
+                                           0.169751,0.179052,
+                                           0.178532,0.198851,
+                                           0.185737,0.239901,
+                                           0.186098};
+  const Int_t nPointsCumulant4th4050ALICE = 
+    sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);   
+  static TGraph alicePointsPt4(nPointsCumulant4th4050ALICE, 
+                              xCumulant4th4050ALICE, 
+                              yCumulant4th4050ALICE);
+
+  const Double_t xCumulant4thTPCrefMultTPConlyAll[] = {1.75,
+                                                      4.225,
+                                                      5.965,
+                                                      7.765,
+                                                      9.215,
+                                                      10.46,
+                                                      11.565,
+                                                      12.575};
+  const Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,
+                                                      0.055818,0.073137,
+                                                      0.083898,0.086690,
+                                                      0.082040,0.077777};
+  const Int_t nPointsCumulant4thTPCrefMultTPConlyAll = 
+    sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
+  TGraph aliceCent(nPointsCumulant4thTPCrefMultTPConlyAll,
+                  xCumulant4thTPCrefMultTPConlyAll,
+                  yCumulant4thTPCrefMultTPConlyAll);
+
+
+  Double_t weight = (20. * gaus.Eval(eta) * (alicePointsPt2.Eval(pt) * 0.5 + 
+                                            alicePointsPt4.Eval(pt) * 0.5) 
+                    * (aliceCent.Eval(b) / aliceCent.Eval(10.46)) 
+                    * 2. * TMath::Cos(2. * (phi - rp)));
+  if      (TMath::Abs(id) == 211)  weight *= 1.3; //pion flow
+  else if (TMath::Abs(id) == 2212) weight *= 1.0;  //proton flow
+  else                             weight *= 0.7;
+  
+  return weight;
+}
 //____________________________________________________________________
 void
 AliFMDMCTrackDensity::Print(Option_t* /*option*/) const 
index 25d28a4..08fa2d3 100644 (file)
@@ -138,7 +138,8 @@ protected:
                     Double_t       vz,
                     UShort_t       nC, 
                     UShort_t       nT,
-                    AliESDFMD&     output) const;
+                    AliESDFMD&     output,
+                     Double_t       w) const;
   /** 
    * Get incident angle of this track reference
    * 
@@ -158,6 +159,21 @@ protected:
    * @return Pointer to mother or null 
    */
   const AliMCParticle* GetMother(Int_t iTr, const AliMCEvent& event) const;
+  /** 
+   * Calculate flow weight 
+   *
+   * @param eta  Pseudo rapidity 
+   * @param pt   Transverse momemtum 
+   * @param b    Impact parameter
+   * @param phi  Azimuthal angle 
+   * @param rp   Reaction plance angle 
+   * @param id   Particle PDG code
+   *
+   * @return Flow weight for the particle
+   */
+  Double_t CalculateWeight(Double_t eta, Double_t pt, Double_t b, 
+                          Double_t phi, Double_t rp, Int_t id) const;
+
   Bool_t   fUseOnlyPrimary;       // Only use primaries 
   UShort_t fMaxConsequtiveStrips; // Max 'cluster' size
   TH1D*    fNr;                   // Number of track-refs per cluster
index e2813fd..532b73a 100644 (file)
@@ -15,7 +15,6 @@
 #include <TList.h>
 #include <iostream>
 #include <TMath.h>
-#include "THnSparse.h"
 #include "TH3D.h"
 #include "TProfile2D.h"
 #include "AliLog.h"
 #include "AliAnalysisManager.h"
 #include "AliAODHandler.h"
 #include "AliAODInputHandler.h"
-#include "AliAODMCParticle.h"
 #include "AliAODForwardMult.h"
+#include "AliAODCentralMult.h"
 #include "AliAODEvent.h"
 
-//
-// Enumeration for adding and retrieving stuff from the histogram
-//
-enum { kW2Two = 1, kW2, kW4Four, kW4, kQnRe, kQnIm, kM,
-       kCosphi1phi2, kSinphi1phi2, kCosphi1phi2phi3m, kSinphi1phi2phi3m, kMm1m2, 
-       kw2two, kw2, kw4four, kw4, kpnRe, kpnIm, kmp, 
-       kCospsi1phi2, kSinpsi1phi2, kCospsi1phi2phi3m, kSinpsi1phi2phi3m,
-       kmpmq, kCospsi1phi2phi3p, kSinpsi1phi2phi3p };
-
 ClassImp(AliForwardFlowTaskQC)
 #if 0
 ; // For emacs 
 #endif
 
 AliForwardFlowTaskQC::AliForwardFlowTaskQC()
-  : fOutputList(0),    // Output list
-    fFlowUtil(0),      // AliForwardFlowUtil
+  : AliAnalysisTaskSE(),
+    fBinsFMD(),         // List with FMD flow histos
+    fBinsSPD(),         // List with SPD flow histos
+    fSumList(0),       // Event sum list
+    fOutputList(0),    // Result output list
     fAOD(0),           // AOD input event
-    fMC(kFALSE),       // MC flag
-    fEtaBins(48),      // # of etaBin bins in histograms
-    fEtaRef(12),        // # of Eta bins for reference flow
-    fAddFlow(0),       // Add flow string
-    fAddType(0),       // Add flow type #
-    fAddOrder(0),      // Add flow order
-    fZvertex(1111),    // Z vertex range
-    fCent(-1)          // Centrality
+    fZvertex(1111),    // Z vertex coordinate
+    fCent(-1),         // Centrality
+    fHistCent(),        // Histo for centrality
+    fHistVertexSel(),   // Histo for selected vertices
+    fHistVertexAll()    // Histo for all vertices
 {
   // 
   // Default constructor
@@ -60,19 +50,18 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC()
   for (Int_t n = 0; n <= 6; n++) fv[n] = kTRUE;
 }
 //_____________________________________________________________________
-AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) :
-  AliAnalysisTaskSE(name),
-  fOutputList(0),      // Output list
-  fFlowUtil(0),                // AliForwardFlowUtil
-  fAOD(0),             // AOD input event
-  fMC(kFALSE),         // MC flag
-  fEtaBins(48),                // # of Eta bins
-  fEtaRef(12),          // # of Eta bins for reference flow
-  fAddFlow(0),         // Add flow string
-  fAddType(0),         // Add flow type #
-  fAddOrder(0),                // Add flow order
-  fZvertex(1111),      // Z vertex range
-  fCent(-1)            // Centrality
+AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) 
+  : AliAnalysisTaskSE(name),
+    fBinsFMD(),         // List with FMD flow histos
+    fBinsSPD(),         // List with SPD flow histos
+    fSumList(0),        // Event sum list           
+    fOutputList(0),     // Result output list       
+    fAOD(0),           // AOD input event          
+    fZvertex(1111),     // Z vertex coordinate      
+    fCent(-1),          // Centrality               
+    fHistCent(),        // Histo for centrality
+    fHistVertexSel(),   // Histo for selected vertices
+    fHistVertexAll()    // Histo for all vertices      
 {
   // 
   // Constructor
@@ -81,22 +70,27 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) :
   //  name: Name of task
   //
   for (Int_t n = 0; n <= 6; n++) fv[n] = kTRUE;
+
+  fHistCent      = new TH1D("hCent", "Centralities", 100, 0, 100);
+  fHistVertexSel = new TH1D("hVertexSel", "Selected vertices", 40, -20, 20);
+  fHistVertexAll = new TH1D("hVertexAll", "All vertices", 40, -20, 20);
+
   DefineOutput(1, TList::Class());
+  DefineOutput(2, TList::Class());
 }
 //_____________________________________________________________________
-AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) :
-  AliAnalysisTaskSE(o),
-  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
-  fEtaRef(o.fEtaRef),           // # of Eta bins for reference flow
-  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(o.fCent)               // Centrality
+AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o)
+  : AliAnalysisTaskSE(o),
+    fBinsFMD(),                        // List with FMD flow histos
+    fBinsSPD(),                        // List with SPD flow histos
+    fSumList(o.fSumList),              // Event sum list           
+    fOutputList(o.fOutputList),        // Result output list       
+    fAOD(o.fAOD),                     // AOD input event          
+    fZvertex(o.fZvertex),              // Z vertex coordinate      
+    fCent(o.fCent),                   // Centrality               
+    fHistCent(o.fHistCent),            // Histo for centrality
+    fHistVertexSel(o.fHistVertexSel),  // Histo for selected vertices
+    fHistVertexAll(o.fHistVertexAll)   // Histo for all vertices      
 {
   // 
   // Copy constructor 
@@ -105,221 +99,557 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) :
   //    o Object to copy from 
   //
   for (Int_t n = 0; n <= 6; n++) fv[n] = o.fv[n];
-  DefineOutput(1, TList::Class());
+}
+//_____________________________________________________________________
+AliForwardFlowTaskQC&
+AliForwardFlowTaskQC::operator=(const AliForwardFlowTaskQC& o)
+{
+  // 
+  // Assignment operator 
+  //
+  fSumList       = o.fSumList;
+  fOutputList    = o.fOutputList;
+  fAOD           = o.fAOD;
+  fZvertex       = o.fZvertex;
+  fCent          = o.fCent;
+  fHistCent      = o.fHistCent;
+  fHistVertexSel = o.fHistVertexSel;
+  fHistVertexAll = o.fHistVertexAll;
+  fHistVertexAll = o.fHistVertexAll;
+
+  for (Int_t n = 0; n <= 6; n++) fv[n] = o.fv[n];
+  return *this;
 }
 //_____________________________________________________________________
 void AliForwardFlowTaskQC::UserCreateOutputObjects()
 {
   //
   // Create output objects
-  // 
-  if (!fOutputList)
-    fOutputList = new TList();
-  fOutputList->SetName("QCumulants");
-  fOutputList->SetOwner();
-  fFlowUtil = new AliForwardFlowUtil(fOutputList);
-  
-  Double_t x[101] = { 0. };
-  // Double_t etaMin = -6;
-  // Double_t etaMax = 6;
-
-  // First we have a number of options for eta binning, also if it is
-  // not really eta, but centrality or pt we want to do flow as a
-  // function of, then this is possible:
-  if (fEtaBins == 5) {
-    x[0] = 0.;
-    x[1] = 1.;
-    x[2] = 2.;
-    x[3] = 3.;
-    x[4] = 4.5;
-    x[5] = 6.0;
-    // etaMin = 0;
-    // etaMax = 6;
-  }
+  //
+  InitVertexBins();
+  InitHists();
 
-  else if (fEtaBins == 100) { 
-    for (Int_t e = 0; e<= 100; e++) {
-      x[e] = e;
-    }
-    // etaMin = 0;
-    // etaMax = 100;
-  }
-  
-  else {
-    if (fEtaBins % 6) fEtaBins = 48;
-    for (Int_t e = 0; e <=fEtaBins; e++) {
-      x[e] = -6. + e*(12./(Double_t)fEtaBins);
+  PostData(1, fSumList);
+  PostData(2, fOutputList);
+
+}
+//_____________________________________________________________________
+void AliForwardFlowTaskQC::InitVertexBins()
+{
+  // 
+  // Init vertexbin objects for FMD and SPD, and add them to the lists
+  //
+  for(Int_t n = 1; n <= 6; n++) {
+    if (!fv[n]) continue;
+    for (Int_t v = -10; v < 10; v++) {
+      fBinsFMD.Add(new VertexBin(v, v+1, n, "FMD"));
+      fBinsSPD.Add(new VertexBin(v, v+1, n, "SPD"));
     }
   }
-  // Reference flow binning
-  Double_t xR[101] = { 0. };
-  for (Int_t r = 0; r <= fEtaRef; r++) {
-    xR[r] = -6 + r*(12./(Double_t)fEtaRef);
-  }
-  // Phi binning
-  Double_t phi[21] =  { 0. };
-  for (Int_t p = 0; p <= 20; p++) {
-    phi[p] = p*2.*TMath::Pi() / 20.;
+
+}
+//_____________________________________________________________________
+void AliForwardFlowTaskQC::InitHists()
+{
+  //
+  // Init histograms and add vertex bin histograms to the sum list
+  //
+  if (!fSumList)
+    fSumList = new TList();
+  fSumList->SetName("Sums");
+  fSumList->SetOwner();
+
+  TList* dList = new TList();
+  dList->SetName("Diagnostics");
+  dList->Add(fHistCent);
+  dList->Add(fHistVertexSel);
+  dList->Add(fHistVertexAll);
+  fSumList->Add(dList);
+
+  TIter nextFMD(&fBinsFMD);
+  VertexBin* bin = 0;
+  while ((bin = static_cast<VertexBin*>(nextFMD()))) {
+    bin->AddOutput(fSumList);
   }
-  Double_t phiMC[201] = { 0. };
-  for (Int_t p = 0; p <= 200; p++) {
-    phiMC[p] = p*2.*TMath::Pi() / 200.;
+  TIter nextSPD(&fBinsSPD);
+  while ((bin = static_cast<VertexBin*>(nextSPD()))) {
+    bin->AddOutput(fSumList);
   }
-
-  // Histograms for cumulants analysis
-
-  // We loop over flow histograms here to add different orders of harmonics
-  TString type = "";
-  for (Int_t loop = 1; loop <= 5; loop++) {
-    
-    if (loop == 1) type = "FMD";
-    if (loop == 2) type = "SPD";
-    if (loop == 3) type = "MC";
-    if (loop == 4) type = "FMDTR";
-    if (loop == 5) type = "SPDTR";
-
-    TList* tList = new TList();
-    tList->SetName(type.Data());
-    fOutputList->Add(tList);
-
-    for (Int_t n = 1; n <= 6; n++) {
-      if (!fv[n]) continue;
-      
-      TList* vList = new TList();
-      vList->SetName(Form("v%d", n));
-      tList->Add(vList);
-      
-      TString tag = TString();
-      for (Int_t c = 0; c < 100; c++) {
-        // Output histograms
-        tag = Form("%d_%d", c, c+1);
-          
-        TH3D* hFlowHist = new TH3D(Form("hQ%dCumuHist%s_%s", n, type.Data(), tag.Data()), 
-                   Form("hQ%dCumuHist%s_%s", n, type.Data(), tag.Data()), fEtaBins, -6, 6, 10, -5, 5, 26, 0.5, 26.5);
-        TAxis* xAxis = hFlowHist->GetXaxis();
-        xAxis->Set(fEtaBins, x);
-        hFlowHist->Sumw2();
-        vList->Add(hFlowHist);
-
-        TProfile* hCumulant2DiffFlow = new TProfile(Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), 
-                     Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
-        hCumulant2DiffFlow->Sumw2();
-        vList->Add(hCumulant2DiffFlow);
-   
-        TProfile* hCumulant4DiffFlow = new TProfile(Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), 
-                     Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
-        hCumulant4DiffFlow->Sumw2();
-        vList->Add(hCumulant4DiffFlow);
-      } // end of centrality loop
-
-    } // end of v_{n} loop
-    // Single Event histograms
-    TH2D* hdNdphiSE = new TH2D(Form("hdNdphiSE%s", type.Data()),
-                 Form("hdNdphiSE%s", type.Data()), fEtaRef, xR, (type.Contains("MC") ? 200 : 20), (type.Contains("MC") ? phiMC : phi));
-    hdNdphiSE->Sumw2();
-    fOutputList->Add(hdNdphiSE);
-
-    TH2D* hdNdetadphiSE = new TH2D(Form("hdNdetadphiSE%s", type.Data()), 
-                 Form("hdNdetadphiSE%s", type.Data()), fEtaBins, x, (type.Contains("MC") ? 200 : 20), (type.Contains("MC") ? phiMC : phi));
-    hdNdetadphiSE->Sumw2();
-    fOutputList->Add(hdNdetadphiSE);
-  } // end of type loop
-
-  TProfile2D* pMCTruth = new TProfile2D("pMCTruth", "pMCTruth", 48, -6, 6, 100, 0, 100);
-  TAxis* xAxis = pMCTruth->GetXaxis();
-  xAxis->Set(fEtaBins, x);
-  pMCTruth->Sumw2();
-  fOutputList->Add(pMCTruth);
-
-  // Monitoring plots
-  TH1D* cent = new TH1D("Centralities", "Centralities", 100, 0, 100);
-  fOutputList->Add(cent);
-
-  TH1D* vertexSel = new TH1D("VertexSelected", "VertexSelected", 50, -10, 10);
-  fOutputList->Add(vertexSel);
-  TH1D* vertexAll = new TH1D("VertexAll", "VertexAll", 50, -10, 10);
-  fOutputList->Add(vertexAll);
-  
-  TH2D* vertex2D = new TH2D("CoverageVsVertex", "CoverageVsVertex", fEtaBins, -6, 6, 20, -10, 10);
-  fOutputList->Add(vertex2D);
-  PostData(1, fOutputList);
 }
 //_____________________________________________________________________
 void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
 {
-  // 
-  // Process each event
+  //
+  // Calls the analyze function - called every event
   //
   // Parameters:
   //  option: Not used
   //
 
+  Analyze();
+  
+  PostData(1, fSumList);
+
+}
+//_____________________________________________________________________
+Bool_t AliForwardFlowTaskQC::Analyze()
+{
+  // 
+  // Load FMD and SPD objects from aod tree and call the cumulants 
+  // calculation for the correct vertexbin
+  //
+
   // Reset data members
   fCent = -1;
   fZvertex = 1111;
 
   // Get input event
   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
-  if (!fAOD) return;
+  if (!fAOD) return kFALSE;
 
-  // Fill histograms
-  if (!fFlowUtil->LoopAODFMD(fAOD)) return;
-  if (!fFlowUtil->LoopAODSPD(fAOD)) return;
+  AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
+  if (!aodfmult) return kFALSE;
+  if (!AODCheck(aodfmult)) return kFALSE;
+  TH2D fmddNdetadphi = aodfmult->GetHistogram();
 
-  // Get centrality and vertex from flow utility and fill monitoring histograms
-  fCent = fFlowUtil->GetCentrality();
-  fZvertex = fFlowUtil->GetVertex();
+  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
+  if (!aodcmult) return kFALSE;
+  TH2D spddNdetadphi = aodcmult->GetHistogram();
 
-  TH1D* cent = (TH1D*)fOutputList->FindObject("Centralities");
-
-  cent->Fill(fCent);
-  TH1D* vertex = (TH1D*)fOutputList->FindObject("VertexSelected");
-  vertex->Fill(fZvertex);
+  // TODO: remove me!
+//  fCent = 0.5;
 
   // Run analysis on FMD and SPD
-  for (Int_t n = 1; n <= 6; n++) {
-    if (fv[n]) {
-      CumulantsMethod("FMD", n);
-      CumulantsMethod("SPD", n);
+  TIter nextFMD(&fBinsFMD);
+  VertexBin* bin = 0;
+  while ((bin = static_cast<VertexBin*>(nextFMD()))) {
+    if (bin->CheckVertex(fZvertex)) {
+      if (!bin->FillHists(fmddNdetadphi)) return kFALSE; 
+      bin->CumulantsAccumulate(fCent);
     }
   }
 
-  // And do analysis if there is.
-  if (fMC) { 
-    ProcessPrimary();
+  TIter nextSPD(&fBinsSPD);
+  while ((bin = static_cast<VertexBin*>(nextSPD()))) {
+    if (bin->CheckVertex(fZvertex)) {
+      if (!bin->FillHists(spddNdetadphi)) return kFALSE;
+      bin->CumulantsAccumulate(fCent);
+    }
   }
 
-  PostData(1, fOutputList);
+  return kTRUE;
 }
 //_____________________________________________________________________
-void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", 
-                                          Int_t harmonic = 2)
+void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
 {
-  // 
-  // Calculate the Q cumulant of order n
+  //
+  // Calls the finalize function, done at the end of the analysis
   //
   // Parameters:
-  //  type: Determines which histograms should be used
-  //        - "FMD" or "SPD" = data histograms
-  //        - "FMDTR" or "SPDTR" = track reference histograms
-  //        - "MC" = MC truth histograms
-  //  harmonic: Which harmonic to calculate
+  //  option: Not used
   //
-  Double_t n = harmonic;
-  
-  TList* tList = (TList*)fOutputList->FindObject(type.Data());
-  TList* vList = (TList*)tList->FindObject(Form("v%d", harmonic));
-  // We get the histograms
-  TH3D* flowHist   = (TH3D*)vList->FindObject(Form("hQ%dCumuHist%s_%d_%d", harmonic, type.Data(), (Int_t)fCent, (Int_t)fCent+1));
-  TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject(Form("hdNdetadphiSE%s", type.Data()));
-  TH2D* dNdphi = (TH2D*)fOutputList->FindObject(Form("hdNdphiSE%s", /*"FMD"*/type.Data()));
+
+  // Make sure pointers are set to the correct lists
+  fSumList = dynamic_cast<TList*> (GetOutputData(1));
+  if(!fSumList) {
+    AliError("Could not retrieve TList fSumList"); 
+    return; 
+  }
+
+  if (!fOutputList)
+    fOutputList = new TList();
+  fOutputList->SetName("Results");
+  fOutputList->SetOwner();
+
+  // Run finalize on VertexBins
+  Finalize();
+
+  // Collect centralities
+  TProfile2D* hist2D = 0;
+  TList* centList = 0;
+  TH1D* hist1D = 0;
+  TIter nextProfile(fOutputList);
+  while ((hist2D = dynamic_cast<TProfile2D*>(nextProfile()))) {
+    for (Int_t cBin = 1; cBin <= 100; ) {
+      Int_t cMin = cBin - 1;
+      Int_t cMax = (cMin < 80 ? (cMin < 20 ? cMin + 5 : cMin + 10) : cMin + 20);
+      TString name = Form("cent_%d-%d", cMin, cMax);
+      centList = (TList*)fOutputList->FindObject(name.Data());
+      if (!centList) { 
+       centList = new TList();
+       centList->SetName(name.Data());
+       fOutputList->Add(centList);
+      }
+      hist1D = hist2D->ProjectionX(Form("%s_%s", hist2D->GetName(), name.Data()), 
+                                  cMin, cMax, "E");
+      hist1D->SetTitle(hist1D->GetName());
+      hist1D->Scale(1./(cMax-cMin));
+      centList->Add(hist1D);
+
+      cBin = cMax+1;
+    }
+  }
+
+  PostData(2, fOutputList);
+
+}
+//_____________________________________________________________________
+void AliForwardFlowTaskQC::Finalize()
+{
+  //
+  // Finalize command, called by Terminate()
+  //
+
+
+  // Reinitiate vertex bins if Terminate is called separately!
+  if (fBinsFMD.GetEntries() == 0) InitVertexBins();
+
+  // Iterate over all vertex bins objects and finalize cumulants
+  // calculations
+  TIter nextFMD(&fBinsFMD);
+  VertexBin* bin = 0;
+  while ((bin = static_cast<VertexBin*>(nextFMD()))) {
+    bin->CumulantsTerminate(fSumList, fOutputList);
+  }
+  TIter nextSPD(&fBinsSPD);
+  while ((bin = static_cast<VertexBin*>(nextSPD()))) {
+    bin->CumulantsTerminate(fSumList, fOutputList);
+  }
+
+} 
+// _____________________________________________________________________
+Bool_t AliForwardFlowTaskQC::AODCheck(const AliAODForwardMult* aodfm) 
+{
+  // 
+  // Function to check that and AOD event meets the cuts
+  //
+  // Parameters: 
+  //  AliAODForwardMult: forward mult object with trigger and vertex info
+  //
+  // Returns false if there is no trigger or if the centrality or vertex
+  // is out of range. Otherwise true.
+  //
+
+  // First check for trigger
+  if (!aodfm->IsTriggerBits(AliAODForwardMult::kOffline)) return kFALSE;
+
+  // Then check for centrality
+  fCent = (Double_t)aodfm->GetCentrality();
+  if (0. >= fCent || fCent >= 100.) return kFALSE;
+  fHistCent->Fill(fCent);
+
+  // And finally check for vertex
+  fZvertex = aodfm->GetIpZ();
+  fHistVertexAll->Fill(fZvertex);
+  if (TMath::Abs(fZvertex) >= 10.) return kFALSE;
+  fHistVertexSel->Fill(fZvertex);
+
+  return kTRUE;
+
+}
+//_____________________________________________________________________
+AliForwardFlowTaskQC::VertexBin::VertexBin()
+  : TNamed(),
+    fMoment(0),      // Flow moment for this vertexbin
+    fVzMin(0),       // Vertex z-coordinate min
+    fVzMax(0),       // Vertex z-coordinate max
+    fType(),         // Data type (FMD/SPD/FMDTR/SPDTR/MC)
+    fCumuRef(),      // Histogram for reference flow
+    fCumuDiff(),     // Histogram for differential flow
+    fCumuHist(),     // Sum histogram for cumulants
+    fHistTwoCorr(),  // Diagnostics histogram for <<2>>
+    fHistW2(),       // Diagnostics histogram for w_<<2>>
+    fHistFourCorr(), // Diagnostics histogram for <<4>>
+    fHistW4(),       // Diagnostics histogram for w_<<4>>      
+    fdNdedpAcc()     // Diagnostics histogram to make acc. maps
+{
+  //
+  // Default constructor
+  //
+}
+//_____________________________________________________________________
+AliForwardFlowTaskQC::VertexBin::VertexBin(const Int_t vLow, const Int_t vHigh, 
+                                           const Int_t moment, const Char_t* name)
+  : TNamed("", ""),
+    fMoment(moment),  // Flow moment for this vertexbin
+    fVzMin(vLow),     // Vertex z-coordinate min
+    fVzMax(vHigh),    // Vertex z-coordinate max
+    fType(name),      // Data type (FMD/SPD/FMDTR/SPDTR/MC)
+    fCumuRef(),       // Histogram for reference flow
+    fCumuDiff(),      // Histogram for differential flow
+    fCumuHist(),      // Sum histogram for cumulants
+    fHistTwoCorr(),   // Diagnostics histogram for <<2>>
+    fHistW2(),        // Diagnostics histogram for w_<<2>>
+    fHistFourCorr(),  // Diagnostics histogram for <<4>>
+    fHistW4(),        // Diagnostics histogram for w_<<4>> 
+    fdNdedpAcc()      // Diagnostics histogram to make acc. maps
+{
+  //
+  // Constructor
+  //
+  // Parameters
+  //  vLow: min z-coordinate
+  //  vHigh: max z-coordinate
+  //  moment: flow moment
+  //  name: data type name (FMD/SPD/FMDTR/SPDTR/MC)
+  //
+  SetName(Form("%svertexBin%d_%d_%d", name, moment, vLow, vHigh));
+  SetTitle(Form("%svertexBin%d_%d_%d", name, moment, vLow, vHigh));
+
+}
+//_____________________________________________________________________
+AliForwardFlowTaskQC::VertexBin&
+AliForwardFlowTaskQC::VertexBin::operator=(const AliForwardFlowTaskQC::VertexBin& o)
+{
+  //
+  // Assignment operator
+  //
+  // Parameters
+  //  o: AliForwardFlowTaskQC::VertexBin
+  //
+  fType         = o.fType;
+  fCumuRef      = o.fCumuRef;
+  fCumuDiff     = o.fCumuDiff;
+  fCumuHist     = o.fCumuHist;
+  fHistTwoCorr  = o.fHistTwoCorr;
+  fHistW2       = o.fHistW2;
+  fHistFourCorr = o.fHistFourCorr;
+  fHistW4       = o.fHistW4;
+  fdNdedpAcc    = o.fdNdedpAcc;
+
+  return *this;
+}
+//_____________________________________________________________________
+void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist)
+{
+  //
+  // Add histograms to outputlist
+  //
+  // Parameters
+  //  outputlist: list of histograms
+  //
+
+  // First we try to find an outputlist for this vertexbin
+  TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d", fType, fVzMin, fVzMax));
+
+  // If it doesn't exist we make one
+  if (!list) {
+    list = new TList();
+    list->SetName(Form("%svertex_%d_%d", fType, fVzMin, fVzMax));
+    outputlist->Add(list);
+  }
+
+  // We initiate the reference histogram according to an acceptance correction map,
+  // so we don't shift the SPD coverage within a reference bin
+  fCumuRef = new TH2D(Form("%s_v%d_%d_%d_ref", fType, fMoment, fVzMin, fVzMax),
+                        Form("%s_v%d_%d_%d_ref", fType, fMoment, fVzMin, fVzMax),
+                        24, -6., 6., 5, 0.5, 5.5);
+  TFile acc("$ALICE_ROOT/PWGLF/FORWARD/corrections/FlowCorrections/FlowAccMap.root", "READ");
+  TH1D* accMap = (TH1D*)acc.Get(Form("%saccVertex_%d_%d", fType, fVzMin, fVzMax));
+  if (accMap) {
+    Int_t nBins = accMap->GetNbinsX();
+    Double_t eta[48] = { 0. };
+    Int_t n = 0;
+    Double_t newOcc[48] = { 0. };
+    Double_t prev = -1;
+    for (Int_t i = 0; i < nBins; i++) {
+      Double_t occ = accMap->GetBinContent(i+1);
+      if (prev != occ && (occ > 0.6 || occ == 0)) {
+        eta[n] = i*0.25-6.;
+        newOcc[n] = occ;
+        n++;
+//        printf("eta: %f \t occ: %f \t Vertex: %d \n", eta[n-1], occ, fVzMin);
+      }
+      prev = occ;
+    }
+    eta[n] = 6.;
+    fCumuRef->GetXaxis()->Set(n, eta);
+  }
+  acc.Close();
+
+  fCumuRef->Sumw2();
+  list->Add(fCumuRef);
+
+  // We initiate the differential histogram
+  fCumuDiff = new TH2D(Form("%s_v%d_%d_%d_diff", fType, fMoment, fVzMin, fVzMax),
+                        Form("%s_v%d_%d_%d_diff", fType, fMoment, fVzMin, fVzMax),
+                        48, -6., 6., 5, 0.5, 5.5);
+  fCumuDiff->Sumw2();
+  list->Add(fCumuDiff);
+
+  // Initiate the cumulant sum histogram
+  fCumuHist = new TH3D(Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax),
+                       Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax),
+                       48, -6., 6., 100, 0., 100., 26, 0.5, 26.5);
+  fCumuHist->Sumw2();
+
+  list->Add(fCumuHist);
+
+  // We check for diagnostics histograms (only done per type and moment, not vertexbin)
+  // If they are not found we create them.
+  TList* dList = (TList*)outputlist->FindObject("Diagnostics");
+  if (!dList) AliFatal("No diagnostics list found, what kind of game are you running here?!?!");
+
+  // Corr. hists are shared over all vertex bins...
+  fHistTwoCorr    = (TH3D*)dList->FindObject(Form("hHistTwoCorr_%s_v%d", fType, fMoment));
+  if (fHistTwoCorr) {
+    fHistW2       = (TH3D*)dList->FindObject(Form("hHistW2_%s_v%d", fType, fMoment));
+    fHistFourCorr = (TH3D*)dList->FindObject(Form("hHistFourCorr_%s_v%d", fType, fMoment));
+    fHistW4       = (TH3D*)dList->FindObject(Form("hHistW4_%s_v%d", fType, fMoment));
+  } else {
+    fHistTwoCorr  = new TH3D(Form("hHistTwoCorr_%s_v%d", fType, fMoment), 
+                            "Two particle correlator: w_{2}<<2>>", 
+                            48, -6., 6., 100, 0., 150000, 100, 0., 100.);
+    fHistTwoCorr->Sumw2();
+    fHistW2       = new TH3D(Form("hHistW2_%s_v%d", fType, fMoment),
+                            "Two particle event weight: w_{2}",
+                            48, -6., 6., 100, 0., 2e+7, 100, 0., 100.);
+    fHistW2->Sumw2();
+    fHistFourCorr = new TH3D(Form("hHistFourCorr_%s_v%d", fType, fMoment),  
+                            "Four particle correlator: w_{4}<<4>>", 
+                            48, -6., 6., 100, 0., 1e+10, 100, 0., 100.);
+    fHistFourCorr->Sumw2();
+    fHistW4       = new TH3D(Form("hHistW4_%s_v%d", fType, fMoment), 
+                            "Four particle event weight: w_{4}",
+                            48, -6., 6., 100, 0., 4e+14, 100, 0., 100.);
+    fHistW4->Sumw2();
+
+    dList->Add(fHistTwoCorr);
+    dList->Add(fHistW2);
+    dList->Add(fHistFourCorr);
+    dList->Add(fHistW4);
+  }
   
-  // We create the coordinate array used to fill the THnSparse, centrality and vertex is set from the beginning.
-  Double_t coord[4] = { 0., fCent, fZvertex, 0. };
+  // Acceptance hists are shared over all moments
+  fdNdedpAcc = (TH2D*)dList->FindObject(Form("h%sdNdedpAcc_%d_%d", fType, fVzMin, fVzMax));
+  if (!fdNdedpAcc) {
+    fdNdedpAcc = new TH2D(Form("h%sdNdedpAcc_%d_%d", fType, fVzMin, fVzMax), 
+                          Form("%s acceptance map for %d cm < v_{z} < %d cm", fType, fVzMin, fVzMax),
+                          48, -6, 6, 20, 0, TMath::TwoPi());
+    fdNdedpAcc->Sumw2();
+    dList->Add(fdNdedpAcc);
+  }
+
+
+}
+//_____________________________________________________________________
+Bool_t AliForwardFlowTaskQC::VertexBin::CheckVertex(Double_t vz)
+{
+  //
+  // We check if this is the correct bin for the current event's vertex
+  //
+  // Parameters:
+  //  vZ: Current event vertex
+  //
+  // Returns false if out of range, true otherwise
+  //
+  if ((Double_t)fVzMin > vz) return kFALSE;
+  if ((Double_t)fVzMax <= vz) return kFALSE;
+
+  return kTRUE; 
+}
+//_____________________________________________________________________
+Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D dNdetadphi) 
+{
+  // 
+  // Fill reference and differential eta-histograms
+  //
+  // Parameters:
+  //  dNdetadphi: 2D histogram with input data
+  //
+
+  if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
+
+  // Fist we reset histograms
+  fCumuRef->Reset();
+  fCumuDiff->Reset();
+
+  // Numbers to cut away bad events and acceptance.
+  Double_t runAvg = 0;
+  Double_t max = 0;
+  Int_t nInAvg = 0;
+  Int_t nBadBins = 0;
+  Int_t nBins = (dNdetadphi.GetNbinsX() * 6) / (fCumuDiff->GetNbinsX() * 5);
+  Int_t nInBin = 0;
+  Int_t nCurBin = 0, nPrevBin = 0;
+
+  // Then we loop over the input and calculate sum cos(k*n*phi)
+  // and fill it in the reference and differential histograms
+  Double_t eta, phi, weight;
+  Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0;
+  for (Int_t etaBin = 1; etaBin <= dNdetadphi.GetNbinsX(); etaBin++) {
+    eta = dNdetadphi.GetXaxis()->GetBinCenter(etaBin);
+    nCurBin = fCumuDiff->GetXaxis()->FindBin(eta);
+    // If we have moved to a new bin in the flow hist, and less than half the eta
+    // region has been covered by it we cut it away.
+    if (nPrevBin && nCurBin != nPrevBin) {
+      if (nInBin <= nBins/2) {
+       for (Int_t pBin = 1; pBin <= fCumuDiff->GetNbinsY(); pBin++) {
+         fCumuDiff->SetBinContent(nPrevBin, pBin, 0);
+         fCumuDiff->SetBinError(nPrevBin, pBin, 0);
+       }
+      }
+      nInBin = 0;
+      nPrevBin = nCurBin;
+      runAvg = 0;
+      nInAvg = 0;
+      max = 0;
+    }
+    Bool_t data = kFALSE;
+    for (Int_t phiBin = 1; phiBin <= dNdetadphi.GetNbinsY(); phiBin++) {
+      phi = dNdetadphi.GetYaxis()->GetBinCenter(phiBin);
+      weight = dNdetadphi.GetBinContent(etaBin, phiBin);
+      if (!weight) continue;
+      if (!data) data = kTRUE;
+      // We calculate the running average Nch per. bin
+      runAvg *= nInAvg;
+      runAvg += weight;
+      nInAvg++;
+      runAvg /= nInAvg;
+      // if the current bin has more than write the avg we count a bad bin
+      if (weight > max) max = weight;
+
+      dQnRe = weight*TMath::Cos(fMoment*phi);
+      dQnIm = weight*TMath::Sin(fMoment*phi);
+      dQ2nRe = weight*TMath::Cos(2.*fMoment*phi);
+      dQ2nIm = weight*TMath::Sin(2.*fMoment*phi);
+
+      fCumuRef->Fill(eta, kHmult, weight);
+      fCumuRef->Fill(eta, kHQnRe, dQnRe);
+      fCumuRef->Fill(eta, kHQnIm, dQnIm);
+      fCumuRef->Fill(eta, kHQ2nRe, dQ2nRe);
+      fCumuRef->Fill(eta, kHQ2nIm, dQ2nIm);
+
+      fCumuDiff->Fill(eta, kHmult, weight);
+      fCumuDiff->Fill(eta, kHQnRe, dQnRe);
+      fCumuDiff->Fill(eta, kHQnIm, dQnIm);
+      fCumuDiff->Fill(eta, kHQ2nRe, dQ2nRe);
+      fCumuDiff->Fill(eta, kHQ2nIm, dQ2nIm);
+
+      // Fill acc. map
+      fdNdedpAcc->Fill(eta, phi, weight);
+    }
+    if (data) {
+      nInBin++;
+      if (max > 35) nBadBins++;
+    }
+    // If there are too many bad bins we throw the event away!
+    if (nBadBins > 2) return kFALSE;
+  }
+  return kTRUE;
+}
+//_____________________________________________________________________
+void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent) 
+{
+  // 
+  // Calculate the Q cumulant of order fMoment
+  //
+  // Parameters:
+  //  cent: Centrality of event
+  //
+
+  if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
+
   // We create the objects needed for the analysis
   Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0, mult = 0;
   Double_t pnRe = 0, p2nRe = 0, qnRe = 0, q2nRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, q2nIm = 0;
@@ -328,16 +658,16 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "",
   Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0;
   Double_t cosPsi1Phi2 = 0, cosPsi1Phi2Phi3m = 0, cosPsi1Phi2Phi3p = 0;
   Double_t sinPsi1Phi2 = 0, sinPsi1Phi2Phi3m = 0, sinPsi1Phi2Phi3p = 0;
-  Double_t phi = 0, eta = 0;
+  Double_t eta = 0;
   Double_t multi = 0, multp = 0, mp = 0, mq = 0;
   Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
   Int_t refEtaBin = 0;
-  Bool_t kEventCount = kFALSE;
-
+  Bool_t eventCount = kFALSE;
+  
   // We loop over the data 1 time!
-  for (Int_t etaBin = 1; etaBin <= dNdetadphi->GetNbinsX(); etaBin++) {
-    eta = dNdetadphi->GetXaxis()->GetBinCenter(etaBin);
-    refEtaBin = dNdphi->GetXaxis()->FindBin(eta);
+  for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
+    eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
+    refEtaBin = fCumuRef->GetXaxis()->FindBin(eta);
     // The values for each individual etaBin bins are reset
     mp = 0;
     pnRe = 0;
@@ -351,53 +681,43 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "",
     dQ2nRe = 0;
     dQ2nIm = 0;
 
-    for (Int_t phiBin = 1; phiBin <= dNdphi->GetNbinsY(); phiBin++) {
-      phi = dNdphi->GetYaxis()->GetBinCenter(phiBin);
-      
-      // Reference flow
-      multi = dNdphi->GetBinContent(refEtaBin, phiBin);
-      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);
-      mult += multi;
-      
-      // For each etaBin bin the necessary values for differential flow
-      // is calculated. Here is the loop over the phi's.
-      multp = dNdetadphi->GetBinContent(etaBin, phiBin);
-      mp += multp;
-      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);
-    }
-    if (mult == 0) continue; 
+    // Reference flow
+    multi = fCumuRef->GetBinContent(refEtaBin, kHmult);
+    dQnRe = fCumuRef->GetBinContent(refEtaBin, kHQnRe);
+    dQnIm = fCumuRef->GetBinContent(refEtaBin, kHQnIm);
+    dQ2nRe = fCumuRef->GetBinContent(refEtaBin, kHQ2nRe);
+    dQ2nIm = fCumuRef->GetBinContent(refEtaBin, kHQ2nIm);
+    mult += multi;
+    
+    // For each etaBin bin the necessary values for differential flow
+    // is calculated. Here is the loop over the phi's.
+    multp = fCumuDiff->GetBinContent(etaBin, kHmult);
+    pnRe = fCumuDiff->GetBinContent(etaBin, kHQnRe);
+    pnIm = fCumuDiff->GetBinContent(etaBin, kHQnIm);
+    p2nRe = fCumuDiff->GetBinContent(etaBin, kHQ2nRe);
+    p2nIm = fCumuDiff->GetBinContent(etaBin, kHQ2nIm);
+    mp += multp;
+    
+    if (mult <= 3) continue; 
 
-    if (!kEventCount) {
+    if (!eventCount) {
      // Count number of events
-      coord[0] = -7;
-      coord[3] = -0.5;
-      flowHist->Fill(coord[0], coord[2], coord[3], 1.);
-      kEventCount = kTRUE;
+      fCumuHist->Fill(-7., cent, -0.5, 1.);
+      eventCount = kTRUE;
     } 
-
+    if (mp == 0) continue; 
     // The reference flow is calculated 
-    coord[0] = eta;
     
     // 2-particle
     w2 = mult * (mult - 1.);
     two = dQnRe*dQnRe + dQnIm*dQnIm - mult;
-    coord[3] = kW2Two;
-    flowHist->Fill(coord[0], coord[2], coord[3], two);
-    coord[3] = kW2;
-    flowHist->Fill(coord[0], coord[2], coord[3], w2);
-
-    coord[3] = kQnRe;
-    flowHist->Fill(coord[0], coord[2], coord[3], dQnRe);
-    coord[3] = kQnIm;
-    flowHist->Fill(coord[0], coord[2], coord[3], dQnIm);
-    coord[3] = kM;
-    flowHist->Fill(coord[0], coord[2], coord[3], mult);
+    
+    fCumuHist->Fill(eta, cent, kW2Two, two);
+    fCumuHist->Fill(eta, cent, kW2, w2);
+
+    fCumuHist->Fill(eta, cent, kQnRe, dQnRe);
+    fCumuHist->Fill(eta, cent, kQnIm, dQnIm);
+    fCumuHist->Fill(eta, cent, kM, mult);
     
     // 4-particle
     w4 = mult * (mult - 1.) * (mult - 2.) * (mult - 3.);
@@ -407,10 +727,8 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "",
              -2.*(TMath::Power(dQnRe,2.)*dQ2nRe+2.*dQnRe*dQnIm*dQ2nIm-TMath::Power(dQnIm,2.)*dQ2nRe)
              +(TMath::Power(dQ2nRe,2.)+TMath::Power(dQ2nIm,2.));
 
-    coord[3] = kW4Four;
-    flowHist->Fill(coord[0], coord[2], coord[3], four);
-    coord[3] = kW4;
-    flowHist->Fill(coord[0], coord[2], coord[3], w4);
+    fCumuHist->Fill(eta, cent, kW4Four, four);
+    fCumuHist->Fill(eta, cent, kW4, w4);
 
     cosPhi1Phi2 = dQnRe*dQnRe - dQnIm*dQnIm - dQ2nRe;
     sinPhi1Phi2 = 2.*dQnRe*dQnIm - dQ2nIm;
@@ -419,18 +737,19 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "",
 
     sinPhi1Phi2Phi3m = -dQnIm*(TMath::Power(dQnRe,2)+TMath::Power(dQnIm,2))+dQnRe*dQ2nIm-dQnIm*dQ2nRe+2.*(mult-1)*dQnIm; 
 
-    coord[3] = kCosphi1phi2;
-    flowHist->Fill(coord[0], coord[2], coord[3], cosPhi1Phi2);
-    coord[3] = kSinphi1phi2;
-    flowHist->Fill(coord[0], coord[2], coord[3], sinPhi1Phi2);
-    coord[3] = kCosphi1phi2phi3m;
-    flowHist->Fill(coord[0], coord[2], coord[3], cosPhi1Phi2Phi3m);
-    coord[3] = kSinphi1phi2phi3m;
-    flowHist->Fill(coord[0], coord[2], coord[3], sinPhi1Phi2Phi3m);
-    coord[3] = kMm1m2;
-    flowHist->Fill(coord[0], coord[2], coord[3], mult*(mult-1.)*(mult-2.));
-
-    // Differential flow calculations for each etaBin bin is done:
+    fCumuHist->Fill(eta, cent, kCosphi1phi2, cosPhi1Phi2);
+    fCumuHist->Fill(eta, cent, kSinphi1phi2, sinPhi1Phi2);
+    fCumuHist->Fill(eta, cent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
+    fCumuHist->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
+    fCumuHist->Fill(eta, cent, kMm1m2, mult*(mult-1.)*(mult-2.));
+
+    // Diagnostics are filled
+    fHistTwoCorr->Fill(eta, two, cent);
+    fHistW2->Fill(eta, w2, cent);
+    fHistFourCorr->Fill(eta, four, cent);
+    fHistW4->Fill(eta, w4, cent);
+
+    // Differential flow calculations for each eta bin bin is done:
     mq = mp;
     qnRe = pnRe;
     qnIm = pnIm;
@@ -441,17 +760,12 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "",
     w2p = mp * mult - mq;
     twoPrime = pnRe*dQnRe + pnIm*dQnIm - mq;
     
-    coord[3] = kw2two;
-    flowHist->Fill(coord[0], coord[2], coord[3], twoPrime);
-    coord[3] = kw2;
-    flowHist->Fill(coord[0], coord[2], coord[3], w2p);
-
-    coord[3] = kpnRe;
-    flowHist->Fill(coord[0], coord[2], coord[3], pnRe);
-    coord[3] = kpnIm;
-    flowHist->Fill(coord[0], coord[2], coord[3], pnIm);
-    coord[3] = kmp;
-    flowHist->Fill(coord[0], coord[2], coord[3], mp);
+    fCumuHist->Fill(eta, cent, kw2two, twoPrime);
+    fCumuHist->Fill(eta, cent, kw2, w2p);
+
+    fCumuHist->Fill(eta, cent, kpnRe, pnRe);
+    fCumuHist->Fill(eta, cent, kpnIm, pnIm);
+    fCumuHist->Fill(eta, cent, kmp, mp);
 
     // 4-particle differential flow
     w4p = (mp * mult - 3.*mq)*(mult - 1.)*(mult - 2.);
@@ -468,11 +782,9 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "",
                       + 2.*(pnRe*dQnRe+pnIm*dQnIm)                       
                       + 2.*mq*mult                      
                       - 6.*mq; 
-   
-    coord[3] = kw4four;
-    flowHist->Fill(coord[0], coord[2], coord[3], fourPrime);
-    coord[3] = kw4;
-    flowHist->Fill(coord[0], coord[2], coord[3], w4p);
+
+    fCumuHist->Fill(eta, cent, kw4four, fourPrime);
+    fCumuHist->Fill(eta, cent, kw4, w4p);
 
     cosPsi1Phi2 = pnRe*dQnRe - pnIm*dQnIm - q2nRe;
     sinPsi1Phi2 = pnRe*dQnIm + pnIm*dQnRe - q2nIm;
@@ -493,43 +805,51 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "",
                           - 1.*(pnIm*dQ2nRe-pnRe*dQ2nIm)
                           + 2.*mq*dQnIm-2.*qnIm;
 
-    coord[3] = kCospsi1phi2;
-    flowHist->Fill(coord[0], coord[2], coord[3], cosPsi1Phi2);
-    coord[3] = kSinpsi1phi2;
-    flowHist->Fill(coord[0], coord[2], coord[3], sinPsi1Phi2);
-    coord[3] = kCospsi1phi2phi3m;
-    flowHist->Fill(coord[0], coord[2], coord[3], cosPsi1Phi2Phi3m);
-    coord[3] = kSinpsi1phi2phi3m;
-    flowHist->Fill(coord[0], coord[2], coord[3], sinPsi1Phi2Phi3m);
-    coord[3] = kmpmq;
-    flowHist->Fill(coord[0], coord[2], coord[3], (mp*mult-2.*mq)*(mult-1.));
-    coord[3] = kCospsi1phi2phi3p;
-    flowHist->Fill(coord[0], coord[2], coord[3], cosPsi1Phi2Phi3p);
-    coord[3] = kSinpsi1phi2phi3p;
-    flowHist->Fill(coord[0], coord[2], coord[3], sinPsi1Phi2Phi3p); 
-
+    fCumuHist->Fill(eta, cent, kCospsi1phi2, cosPsi1Phi2);
+    fCumuHist->Fill(eta, cent, kSinpsi1phi2, sinPsi1Phi2);
+    fCumuHist->Fill(eta, cent, kCospsi1phi2phi3m, cosPsi1Phi2Phi3m);
+    fCumuHist->Fill(eta, cent, kSinpsi1phi2phi3m, sinPsi1Phi2Phi3m);
+    fCumuHist->Fill(eta, cent, kmpmq, (mp*mult-2.*mq)*(mult-1.));
+    fCumuHist->Fill(eta, cent, kCospsi1phi2phi3p, cosPsi1Phi2Phi3p);
+    fCumuHist->Fill(eta, cent, kSinpsi1phi2phi3p, sinPsi1Phi2Phi3p); 
   }
 
 }
 //_____________________________________________________________________
-void AliForwardFlowTaskQC::Terminate(Option_t */*option*/) 
+void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* outlist) 
 {
   // 
-  //  End of job
   //  Finalizes the Q cumulant calculations
   // 
   //  Parameters:
-  //    option Not used 
+  //   inlist: input sumlist
+  //   outlist: output result list 
   //
 
-  TH3D* cumulantsHist = 0;
-  TProfile* cumulant2diffHist = 0;
-  TProfile* cumulant4diffHist = 0;
-  TList* tList = 0;
-  TList* vList = 0;
+  // Re-find cumulants hist if Terminate is called separately
+  if (!fCumuHist) {
+    TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d", fType, fVzMin, fVzMax));
+    fCumuHist = (TH3D*)list->FindObject(Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax));
+  }
+
+  // Create result profiles
+  TProfile2D* cumu2 = (TProfile2D*)outlist->FindObject(Form("%sQC2_v%d_unCorr", fType, fMoment)); 
+  TProfile2D* cumu4 = (TProfile2D*)outlist->FindObject(Form("%sQC4_v%d_unCorr", fType, fMoment)); 
+  if (!cumu2) {
+    cumu2 = new TProfile2D(Form("%sQC2_v%d_unCorr", fType, fMoment),
+                           Form("%sQC2_v%d_unCorr", fType, fMoment),
+                           48, -6., 6., 100, 0., 100);
+    outlist->Add(cumu2);
+  }
+  if (!cumu4) {
+    cumu4 = new TProfile2D(Form("%sQC4_v%d_unCorr", fType, fMoment),
+                           Form("%sQC4_v%d_unCorr", fType, fMoment),
+                           48, -6., 6., 100, 0., 100);
+    outlist->Add(cumu4);
+  }
 
   // For flow calculations
-  Double_t two = 0, qc2 = 0, /* vnTwo = 0, */ four = 0, qc4 = 0 /*, vnFour = 0*/; 
+  Double_t two = 0, qc2 = 0, vnTwo = 0, four = 0, qc4 = 0, vnFour = 0; 
   Double_t twoPrime = 0, qc2Prime = 0, vnTwoDiff = 0, fourPrime = 0, qc4Prime = 0, vnFourDiff = 0;
   Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
   Double_t w2Two = 0, w2pTwoPrime = 0, w4Four = 0, w4pFourPrime = 0;
@@ -539,191 +859,126 @@ void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
   Double_t cosP1nPsi1M1nPhi2M1nPhi3 = 0, sinP1nPsi1M1nPhi2M1nPhi3 = 0, mpqMult = 0;
   Double_t cosP1nPsi1P1nPhi2M1nPhi3 = 0, sinP1nPsi1P1nPhi2M1nPhi3 = 0;
 
-  Int_t nLoops = (fMC ? 5 : 2);
-  TString type = "";
-
-  // Do a loop over the difference analysis types, calculating flow
-  // 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++) {
-
-    if (loop == 1) type = "FMD";
-    if (loop == 2) type = "SPD";
-    if (loop == 3) type = "MC";
-    if (loop == 4) type = "FMDTR";
-    if (loop == 5) type = "SPDTR";
-    
-    for (Int_t n = 1; n <= 6; n++) {
-      if (!fv[n]) continue;
-     
-      tList = (TList*)fOutputList->FindObject(type.Data());
-      vList = (TList*)tList->FindObject(Form("v%d", n));
-      TString tag = "";
-      // Centrality loop
-      for (Int_t c = 0; c < 100; c++) {
-
-        Double_t nEv = 0;
-        tag = Form("%d_%d", c, c+1);
-        cumulantsHist = (TH3D*)vList->FindObject(Form("hQ%dCumuHist%s_%s", n, type.Data(), tag.Data()));
-        
-        cumulant2diffHist = (TProfile*)vList->FindObject(Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()));
-        cumulant4diffHist = (TProfile*)vList->FindObject(Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()));
+  // Loop over cumulant histogram for final calculations   
+  // Centrality loop
+  for (Int_t c = 1; c <= 80; c++) {
+    Double_t nEv = 0;
+    // Eta loop
+    for (Int_t etaBin = 1; etaBin <= fCumuHist->GetNbinsX(); etaBin++) {
+      Double_t eta = fCumuHist->GetXaxis()->GetBinCenter(etaBin);
+      // 2-particle reference flow
+      w2Two = fCumuHist->GetBinContent(etaBin, c, kW2Two);
+      w2 = fCumuHist->GetBinContent(etaBin, c, kW2);
+      mult = fCumuHist->GetBinContent(etaBin, c, kM);
+      if (!w2 || !mult) continue;
+      cosP1nPhi = fCumuHist->GetBinContent(etaBin, c, kQnRe);
+      sinP1nPhi = fCumuHist->GetBinContent(etaBin, c, kQnIm);
         
-        for (Int_t vertexBin = 1; vertexBin <= cumulantsHist->GetNbinsY(); vertexBin++) {
-        
-          for (Int_t etaBin = 1; etaBin <= cumulantsHist->GetNbinsX(); etaBin++) {
-            Double_t eta = cumulantsHist->GetXaxis()->GetBinCenter(etaBin);
-            // 2-particle reference flow
-            w2Two = cumulantsHist->GetBinContent(etaBin, vertexBin, kW2Two);
-            if (!w2Two) continue;
-            w2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kW2);
-            cosP1nPhi = cumulantsHist->GetBinContent(etaBin, vertexBin, kQnRe);
-            sinP1nPhi = cumulantsHist->GetBinContent(etaBin, vertexBin, kQnIm);
-            mult = cumulantsHist->GetBinContent(etaBin, vertexBin, kM);
-              
-            cosP1nPhi /= mult;
-            sinP1nPhi /= mult;
-            two = w2Two / w2;
-            qc2 = two - TMath::Power(cosP1nPhi, 2) - TMath::Power(sinP1nPhi, 2);
-            if (qc2 <= 0) continue;
-            // vnTwo = TMath::Sqrt(qc2);
-       //     if (!TMath::IsNaN(vnTwo*mult)) 
-       //       cumulant2diffHist->Fill(eta, vnTwo, cumulantsHist->GetBinContent(0,vertexBin,0)); 
-
-            // 4-particle reference flow
-            w4Four = cumulantsHist->GetBinContent(etaBin, vertexBin, kW4Four);
-            w4 = cumulantsHist->GetBinContent(etaBin, vertexBin, kW4);
-            cosP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCosphi1phi2);
-            sinP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinphi1phi2);
-            cosP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCosphi1phi2phi3m);
-            sinP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinphi1phi2phi3m);
-            multm1m2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kMm1m2);
-
-            cosP1nPhi1P1nPhi2 /= w2;
-            sinP1nPhi1P1nPhi2 /= w2;
-            cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
-            sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
-            four = w4Four / w4;
-            qc4 = four-2.*TMath::Power(two,2.)
-               - 4.*cosP1nPhi*cosP1nPhi1M1nPhi2M1nPhi3
-               + 4.*sinP1nPhi*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
-               + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
-               + 8.*sinP1nPhi1P1nPhi2*sinP1nPhi*cosP1nPhi
-               + 8.*two*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.))
-               - 6.*TMath::Power((TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)),2.);
-            
-            if (qc4 >= 0) continue;
-            // vnFour = TMath::Power(-qc4, 0.25);
-       //     if (!TMath::IsNaN(vnFour*mult)) 
-       //         cumulant4diffHist->Fill(eta, vnFour, cumulantsHist->GetBinContent(0,vertexBin,0));
-  
-            // 2-particle differential flow
-            w2pTwoPrime = cumulantsHist->GetBinContent(etaBin, vertexBin, kw2two);
-            if (!w2pTwoPrime) continue;
-            w2p = cumulantsHist->GetBinContent(etaBin, vertexBin, kw2);
-            cosP1nPsi = cumulantsHist->GetBinContent(etaBin, vertexBin, kpnRe);
-            sinP1nPsi = cumulantsHist->GetBinContent(etaBin, vertexBin, kpnIm);
-            mp = cumulantsHist->GetBinContent(etaBin, vertexBin, kmp);
-
-            cosP1nPsi /= mp;
-            sinP1nPsi /= mp;
-            twoPrime = w2pTwoPrime / w2p;
-            qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi;
-
-            vnTwoDiff = qc2Prime / TMath::Sqrt(qc2);
-            if (!TMath::IsNaN(vnTwoDiff*mp) && vnTwoDiff > 0) cumulant2diffHist->Fill(eta, vnTwoDiff, cumulantsHist->GetBinContent(0,vertexBin,0));
-
-            // 4-particle differential flow
-            w4pFourPrime = cumulantsHist->GetBinContent(etaBin, vertexBin, kw4four);
-            w4p = cumulantsHist->GetBinContent(etaBin, vertexBin, kw4);
-            cosP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCospsi1phi2);
-            sinP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinpsi1phi2);
-            cosP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCospsi1phi2phi3m);
-            sinP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinpsi1phi2phi3m);
-            cosP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCospsi1phi2phi3p);
-            sinP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinpsi1phi2phi3p); 
-            mpqMult = cumulantsHist->GetBinContent(etaBin, vertexBin, kmpmq);
-
-            cosP1nPsi1P1nPhi2 /= w2p;
-            sinP1nPsi1P1nPhi2 /= w2p;
-            cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
-            sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
-            cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
-            sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
-            fourPrime = w4pFourPrime / 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*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
-                      + 4.*sinP1nPsi1P1nPhi2*cosP1nPhi*sinP1nPhi
-                      + 4.*twoPrime*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.))
-                      - 6.*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.)) 
-                      * (cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi)
-                      - 12.*cosP1nPhi*sinP1nPhi
-                      * (sinP1nPsi*cosP1nPhi+cosP1nPsi*sinP1nPhi);
-  
-            vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
-            if (!TMath::IsNaN(vnFourDiff*mp) && vnFourDiff > 0) cumulant4diffHist->Fill(eta, vnFourDiff, cumulantsHist->GetBinContent(0,vertexBin,0));
-          } // End of etaBin loop
-          nEv += cumulantsHist->GetBinContent(0,vertexBin,0);
-        } // End of vertexBin loop
-        // Number of events:
-        cumulant2diffHist->Fill(7., nEv);
-        cumulant4diffHist->Fill(7., nEv);
-  
-      } // End of centrality loop
-    } // End of harmonics loop
-  }  // End of type loop
-
-  PostData(1, fOutputList);
+      cosP1nPhi /= mult;
+      sinP1nPhi /= mult;
+      two = w2Two / w2;
+      qc2 = two - TMath::Power(cosP1nPhi, 2) - TMath::Power(sinP1nPhi, 2);
+      if (qc2 <= 0) continue;
+      vnTwo = TMath::Sqrt(qc2);
+ //     if (!TMath::IsNaN(vnTwo*mult)) 
+ //       cumu2->Fill(eta, vnTwo, fCumuHist->GetBinContent(0,c,0)); 
+
+      // 4-particle reference flow
+      w4Four = fCumuHist->GetBinContent(etaBin, c, kW4Four);
+      w4 = fCumuHist->GetBinContent(etaBin, c, kW4);
+      multm1m2 = fCumuHist->GetBinContent(etaBin, c, kMm1m2);
+      if (!w4 || !multm1m2) continue;
+      cosP1nPhi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kCosphi1phi2);
+      sinP1nPhi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kSinphi1phi2);
+      cosP1nPhi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kCosphi1phi2phi3m);
+      sinP1nPhi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kSinphi1phi2phi3m);
+
+      cosP1nPhi1P1nPhi2 /= w2;
+      sinP1nPhi1P1nPhi2 /= w2;
+      cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
+      sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
+      four = w4Four / w4;
+      qc4 = four-2.*TMath::Power(two,2.)
+         - 4.*cosP1nPhi*cosP1nPhi1M1nPhi2M1nPhi3
+         + 4.*sinP1nPhi*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
+         + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
+         + 8.*sinP1nPhi1P1nPhi2*sinP1nPhi*cosP1nPhi
+         + 8.*two*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.))
+         - 6.*TMath::Power((TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)),2.);
+      
+      if (qc4 >= 0) continue;
+      vnFour = TMath::Power(-qc4, 0.25);
+ //     if (!TMath::IsNaN(vnFour*mult)) 
+ //         cumu4->Fill(eta, vnFour, fCumuHist->GetBinContent(0,c,0));
+
+      // 2-particle differential flow
+      w2pTwoPrime = fCumuHist->GetBinContent(etaBin, c, kw2two);
+      w2p = fCumuHist->GetBinContent(etaBin, c, kw2);
+      mp = fCumuHist->GetBinContent(etaBin, c, kmp);
+      if (!w2p || !mp) continue;
+      cosP1nPsi = fCumuHist->GetBinContent(etaBin, c, kpnRe);
+      sinP1nPsi = fCumuHist->GetBinContent(etaBin, c, kpnIm);
+
+      cosP1nPsi /= mp;
+      sinP1nPsi /= mp;
+      twoPrime = w2pTwoPrime / w2p;
+      qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi;
+
+      vnTwoDiff = qc2Prime / TMath::Sqrt(qc2);
+      if (!TMath::IsNaN(vnTwoDiff*mp)) cumu2->Fill(eta, (Double_t)c-1., vnTwoDiff, fCumuHist->GetBinContent(0,c,0));
+
+      // 4-particle differential flow
+      w4pFourPrime = fCumuHist->GetBinContent(etaBin, c, kw4four);
+      w4p = fCumuHist->GetBinContent(etaBin, c, kw4);
+      mpqMult = fCumuHist->GetBinContent(etaBin, c, kmpmq);
+      if (!w4p || !mpqMult) continue;
+      cosP1nPsi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kCospsi1phi2);
+      sinP1nPsi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kSinpsi1phi2);
+      cosP1nPsi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kCospsi1phi2phi3m);
+      sinP1nPsi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kSinpsi1phi2phi3m);
+      cosP1nPsi1P1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kCospsi1phi2phi3p);
+      sinP1nPsi1P1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kSinpsi1phi2phi3p); 
+      
+      cosP1nPsi1P1nPhi2 /= w2p;
+      sinP1nPsi1P1nPhi2 /= w2p;
+      cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
+      sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
+      cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
+      sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
+
+      fourPrime = w4pFourPrime / 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*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
+                + 4.*sinP1nPsi1P1nPhi2*cosP1nPhi*sinP1nPhi
+                + 4.*twoPrime*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.))
+                - 6.*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.)) 
+                * (cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi)
+                - 12.*cosP1nPhi*sinP1nPhi
+                * (sinP1nPsi*cosP1nPhi+cosP1nPsi*sinP1nPhi);
+
+      vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
+      if (!TMath::IsNaN(vnFourDiff*mp)) cumu4->Fill(eta, (Double_t)c-1., vnFourDiff, fCumuHist->GetBinContent(0,c,0));
+    } // End of eta loop
+    // Number of events:
+    nEv += fCumuHist->GetBinContent(0,c,0);
+    cumu2->Fill(7., (Double_t)c-1., nEv);
+    cumu4->Fill(7., (Double_t)c-1., nEv);
+  } // End of centrality loop
   
   return;
 }
-// _____________________________________________________________________
-void AliForwardFlowTaskQC::ProcessPrimary() 
-{
-  //
-  // If fMC == kTRUE this function takes care of organizing the input 
-  // Monte Carlo data and histograms so AliForwardFlowTaskQC::QCumulants 
-  // can be run on it.
-  //
-
-  if (fFlowUtil->LoopAODFMDtrrefHits(fAOD)) {
-   // Run analysis on TrackRefs
-    for (Int_t n = 1; n <= 6; n++) {
-      if (fv[n])
-        CumulantsMethod("FMDTR", n);
-    }
-  }
-
-  if (fFlowUtil->LoopAODSPDtrrefHits(fAOD)) {
-    // Run analysis on SPD TrackRefs
-    for (Int_t n = 1; n <= 6; n++) {
-      if (fv[n])
-        CumulantsMethod("SPDTR", n);
-    }
-  }
-
-  if (fFlowUtil->LoopAODmc(fAOD, fAddFlow, fAddType, fAddOrder)) {
-    // Run analysis on MC truth
-    for (Int_t n = 1; n <= 6; n++) {
-      if (fv[n])
-        CumulantsMethod("MC", n);
-    }
-  }
-
-}
 //_____________________________________________________________________
 //
 //
index b734aa0..09d652c 100644 (file)
@@ -4,18 +4,20 @@
 #ifndef ALIFORWARDFLOWTASKQC_H
 #define ALIFORWARDFLOWTASKQC_H
 /**
- * @file   AliForwardFlowTaskQC.h
- * @author Alexander Hansen alexander.hansen@cern.ch
- * @date   Fri Mar 25 13:53:00 2011
+ * @file AliForwardFlowTaskQC.h
+ * @author Alexander Hansen
+ * @date   Tue Feb 14 2012
  * 
- * @brief  
+ * @brief
  * 
  * 
  * @ingroup pwglf_forward_flow
  */
 #include "AliAnalysisTaskSE.h"
-#include "AliForwardFlowUtil.h"
-class AliAODEvent;
+class AliAODForwardMult;
+class TH1D;
+class TH2D;
+class TH3D;
 
  /**
  * @defgroup pwglf_forward_tasks_flow Flow tasks 
@@ -33,7 +35,6 @@ class AliAODEvent;
  * @ingroup pwglf_forward_tasks_flow
  * @ingroup pwglf_forward_flow
  *
- * @todo Add centrality stuff
  *
  */
 class AliForwardFlowTaskQC : public AliAnalysisTaskSE
@@ -59,14 +60,12 @@ public:
    */
   /** 
    * Create output objects 
-   * 
    */
   virtual void UserCreateOutputObjects();
   /**
    * Initialize the task
-   *
    */
-  virtual void Init() {}
+  virtual void Init() {} 
   /** 
    * Process each event 
    *
@@ -80,60 +79,152 @@ public:
    */
   virtual void Terminate(Option_t *option);
   /* @} */
-  /*
+  /**
    * Returns the outputlist
-   *
-   */
+   * 
+   * @return TList* 
+  */
   TList* GetOutputList() { return fOutputList; }
-  /* 
-   * Set Number of @f$ \eta@f$ bins to be used in flow analysis
-   *
-   */
-  void SetUseNEtaBins(Int_t nbins) { fEtaBins = nbins; }
-  /*
+ /**
+  * Check AODForwardMult object for trigger, vertex and centrality
+  * returns true if event is OK
+  * 
+  * @param const aodfm
+  * 
+  * @return Bool_t 
+  */
+  Bool_t AODCheck(const AliAODForwardMult* aodfm);
+   /*
    * Set which harmonics to calculate. @f$ v_{1}@f$ to @f$ v_{4}@f$ is
    * available and calculated as default
-   *
-   * @param v1  Do @f$ v_{1}@f$ 
-   * @param v2  Do @f$ v_{2}@f$ 
-   * @param v3  Do @f$ v_{3}@f$ 
-   * @param v4  Do @f$ v_{4}@f$ 
+   * 
+   * @param  v1 Do @f$ v_{1}$f$
+   * @param  v2 Do @f$ v_{2}$f$
+   * @param  v3 Do @f$ v_{3}$f$
+   * @param  v4 Do @f$ v_{4}$f$
+   * @param  v5 Do @f$ v_{5}$f$
+   * @param  v6 Do @f$ v_{6}$f$
+   * 
+   * @return void 
    */
   void SetDoHarmonics(Bool_t v1 = kTRUE, Bool_t v2 = kTRUE, 
                      Bool_t v3 = kTRUE, Bool_t v4 = kTRUE,
                      Bool_t v5 = kTRUE, Bool_t v6 = kTRUE) { 
     fv[1] = v1; fv[2] = v2; fv[3] = v3; fv[4] = v4; fv[5] = v5; fv[6] = v6;}
-  /*
-   * Set string to add flow to MC truth particles
-   *
-   * @param type String
-   */
-  void AddFlow(TString type = "") { fAddFlow = type; }
-  /*
-   * Set which function fAddFlow should use
-   *
-   * @param type of AddFlow 
-   */
-  void AddFlowType(Int_t number = 0) { fAddType = number; }
-  /*
-   * Set which order of flow to add
-   *
-   * @param order Flow order 
-   */
-  void AddFlowOrder(Int_t order = 2) { fAddOrder = order; }
-  /*
-   *
-   * Set MC input flag
-   *
-   * @param mc MC input
-   */
-  void SetMCinput(Bool_t mc = kTRUE) { fMC = mc; }
-  /*
-   * Set number of eta bins to be used in reference flow
-   *
-   * @param bins Ref Eta Bins
-   */
-  void SetRefEtaBins(Int_t bins) { fEtaRef = bins; } 
+ /**
+  * Nested class to handle cumulant calculations in vertex bins
+  */
+  class VertexBin : public TNamed
+  {
+    public:
+    /*
+     * Constructor
+     */
+    VertexBin();
+   /**
+    * Constructor
+    * 
+    * @param vLow Min vertex z-coordinate
+    * @param vHigh Max vertex z-coordinate
+    * @param moment Flow moment
+    * @param type Data type (FMD/SPD/FMDTR/SPDTR/MC)
+    */
+    VertexBin(const Int_t vLow, const Int_t vHigh, 
+              const Int_t moment, const Char_t* type);
+   /**
+    * Copy constructor 
+    * 
+    * @param o Object to copy from
+    * 
+    * @return VertexBin
+    */
+    VertexBin(const VertexBin& o);
+   /**
+    * Assignment operator 
+    * 
+    * @param VertexBin&
+    * 
+    * @return VertexBin& 
+    */
+    VertexBin& operator=(const VertexBin&);
+   /**
+    * Destructor 
+    */
+    ~VertexBin(){}
+   /**
+    * Add vertex bin output to list
+    * 
+    * @param list Histograms are added to this list
+    * 
+    * @return void 
+    */
+    virtual void AddOutput(TList* list);
+   /**
+    * Check if vertex vZ is in range for this bin
+    * 
+    * @param vZ z-coordinate of vertex to check
+    * 
+    * @return Bool_t 
+    */
+    Bool_t CheckVertex(Double_t vZ);
+  /**
+    * Fill reference and differential flow histograms for analysis
+    *
+    * @param dNdetadphi 2D data histogram
+    *
+    * @return false if bad event (det. hotspot)
+    */
+    Bool_t FillHists(TH2D dNdetadphi);
+   /**
+    * Do cumulants calculations for current event with 
+    * centrality cent
+    * 
+    * @param cent Event centrality
+    * 
+    * @return void 
+    */
+    void CumulantsAccumulate(Double_t cent);
+   /**
+    * Finish cumulants calculations. Takes input and
+    * output lists in case Terminate is called separately
+    * 
+    * @param inlist List with input histograms
+    * @param outlist List with output histograms
+    * 
+    * @return void 
+    */
+    void CumulantsTerminate(TList* inlist, TList* outlist);
+
+    protected:
+   /*
+    * Enumeration for ref/diff histograms
+    */
+    enum { kHmult = 1, kHQnRe, kHQnIm, kHQ2nRe, kHQ2nIm };
+   /*
+    * Enumeration for cumulant histogram
+    */
+    enum { kW2Two = 1, kW2, kW4Four, kW4, kQnRe, kQnIm, kM,
+       kCosphi1phi2, kSinphi1phi2, kCosphi1phi2phi3m, kSinphi1phi2phi3m, kMm1m2, 
+       kw2two, kw2, kw4four, kw4, kpnRe, kpnIm, kmp, 
+       kCospsi1phi2, kSinpsi1phi2, kCospsi1phi2phi3m, kSinpsi1phi2phi3m,
+       kmpmq, kCospsi1phi2phi3p, kSinpsi1phi2phi3p };
+
+    const Int_t   fMoment;        // flow moment 
+    const Int_t   fVzMin;         // z-vertex min
+    const Int_t   fVzMax;         // z-vertex max
+    const Char_t* fType;          // data type
+    TH2D*         fCumuRef;       // histogram for reference flow
+    TH2D*         fCumuDiff;      // histogram for differential flow
+    TH3D*         fCumuHist;      // histogram for cumulants calculations
+    TH3D*         fHistTwoCorr;   // Diagnostics histogram for <2>
+    TH3D*         fHistW2;        // Diagnostics histogram for w_<2>
+    TH3D*         fHistFourCorr;  // Diagnostics histogram for <4>
+    TH3D*         fHistW4;        // Diagnostics histogram for w_<4>
+    TH2D*         fdNdedpAcc;     // Diagnostics histogram to make acc. maps
+
+    ClassDef(VertexBin, 1); // object for cumulants ananlysis in FMD
+  };
+
 protected:
   /** 
    * Copy constructor 
@@ -146,40 +237,37 @@ protected:
    * 
    * @return Reference to this object 
    */
-  AliForwardFlowTaskQC& operator=(const AliForwardFlowTaskQC&) { return *this; }
+  AliForwardFlowTaskQC& operator=(const AliForwardFlowTaskQC&);
   /*
-   * if MC information is available do analysis on Monte Carlo truth
-   * and track references
-   *
+   * Initiate vertex bin objects
    */
-  void ProcessPrimary();
-  /**
-   * Calculate Q cumulant
-   * 
-   * @param type     Determines which histograms should be used
-   *                 - "FMD" = FMD data histograms
-   *                 - "SPD" = SPD data histograms
-   *                 - "TR" = track reference histograms
-   *                 - "MC" = MC truth histograms
-   * @param harmonic Which harmonic to calculate
+  virtual void InitVertexBins();
+  /*
+   * Initiate diagnostics histograms
+   */
+  virtual void InitHists();
+  /*
+   * Analyze event
+   */
+  virtual Bool_t Analyze();
+  /*
+   * Finalize analysis
    */
-  void CumulantsMethod(TString type, Int_t harmonic);
+  virtual void Finalize();
 
-  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
-  Int_t          fEtaRef;       //  Number of eta bins for reference flow
-  Bool_t         fv[7];         //  Calculate v_{n} flag
-  TString        fAddFlow;     //  Add flow string
-  Int_t          fAddType;     //  Add flow type #
-  Int_t          fAddOrder;    //  Add flow order
-  Float_t       fZvertex;      //  Z vertex bin
-  Double_t       fCent;         //  Centrality
-  
+  TList          fBinsFMD;        //  list with FMD VertexBin objects 
+  TList          fBinsSPD;        //  list with SPD VertexBin objects
+  TList*         fSumList;        //  sum list
+  TList*         fOutputList;     //  Output list
+  AliAODEvent*   fAOD;            //  AOD event
+  Bool_t         fv[7];           //  Calculate v_{n} flag
+  Float_t       fZvertex;        //  Z vertex bin
+  Double_t       fCent;           //  Centrality
+  TH1D*          fHistCent;       //  Diagnostics hist for centrality
+  TH1D*          fHistVertexSel;  //  Diagnostics hist for selected vertices
+  TH1D*          fHistVertexAll;  //  Diagnostics hist for all vertices
 
-  ClassDef(AliForwardFlowTaskQC, 2); // Analysis task for FMD analysis
+  ClassDef(AliForwardFlowTaskQC, 1); // Analysis task for FMD analysis
 };
  
 #endif
diff --git a/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.cxx b/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.cxx
deleted file mode 100644 (file)
index c991f03..0000000
+++ /dev/null
@@ -1,433 +0,0 @@
-//
-// 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"
-#include "TProfile2D.h"
-#include "TGraph.h"
-#include "TF1.h"
-
-ClassImp(AliForwardFlowUtil)
-
-//_____________________________________________________________________
-AliForwardFlowUtil::AliForwardFlowUtil() : 
-  fList(0),
-  fCent(-1),
-  fVertex(1111),
-  fAliceCent4th(),
-  fAlicePt2nd4050(),
-  fAlicePt4th3040(),
-  fAlicePt4th4050(),
-  fImpactParToCent()
-   {} 
-  //
-  // Default Constructor
-  //
-//_____________________________________________________________________
-AliForwardFlowUtil::AliForwardFlowUtil(TList* outputList) :
-  fList(0),
-  fCent(-1),
-  fVertex(1111),
-  fAliceCent4th(),
-  fAlicePt2nd4050(),
-  fAlicePt4th3040(),
-  fAlicePt4th4050(),
-  fImpactParToCent()
-{ 
-  // 
-  // Constructor
-  //
-  // Parameters:
-  // TList: list containing histograms for flow analysis
-  //
-  fList = outputList;
-
-  Double_t xCumulant4thTPCrefMultTPConlyAll[] = {2.5,7.5,15,25,35,45,55,65};
-  Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,0.055818,0.073137,0.083898,0.086690,0.082040,0.077777};
-  Int_t nPointsCumulant4thTPCrefMultTPConlyAll = sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
-  fAliceCent4th = new TGraph(nPointsCumulant4thTPCrefMultTPConlyAll,xCumulant4thTPCrefMultTPConlyAll,
-                                                        yCumulant4thTPCrefMultTPConlyAll);
-  Double_t xCumulant2nd4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
-  1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
-  Double_t yCumulant2nd4050ALICE[] = {0.000000,0.043400,0.059911,0.073516,0.089756,0.105486,0.117391,0.128199,0.138013,
-  0.158271,0.177726,0.196383,0.208277,0.216648,0.242954,0.249961,0.240131,0.269006,0.207796};
-  Int_t nPointsCumulant2nd4050ALICE = sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);                                      
-  fAlicePt2nd4050 = new TGraph(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
-
-  Double_t xCumulant4th3040ALICE[] = {0.00000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
-  1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000,
-  5.500000,7.000000,9.000000};
-  Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,0.048566,0.061083,0.070910,0.078831,0.091396,0.102026,0.109691,
-  0.124449,0.139819,0.155561,0.165701,0.173678,0.191149,0.202015,0.204540,0.212560,0.195885,
-  0.000000,0.000000,0.000000};
-  Int_t nPointsCumulant4th3040ALICE = sizeof(xCumulant4th3040ALICE)/sizeof(Double_t);                                      
-  fAlicePt4th3040 = new TGraph(nPointsCumulant4th3040ALICE,xCumulant4th3040ALICE,yCumulant4th3040ALICE);
-
-  Double_t xCumulant4th4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
-  1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
-  Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,0.049824,0.066662,0.075856,0.081583,0.099778,0.104674,0.118545,
-  0.131874,0.152959,0.155348,0.169751,0.179052,0.178532,0.198851,0.185737,0.239901,0.186098};
-  Int_t nPointsCumulant4th4050ALICE = sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);   
-  fAlicePt4th4050 = new TGraph(nPointsCumulant4th4050ALICE, xCumulant4th4050ALICE, yCumulant4th4050ALICE);
-
-  Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,11.565,12.575,13.515,16.679};
-  Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
-
-  Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
-  fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
-
-}
-//_____________________________________________________________________
-Bool_t AliForwardFlowUtil::AODCheck(const AliAODForwardMult* aodfm) 
-{
-  // 
-  // Function to check that and AOD event meets the cuts
-  //
-  // Parameters: 
-  //  AliAODForwardMult: forward mult object with trigger and vertex info
-  //
-  fCent = -1;
-  fVertex = 1111;
-
-  if (!aodfm->IsTriggerBits(AliAODForwardMult::kOffline)) return kFALSE;
-  fCent = (Double_t)aodfm->GetCentrality();
-  if (0. >= fCent || fCent >= 100.) return kFALSE;
-  TH1D* vertex = (TH1D*)fList->FindObject("VertexAll");
-  fVertex = aodfm->GetIpZ();
-  vertex->Fill(fVertex);
-  if (TMath::Abs(fVertex) >= 5.) return kFALSE;
-  return kTRUE;
-}
-//_____________________________________________________________________
-Bool_t AliForwardFlowUtil::LoopAODFMD(const AliAODEvent* aodevent)
-{
-  //
-  // 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();
-  TH2D* vertex = (TH2D*)fList->FindObject("CoverageVsVertex");
-  TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSEFMD");
-  TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEFMD");
-  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);
-    if (TMath::Abs(eta) < 1.75) continue;
-    for (Int_t phiBin = 0; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
-      phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
-      weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
-      if (phiBin == 0) {
-        vertex->Fill(eta, fVertex, weight);
-        continue;
-      }
-      dNdetadphi->Fill(eta, phi, weight);
-      dNdphi->Fill(eta, phi, weight);
-      dNdphi->Fill(-1.*eta, phi, weight);
-      mult += weight;
-    }
-  }
-//  fCent = GetCentFromMC(aodevent);
-//  fCent = 0.5;
-  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* vertex = (TH2D*)fList->FindObject("CoverageVsVertex");
-  TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSESPD");
-  TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSESPD");
-  dNdphi->Reset();
-  dNdetadphi->Reset();
-
-  Double_t eta = 0;
-  Double_t phi = 0;
-  Double_t weight = 0;
-  Double_t mult = 0;
-  for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
-    eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
-    if (TMath::Abs(eta) > 1.75) continue;
-    for (Int_t phiBin = 0; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
-      phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
-      weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
-      if (phiBin == 0) {
-        vertex->Fill(eta, fVertex, weight);
-        continue;
-      }
-      dNdetadphi->Fill(eta, phi, weight);
-      if (TMath::Abs(eta) < 0.5) continue;
-      dNdphi->Fill(eta, phi, weight);
-//      dNdphi->Fill(-1.*eta, phi, weight);
-      mult += weight;
-    }
-  }
-
-  return kTRUE;
-}
-//_____________________________________________________________________
-Bool_t AliForwardFlowUtil::LoopAODFMDtrrefHits(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;
-
-  TH2D hdNdetadphi = aodfmult->GetHistogram();
-  TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSEFMDTR");
-  TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEFMDTR");
-  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);
-    if (TMath::Abs(eta) < 1.75) continue;
-    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(eta, phi, weight);
-      dNdphi->Fill(-1.*eta, phi, weight);
-      mult += weight;
-    }
-  }
-
-  return kTRUE;
-}
-//_____________________________________________________________________
-Bool_t AliForwardFlowUtil::LoopAODSPDtrrefHits(const AliAODEvent* aodevent) const
-{
-  // 
-  // Loop over AliAODCentralMult object and fill histograms
-  // provided by flow task
-  //
-  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(aodevent->FindListObject("CentralClustersMC"));
-  if (!aodcmult) return kFALSE;
-
-  TH2D hdNdetadphi = aodcmult->GetHistogram();
-  TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSESPDTR");
-  TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSESPDTR");
-  dNdphi->Reset();
-  dNdetadphi->Reset();
-
-  Double_t eta = 0;
-  Double_t phi = 0;
-  Double_t weight = 0;
-  Double_t mult = 0;
-  for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
-    eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
-    if (TMath::Abs(eta) > 1.75) continue;
-    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(eta, phi, weight);
-      dNdphi->Fill(-1.*eta, phi, weight);
-      mult += weight;
-    }
-  }
-
-  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;
-  AliAODMCHeader* header = dynamic_cast<AliAODMCHeader*>(aodevent->FindListObject(AliAODMCHeader::StdBranchName()));
-  if (!header) {
-    AliWarning("No header file found.");
-  }
-  else {
-    rp = header->GetReactionPlaneAngle();
-  }
-
-  TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEMC");
-  TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSEMC");
-  TProfile2D* mcTruth = (TProfile2D*)fList->FindObject("pMCTruth");
-  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 = 1;
-    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;
-    Double_t pT = particle->Pt();
-//    if (pT <= 0.2 || pT >= 5.) continue;
-    Double_t eta = particle->Eta();
-    Double_t phi = particle->Phi();
-    if (TMath::Abs(eta) < 6.) {
-      // Add flow if it is in the argument
-      if (addFlow.Length() > 1) {
-        if (addFlow.Contains("pt"))
-          weight *= AddptFlow(pT, type);
-//          weight *= AddptFlow(pT, 1);
-        if (addFlow.Contains("pid"))
-          weight *= AddpidFlow(particle->PdgCode(), 1);
-        if (addFlow.Contains("eta"))
-          weight *= AddetaFlow(eta, type);
-//          weight *= AddetaFlow(eta, 2);
-        if (addFlow.Contains("cent")) 
-          weight *= fAliceCent4th->Eval(fCent)/fAliceCent4th->Eval(45);
-        
-        weight *= 20*2.*TMath::Cos((Double_t)order*(phi-rp)); 
-        weight += 1;
-      }
-      dNdphi->Fill(eta, phi, weight);
-      dNdphi->Fill(-1.*eta, phi, weight);
-      dNdetadphi->Fill(eta, phi, weight);
-      mcTruth->Fill(eta, fCent, weight*TMath::Cos(2.*(phi-rp)));
-      mult += weight;
-    }
-  }
-
-  return kTRUE;
-}
-//_____________________________________________________________________
-Double_t AliForwardFlowUtil::AddptFlow(Double_t pt = 0, Int_t type = 0) const 
-{
-  //
-  // Add pt dependent flow factor
-  //
-  Double_t weight = 0;
-
-  switch(type) 
-  {
-    case 1: weight = fAlicePt2nd4050->Eval(pt)*0.5+fAlicePt4th4050->Eval(pt)*0.5;
-            break;
-    case 2: weight = fAlicePt2nd4050->Eval(pt);
-            break;
-    case 3: weight = fAlicePt4th3040->Eval(pt);
-            break;
-    case 4: weight = fAlicePt4th4050->Eval(pt);
-            break;
-  } 
-  
-  return weight;
-}
-//_____________________________________________________________________
-Double_t AliForwardFlowUtil::AddpidFlow(Int_t id = 0, Int_t type = 0) const 
-{
-  //
-  // Add pid dependent flow factor 
-  //
-  Double_t weight = 0;
-
-  switch(type)
-  {
-    case 1: if (TMath::Abs(id) ==  211) // pion flow
-              weight = 1.;
-            else if (TMath::Abs(id) ==  2212) // proton flow
-              weight = 1.3;
-            else 
-              weight = 0.7;
-            break;
-    case 2: weight = 1.207;
-            break;
-  }
-
-  return weight;
-}
-//_____________________________________________________________________
-Double_t AliForwardFlowUtil::AddetaFlow(Double_t eta = 0, Int_t type = 0) const 
-{
-  //
-  // Add eta dependent flow factor 
-  //
-  Double_t weight = 0;
-
-  TF1 gaus = TF1("gaus", "gaus", -6, 6);
-
-  switch(type)
-  {
-     case 1: gaus.SetParameters(0.1, 0., 9);
-             break;
-     case 2: gaus.SetParameters(0.1, 0., 3);
-             break;
-     case 3: gaus.SetParameters(0.1, 0., 15);
-             break;
-  }
-
-  weight = gaus.Eval(eta);
-
-  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();
-
-  cent = fImpactParToCent->Eval(b);
-
-  return cent;
-}
-//_____________________________________________________________________
-//
-//
-// EOF
-
diff --git a/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.h b/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.h
deleted file mode 100644 (file)
index f69c26a..0000000
+++ /dev/null
@@ -1,176 +0,0 @@
-//
-// 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
-/**
- * @file   AliForwardFlowUtil.h
- * @author Alexander Hansen alexander.hansen@cern.ch
- * @date   Fri Mar 25 13:15:40 2011
- * 
- * @brief  
- * 
- * 
- * @ingroup pwglf_forward_flow
- */
-#include "TNamed.h"
-class AliAODForwardMult;
-class AliAODEvent;
-class TList;
-class TGraph;
-
-/**
- * 
- * Class used to handle the input from AODs and put it into histograms
- * the Forward Flow tasks can run on.
- *
- * @ingroup pwglf_forward_tasks_flow
- * @ingroup pwglf_forward_flow
- */
-class AliForwardFlowUtil : public TNamed
-{
-public:
-  /**
-   * Constructor
-   */
-  AliForwardFlowUtil();
-  /*
-   * Constructor
-   *
-   * @param l list of histograms for flow analysis
-   */
-  AliForwardFlowUtil(TList* l);
-  /**
-   * Check that AOD event meet trigger requirements
-   * 
-   * @param aodfm Forward multplicity AOD event structure 
-   * 
-   * @return true on success 
-   */
-  Bool_t AODCheck(const AliAODForwardMult* aodfm);
-  /**
-   * Loop over AliAODForwardMult object and fill flow histograms
-   * 
-   * @param AODevent AOD event structure 
-   * 
-   * @return true on success
-   */
-  Bool_t LoopAODFMD(const AliAODEvent* AODevent);
-  /*
-   * Loop over AliAODCentralMult object and fill flow histograms
-   * 
-   * @param AODevent AOD event structure 
-   * 
-   * @return true on success
-   */
-  Bool_t LoopAODSPD(const AliAODEvent* AODevent) const;
-  /**
-   * Loop over AliAODForwardMult object and fill flow histograms from
-   * track refs
-   * 
-   * @param AODevent AOD event structure 
-   * 
-   * @return true on success
-   */
-  Bool_t LoopAODFMDtrrefHits(const AliAODEvent* AODevent) const;
-  /**
-   * Loop over AliAODCentralMult object and fill flow histograms from
-   * track refs
-   * 
-   * @param AODevent AOD event structure 
-   * 
-   * @return true on success
-   */
-  Bool_t LoopAODSPDtrrefHits(const AliAODEvent* AODevent) const;
-  /**
-   * Loop over AliAODMCParticle branch object and fill flow histograms
-   * add flow if arguments are set
-   * 
-   * @param AODevent AOD event structure 
-   * @param addFlow  What to add flow to 
-   * @param type     Type of flow 
-   * @param order    Order of added flow 
-   * 
-   * @return true on success
-   */
-  Bool_t LoopAODmc(const AliAODEvent* AODevent, TString addFlow, 
-                  Int_t type, Int_t order) const;
- /**
-   * Get centrality from newest processed event
-   */
-  Double_t GetCentrality() const { return fCent; }
- /**
-   * Get z vertex coordinate from newest processed event
-   */
-  Float_t GetVertex() const { return fVertex; }
- /**
-   * Parametrize ALICE data points
-   */
-   
-protected:
-  /*
-   * Copy constructor
-   *
-   * @param o Object to copy from
-   */
-  AliForwardFlowUtil(const AliForwardFlowUtil& o) : TNamed(),
-                                                   fList(o.fList),
-                                                   fCent(o.fCent),
-                                                   fVertex(o.fVertex),
-                                                    fAliceCent4th(o.fAliceCent4th),
-                                                    fAlicePt2nd4050(o.fAlicePt2nd4050),
-                                                    fAlicePt4th3040(o.fAlicePt4th3040),
-                                                    fAlicePt4th4050(o.fAlicePt4th4050),
-                                                    fImpactParToCent(o.fImpactParToCent)
-                                                    {}
-  /** 
-   * Assignment operator 
-   * 
-   * @return Reference to this object 
-   */
-  AliForwardFlowUtil& operator=(const AliForwardFlowUtil&) { return *this; }
-  /**
-   * Add pt dependent flow factor
-   *
-   * @param Pt   @f$ p_T@f$
-   * @param type Type of flow 
-   */
-  Double_t AddptFlow(Double_t Pt, Int_t type) const;
-  /**
-   * Add pid dependent flow factor
-   *
-   * @param ID   Particle ID 
-   * @param type Type of flow
-   */
-  Double_t AddpidFlow(Int_t ID, Int_t type) const;
-  /**
-   * Add eta dependent flow factor
-   * 
-   * @param Eta  @f$\eta@f$ 
-   * @param type Type of flow 
-   */
-  Double_t AddetaFlow(Double_t Eta, Int_t type) const;
-  /**
-   * Get centrality form MC impact parameter
-   * 
-   * @param AODevent AOD event structure 
-   */
-  Double_t GetCentFromMC(const AliAODEvent* AODevent) const;
-
-  TList*       fList;   // List of flow histograms
-  Double_t     fCent;   // centrality
-  Float_t       fVertex; // z vertex coordinate
-  TGraph*       fAliceCent4th;
-  TGraph*       fAlicePt2nd4050;
-  TGraph*       fAlicePt4th3040;
-  TGraph*       fAlicePt4th4050;
-  TGraph*       fImpactParToCent;
-
-  ClassDef(AliForwardFlowUtil, 1); 
-};
-#endif
-// Local Variables:
-//   mode: C++ 
-// End:
diff --git a/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx b/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx
new file mode 100644 (file)
index 0000000..bfc407e
--- /dev/null
@@ -0,0 +1,441 @@
+//
+// Calculate flow on MC data in the forward and central regions using the Q cumulants method.
+//
+// Inputs:
+//  - AliAODEvent
+//
+// Outputs:
+//  - AnalysisResults.root
+//
+#include "AliForwardMCFlowTaskQC.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "TGraph.h"
+#include "TF1.h"
+#include "TProfile2D.h"
+#include "AliAODEvent.h"
+#include "AliAODForwardMult.h"
+#include "AliAODCentralMult.h"
+
+ClassImp(AliForwardMCFlowTaskQC)
+
+//_____________________________________________________________________
+AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC() :
+  AliForwardFlowTaskQC(), 
+  fBinsFMDTR(),          // List of FMDTR analysis objects
+  fBinsSPDTR(),          // List of SPDTR analysis objects
+  fBinsMC(),             // List of MC truth analysis objects
+  fdNdedpMC(),           // MC truth d^2N/detadphi histogram
+  fAliceCent4th(),       // Alice QC4 vs. centrality data points
+  fAlicePt2nd4050(),     // Alice QC2 vs. pT data points
+  fAlicePt4th3040(),     // Alice QC4 vs. pT data points
+  fAlicePt4th4050(),     // Alice QC4 vs. pT data points
+  fImpactParToCent(),    // Impact parameter to centrality graph
+  fAddFlow(0),           // Add flow to MC truth
+  fAddType(0),           // Add type of flow to MC truth
+  fAddOrder(0)           // Add order of flow to MC truth        
+   {} 
+  //
+  // Default Constructor
+  //
+//_____________________________________________________________________
+AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name) :
+  AliForwardFlowTaskQC(name),
+  fBinsFMDTR(),          // List of FMDTR analysis objects
+  fBinsSPDTR(),          // List of SPDTR analysis objects
+  fBinsMC(),             // List of MC truth analysis objects
+  fdNdedpMC(),           // MC truth d^2N/detadphi histogram
+  fAliceCent4th(),       // Alice QC4 vs. centrality data points
+  fAlicePt2nd4050(),     // Alice QC2 vs. pT data points
+  fAlicePt4th3040(),     // Alice QC4 vs. pT data points
+  fAlicePt4th4050(),     // Alice QC4 vs. pT data points
+  fImpactParToCent(),    // Impact parameter to centrality graph
+  fAddFlow(0),           // Add flow to MC truth
+  fAddType(0),           // Add type of flow to MC truth
+  fAddOrder(0)           // Add order of flow to MC truth        
+{ 
+  // 
+  // Constructor
+  // Parameters:
+  //  name: Name of task
+  //
+  fdNdedpMC = TH2D("fdNdedpMC", "fdNdedpMC", 48, -6., 6., 200, 0., 2.*TMath::Pi());
+  fdNdedpMC.Sumw2();
+  
+  // Add parametrizations of ALICE data
+  Double_t xCumulant4thTPCrefMultTPConlyAll[] = {2.5,7.5,15,25,35,45,55,65};
+  Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,0.055818,0.073137,0.083898,0.086690,0.082040,0.077777};
+  Int_t nPointsCumulant4thTPCrefMultTPConlyAll = sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
+  fAliceCent4th = new TGraph(nPointsCumulant4thTPCrefMultTPConlyAll,xCumulant4thTPCrefMultTPConlyAll,
+                                                        yCumulant4thTPCrefMultTPConlyAll);
+  Double_t xCumulant2nd4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
+  1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
+  Double_t yCumulant2nd4050ALICE[] = {0.000000,0.043400,0.059911,0.073516,0.089756,0.105486,0.117391,0.128199,0.138013,
+  0.158271,0.177726,0.196383,0.208277,0.216648,0.242954,0.249961,0.240131,0.269006,0.207796};
+  Int_t nPointsCumulant2nd4050ALICE = sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);                                      
+  fAlicePt2nd4050 = new TGraph(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
+
+  Double_t xCumulant4th3040ALICE[] = {0.00000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
+  1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000,
+  5.500000,7.000000,9.000000};
+  Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,0.048566,0.061083,0.070910,0.078831,0.091396,0.102026,0.109691,
+  0.124449,0.139819,0.155561,0.165701,0.173678,0.191149,0.202015,0.204540,0.212560,0.195885,
+  0.000000,0.000000,0.000000};
+  Int_t nPointsCumulant4th3040ALICE = sizeof(xCumulant4th3040ALICE)/sizeof(Double_t);                                      
+  fAlicePt4th3040 = new TGraph(nPointsCumulant4th3040ALICE,xCumulant4th3040ALICE,yCumulant4th3040ALICE);
+
+  Double_t xCumulant4th4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
+  1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
+  Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,0.049824,0.066662,0.075856,0.081583,0.099778,0.104674,0.118545,
+  0.131874,0.152959,0.155348,0.169751,0.179052,0.178532,0.198851,0.185737,0.239901,0.186098};
+  Int_t nPointsCumulant4th4050ALICE = sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);   
+  fAlicePt4th4050 = new TGraph(nPointsCumulant4th4050ALICE, xCumulant4th4050ALICE, yCumulant4th4050ALICE);
+
+  Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,11.565,12.575,13.515,16.679};
+  Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
+
+  Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
+  fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
+
+}
+//_____________________________________________________________________
+AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o) :
+  AliForwardFlowTaskQC(o), 
+  fBinsFMDTR(),                          // List of FMDTR analysis objects
+  fBinsSPDTR(),                          // List of SPDTR analysis objects
+  fBinsMC(),                             // List of MC truth analysis objects
+  fdNdedpMC(o.fdNdedpMC),                // MC truth d^2N/detadphi histogram
+  fAliceCent4th(o.fAliceCent4th),        // Alice QC4 vs. centrality data points
+  fAlicePt2nd4050(o.fAlicePt2nd4050),    // Alice QC2 vs. pT data points
+  fAlicePt4th3040(o.fAlicePt4th3040),    // Alice QC4 vs. pT data points
+  fAlicePt4th4050(o.fAlicePt4th4050),    // Alice QC4 vs. pT data points
+  fImpactParToCent(o.fImpactParToCent),  // Impact parameter to centrality graph
+  fAddFlow(o.fAddFlow),                  // Add flow to MC truth
+  fAddType(o.fAddType),                  // Add type of flow to MC truth
+  fAddOrder(o.fAddOrder)                 // Add order of flow to MC truth        
+   {} 
+  //
+  // Copy Constructor
+  //
+//_____________________________________________________________________
+AliForwardMCFlowTaskQC&
+AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
+{
+  // 
+  // Assignment operator
+  // Parameters:
+  //    o Object to copy from 
+  //
+  fdNdedpMC        = o.fdNdedpMC;
+  fAliceCent4th    = o.fAliceCent4th;
+  fAlicePt2nd4050  = o.fAlicePt2nd4050;
+  fAlicePt4th3040  = o.fAlicePt4th3040;
+  fAlicePt4th4050  = o.fAlicePt4th4050;
+  fImpactParToCent = o.fImpactParToCent;
+  fAddFlow         = o.fAddFlow;
+  fAddType         = o.fAddType;
+  fAddOrder        = o.fAddOrder;
+  return *this;
+}
+//_____________________________________________________________________
+void AliForwardMCFlowTaskQC::InitVertexBins()
+{
+  //
+  // Initiate VertexBin lists
+  //
+  AliForwardFlowTaskQC::InitVertexBins();
+
+  for(Int_t n = 1; n <= 6; n++) {
+    if (!fv[n]) continue;
+    for (Int_t v = -10; v < 10; v++) {
+      fBinsFMDTR.Add(new VertexBin(v, v+1, n, "FMDTR"));
+      fBinsSPDTR.Add(new VertexBin(v, v+1, n, "SPDTR"));
+      fBinsMC.Add(new VertexBin(v, v+1, n, "MC"));
+    }
+  }
+
+}
+//_____________________________________________________________________
+void AliForwardMCFlowTaskQC::InitHists()
+{
+  //
+  // Initiate diagnostics hists and add to outputlist
+  //
+  AliForwardFlowTaskQC::InitHists();
+
+  TIter nextFMDTR(&fBinsFMDTR);
+  VertexBin* bin = 0;
+  while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
+    bin->AddOutput(fSumList);
+  }
+  TIter nextSPDTR(&fBinsSPDTR);
+  while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
+    bin->AddOutput(fSumList);
+  }
+  TIter nextMC(&fBinsMC);
+  while ((bin = static_cast<VertexBin*>(nextMC()))) {
+    bin->AddOutput(fSumList);
+  }
+
+}
+//_____________________________________________________________________
+Bool_t AliForwardMCFlowTaskQC::Analyze() 
+{
+  // 
+  // Load FMD and SPD MC objects from aod tree and call the cumulants 
+  // calculation for the correct vertexbin
+  //
+  if (!AliForwardFlowTaskQC::Analyze()) return kFALSE;
+
+  // Run analysis on trackrefs from FMD and SPD
+  AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
+  if (!aodfmult) return kFALSE;
+  TH2D fmdTRdNdetadphi = aodfmult->GetHistogram();
+
+  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
+  if (!aodcmult) return kFALSE;
+  TH2D spdTRdNdetadphi = aodcmult->GetHistogram();
+  
+  TIter nextFMDTR(&fBinsFMDTR);
+  VertexBin* bin = 0;
+  while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
+    if (bin->CheckVertex(fZvertex)) {
+      if (!bin->FillHists(fmdTRdNdetadphi)) return kFALSE;
+      bin->CumulantsAccumulate(fCent);
+    }
+  }
+
+  TIter nextSPDTR(&fBinsSPDTR);
+  while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
+    if (bin->CheckVertex(fZvertex)) {
+      if (!bin->FillHists(spdTRdNdetadphi)) return kFALSE;
+      bin->CumulantsAccumulate(fCent);
+    }
+  }
+
+  // Run analysis on MC branch
+  if (!LoopAODMC()) return kFALSE;
+
+  TIter nextMC(&fBinsMC);
+  while ((bin = static_cast<VertexBin*>(nextMC()))) {
+    if (bin->CheckVertex(fZvertex)) {
+      if (!bin->FillHists(fdNdedpMC)) return kFALSE;
+      bin->CumulantsAccumulate(fCent);
+    }
+  }
+
+  return kTRUE;
+}
+//_____________________________________________________________________
+void AliForwardMCFlowTaskQC::Finalize() 
+{
+  //
+  // Finalize command, called by Terminate()
+  //
+  AliForwardFlowTaskQC::Finalize();
+
+  TIter nextFMDTR(&fBinsFMDTR);
+  VertexBin* bin = 0;
+  while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
+    bin->CumulantsTerminate(fSumList, fOutputList);
+  }
+  TIter nextSPDTR(&fBinsSPDTR);
+  while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
+    bin->CumulantsTerminate(fSumList, fOutputList);
+  }
+  TIter nextMC(&fBinsMC);
+  while ((bin = static_cast<VertexBin*>(nextMC()))) {
+    bin->CumulantsTerminate(fSumList, fOutputList);
+  }
+
+  TProfile2D* fmdHist = 0;
+  TProfile2D* spdHist = 0;
+  TProfile2D* mcHist = 0;
+  TList* correctionList = new TList();
+  correctionList->SetName("CorrectionList");
+//  fOutputList->Add(correctionList);
+
+  for (Int_t i = 2; i <= 4; i += 2) {
+    for (Int_t n = 1; n <= 6; n++) {
+      if (!fv[n]) continue;
+      fmdHist = (TProfile2D*)fOutputList->FindObject(Form("FMDQC%d_v%d_unCorr", i, n))
+                 ->Clone(Form("FMDQC%d_v%d_Correction", i, n));
+      spdHist = (TProfile2D*)fOutputList->FindObject(Form("SPDQC%d_v%d_unCorr", i, n))
+                 ->Clone(Form("SPDQC%d_v%d_Correction", i, n));
+      mcHist = (TProfile2D*)fOutputList->FindObject(Form("MCQC%d_v%d_unCorr", i, n));
+     
+      if (!fmdHist || !spdHist || !mcHist) {
+       AliError(Form("Histogram missing, correction object not created for v%d", n));
+       continue;
+      }
+
+      fmdHist->Divide(mcHist);
+      spdHist->Divide(mcHist);
+      fmdHist->SetTitle(Form("FMD QC{%d} v_{%d} Correction Object", i, n));
+      fmdHist->SetTitle(Form("SPD QC{%d} v_{%d} Correction Object", i, n));
+
+//      correctionList->Add(fmdHist);
+//      correctionList->Add(spdHist);
+    }
+  }
+
+}
+//_____________________________________________________________________
+Bool_t AliForwardMCFlowTaskQC::LoopAODMC()  
+{
+  // 
+  // Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
+  // Add flow if set to do so in AddTask function
+  fdNdedpMC.Reset();
+
+  //retreive MC particles from event
+  TClonesArray* mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
+  if(!mcArray){
+//    AliWarning("No MC array found in AOD. Try making it again.");
+    return kFALSE;
+  }
+
+  Double_t rp = 0;
+  AliAODMCHeader* header = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::StdBranchName()));
+  if (!header) {
+    AliWarning("No header file found.");
+  }
+  else {
+    rp = header->GetReactionPlaneAngle();
+  }
+
+  Int_t ntracks = mcArray->GetEntriesFast();
+
+  // Track loop: chek how many particles will be accepted
+  Double_t weight = 0;
+  for (Int_t it = 0; it < ntracks; it++) {
+    weight = 1;
+    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;
+    Double_t pT = particle->Pt();
+    Double_t eta = particle->Eta();
+    Double_t phi = particle->Phi();
+    if (TMath::Abs(eta) < 6.) {
+      // Add flow if it is in the argument
+      if (fAddFlow.Length() > 1) {
+        if (fAddFlow.Contains("pt"))
+          weight *= AddptFlow(pT);
+        if (fAddFlow.Contains("pid"))
+          weight *= AddpidFlow(particle->PdgCode());
+        if (fAddFlow.Contains("eta"))
+          weight *= AddetaFlow(eta);
+        if (fAddFlow.Contains("cent")) 
+          weight *= fAliceCent4th->Eval(fCent)/fAliceCent4th->Eval(45);
+        
+        weight *= 20*2.*TMath::Cos((Double_t)fAddOrder*(phi-rp)); 
+        weight += 1;
+      }
+      fdNdedpMC.Fill(eta, phi, weight);
+    }
+  }
+
+  return kTRUE;
+}
+//_____________________________________________________________________
+Double_t AliForwardMCFlowTaskQC::AddptFlow(Double_t pt = 0) const 
+{
+  //
+  // Add pt dependent flow factor
+  // Parameters:
+  //  pt: pT parametrization to use
+  //
+  Double_t weight = 0;
+
+  switch(fAddType) 
+  {
+    case 1: weight = fAlicePt2nd4050->Eval(pt)*0.5+fAlicePt4th4050->Eval(pt)*0.5;
+            break;
+    case 2: weight = fAlicePt2nd4050->Eval(pt);
+            break;
+    case 3: weight = fAlicePt4th3040->Eval(pt);
+            break;
+    case 4: weight = fAlicePt4th4050->Eval(pt);
+            break;
+  } 
+  
+  return weight;
+}
+//_____________________________________________________________________
+Double_t AliForwardMCFlowTaskQC::AddpidFlow(Int_t id = 0) const 
+{
+  //
+  // Add pid dependent flow factor 
+  // Parameters:
+  //  id: choose PID dependent setup
+  //
+  Double_t weight = 0;
+
+  switch(fAddType)
+  {
+    case 1: if (TMath::Abs(id) ==  211) // pion flow
+              weight = 1.3;
+            else if (TMath::Abs(id) ==  2212) // proton flow
+              weight = 1.;
+            else 
+              weight = 0.7;
+            break;
+    case 2: weight = 1.207;
+            break;
+  }
+
+  return weight;
+}
+//_____________________________________________________________________
+Double_t AliForwardMCFlowTaskQC::AddetaFlow(Double_t eta = 0) const 
+{
+  //
+  // Add eta dependent flow factor 
+  // Parameters:
+  //  eta: choose v_n(eta) shape
+  //
+  Double_t weight = 0;
+
+  TF1 gaus = TF1("gaus", "gaus", -6, 6);
+
+  switch(fAddType)
+  {
+     case 1: gaus.SetParameters(0.1, 0., 9);
+             break;
+     case 2: gaus.SetParameters(0.1, 0., 3);
+             break;
+     case 3: gaus.SetParameters(0.1, 0., 15);
+             break;
+  }
+
+  weight = gaus.Eval(eta);
+
+  return weight;
+}
+//_____________________________________________________________________
+Double_t AliForwardMCFlowTaskQC::GetCentFromB() 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*)fAOD->FindListObject(AliAODMCHeader::StdBranchName());
+  if (!header) return cent;
+  b = header->GetImpactParameter();
+
+  cent = fImpactParToCent->Eval(b);
+
+  return cent;
+}
+//_____________________________________________________________________
+//
+//
+// EOF
+
diff --git a/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.h b/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.h
new file mode 100644 (file)
index 0000000..84cf949
--- /dev/null
@@ -0,0 +1,159 @@
+//
+// Calculate the flow in the forward regions using the Q cumulants method
+//
+#ifndef ALIFORWARDMCFLOWTASKQC_H
+#define ALIFORWARDMCFLOWTASKQC_H
+/**
+ * @file   AliForwardMCFlowTaskQC.h
+ * @author Alexander Hansen alexander.hansen@cern.ch
+ * @date   Tue Feb 14 2012
+ * 
+ * @brief  
+ * 
+ * 
+ * @ingroup pwglf_forward_flow
+ */
+#include "AliForwardFlowTaskQC.h"
+#include <TH2D.h>
+class TGraph;
+
+ /**
+ * @defgroup pwg2_forward_tasks_flow Flow tasks 
+ * @ingroup pwg2_forward_tasks
+ */
+/**
+ * Calculate the flow in the forward regions using the Q cumulants method
+ *
+ * @par Inputs:
+ *   - AliAODEvent
+ *
+ * Outputs:
+ *   - AnalysisResults.root
+ *
+ * @ingroup pwg2_forward_tasks_flow
+ * @ingroup pwg2_forward_flow
+ *
+ *
+ */
+class AliForwardMCFlowTaskQC : public AliForwardFlowTaskQC
+{
+public:
+  /**
+   * Constructor
+   */
+  AliForwardMCFlowTaskQC();
+  /*
+   * Constructor
+   *
+   * @param name Name of task
+   */
+  AliForwardMCFlowTaskQC(const char* name);
+  /**
+   * Destructor 
+   */
+  virtual ~AliForwardMCFlowTaskQC() {}
+  /** 
+   * @{ 
+   * @name Task interface methods 
+   */
+  /**
+   * Loop over AliAODMCParticle branch object and fill d^2N/detadphi histograms
+   * add flow if arguments are set
+   * 
+   * @return true on success
+   */
+  Bool_t LoopAODMC();
+    /*
+   * Set string to add flow to MC truth particles
+   *
+   * @param type String
+   */
+  void AddFlow(TString type = "") { fAddFlow = type; }
+  /*
+   * Set which function fAddFlow should use
+   *
+   * @param type of AddFlow 
+   */
+  void AddFlowType(Int_t number = 0) { fAddType = number; }
+  /*
+   * Set which order of flow to add
+   *
+   * @param order Flow order 
+   */
+  void AddFlowOrder(Int_t order = 2) { fAddOrder = order; }
+protected:
+  /*
+   * Copy constructor
+   *
+   * @param o Object to copy from
+   */
+  AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o);
+  /** 
+   * Assignment operator 
+   * 
+   * @return Reference to this object 
+   */
+  AliForwardMCFlowTaskQC& operator=(const AliForwardMCFlowTaskQC& o);
+  /*
+   * Initiate vertex bin objects
+   */
+  void InitVertexBins();
+   /*
+   * Initiate diagnostics histograms
+   */
+  void InitHists();
+  /*
+   * Analyze event
+   */
+  Bool_t Analyze();
+  /*
+   * Finalize analysis
+   */
+  void Finalize();
+  /**
+   * Add pt dependent flow factor
+   *
+   * @param Pt   @f$ p_T@f$
+   * @param type Type of flow 
+   */
+  Double_t AddptFlow(Double_t pt) const;
+  /**
+   * Add pid dependent flow factor
+   *
+   * @param ID   Particle ID 
+   * @param type Type of flow
+   */
+  Double_t AddpidFlow(Int_t id) const;
+  /**
+   * Add eta dependent flow factor
+   * 
+   * @param Eta  @f$\eta@f$ 
+   * @param type Type of flow 
+   */
+  Double_t AddetaFlow(Double_t eta) const;
+  /**
+   * Get centrality form MC impact parameter
+   */
+  Double_t GetCentFromB() const;
+  
+  TList         fBinsFMDTR;         //  List with FMDTR VertexBin objects
+  TList         fBinsSPDTR;         //  List with SPDTR VertexBin objects
+  TList         fBinsMC;            //  List with MC VertexBin objects
+  TH2D          fdNdedpMC;          //  d^2N/detadphi MC truth histogram
+  TGraph*       fAliceCent4th;      //  Parametrization of ALICE QC4 vs. cent. data
+  TGraph*       fAlicePt2nd4050;    //  Parametrization of ALICE QC2 vs. pT data
+  TGraph*       fAlicePt4th3040;    //  Parametrization of ALICE QC4 vs. pT data
+  TGraph*       fAlicePt4th4050;    //  Parametrization of ALICE QC4 vs. pT data
+  TGraph*       fImpactParToCent;   //  Parametrization of b to centrality datapoints
+  TString       fAddFlow;           //  Add flow string
+  Int_t         fAddType;           //  Add flow type #
+  Int_t         fAddOrder;          //  Add flow order
+
+  ClassDef(AliForwardMCFlowTaskQC, 1); // FMD MC analysis task 
+};
+#endif
+// Local Variables:
+//   mode: C++ 
+// End:
index 1fcad32..7b1e6bc 100644 (file)
@@ -35,6 +35,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
     fESDFMD(),
     fHistos(),
     fAODFMD(),
+    fAODEP(),
     fMCESDFMD(),
     fMCHistos(),
     fMCAODFMD(),
@@ -46,6 +47,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
     fDensityCalculator(),
     fCorrections(),
     fHistCollector(),
+    fEventPlaneFinder(),
     fList(0)
 {
   // 
@@ -60,6 +62,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
     fESDFMD(),
     fHistos(),
     fAODFMD(kFALSE),
+    fAODEP(kFALSE),
     fMCESDFMD(),
     fMCHistos(),
     fMCAODFMD(kTRUE),
@@ -71,6 +74,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
     fDensityCalculator("density"),
     fCorrections("corrections"),
     fHistCollector("collector"),
+    fEventPlaneFinder("eventplane"),
     fList(0)
 {
   // 
@@ -89,6 +93,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMul
     fESDFMD(o.fESDFMD),
     fHistos(o.fHistos),
     fAODFMD(o.fAODFMD),
+    fAODEP(o.fAODEP),
     fMCESDFMD(o.fMCESDFMD),
     fMCHistos(o.fMCHistos),
     fMCAODFMD(o.fMCAODFMD),
@@ -100,6 +105,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMul
     fDensityCalculator(o.fDensityCalculator),
     fCorrections(o.fCorrections),
     fHistCollector(o.fHistCollector),
+    fEventPlaneFinder(o.fEventPlaneFinder),
     fList(o.fList) 
 {
   // 
@@ -133,8 +139,10 @@ AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
   fDensityCalculator = o.fDensityCalculator;
   fCorrections       = o.fCorrections;
   fHistCollector     = o.fHistCollector;
+  fEventPlaneFinder  = o.fEventPlaneFinder;
   fHistos            = o.fHistos;
   fAODFMD            = o.fAODFMD;
+  fAODEP             = o.fAODEP;
   fMCHistos          = o.fMCHistos;
   fMCAODFMD          = o.fMCAODFMD;
   fRingSums          = o.fRingSums;
@@ -160,6 +168,7 @@ AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
   fDensityCalculator.SetDebug(dbg);
   fCorrections.SetDebug(dbg);
   fHistCollector.SetDebug(dbg);
+  fEventPlaneFinder.SetDebug(dbg);
 }
 //____________________________________________________________________
 void
@@ -184,6 +193,7 @@ AliForwardMCMultiplicityTask::InitializeSubs()
 
   fHistos.Init(*pe);
   fAODFMD.Init(*pe);
+  fAODEP.Init(*pe);
   fMCHistos.Init(*pe);
   fMCAODFMD.Init(*pe);
   fRingSums.Init(*pe);
@@ -234,6 +244,7 @@ AliForwardMCMultiplicityTask::InitializeSubs()
   fDensityCalculator.Init(*pe);
   fCorrections.Init(*pe);
   fHistCollector.Init(*pv,*pe);
+  fEventPlaneFinder.Init(*pe);
 
   this->Print();
 }
@@ -261,6 +272,9 @@ AliForwardMCMultiplicityTask::UserCreateOutputObjects()
   TObject* mcobj = &fMCAODFMD;
   ah->AddBranch("AliAODForwardMult", &mcobj);
 
+  TObject* epobj = &fAODFMD;
+  ah->AddBranch("AliAODForwardEP", &epobj);
+  
   fPrimary = new TH2D("primary", "MC Primaries", 
                      200, -4, 6, 20, 0, 2*TMath::Pi());
   fPrimary->SetXTitle("#eta");
@@ -276,6 +290,7 @@ AliForwardMCMultiplicityTask::UserCreateOutputObjects()
   fDensityCalculator.DefineOutput(fList);
   fCorrections.DefineOutput(fList);
   fHistCollector.DefineOutput(fList);
+  fEventPlaneFinder.DefineOutput(fList);
 
   PostData(1, fList);
 }
@@ -302,6 +317,7 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   fHistos.Clear();
   fESDFMD.Clear();
   fAODFMD.Clear();
+  fAODEP.Clear();
   fMCHistos.Clear();
   fMCESDFMD.Clear();
   fMCAODFMD.Clear();
@@ -396,6 +412,11 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   }
   fDensityCalculator.CompareResults(fHistos, fMCHistos);
   
+  if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
+    if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
+      AliWarning("Eventplane finder failed!");
+  }
+
   // Do the secondary and other corrections. 
   if (!fCorrections.Correct(fHistos, ivz)) { 
     AliWarning("Corrections failed");
index 65f21aa..ca2d20a 100644 (file)
@@ -21,7 +21,9 @@
 #include "AliFMDMCCorrector.h"
 #include "AliFMDHistCollector.h"
 #include "AliAODForwardMult.h"
+#include "AliAODForwardEP.h"
 #include "AliFMDEnergyFitter.h"
+#include "AliFMDEventPlaneFinder.h"
 #include <AliESDFMD.h>
 class AliESDEvent;
 class TH2D;
@@ -169,6 +171,18 @@ public:
    * @return Reference to AliFMDHistCollector object 
    */
   const AliFMDHistCollector& GetHistCollector() const { return fHistCollector; }
+   /**
+   * Get reference to the EventPlaneFinder algorithm 
+   * 
+   * @return Reference to AliFMDEventPlaneFinder object 
+   */
+  AliFMDEventPlaneFinder& GetEventPlaneFinder() { return fEventPlaneFinder; }
+  /**
+   * Get reference to the EventPlaneFinder algorithm 
+   * 
+   * @return Reference to AliFMDEventPlaneFinder object 
+   */
+  const AliFMDEventPlaneFinder& GetEventPlaneFinder() const { return fEventPlaneFinder; }
   /** 
    * @} 
    */
@@ -189,6 +203,7 @@ protected:
   AliESDFMD              fESDFMD;       // Sharing corrected ESD object
   AliForwardUtil::Histos fHistos;       // Cache histograms 
   AliAODForwardMult      fAODFMD;       // Output object
+  AliAODForwardEP        fAODEP;       // Output object
   AliESDFMD              fMCESDFMD;     // MC 'Sharing corrected' ESD object
   AliForwardUtil::Histos fMCHistos;     // MC Cache histograms 
   AliAODForwardMult      fMCAODFMD;     // MC Output object
@@ -201,6 +216,7 @@ protected:
   AliFMDMCDensityCalculator fDensityCalculator; // Algorithm
   AliFMDMCCorrector         fCorrections;       // Algorithm
   AliFMDHistCollector       fHistCollector;     // Algorithm
+  AliFMDEventPlaneFinder    fEventPlaneFinder;  // Algorithm
 
   TList* fList; // Output list 
 
index 3421855..3f2cdf7 100644 (file)
@@ -24,6 +24,7 @@
 #include "AliFMDDensityCalculator.h"
 #include "AliFMDCorrector.h"
 #include "AliFMDHistCollector.h"
+#include "AliFMDEventPlaneFinder.h"
 #include "AliESDEvent.h"
 #include <TROOT.h>
 #include <TAxis.h>
@@ -191,6 +192,7 @@ AliForwardMultiplicityBase::GetESDEvent()
 
     fFirstEvent = false;
 
+    GetEventPlaneFinder().SetRunNumber(esd->GetRunNumber());
     InitializeSubs();
   }
   return esd;
@@ -315,6 +317,7 @@ AliForwardMultiplicityBase::Print(Option_t* option) const
   GetDensityCalculator().Print(option);
   GetCorrections()      .Print(option);
   GetHistCollector()    .Print(option);
+  GetEventPlaneFinder() .Print(option);
   gROOT->DecreaseDirLevel();
 }
 
index e78a615..4a25f5e 100644 (file)
@@ -22,6 +22,7 @@ class AliFMDDensityCalculator;
 class AliFMDCorrector;
 class AliFMDHistCollector;
 class AliForwardCorrectionManager;
+class AliFMDEventPlaneFinder;
 class AliESDEvent;
 class TH2D;
 class TList;
@@ -234,6 +235,19 @@ public:
    * @return Reference to AliFMDHistCollector object 
    */
   virtual const AliFMDHistCollector& GetHistCollector() const = 0;
+   /**
+   * Get reference to the EventPlaneFinder algorithm 
+   * 
+   * @return Reference to AliFMDEventPlaneFinder object 
+   */
+  virtual AliFMDEventPlaneFinder& GetEventPlaneFinder() = 0;
+  /**
+   * Get reference to the EventPlaneFinder algorithm 
+   * 
+   * @return Reference to AliFMDEventPlaneFinder object 
+   */
+  virtual const AliFMDEventPlaneFinder& GetEventPlaneFinder() const = 0;
+
   /** 
    * @} 
    */
index b7508e6..20b3052 100644 (file)
@@ -34,12 +34,14 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask()
     fESDFMD(),
     fHistos(),
     fAODFMD(),
+    fAODEP(),
     fRingSums(),
     fEventInspector(),
     fSharingFilter(),
     fDensityCalculator(),
     fCorrections(),
     fHistCollector(),
+    fEventPlaneFinder(),
     fList(0)
 {
   // 
@@ -54,12 +56,14 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
     fESDFMD(),
     fHistos(),
     fAODFMD(false),
+    fAODEP(false),
     fRingSums(),
     fEventInspector("event"),
     fSharingFilter("sharing"), 
     fDensityCalculator("density"),
     fCorrections("corrections"),
     fHistCollector("collector"),
+    fEventPlaneFinder("eventplane"),
     fList(0)
 {
   // 
@@ -78,12 +82,14 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplic
     fESDFMD(o.fESDFMD),
     fHistos(o.fHistos),
     fAODFMD(o.fAODFMD),
+    fAODEP(o.fAODEP),
     fRingSums(o.fRingSums),
     fEventInspector(o.fEventInspector),
     fSharingFilter(o.fSharingFilter),
     fDensityCalculator(o.fDensityCalculator),
     fCorrections(o.fCorrections),
     fHistCollector(o.fHistCollector),
+    fEventPlaneFinder(o.fEventPlaneFinder),
     fList(o.fList) 
 {
   // 
@@ -117,8 +123,10 @@ AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
   fDensityCalculator = o.fDensityCalculator;
   fCorrections       = o.fCorrections;
   fHistCollector     = o.fHistCollector;
+  fEventPlaneFinder  = o.fEventPlaneFinder;
   fHistos            = o.fHistos;
   fAODFMD            = o.fAODFMD;
+  fAODEP             = o.fAODEP;
   fRingSums          = o.fRingSums;
   fList              = o.fList;
 
@@ -140,6 +148,7 @@ AliForwardMultiplicityTask::SetDebug(Int_t dbg)
   fDensityCalculator.SetDebug(dbg);
   fCorrections.SetDebug(dbg);
   fHistCollector.SetDebug(dbg);
+  fEventPlaneFinder.SetDebug(dbg);
 }
 
 //____________________________________________________________________
@@ -157,6 +166,7 @@ AliForwardMultiplicityTask::InitializeSubs()
 
   fHistos.Init(*pe);
   fAODFMD.Init(*pe);
+  fAODEP.Init(*pe);
   fRingSums.Init(*pe);
 
   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
@@ -185,6 +195,7 @@ AliForwardMultiplicityTask::InitializeSubs()
   fDensityCalculator.Init(*pe);
   fCorrections.Init(*pe);
   fHistCollector.Init(*pv,*pe);
+  fEventPlaneFinder.Init(*pe);
 
   this->Print();
 }
@@ -208,12 +219,15 @@ AliForwardMultiplicityTask::UserCreateOutputObjects()
     
   TObject* obj = &fAODFMD;
   ah->AddBranch("AliAODForwardMult", &obj);
+  TObject* epobj = &fAODEP;
+  ah->AddBranch("AliAODForwardEP", &epobj);
 
   fEventInspector.DefineOutput(fList);
   fSharingFilter.DefineOutput(fList);
   fDensityCalculator.DefineOutput(fList);
   fCorrections.DefineOutput(fList);
   fHistCollector.DefineOutput(fList);
+  fEventPlaneFinder.DefineOutput(fList);
 
   PostData(1, fList);
 }
@@ -237,6 +251,7 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
   fHistos.Clear();
   fESDFMD.Clear();
   fAODFMD.Clear();
+  fAODEP.Clear();
   
   Bool_t   lowFlux   = kFALSE;
   UInt_t   triggers  = 0;
@@ -257,7 +272,7 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
   fAODFMD.SetCentrality(cent);
   fAODFMD.SetNClusters(nClusters);
   MarkEventForStore();
-  
   if (found & AliFMDEventInspector::kNoSPD)      return;
   if (found & AliFMDEventInspector::kNoFMD)      return;
   if (found & AliFMDEventInspector::kNoVertex)   return;
@@ -278,13 +293,18 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
     AliWarning("Sharing filter failed!");
     return;
   }
-
+  
   // Calculate the inclusive charged particle density 
   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent)) { 
     // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) { 
     AliWarning("Density calculator failed!");
     return;
   }
+
+  if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
+    if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()), &fHistos))
+      AliWarning("Eventplane finder failed!");
+  }
   
   // Do the secondary and other corrections. 
   if (!fCorrections.Correct(fHistos, ivz)) { 
index 4792d2f..bcef960 100644 (file)
@@ -21,7 +21,9 @@
 #include "AliFMDCorrector.h"
 #include "AliFMDHistCollector.h"
 #include "AliAODForwardMult.h"
+#include "AliAODForwardEP.h"
 #include "AliFMDEnergyFitter.h"
+#include "AliFMDEventPlaneFinder.h"
 #include <AliESDFMD.h>
 class AliESDEvent;
 class TH2D;
@@ -159,9 +161,22 @@ public:
    * @return Reference to AliFMDHistCollector object 
    */
   const AliFMDHistCollector& GetHistCollector() const { return fHistCollector; }
+ /**
+   * Get reference to the EventPlaneFinder algorithm 
+   * 
+   * @return Reference to AliFMDEventPlaneFinder object 
+   */
+  AliFMDEventPlaneFinder& GetEventPlaneFinder() { return fEventPlaneFinder; }
+  /**
+   * Get reference to the EventPlaneFinder algorithm 
+   * 
+   * @return Reference to AliFMDEventPlaneFinder object 
+   */
+  const AliFMDEventPlaneFinder& GetEventPlaneFinder() const { return fEventPlaneFinder; }
   /** 
    * @} 
    */
+
   /** 
    * Set debug level 
    * 
@@ -179,6 +194,7 @@ protected:
   AliESDFMD              fESDFMD;       // Sharing corrected ESD object
   AliForwardUtil::Histos fHistos;       // Cache histograms 
   AliAODForwardMult      fAODFMD;       // Output object
+  AliAODForwardEP        fAODEP;        // Output object
   AliForwardUtil::Histos fRingSums;     // Cache histograms 
 
   AliFMDEventInspector    fEventInspector;    // Algorithm
@@ -186,6 +202,7 @@ protected:
   AliFMDDensityCalculator fDensityCalculator; // Algorithm
   AliFMDCorrector         fCorrections;       // Algorithm
   AliFMDHistCollector     fHistCollector;     // Algorithm
+  AliFMDEventPlaneFinder  fEventPlaneFinder;  // Algorithm
 
   TList* fList; // Output list 
 
index ad072af..e398497 100644 (file)
@@ -10,6 +10,9 @@
 #include <TList.h>
 #include <TROOT.h>
 #include <iostream>
+#include "AliGenHijingEventHeader.h"
+#include <TF1.h>
+#include <TGraph.h>
 
 //____________________________________________________________________
 AliSPDMCTrackDensity::AliSPDMCTrackDensity()
@@ -117,7 +120,8 @@ AliSPDMCTrackDensity::StoreParticle(AliMCParticle* particle,
                                    const AliMCParticle* mother, 
                                    Int_t          refNo,
                                    Double_t       vz, 
-                                   TH2D&     output) const
+                                   TH2D&     output,
+                                    Double_t  w) const
 {
   // Store a particle. 
   if (refNo < 0) return;
@@ -136,14 +140,14 @@ AliSPDMCTrackDensity::StoreParticle(AliMCParticle* particle,
   Double_t et = -TMath::Log(TMath::Tan(th/2));
   Double_t ph = TMath::ATan2(y,x);
   if (ph < 0) ph += 2*TMath::Pi();
-  output.Fill(et,ph);
+  output.Fill(et,ph,w);
 
   const AliMCParticle* mp = (mother ? mother : particle);
   Double_t dEta = mp->Eta() - et;
   Double_t dPhi = (mp->Phi() - ph) * 180 / TMath::Pi();
   if (dPhi >  180) dPhi -= 360;
   if (dPhi < -180) dPhi += 360;
-  fBinFlow->Fill(dEta, dPhi);
+  fBinFlow->Fill(dEta, dPhi,w);
 }
 
 
@@ -166,6 +170,7 @@ AliSPDMCTrackDensity::GetMother(Int_t     iTr,
   return 0;
 }  
 
+// #define USE_FLOW_WEIGHTS
 //____________________________________________________________________
 Bool_t
 AliSPDMCTrackDensity::Calculate(const AliMCEvent& event, 
@@ -188,6 +193,13 @@ AliSPDMCTrackDensity::Calculate(const AliMCEvent& event,
   //    True on succes, false otherwise 
   //
 
+#ifdef USE_FLOW_WEIGHTS
+  AliGenHijingEventHeader* hd = dynamic_cast<AliGenHijingEventHeader*>
+                (event.GenEventHeader());
+  Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.);
+  Double_t b = (hd ? hd->ImpactParameter() : -1 );
+#endif
+
   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
   Int_t nTracks   = stack->GetNtrack();
   Int_t nPrim     = stack->GetNtrack();
@@ -243,13 +255,127 @@ AliSPDMCTrackDensity::Calculate(const AliMCEvent& event,
       // Only fill first reference 
       if (nRef == 1) { 
        const AliMCParticle* mother = GetMother(iTr, event);
-       StoreParticle(particle, mother, iTrRef, vz, output);
+#ifdef USE_FLOW_WEIGHTS
+        Double_t phi = (mother ? mother->Phi() : particle->Phi());
+        Double_t eta = (mother ? mother->Eta() : particle->Eta());
+        Double_t pt  = (mother ? mother->Pt() : particle->Pt());
+        Int_t    id  = (mother ? mother->PdgCode() : 2212);
+       Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
+#else
+        Double_t weight = 1.; 
+#endif
+       StoreParticle(particle, mother, iTrRef, vz, output, weight);
       }
     }
     fNRefs->Fill(nRef);
   }
   return kTRUE;
 }
+
+//____________________________________________________________________
+Double_t
+AliSPDMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, Double_t b, 
+                                     Double_t phi, Double_t rp, Int_t id) const
+{
+  static TF1 gaus = TF1("gaus", "gaus", -6, 6);
+  gaus.SetParameters(0.1, 0., 9);
+  //  gaus.SetParameters(0.1, 0., 3);
+  //  gaus.SetParameters(0.1, 0., 15);
+  
+  const Double_t xCumulant2nd4050ALICE[] = {0.00, 0.25, 0.350,
+                                           0.45, 0.55, 0.650, 
+                                           0.75, 0.85, 0.950,
+                                           1.10, 1.30, 1.500,
+                                           1.70, 1.90, 2.250,
+                                           2.75, 3.25, 3.750,
+                                           4.50};
+  const Double_t yCumulant2nd4050ALICE[] = {0.00000, 0.043400,
+                                           0.059911,0.073516,
+                                           0.089756,0.105486,
+                                           0.117391,0.128199,
+                                           0.138013,0.158271,
+                                           0.177726,0.196383,
+                                           0.208277,0.216648,
+                                           0.242954,0.249961,
+                                           0.240131,0.269006,
+                                           0.207796};
+  const Int_t nPointsCumulant2nd4050ALICE = 
+    sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);                                      
+  static TGraph alicePointsPt2(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
+#if 0
+  const Double_t xCumulant4th3040ALICE[] = {0.00,0.250,0.35,
+                                           0.45,0.550,0.65,
+                                           0.75,0.850,0.95,
+                                           1.10,1.300,1.50,
+                                           1.70,1.900,2.25,
+                                           2.75,3.250,3.75,
+                                           4.50,5.500,7.00,
+                                           9.000000};
+  const Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,
+                                           0.048566,0.061083,
+                                           0.070910,0.078831,
+                                           0.091396,0.102026,
+                                           0.109691,0.124449,
+                                           0.139819,0.155561,
+                                           0.165701,0.173678,
+                                           0.191149,0.202015,
+                                           0.204540,0.212560,
+                                           0.195885,0.000000,
+                                           0.000000,0.000000};
+#endif
+  const Double_t xCumulant4th4050ALICE[] = {0.00,0.25,0.350,
+                                           0.45,0.55,0.650,
+                                           0.75,0.85,0.950,
+                                           1.10,1.30,1.500,
+                                           1.70,1.90,2.250,
+                                           2.75,3.25,3.750,
+                                           4.50};
+  const Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,
+                                           0.049824,0.066662,
+                                           0.075856,0.081583,
+                                           0.099778,0.104674,
+                                           0.118545,0.131874,
+                                           0.152959,0.155348,
+                                           0.169751,0.179052,
+                                           0.178532,0.198851,
+                                           0.185737,0.239901,
+                                           0.186098};
+  const Int_t nPointsCumulant4th4050ALICE = 
+    sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);   
+  static TGraph alicePointsPt4(nPointsCumulant4th4050ALICE, 
+                              xCumulant4th4050ALICE, 
+                              yCumulant4th4050ALICE);
+
+  const Double_t xCumulant4thTPCrefMultTPConlyAll[] = {1.75,
+                                                      4.225,
+                                                      5.965,
+                                                      7.765,
+                                                      9.215,
+                                                      10.46,
+                                                      11.565,
+                                                      12.575};
+  const Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,
+                                                      0.055818,0.073137,
+                                                      0.083898,0.086690,
+                                                      0.082040,0.077777};
+  const Int_t nPointsCumulant4thTPCrefMultTPConlyAll = 
+    sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
+  TGraph aliceCent(nPointsCumulant4thTPCrefMultTPConlyAll,
+                  xCumulant4thTPCrefMultTPConlyAll,
+                  yCumulant4thTPCrefMultTPConlyAll);
+
+
+  Double_t weight = (20. * gaus.Eval(eta) * (alicePointsPt2.Eval(pt) * 0.5 + 
+                                            alicePointsPt4.Eval(pt) * 0.5) 
+                    * (aliceCent.Eval(b) / aliceCent.Eval(10.46)) 
+                    * 2. * TMath::Cos(2. * (phi - rp)));
+  if      (TMath::Abs(id) == 211)  weight *= 1.3; //pion flow
+  else if (TMath::Abs(id) == 2212) weight *= 1.0;  //proton flow
+  else                             weight *= 0.7;
+  
+  return weight;
+}
+
 //____________________________________________________________________
 void
 AliSPDMCTrackDensity::Print(Option_t* /*option*/) const 
index a0e5f64..36af384 100644 (file)
@@ -117,7 +117,8 @@ protected:
                     const AliMCParticle* mother,
                     Int_t                refNo,
                     Double_t             vz, 
-                    TH2D&                output) const;
+                    TH2D&                output,
+                     Double_t             w) const;
   /** 
    * Get ultimate mother of a track 
    * 
@@ -137,6 +138,21 @@ protected:
    */
   Double_t GetTrackRefTheta(const AliTrackReference* ref,
                            Double_t vz) const;
+  /** 
+   * Calculate flow weight 
+   *
+   * @param eta  Pseudo rapidity 
+   * @param pt   Transverse momemtum 
+   * @param b    Impact parameter
+   * @param phi  Azimuthal angle 
+   * @param rp   Reaction plance angle 
+   * @param id   Particle PDG code
+   *
+   * @return Flow weight for the particle
+   */
+  Double_t CalculateWeight(Double_t eta, Double_t pt, Double_t b, 
+                          Double_t phi, Double_t rp, Int_t id) const;
+
   Bool_t   fUseOnlyPrimary;       // Only use primaries 
   Double_t fMinR;             // Min radius 
   Double_t fMaxR;             // Max radius 
index aed6bf9..6d2edd3 100644 (file)
@@ -151,6 +151,9 @@ ForwardAODConfig(AliForwardMultiplicityBase* task)
   // Set the debug level of a single algorithm 
   task->GetSharingFilter().SetDebug(0);
 
+  // --- Eventplane Finder -------------------------------------------
+  task->GetEventPlaneFinder().SetUsePhiWeights(false);
+
   // --- Set limits on fits the energy -------------------------------
   // Maximum relative error on parameters 
   AliFMDCorrELossFit::ELossFit::fgMaxRelError = .12;
index ba7d7bc..424847a 100644 (file)
@@ -5,7 +5,7 @@
  * 
  * @brief  
  * 
- * @ingroup pwglf_forward_scripts_makers
+ * @ingroup pwg2_forward_scripts_makers
  * 
  */
 /**
  * @par Outputs: 
  * - 
  *
- * @ingroup pwglf_forward_flow
+ * @ingroup pwg2_forward_flow
  */
 void MakeFlow(TString data      = "", 
-             Int_t   nevents   = 0, 
+             Int_t   nEvents   = 0, 
              TString type      = "", 
-             Int_t   etabins   = 48,
               Bool_t  mc        = kFALSE,
+             const char* name  = 0,
+             Int_t   proof     = 0,
              TString addFlow   = "",
               Int_t   addFType  = 0,
               Int_t   addFOrder = 0,
-             Bool_t  proof     = kFALSE)
+              Bool_t  gdb       = kFALSE)
 {
-  Bool_t proof = kFALSE;
-
   // --- Load libs ---------------------------------------------------
   gROOT->Macro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");
 
-  // --- Check for proof mode, and possibly upload pars --------------
-  if (proof> 0) { 
-    gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadPars.C");
-    if (!LoadPars(proof)) { 
-      AliError("MakeFlow", "Failed to load PARs");
-      return;
-    }
+  // --- Possibly use plug-in for this -------------------------------
+  if ((name && name[0] != '\0') && gSystem->Load("libRAliEn") >= 0) {
+    gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/trains/TrainSetup.C+");
+    gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C+");
+
+    MakeFlowTrain t(name, type.Data(), mc, addFlow.Data(), addFType, addFOrder, false);
+    t.SetDataDir(data.Data());
+    t.SetDataSet("");
+    t.SetProofServer(Form("workers=%d", proof));
+    t.SetUseGDB(gdb);
+    t.Run(proof > 0 ? "proof" : "local", "full", nEvents, proof > 0);
+    return;
   }
 
   // --- Set the macro path ------------------------------------------
@@ -57,19 +62,10 @@ void MakeFlow(TString data      = "",
     AliError("You didn't add a data file");
     return;
   }
-  TChain* chain = new TChain("aodTree");
-
-  if (data.Contains(".txt")) MakeChain(data, chain);
-
-  if (data.Contains(".root")) {
-    TFile* test = TFile::Open(data.Data());
-    if (!test) {
-      AliError(Form("AOD file %s not found", data.Data()));
-      return;
-    }
-    test->Close(); // Remember to close!
-    chain->Add(data.Data());
-  }
+  gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/MakeChain.C");
+  TChain* chain = MakeChain("AOD", data.Data(), true);
+  // If 0 or less events is select, choose all 
+  if (nEvents <= 0) nEvents = chain->GetEntries();
 
   // --- Initiate the event handlers --------------------------------
   AliAnalysisManager *mgr  = new AliAnalysisManager("Forward Flow", 
@@ -81,7 +77,7 @@ void MakeFlow(TString data      = "",
 
   // --- Add the tasks ---------------------------------------------
   gROOT->LoadMacro("AddTaskForwardFlow.C");
-  AddTaskForwardFlow(type, etabins, mc, addFlow, addFType, addFOrder);
+  AddTaskForwardFlow(type, mc, addFlow, addFType, addFOrder);
 
   // --- Run the analysis --------------------------------------------
   TStopwatch t;
@@ -91,51 +87,21 @@ void MakeFlow(TString data      = "",
   }
   mgr->PrintStatus();
   Printf("****************************************");
-  Printf("Doing flow analysis on %d Events", nevents == 0 ? chain->GetEntries() : nevents);
+  Printf("Doing flow analysis on %d Events", nEvents);
   Printf("****************************************");
   // 
-  if (proof) mgr->SetDebugLevel(3);
-  if (mgr->GetDebugLevel() < 1 && !proof) 
-    mgr->SetUseProgressBar(kTRUE,chain->GetEntries() < 10000 ? 100 : 1000);
+  mgr->SetDebugLevel(0);
+  if (mgr->GetDebugLevel() < 1) 
+    mgr->SetUseProgressBar(kTRUE, nEvents < 10000 ? 100 : 1000);
 
 //  mgr->SetSkipTerminate(true);
 
   t.Start();
-  if (nevents == 0) mgr->StartAnalysis("local", chain);
-  if (nevents != 0) mgr->StartAnalysis("local", chain, nevents);
+  mgr->StartAnalysis("local", chain, nEvents);
   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;
-  TFile* file;
-  while(in.good()) 
-  {
-    in >> line;
-      
-    if (line.Length() == 0)
-      continue;      
-    if (!(file = TFile::Open(line))) 
-      gROOT->ProcessLine(Form(".!rm %s", line.Data()));
-    else {
-      chain->Add(line);
-      file->Close();
-    }
-  }
-
-  in.close();
-
-  return;
-}
 //
 // EOF
 //
index 3ac1c70..9f0c7d0 100644 (file)
@@ -1,3 +1,5 @@
+#include "TrainSetup.C"
+
 /**
  * @file   MakeAODTrain.C
  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
@@ -8,7 +10,6 @@
  * @ingroup pwglf_forward_trains
  */
 //====================================================================
-#include "TrainSetup.C"
 
 /**
  * Analysis train to make Forward and Central multiplicity
@@ -117,10 +118,15 @@ protected:
     // --- Set load path ---------------------------------------------
     gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2",
                             gROOT->GetMacroPath()));
+    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/ANALYSIS/macros",
+                            gROOT->GetMacroPath()));
 
     // --- Check if this is MC ---------------------------------------
     Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
     
+    // --- Add TPC eventplane task
+    gROOT->Macro("AddTaskEventplane.C");
+
     // --- Task to copy header information ---------------------------
     gROOT->Macro("AddTaskCopyHeader.C");
 
@@ -133,6 +139,10 @@ protected:
     gROOT->Macro(Form("AddTaskCentralMult.C(%d,%d,%d,%d)", 
                      mc, fSys, fSNN, fField));
     AddExtraFile(gSystem->Which(gROOT->GetMacroPath(), "CentralAODConfig.C"));
+
+    // --- Add MC particle task --------------------------------------
+    if (mc) gROOT->Macro("AddTaskMCParticleFilter.C");
+
   }
   //__________________________________________________________________
   /** 
@@ -161,6 +171,41 @@ protected:
     // --- Ignore trigger class when selecting events.  This means ---
     // --- that we get offline+(A,C,E) events too --------------------
     // ps->SetSkipTriggerClassSelection(true);
+/*
+    if (mc) {
+      AliOADBPhysicsSelection * oadbDefaultPbPb = new AliOADBPhysicsSelection("oadbDefaultPbPb");
+      oadbDefaultPbPb->AddCollisionTriggerClass   ( AliVEvent::kHighMult,"+C0SMH-B-NOPF-ALLNOTRD","B",0);
+      oadbDefaultPbPb->AddBGTriggerClass          ( AliVEvent::kHighMult,"+C0SMH-A-NOPF-ALLNOTRD","A",0);
+      oadbDefaultPbPb->AddBGTriggerClass          ( AliVEvent::kHighMult,"+C0SMH-C-NOPF-ALLNOTRD","C",0);
+      oadbDefaultPbPb->AddBGTriggerClass          ( AliVEvent::kHighMult,"+C0SMH-E-NOPF-ALLNOTRD","E",0);  
+      oadbDefaultPbPb->SetHardwareTrigger         ( 0,"SPDGFO >= 100");
+      oadbDefaultPbPb->SetOfflineTrigger          ( 0,"SPDGFO >= 100 && !V0ABG && !V0CBG");
+  
+      oadbDefaultPbPb->AddCollisionTriggerClass   ( AliVEvent::kMB,"+CMBACS2-B-NOPF-ALLNOTRD","B",1);
+      oadbDefaultPbPb->AddBGTriggerClass          ( AliVEvent::kMB,"+CMBACS2-A-NOPF-ALLNOTRD","A",1);
+      oadbDefaultPbPb->AddBGTriggerClass          ( AliVEvent::kMB,"+CMBACS2-C-NOPF-ALLNOTRD","C",1);
+      oadbDefaultPbPb->AddBGTriggerClass          ( AliVEvent::kMB,"+CMBACS2-E-NOPF-ALLNOTRD","E",1);
+      oadbDefaultPbPb->SetHardwareTrigger         ( 1,"(V0A && V0C && SPDGFO > 1)");
+      oadbDefaultPbPb->SetOfflineTrigger          ( 1,"(V0A && V0C && SPDGFO > 1) && !V0ABG && !V0CBG");
+    // LHC10h8
+      AliOADBPhysicsSelection * oadbLHC10h8 = new AliOADBPhysicsSelection("oadbLHC10h8");
+      oadbLHC10h8->AddCollisionTriggerClass   ( AliVEvent::kHighMult,"+C0SMH-B-NOPF-ALL","B",0);
+      oadbLHC10h8->AddBGTriggerClass          ( AliVEvent::kHighMult,"+C0SMH-A-NOPF-ALL","A",0);
+      oadbLHC10h8->AddBGTriggerClass          ( AliVEvent::kHighMult,"+C0SMH-C-NOPF-ALL","C",0);
+      oadbLHC10h8->AddBGTriggerClass          ( AliVEvent::kHighMult,"+C0SMH-E-NOPF-ALL","E",0);  
+      oadbLHC10h8->SetHardwareTrigger         ( 0,"SPDGFO >= 100");
+      oadbLHC10h8->SetOfflineTrigger          ( 0,"SPDGFO >= 100 && !V0ABG && !V0CBG");
+  
+      oadbLHC10h8->AddCollisionTriggerClass   ( AliVEvent::kMB,"+CMBACS2-B-NOPF-ALL","B",1);
+      oadbLHC10h8->AddBGTriggerClass          ( AliVEvent::kMB,"+CMBACS2-A-NOPF-ALL","A",1);
+      oadbLHC10h8->AddBGTriggerClass          ( AliVEvent::kMB,"+CMBACS2-C-NOPF-ALL","C",1);
+      oadbLHC10h8->AddBGTriggerClass          ( AliVEvent::kMB,"+CMBACS2-E-NOPF-ALL","E",1);
+      oadbLHC10h8->SetHardwareTrigger         ( 1,"(V0A && V0C && SPDGFOL1 > 1)");
+      oadbLHC10h8->SetOfflineTrigger          ( 1,"(V0A && V0C && SPDGFOL1 > 1) && !V0ABG && !V0CBG");
+
+      ps->SetCustomOADBObjects(oadbDefaultPbPb, 0, 0);
+      ps->SetCustomOADBObjects(oadbLHC10h8, 0, 0);
+    }*/
   }
   //__________________________________________________________________
   /** 
diff --git a/PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C b/PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C
new file mode 100644 (file)
index 0000000..2f92142
--- /dev/null
@@ -0,0 +1,124 @@
+#include "TrainSetup.C"
+
+//====================================================================
+/**
+ * Analysis train to make @f$ flow@f$
+ * 
+ * To run, do 
+ * @code 
+ * gROOT->LoadMacro("TrainSetup.C");
+ * // Make train 
+ * MakeFlowTrain t("My Analysis");
+ * // Set variaous parameters on the train 
+ * t.SetDataDir("/home/of/data");
+ * t.AddRun(118506)
+ * // Run it 
+ * t.Run("LOCAL", "FULL", -1, false, false);
+ * @endcode 
+ *
+ * @ingroup pwg2_forward_flow
+ * @ingroup pwg2_forward_trains
+ */
+class MakeFlowTrain : public TrainSetup
+{
+public:
+  /** 
+   * Constructor.  Date and time must be specified when running this
+   * in Termiante mode on Grid
+   * 
+   * @param name     Name of train (free form)
+   * @param dateTime Append date and time to name 
+   * @param year     Year     - if not specified, current year
+   * @param month    Month    - if not specified, current month
+   * @param day      Day      - if not specified, current day
+   * @param hour     Hour     - if not specified, current hour
+   * @param min      Minutes  - if not specified, current minutes
+   */
+  MakeFlowTrain(const char* name, 
+               char*       type="", 
+               Bool_t      mc = false,
+               char*       addFlow = "",
+               Int_t       addFType = 0,
+               Int_t       addFOrder = 0,
+               Bool_t      dateTime=false,
+               UShort_t    year  = 0, 
+               UShort_t    month = 0, 
+               UShort_t    day   = 0, 
+               UShort_t    hour  = 0, 
+               UShort_t    min   = 0) 
+    : TrainSetup(name, dateTime, year, month, day, hour, min),
+      fType(type), 
+      fMC(mc),
+      fAddFlow(addFlow),
+      fAddFType(addFType),
+      fAddFOrder(addFOrder)
+  {
+  }
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode - see TrainSetup::EMode
+   * @param oper     Operation - see TrainSetup::EOperation
+   * @param nEvents  Number of events (negative means all)
+   * @param usePar   If true, use PARs 
+   */
+  void Run(const char* mode, const char* oper, 
+          Int_t nEvents=-1, Bool_t usePar=false)
+  {
+    Exec("AOD", mode, oper, nEvents, false, usePar);
+  }
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode - see TrainSetup::EMode
+   * @param oper     Operation - see TrainSetup::EOperation
+   * @param nEvents  Number of events (negative means all)
+   * @param usePar   If true, use PARs 
+   */
+  void Run(EMode mode, EOper oper, Int_t nEvents=-1, 
+          Bool_t usePar=false)
+  {
+    Exec(kAOD, mode, oper, nEvents, false, usePar);
+  }
+protected:
+  /** 
+   * Create the tasks 
+   * 
+   * @param mode Processing mode
+   * @param par  Whether to use par files 
+   */
+  void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager*)
+  {
+    // --- Output file name ------------------------------------------
+    AliAnalysisManager::SetCommonFileName("AnalysisResults.root");
+
+    // --- Load libraries/pars ---------------------------------------
+    LoadLibrary("PWGLFforward2", mode, par, true);
+    
+    // --- Set load path ---------------------------------------------
+    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2",
+                            gROOT->GetMacroPath()));
+
+    // --- Add the task ----------------------------------------------
+    gROOT->Macro(Form("AddTaskForwardFlow.C(\"%s\",%d,\"%s\",%d,%d)",
+                     fType, fMC, fAddFlow, fAddFType, fAddFOrder));
+  }
+  /** 
+   * Do not the centrality selection
+   */
+  void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {}
+  /** 
+   * Crete output handler - we don't want one here. 
+   * 
+   * @return 0
+   */
+  AliVEventHandler* CreateOutputHandler(EType) { return 0; }
+  char* fType;
+  Bool_t fMC;
+  char* fAddFlow;
+  Int_t fAddFType;
+  Int_t fAddFOrder;
+};
+//
+// EOF
+//
index 99514c0..9e30e89 100644 (file)
@@ -34,6 +34,7 @@
 #pragma link C++ class AliForwardUtil::RingHistos+;
 #pragma link C++ class AliFMDEventInspector+;
 #pragma link C++ class AliFMDMCEventInspector+;
+#pragma link C++ class AliFMDEventPlaneFinder+;
 #pragma link C++ class AliFMDSharingFilter+;
 #pragma link C++ class AliFMDSharingFilter::RingHistos+;
 #pragma link C++ class AliFMDMCSharingFilter+;
@@ -57,7 +58,8 @@
 #pragma link C++ class AliFMDCorrDoubleHit+;
 #pragma link C++ class AliFMDCorrVertexBias+;
 #pragma link C++ class AliFMDCorrMergingEfficiency+;
-#pragma link C++ class AliAODForwardMult+;
+#pragma link C++ class AliAODForwardMult+; 
+#pragma link C++ class AliAODForwardEP+;
 #pragma link C++ class AliForwardMultiplicityBase+;
 #pragma link C++ class AliForwardMultiplicityTask+;
 #pragma link C++ class AliForwardMCMultiplicityTask+;
@@ -81,8 +83,9 @@
 #pragma link C++ class AliCentralCorrSecondaryMap+;
 #pragma link C++ class AliCentralCorrAcceptance+;
 #pragma link C++ class AliCentraldNdetaTask+;
-#pragma link C++ class AliForwardFlowUtil+;
 #pragma link C++ class AliForwardFlowTaskQC+;
+#pragma link C++ class AliForwardFlowTaskQC::VertexBin+;
+#pragma link C++ class AliForwardMCFlowTaskQC+;
 #pragma link C++ class AliSPDMCTrackDensity+;
 #pragma link C++ class AliFMDMultCuts+;
 #pragma link C++ class AliPoissonCalculator+;