make cumulants work on caf (not using tprof3d for the time being)
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Aug 2008 11:17:56 +0000 (11:17 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Aug 2008 11:17:56 +0000 (11:17 +0000)
PWG2/FLOW/AliAnalysisTaskCumulants.cxx
PWG2/FLOW/AliAnalysisTaskCumulants.h
PWG2/FLOW/AliCumulantsFunctions.cxx
PWG2/FLOW/AliCumulantsFunctions.h
PWG2/FLOW/AliFlowAnalysisWithCumulants.cxx
PWG2/FLOW/AliFlowAnalysisWithCumulants.h
PWG2/FLOW/macros/runAliAnalysisTaskCumulants.C
PWG2/FLOW/macros/runProofCumulants.C

index 8ac3914..5808406 100644 (file)
 * provided "as is" without express or implied warranty.                  * 
 **************************************************************************/
 
-#include "Riostream.h" //needed as include
+/**************************************
+ * analysis task for cumulant method  * 
+ *                                    * 
+ * authors: Naomi van der Kolk        *
+ *           (kolk@nikhef.nl)         *  
+ *          Raimond Snellings         *
+ *           (snelling@nikhef.nl)     * 
+ *          Ante Bilandzic            *
+ *           (anteb@nikhef.nl)        * 
+ * ***********************************/
+#include "Riostream.h"
 #include "TChain.h"
 #include "TTree.h"
 #include "TFile.h"
 #include "TList.h"
 #include "TH1.h"
-#include "TH3D.h"
 #include "TProfile.h"
 #include "TProfile2D.h"
 #include "TProfile3D.h"
 
-#include "AliCumulantsFunctions.h"
-
 #include "AliAnalysisTask.h"
 #include "AliAnalysisDataSlot.h"
 #include "AliAnalysisDataContainer.h"
-
-//class AliAnalysisTask;
 #include "AliAnalysisManager.h"
 
 #include "AliESDEvent.h"
 #include "AliFlowCumuConstants.h"
 #include "AliFlowCommonConstants.h"
 #include "AliFlowCommonHistResults.h"
-
 #include "AliCumulantsFunctions.h"
 
-// AliAnalysisTaskCumulants:
-// analysis task for 
-// Cumulant method
-// with many authors (N.K. R.S. A.B.)
-// who do something
-
 ClassImp(AliAnalysisTaskCumulants)
 
-//________________________________________________________________________
-AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name) : 
-  AliAnalysisTask(name, ""), 
-  fESD(NULL),
-  fAOD(NULL),
-  fMyCumuAnalysis(NULL),
-  fEventMaker(NULL),
-  fAnalysisType("ESD"), 
-  fCFManager1(NULL),
-  fCFManager2(NULL),
-  fListHistos(NULL)
+//================================================================================================================
+
+AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name): 
+ AliAnalysisTask(name,""), 
+ fESD(NULL),
+ fAOD(NULL),
+ fCA(NULL),//Cumulant Analysis (CA) object
+ fEventMaker(NULL),
+ fAnalysisType("ESD"), 
+ fCFManager1(NULL),
+ fCFManager2(NULL),
+ fListHistos(NULL)
 {
-  // Constructor
-  cout<<"AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name)"<<endl;
-
-  // Define input and output slots here
-  // Input slot #0 works with a TChain
-  DefineInput(0, TChain::Class());
-
-  // Output slot #0 writes into a TList container
-  DefineOutput(0, TList::Class());  
+ //constructor
+ cout<<"AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name)"<<endl;
+ //input and output slots:
+ //input slot #0 works with a TChain
+ DefineInput(0, TChain::Class());
+ //output slot #0 writes into a TList container
+ DefineOutput(0, TList::Class());  
 }
-//________________________________________________________________________
-AliAnalysisTaskCumulants::AliAnalysisTaskCumulants() : 
-  fESD(NULL),
-  fAOD(NULL), 
-  fMyCumuAnalysis(NULL),
-  fEventMaker(NULL),
-  fAnalysisType("ESD"),
-  fCFManager1(NULL),
-  fCFManager2(NULL),
-  fListHistos(NULL)
+
+AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(): 
+ fESD(NULL),
+ fAOD(NULL), 
+ fCA(NULL),//Cumulant Analysis (CA) object
+ fEventMaker(NULL),
+ fAnalysisType("ESD"),
+ fCFManager1(NULL),
+ fCFManager2(NULL),
+ fListHistos(NULL)
 {
-  // Constructor
-  cout<<"AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name)"<<endl;
-  }
+ //dummy constructor
+ cout<<"AliAnalysisTaskCumulants::AliAnalysisTaskCumulants()"<<endl;
+}
+
+//================================================================================================================
 
-//________________________________________________________________________
 void AliAnalysisTaskCumulants::ConnectInputData(Option_t *) 
 {
-  // Connect ESD or AOD here
-  // Called once
-  cout<<"AliAnalysisTaskCumulants::ConnectInputData(Option_t *)"<<endl;
-
-  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
-  if (!tree) {
-    Printf("ERROR: Could not read chain from input slot 0");
-  } 
-  else {
-    // Disable all branches and enable only the needed ones
-    if (fAnalysisType == "MC") {
-      // we want to process only MC
+ //connect ESD or AOD (called once)
+ cout<<"AliAnalysisTaskCumulants::ConnectInputData(Option_t *)"<<endl;
+
+ TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+ if (!tree) 
+ {
+  Printf("ERROR: Could not read chain from input slot 0");
+ } 
+ else 
+ {
+ //disable all branches and enable only the needed ones
+  if (fAnalysisType == "MC") {
+     // we want to process only MC
       tree->SetBranchStatus("*", kFALSE);
 
       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
@@ -151,38 +151,50 @@ void AliAnalysisTaskCumulants::ConnectInputData(Option_t *)
   }
 }
 
-//________________________________________________________________________
+//================================================================================================================
+
 void AliAnalysisTaskCumulants::CreateOutputObjects() 
 {
-  // Called at every worker node to initialize
-  cout<<"AliAnalysisTaskCumulants::CreateOutputObjects()"<<endl;
+ //called at every worker node to initialize
+ cout<<"AliAnalysisTaskCumulants::CreateOutputObjects()"<<endl;
 
-  if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMC0"  || fAnalysisType == "ESDMC1" || fAnalysisType == "MC")) {
-    cout<<"WRONG ANALYSIS TYPE! only ESD, ESDMC0, ESDMC1, AOD and MC are allowed."<<endl;
-    exit(1);
-  }
  
-  //event maker
-  fEventMaker = new AliFlowEventSimpleMaker();
-  //analyser
-  fMyCumuAnalysis = new AliFlowAnalysisWithCumulants() ;
-   
-  fMyCumuAnalysis->CreateOutputObjects();
-
-  if (fMyCumuAnalysis->GetHistList()) {
-    // fSP->GetHistList()->Print();
-    fListHistos = fMyCumuAnalysis->GetHistList();
-    // fListHistos->Print();
-  }
-  else {Printf("ERROR: Could not retrieve histogram list"); }
+ //OpenFile(0);
+
+ if(!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1" || fAnalysisType == "MC")) 
+ {
+  cout<<"WRONG ANALYSIS TYPE! only ESD, ESDMC0, ESDMC1, AOD and MC are allowed."<<endl;
+  exit(1);
+ }
+ //event maker
+ fEventMaker = new AliFlowEventSimpleMaker();
+  
+ //analyser
+ fCA = new AliFlowAnalysisWithCumulants();
+ fCA->CreateOutputObjects();
+
+ if(fCA->GetHistList()) 
+ {
+  fListHistos = fCA->GetHistList();
+  //fListHistos->Print();
  }
+ else 
+ {
+  Printf("ERROR: Could not retrieve histogram list"); 
+ }
+ //PostData(0,fListHistos);
+}
+
+//================================================================================================================
 
-//________________________________________________________________________
 void AliAnalysisTaskCumulants::Exec(Option_t *) 
 {
-  // Main loop
-  // Called for each event
-  if (fAnalysisType == "MC") {
+ //main loop (called for each event)
+ if (fAnalysisType == "MC") {
     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
     // This handler can return the current MC event
 
@@ -204,7 +216,7 @@ void AliAnalysisTaskCumulants::Exec(Option_t *)
 
     //cumulant analysis 
     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(mcEvent,fCFManager1,fCFManager2);
-    fMyCumuAnalysis->Make(fEvent);
+    fCA->Make(fEvent);
     delete fEvent;
   }
   else if (fAnalysisType == "ESD") {
@@ -214,9 +226,9 @@ void AliAnalysisTaskCumulants::Exec(Option_t *)
     }
     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
     
-    //Cumulant analysis 
+    //cumulant analysis 
     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD);
-    fMyCumuAnalysis->Make(fEvent);
+    fCA->Make(fEvent);
     delete fEvent;
   }
   else if (fAnalysisType == "ESDMC0") {
@@ -241,14 +253,14 @@ void AliAnalysisTaskCumulants::Exec(Option_t *)
     fCFManager1->SetEventInfo(mcEvent);
     fCFManager2->SetEventInfo(mcEvent);
 
-    //Cumulant analysis 
+    //cumulant analysis 
     AliFlowEventSimple* fEvent=NULL;
     if (fAnalysisType == "ESDMC0") { 
       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 0); //0 = kine from ESD, 1 = kine from MC
     } else if (fAnalysisType == "ESDMC1") {
       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 1); //0 = kine from ESD, 1 = kine from MC
     }
-    fMyCumuAnalysis->Make(fEvent);
+    fCA->Make(fEvent);
     delete fEvent;
     //delete mcEvent;
   }
@@ -263,7 +275,7 @@ void AliAnalysisTaskCumulants::Exec(Option_t *)
     // analysis 
     //For the moment don't use CF //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD,fCFManager1,fCFManager2);
     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD);
-    fMyCumuAnalysis->Make(fEvent);
+    fCA->Make(fEvent);
     delete fEvent;
   }
 
@@ -271,41 +283,65 @@ void AliAnalysisTaskCumulants::Exec(Option_t *)
   PostData(0,fListHistos); 
 }
 
-//________________________________________________________________________
+//================================================================================================================
+
 void AliAnalysisTaskCumulants::Terminate(Option_t *) 
 {  
-  //=====================================================================================================
-  // Accessing the output file which contains the merged results from all workers;
-        fListHistos = (TList*)GetOutputData(0);
-        //fListHistos->Print();//printing the list of stored histos
-  //=====================================================================================================
-  
-     if(fListHistos)
-     { 
-      // Profiles with values of generating functions
-      TProfile2D *fIntFlowGenFun=dynamic_cast<TProfile2D*>(fListHistos->FindObject("fIntFlowGenFun")); 
-      TProfile3D *fDiffFlowGenFunRe=dynamic_cast<TProfile3D*>(fListHistos->FindObject("fDiffFlowGenFunRe"));
-      TProfile3D *fDiffFlowGenFunIm=dynamic_cast<TProfile3D*>(fListHistos->FindObject("fDiffFlowGenFunIm"));
+ //accessing the output list which contains the merged 2D and 3D profiles from all worker nodes
+ fListHistos = (TList*)GetOutputData(0);
+ //fListHistos->Print();
+ if(fListHistos)
+ {         
+  //profiles with avarage values of generating functions for int. and diff. flow
+  TProfile2D *intFlowGenFun = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fIntFlowGenFun")); 
+  //TProfile3D *diffFlowGenFunRe=dynamic_cast<TProfile3D*>(fListHistos->FindObject("fDiffFlowGenFunRe"));
+  //TProfile3D *diffFlowGenFunIm=dynamic_cast<TProfile3D*>(fListHistos->FindObject("fDiffFlowGenFunIm"));
       
-      // Histograms to store final results
-      TH1D *fIntFlowResults=dynamic_cast<TH1D*>(fListHistos->FindObject("fIntFlowResults"));
-      TH1D *fDiffFlowResults2=dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults2"));
-      TH1D *fDiffFlowResults4=dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults4"));
-      TH1D *fDiffFlowResults6=dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults6"));
-      TH1D *fDiffFlowResults8=dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults8"));
-     
-      // Avarage multiplicity 
-      TProfile *AvMult=dynamic_cast<TProfile*>(fListHistos->FindObject("fHistProAvM"));
-      Double_t BvM=AvMult->GetBinContent(1);//avarage multiplicity
-   
-      // Calling the function which from values of generating functions calculate integrated and differential flow and store the final results into histograms 
-      AliCumulantsFunctions FinalResults(fIntFlowGenFun,fDiffFlowGenFunRe,fDiffFlowGenFunIm,fIntFlowResults,fDiffFlowResults2,fDiffFlowResults4,fDiffFlowResults6,fDiffFlowResults8,BvM);    
-      FinalResults.Calculate();
-    }
-    else
-    {
-     cout<<"histogram list pointer is empty"<<endl;
-    }
+  TProfile2D *diffFlowGenFunRe0 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe0"));
+  TProfile2D *diffFlowGenFunRe1 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe1")); 
+  TProfile2D *diffFlowGenFunRe2 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe2")); 
+  TProfile2D *diffFlowGenFunRe3 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe3")); 
+  TProfile2D *diffFlowGenFunRe4 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe4")); 
+  TProfile2D *diffFlowGenFunRe5 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe5")); 
+  TProfile2D *diffFlowGenFunRe6 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe6")); 
+  TProfile2D *diffFlowGenFunRe7 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe7")); 
+  TProfile2D *diffFlowGenFunIm0 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm0")); 
+  TProfile2D *diffFlowGenFunIm1 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm1")); 
+  TProfile2D *diffFlowGenFunIm2 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm2")); 
+  TProfile2D *diffFlowGenFunIm3 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm3")); 
+  TProfile2D *diffFlowGenFunIm4 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm4")); 
+  TProfile2D *diffFlowGenFunIm5 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm5")); 
+  TProfile2D *diffFlowGenFunIm6 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm6")); 
+  TProfile2D *diffFlowGenFunIm7 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm7")); 
+    
+  //histograms to store the final results
+  TH1D *intFlowResults   = dynamic_cast<TH1D*>(fListHistos->FindObject("fIntFlowResults"));
+  TH1D *diffFlowResults2 = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults2"));
+  TH1D *diffFlowResults4 = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults4"));
+  TH1D *diffFlowResults6 = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults6"));
+  TH1D *diffFlowResults8 = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults8"));
+
+  //profile with avarage selected multiplicity for int. flow 
+  TProfile *avMult = dynamic_cast<TProfile*>(fListHistos->FindObject("fAvMultIntFlow"));
+  
+  //profile with avarage values of Q-vector components (1st bin: <Q_x>, 2nd bin: <Q_y>, 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>) 
+  TProfile *QVectorComponents = dynamic_cast<TProfile*>(fListHistos->FindObject("fQVectorComponents"));
+  
+  //q-distribution
+  TH1D *qDist = dynamic_cast<TH1D*>(fListHistos->FindObject("fQDist"));
+  
+  AliCumulantsFunctions finalResults(intFlowGenFun,NULL,NULL, intFlowResults,diffFlowResults2,diffFlowResults4,diffFlowResults6,diffFlowResults8,avMult,QVectorComponents,qDist,diffFlowGenFunRe0,diffFlowGenFunRe1,diffFlowGenFunRe2,
+diffFlowGenFunRe3,diffFlowGenFunRe4,diffFlowGenFunRe5,diffFlowGenFunRe6,diffFlowGenFunRe7,diffFlowGenFunIm0,diffFlowGenFunIm1,
+diffFlowGenFunIm2,diffFlowGenFunIm3,diffFlowGenFunIm4,diffFlowGenFunIm5,diffFlowGenFunIm6,diffFlowGenFunIm7);
+         
+  finalResults.Calculate();  
+ }
+ else
+ {
+  cout<<"histogram list pointer is empty"<<endl;
+ }
 }
 
 
index 2f51e03..3f9c651 100644 (file)
@@ -1,16 +1,24 @@
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
-* See cxx source for full Copyright notice */
-/* $Id$ */
-
+/* 
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. 
+ * See cxx source for full Copyright notice 
+ * $Id$ 
+ */
+
+/**************************************
+ * analysis task for cumulant method  * 
+ *                                    * 
+ * authors: Naomi van der Kolk        *
+ *           (kolk@nikhef.nl)         *  
+ *          Raimond Snellings         *
+ *           (snelling@nikhef.nl)     * 
+ *          Ante Bilandzic            *
+ *           (anteb@nikhef.nl)        * 
+ * ***********************************/
 
 #ifndef AliAnalysisTaskCumulants_H
 #define AliAnalysisTaskCumulants_H
 
-// AliAnalysisTaskCumulants:
-// analysis task for 
-// Cumulant method
-// with many authors (N.K. R.S. A.B.)
-// who do something
+#include "AliAnalysisTask.h"
 
 class AliESDEvent;
 class AliAODEvent;
@@ -19,45 +27,45 @@ class AliFlowAnalysisWithCumulants;
 class AliFlowEventSimpleMaker;
 class TFile;
 
-#include "AliAnalysisTask.h"
+//================================================================================================================
 
-class AliAnalysisTaskCumulants : public AliAnalysisTask {
+class AliAnalysisTaskCumulants : public AliAnalysisTask{
  public:
-  //AliAnalysisTaskCumulants(const char *name = "AliAnalysisTaskCumulants");
   AliAnalysisTaskCumulants();
   AliAnalysisTaskCumulants(const char *name);
-  virtual ~AliAnalysisTaskCumulants() {}
+  virtual ~AliAnalysisTaskCumulants(){}; 
   
   virtual void   ConnectInputData(Option_t *);
   virtual void   CreateOutputObjects();
   virtual void   Exec(Option_t *option);
   virtual void   Terminate(Option_t *);
-
+  
   void           SetAnalysisType(TString type) {this->fAnalysisType = type;}
-  TString GetAnalysisType() const    { return this->fAnalysisType; }
-
-  void SetCFManager1(AliCFManager* cfmgr) {this->fCFManager1 = cfmgr; } 
-  AliCFManager* GetCFManager1()           {return this->fCFManager1; }
-  void SetCFManager2(AliCFManager* cfmgr) {this->fCFManager2 = cfmgr; } 
-  AliCFManager* GetCFManager2()           {return this->fCFManager2; } 
+  TString GetAnalysisType() const {return this->fAnalysisType;}
+  
+  void SetCFManager1(AliCFManager* cfmgr) {this->fCFManager1 = cfmgr;} 
+  AliCFManager* GetCFManager1()           {return this->fCFManager1;}
+  void SetCFManager2(AliCFManager* cfmgr) {this->fCFManager2 = cfmgr;} 
+  AliCFManager* GetCFManager2()           {return this->fCFManager2;} 
  
  private:
-  AliAnalysisTaskCumulants(const AliAnalysisTaskCumulants& aAnalysis);
-  AliAnalysisTaskCumulants& operator=(const AliAnalysisTaskCumulants& aAnalysis);
-
-  AliESDEvent *fESD;                      //ESD object
-  AliAODEvent* fAOD;                      //AOD object
-  AliFlowAnalysisWithCumulants* fMyCumuAnalysis;  //Cumulant analysis object
-  AliFlowEventSimpleMaker* fEventMaker;   //FlowEventSimple maker object
-  TString fAnalysisType;                  //string to select which kind of input to analyse: ESD, AOD or MC
-  AliCFManager* fCFManager1;              // correction framework manager
-  AliCFManager* fCFManager2;              // correction framework manager
-  TList  *fListHistos;                    //collection of output //NEW
-   
-  ClassDef(AliAnalysisTaskCumulants, 1); // example of analysis
+  AliAnalysisTaskCumulants(const AliAnalysisTaskCumulants& aatc);
+  AliAnalysisTaskCumulants& operator=(const AliAnalysisTaskCumulants& aatc);
+
+  AliESDEvent *fESD;                             //ESD object
+  AliAODEvent* fAOD;                             //AOD object
+  AliFlowAnalysisWithCumulants* fCA;             //Cumulant Analysis (CA) object
+  AliFlowEventSimpleMaker* fEventMaker;          //FlowEventSimple maker object
+  TString fAnalysisType;                         //string to select which kind of input to analyse (ESD, AOD or MC)
+  AliCFManager* fCFManager1;                     //correction framework manager
+  AliCFManager* fCFManager2;                     //correction framework manager
+  TList  *fListHistos;                           //collection of output 
+     
+  ClassDef(AliAnalysisTaskCumulants, 1); 
 };
 
+//================================================================================================================
+
 #endif
 
 
index 2949985..be3e41e 100644 (file)
@@ -1,8 +1,26 @@
-//******************************* 
-// flow analysis with cumulants *   
-// author: Ante Bilandzic       * 
-// email: anteb@nikhef.nl       *
-//******************************* 
+/*************************************************************************
+* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  * 
+**************************************************************************/
+
+/********************************** 
+ * functions and equations needed * 
+ * for calculation of cumulants   *
+ * and final flow estimates       *
+ *                                *   
+ * author: Ante Bilandzic         * 
+ *          (anteb@nikhef.nl)     *
+ *********************************/ 
 
 #define AliCumulantsFunctions_cxx
 
 #include "AliFlowCommonConstants.h"
 #include "TChain.h"
 #include "TFile.h"
-#include "TList.h" //NEW
+#include "TList.h"
 #include "TParticle.h"
 
 #include "TProfile.h"
 #include "TProfile2D.h" 
 #include "TProfile3D.h"
+#include "TF1.h"//VOLOSHIN
+#include "TAxis.h"//VOLOSHIN
+#include "TH1.h"
+#include "TH1D.h"
+
 
 #include "AliFlowEventSimple.h"
 #include "AliFlowTrackSimple.h"
 #include "AliFlowCommonConstants.h"
 #include "AliCumulantsFunctions.h"
 
-
 ClassImp(AliCumulantsFunctions)
 
-//________________________________________________________________________
+//================================================================================================================_
 
 AliCumulantsFunctions::AliCumulantsFunctions():  
-  fIG(NULL),
-  fDGRe(NULL),
-  fDGIm(NULL),
-  fifr(NULL),
-  fdfr2(NULL), 
-  fdfr4(NULL), 
-  fdfr6(NULL),
-  fdfr8(NULL),
-  fAvMult(0)
- {
-   //default constructor 
- }
+ fIntGenFun(NULL),
+ fDiffGenFunRe(NULL),
+ fDiffGenFunIm(NULL),
+ fifr(NULL),
+ fdfr2(NULL), 
+ fdfr4(NULL), 
+ fdfr6(NULL),
+ fdfr8(NULL),
+ fAvMult(NULL),
+ fQVector(NULL),
+ fQDistrib(NULL),//q-distribution
+ fdRe0(NULL),
+ fdRe1(NULL),
+ fdRe2(NULL),
+ fdRe3(NULL),
+ fdRe4(NULL),
+ fdRe5(NULL),
+ fdRe6(NULL),
+ fdRe7(NULL),
+ fdIm0(NULL),
+ fdIm1(NULL),
+ fdIm2(NULL),
+ fdIm3(NULL),
+ fdIm4(NULL),
+ fdIm5(NULL),
+ fdIm6(NULL),
+ fdIm7(NULL)
+{
+ //default constructor 
+}
 
-AliCumulantsFunctions::~AliCumulantsFunctions(){
-  //desctructor
+AliCumulantsFunctions::~AliCumulantsFunctions()
+{
+ //destructor
 }
 
-AliCumulantsFunctions::AliCumulantsFunctions(TProfile2D *IntGen, TProfile3D *DiffGenRe, TProfile3D *DiffGenIm, TH1D *ifr, TH1D *dfr2, TH1D *dfr4, TH1D *dfr6, TH1D *dfr8, Double_t CvM):
-  fIG(IntGen),
-  fDGRe(DiffGenRe),
-  fDGIm(DiffGenIm),
-  fifr(ifr),
-  fdfr2(dfr2), 
-  fdfr4(dfr4), 
-  fdfr6(dfr6),
-  fdfr8(dfr8),
-  fAvMult(CvM)
- {
-   //custom constructor 
- }
+AliCumulantsFunctions::AliCumulantsFunctions(TProfile2D *IntGenFun, TProfile3D *DiffGenFunRe, TProfile3D *DiffGenFunIm, TH1D *ifr, TH1D *dfr2, TH1D *dfr4, TH1D *dfr6, TH1D *dfr8, TProfile *AvMult, TProfile *QVector, TH1D *QDistrib, TProfile2D *dRe0, TProfile2D *dRe1, TProfile2D *dRe2, TProfile2D *dRe3, TProfile2D *dRe4, TProfile2D *dRe5, TProfile2D *dRe6, TProfile2D *dRe7, TProfile2D *dIm0, TProfile2D *dIm1, TProfile2D *dIm2, TProfile2D *dIm3, TProfile2D *dIm4, TProfile2D *dIm5, TProfile2D *dIm6, TProfile2D *dIm7):
+ fIntGenFun(IntGenFun),
+ fDiffGenFunRe(DiffGenFunRe),
+ fDiffGenFunIm(DiffGenFunIm),
+ fifr(ifr),
+ fdfr2(dfr2), 
+ fdfr4(dfr4), 
+ fdfr6(dfr6),
+ fdfr8(dfr8),
+ fAvMult(AvMult),
+ fQVector(QVector),
+ fQDistrib(QDistrib),//q-distribution
+ fdRe0(dRe0),
+ fdRe1(dRe1),
+ fdRe2(dRe2),
+ fdRe3(dRe3),
+ fdRe4(dRe4),
+ fdRe5(dRe5),
+ fdRe6(dRe6),
+ fdRe7(dRe7),
+ fdIm0(dIm0),
+ fdIm1(dIm1),
+ fdIm2(dIm2),
+ fdIm3(dIm3),
+ fdIm4(dIm4),
+ fdIm5(dIm5),
+ fdIm6(dIm6),
+ fdIm7(dIm7)
+{
+ //custom constructor 
+}
   
-//___________________________________________________________________________
-void AliCumulantsFunctions::Calculate(){
- //calculate final flow estimates
-     static const Int_t fgkQmax=AliFlowCumuConstants::kQmax;//needed for numerics
-     static const Int_t fgkPmax=AliFlowCumuConstants::kPmax;//needed for numerics  
-     static const Int_t fgkFlow=AliFlowCumuConstants::kFlow;//integrated flow coefficient to be calculated
-     static const Int_t fgkMltpl=AliFlowCumuConstants::kMltpl;//the multiple in p=m*n (diff. flow) 
-     static const Int_t fgknBins=100;//number of pt bins //to be improved
-     Double_t fR0=AliFlowCumuConstants::fgR0;//needed for numerics
-     Double_t fPtMax=AliFlowCommonConstants::GetPtMax();//maximum pt
-     Double_t fPtMin=AliFlowCommonConstants::GetPtMin();//minimum pt
-     Double_t fBinWidth=(fPtMax-fPtMin)/fgknBins;//width of pt bin (in GeV)   
+//================================================================================================================
+
+void AliCumulantsFunctions::Calculate()
+{
+ //calculate cumulants and final integrated and differential flow estimates and store the results into output histograms
+ static const Int_t fgkQmax=AliFlowCumuConstants::kQmax;   //needed for numerics
+ static const Int_t fgkPmax=AliFlowCumuConstants::kPmax;   //needed for numerics  
+ static const Int_t fgkFlow=AliFlowCumuConstants::kFlow;   //integrated flow coefficient to be calculated
+ static const Int_t fgkMltpl=AliFlowCumuConstants::kMltpl; //the multiple in p=m*n (diff. flow) 
+ static const Int_t fgknBins=100;                          //number of pt bins //to be improved
+ Double_t fR0=AliFlowCumuConstants::fgR0;                  //needed for numerics
+ //Double_t fPtMax=AliFlowCommonConstants::GetPtMax();       //maximum pt
+ //Double_t fPtMin=AliFlowCommonConstants::GetPtMin();       //minimum pt
+ //Double_t fBinWidth=(fPtMax-fPtMin)/fgknBins;              //width of pt bin (in GeV)   
      
-  Double_t fBvG[fgkPmax][fgkQmax]={0.};
-        
-  for(Int_t p=0;p<fgkPmax;p++){
-   for(Int_t q=0;q<fgkQmax;q++){ 
-    fBvG[p][q]=fIG->GetBinContent(p+1,q+1);
-   }
-  } 
-      
-  /////////////////////////////////////////////////////////////////////////////      
-  //////////////////gen. function for the cumulants////////////////////////////
-  /////////////////////////////////////////////////////////////////////////////
-  
-  Double_t fC[fgkPmax][fgkQmax]={0.};
-  for (Int_t p=0;p<fgkPmax;p++){
-    for (Int_t q=0;q<fgkQmax;q++){
-      fC[p][q]=1.*fAvMult*(pow(fBvG[p][q],(1./fAvMult))-1.); //to be improved
-    }
+ //<G[p][q]>
+ Double_t AvG[fgkPmax][fgkQmax]={0.};                                  
+ for(Int_t p=0;p<fgkPmax;p++){
+  for(Int_t q=0;q<fgkQmax;q++){ 
+   AvG[p][q]=fIntGenFun->GetBinContent(p+1,q+1);
   }
+ } 
+
+ //avarage selected multiplicity
+ Double_t AvM = fAvMult->GetBinContent(1);
+ //mumber of events
+ Int_t nEvents = (Int_t)(fAvMult->GetBinEntries(1));
  
-  /////////////////////////////////////////////////////////////////////////////
-  ///////avaraging the gen. function for the cumulants over azimuth////////////
-  /////////////////////////////////////////////////////////////////////////////
+ //<Q-vector stuff>
+ Double_t AvQx  = fQVector->GetBinContent(1); //<Q_x>
+ Double_t AvQy  = fQVector->GetBinContent(2); //<Q_y>
+ Double_t AvQ2x = fQVector->GetBinContent(3); //<(Q_x)^2>
+ Double_t AvQ2y = fQVector->GetBinContent(4); //<(Q_y)^2>
   
-  Double_t fCAv[fgkPmax]={0.}; 
-  for (Int_t p=0;p<fgkPmax;p++){
-    Double_t fTempHere=0.; 
-    for (Int_t q=0;q<fgkQmax;q++){
-      fTempHere+=1.*fC[p][q];
-    } 
-    fCAv[p]=1.*fTempHere/fgkQmax;
-  }
+ /////////////////////////////////////////////////////////////////////////////      
+ //////////////////gen. function for the cumulants////////////////////////////
+ /////////////////////////////////////////////////////////////////////////////
   
-  /////////////////////////////////////////////////////////////////////////////
-  //////////////////////////////////final results//////////////////////////////
-  /////////////////////////////////////////////////////////////////////////////
+ Double_t C[fgkPmax][fgkQmax]={0.};//C[p][q]
+ for (Int_t p=0;p<fgkPmax;p++){
+  for (Int_t q=0;q<fgkQmax;q++){
+   C[p][q]=1.*AvM*(pow(AvG[p][q],(1./AvM))-1.); 
+  }
+ }
+ /////////////////////////////////////////////////////////////////////////////
+ ///////avaraging the gen. function for the cumulants over azimuth////////////
+ /////////////////////////////////////////////////////////////////////////////
   
-  Double_t fCumulant[fgkPmax];//c_{iFlow}\{2(p+1)\}
+ Double_t AvC[fgkPmax]={0.};//<C[p][q]>
+ for (Int_t p=0;p<fgkPmax;p++){
+  Double_t tempHere=0.; 
+  for (Int_t q=0;q<fgkQmax;q++){
+   tempHere+=1.*C[p][q];
+  } 
+  AvC[p]=1.*tempHere/fgkQmax;
+ }
+ /////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////final results//////////////////////////////
+ /////////////////////////////////////////////////////////////////////////////
+ Double_t cumulant[fgkPmax];//array to store various order cumulants
   
-  fCumulant[0]=(1./(fR0*fR0)) * (8.*fCAv[0] - 14.*fCAv[1] + (56./3.)*fCAv[2] - (35./2.)*fCAv[3] + 
-                             (56./5.)*fCAv[4] - (14./3.)*fCAv[5] + (8./7.)*fCAv[6] - (1./8.)*fCAv[7]);
+ cumulant[0] = (1./(fR0*fR0)) * (8.*AvC[0] - 14.*AvC[1] + (56./3.)*AvC[2] - (35./2.)*AvC[3] + 
+                             (56./5.)*AvC[4] - (14./3.)*AvC[5] + (8./7.)*AvC[6] - (1./8.)*AvC[7]);
 
-  fCumulant[1]=(1./pow(fR0,4.)) * ((-1924./35.)*fCAv[0] + (621./5.)*fCAv[1] - (8012./45.)*fCAv[2] + 
-                                (691./4.)*fCAv[3] - (564./5.)*fCAv[4] + (2143./45.)*fCAv[5] - 
-                                (412./35.)*fCAv[6] + (363./280.)*fCAv[7]);
+ cumulant[1] = (1./pow(fR0,4.)) * ((-1924./35.)*AvC[0] + (621./5.)*AvC[1] - (8012./45.)*AvC[2] + 
+                                (691./4.)*AvC[3] - (564./5.)*AvC[4] + (2143./45.)*AvC[5] - 
+                                (412./35.)*AvC[6] + (363./280.)*AvC[7]);
 
-  fCumulant[2]=(1./pow(fR0,6.)) * (349.*fCAv[0] - (18353./20.)*fCAv[1] + (7173./5.)*fCAv[2] - 
-                                1457.*fCAv[3] + (4891./5.)*fCAv[4] - (1683./4.)*fCAv[5] + 
-                                (527./5.)*fCAv[6] - (469./40.)*fCAv[7]);
+ cumulant[2] = (1./pow(fR0,6.)) * (349.*AvC[0] - (18353./20.)*AvC[1] + (7173./5.)*AvC[2] - 
+                                1457.*AvC[3] + (4891./5.)*AvC[4] - (1683./4.)*AvC[5] + 
+                                (527./5.)*AvC[6] - (469./40.)*AvC[7]);
 
-  fCumulant[3]=(1./pow(fR0,8.)) * ((-10528./5.)*fCAv[0] + (30578./5.)*fCAv[1] - (51456./5.)*fCAv[2] + 
-                                10993.*fCAv[3] - (38176./5.)*fCAv[4] + (16818./5.)*fCAv[5] - 
-                                (4288./5.)*fCAv[6] + (967./10.)*fCAv[7]);
+ cumulant[3] = (1./pow(fR0,8.)) * ((-10528./5.)*AvC[0] + (30578./5.)*AvC[1] - (51456./5.)*AvC[2] + 
+                                10993.*AvC[3] - (38176./5.)*AvC[4] + (16818./5.)*AvC[5] - 
+                                (4288./5.)*AvC[6] + (967./10.)*AvC[7]);
 
-  fCumulant[4]=(1./pow(fR0,10.)) * (11500.*fCAv[0] - 35800.*fCAv[1] + 63900.*fCAv[2] - 71600.*fCAv[3] + 
-                                 51620.*fCAv[4] - 23400.*fCAv[5] + 6100.*fCAv[6] - 700.*fCAv[7]);
+ cumulant[4] = (1./pow(fR0,10.)) * (11500.*AvC[0] - 35800.*AvC[1] + 63900.*AvC[2] - 71600.*AvC[3] + 
+                                 51620.*AvC[4] - 23400.*AvC[5] + 6100.*AvC[6] - 700.*AvC[7]);
 
-  fCumulant[5]=(1./pow(fR0,12.)) * (-52560.*fCAv[0] + 172080.*fCAv[1] - 321840.*fCAv[2] + 376200.*fCAv[3] - 
-                                 281520.*fCAv[4] + 131760.*fCAv[5] - 35280.*fCAv[6] + 4140.*fCAv[7]);
+ cumulant[5] = (1./pow(fR0,12.)) * (-52560.*AvC[0] + 172080.*AvC[1] - 321840.*AvC[2] + 376200.*AvC[3] - 
+                                 281520.*AvC[4] + 131760.*AvC[5] - 35280.*AvC[6] + 4140.*AvC[7]);
 
-  fCumulant[6]=(1./pow(fR0,14.)) * (176400.*fCAv[0] - 599760.*fCAv[1] + 1164240.*fCAv[2] - 1411200.*fCAv[3] + 
-                                 1093680.*fCAv[4] - 529200.*fCAv[5] + 146160.*fCAv[6] - 17640.*fCAv[7]);
+ cumulant[6] = (1./pow(fR0,14.)) * (176400.*AvC[0] - 599760.*AvC[1] + 1164240.*AvC[2] - 1411200.*AvC[3] + 
+                                 1093680.*AvC[4] - 529200.*AvC[5] + 146160.*AvC[6] - 17640.*AvC[7]);
 
-  fCumulant[7]=(1./pow(fR0,16.)) * (-322560*fCAv[0] + 1128960.*fCAv[1] - 2257920.*fCAv[2] + 2822400.*fCAv[3] - 
-                                 2257920.*fCAv[4] + 1128960.*fCAv[5] - 322560.*fCAv[6] + 40320.*fCAv[7]);
-  
+ cumulant[7] = (1./pow(fR0,16.)) * (-322560*AvC[0] + 1128960.*AvC[1] - 2257920.*AvC[2] + 2822400.*AvC[3] - 
+                                 2257920.*AvC[4] + 1128960.*AvC[5] - 322560.*AvC[6] + 40320.*AvC[7]);
+    
+ cout<<""<<endl;
+ cout<<"***************************"<<endl;
+ cout<<"cumulants:"<<endl;
   
-  cout<<""<<endl;
-  cout<<"***************************"<<endl;
-  cout<<"cumulants:"<<endl;
+ cout<<" c_"<<fgkFlow<<"{2} = "<<cumulant[0]<<endl; 
+ cout<<" c_"<<fgkFlow<<"{4} = "<<cumulant[1]<<endl;
+ cout<<" c_"<<fgkFlow<<"{6} = "<<cumulant[2]<<endl;
+ cout<<" c_"<<fgkFlow<<"{8} = "<<cumulant[3]<<endl; 
+ cout<<"c_"<<fgkFlow<<"{10} = "<<cumulant[4]<<endl; 
+ cout<<"c_"<<fgkFlow<<"{12} = "<<cumulant[5]<<endl;
+ cout<<"c_"<<fgkFlow<<"{14} = "<<cumulant[6]<<endl; 
+ cout<<"c_"<<fgkFlow<<"{16} = "<<cumulant[7]<<endl; 
   
-  cout<<" c_"<<fgkFlow<<"{2} = "<<fCumulant[0]<<endl; 
-  cout<<" c_"<<fgkFlow<<"{4} = "<<fCumulant[1]<<endl;
-  cout<<" c_"<<fgkFlow<<"{6} = "<<fCumulant[2]<<endl;
-  cout<<" c_"<<fgkFlow<<"{8} = "<<fCumulant[3]<<endl; 
-  cout<<"c_"<<fgkFlow<<"{10} = "<<fCumulant[4]<<endl; 
-  cout<<"c_"<<fgkFlow<<"{12} = "<<fCumulant[5]<<endl;
-  cout<<"c_"<<fgkFlow<<"{14} = "<<fCumulant[6]<<endl; 
-  cout<<"c_"<<fgkFlow<<"{16} = "<<fCumulant[7]<<endl; 
-  
-  cout<<""<<endl;
-  cout<<"integrated flow: "<<endl;
+ cout<<""<<endl;
+ cout<<"integrated flow: "<<endl;
+ Double_t V2=0.,V4=0.,V6=0.,V8=0.,V10=0.,V12=0.,V14=0.,V16=0.;
+ if(cumulant[0]>=0.){
+  V2 = pow(cumulant[0],(1./2.));
+ }
+ if(cumulant[1]<=0.){
+  V4 = pow(-cumulant[1],(1./4.));
+ }
+ if(cumulant[2]>=0.){
+  V6 = pow((1./4.)*cumulant[2],(1./6.));
+ }
+ if(cumulant[3]<=0.){
+  V8 = pow(-(1./33.)*cumulant[3],(1./8.));
+ }
+ if(cumulant[4]>=0.){
+  V10 = pow((1./456.)*cumulant[4],(1./10.));
+ }
+ if(cumulant[5]<=0.){
+  V12 = pow(-(1./9460.)*cumulant[5],(1./12.));
+ }
+ if(cumulant[6]>=0.){
+  V14 = pow((1./274800.)*cumulant[6],(1./14.));
+ }
+ if(cumulant[7]<=0.){
+  V16 = pow(-(1./10643745.)*cumulant[7],(1./16.));
+ }
+    
+ Double_t SdQ[4]={0.};
+ Double_t ChiQ[4]={0.};
+   
+   //v_2{2}
+   if((cumulant[0]>=0.) && (AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(pow(cumulant[0],(1./2.))*AvM,2.)>0.))
+   {        
+    ChiQ[0]=AvM*V2/pow(AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(V2*AvM,2.),0.5);
+    SdQ[0]=pow(((1./(2.*AvM*nEvents))*((1.+1.*pow(ChiQ[0],2))/(1.*pow(ChiQ[0],2)))),0.5);
+    cout<<" v_"<<fgkFlow<<"{2} = "<<100.*V2<<"%, chi{2} = "<<ChiQ[0]<<", sd{2} = "<<100.*SdQ[0]<<"%"<<endl;
+    fifr->SetBinContent(1,100.*V2);
+    fifr->SetBinError(1,100.*SdQ[0]);
+   } 
+   else 
+   {
+    cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl; 
+   }
+   
+   //v_2{4}   
+   if((cumulant[1]<=0.) && (AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(pow(-cumulant[1],(1./4.))*AvM,2.)>0.))
+   {
+    ChiQ[1]=AvM*V4/pow(AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(V4*AvM,2.),0.5);
+    SdQ[1]=(1./(pow(2.*AvM*nEvents,.5)))*pow((1.+2.*pow(ChiQ[1],2)+(1./4.)*pow(ChiQ[1],4.)+(1./4.)*pow(ChiQ[1],6.))/((1./4.)*pow(ChiQ[1],6.)),.5);
+    cout<<" v_"<<fgkFlow<<"{4} = "<<100.*V4<<"%, chi{4} = "<<ChiQ[1]<<", sd{4} = "<<100.*SdQ[1]<<"%"<<endl;
+    fifr->SetBinContent(2,100.*V4);
+    fifr->SetBinError(2,100.*SdQ[1]);
+   } 
+   else 
+   {
+    cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;  
+   } 
   
+  //v_2{6}
+  if((cumulant[2]>=0.) && (AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(pow((1./4.)*cumulant[2],(1./6.))*AvM,2.)>0.))
+  {
+   ChiQ[2]=AvM*V6/pow(AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(V6*AvM,2.),0.5);
+   SdQ[2]=(1./(pow(2.*AvM*nEvents,.5)))*pow((3.+18.*pow(ChiQ[2],2)+9.*pow(ChiQ[2],4.)+28.*pow(ChiQ[2],6.)+12.*pow(ChiQ[2],8.)+24.*pow(ChiQ[2],10.))/(24.*pow(ChiQ[2],10.)),.5);
+   cout<<" v_"<<fgkFlow<<"{6} = "<<100.*V6<<"%, chi{6} = "<<ChiQ[2]<<", sd{6} = "<<100.*SdQ[2]<<"%"<<endl;
+   fifr->SetBinContent(3,100.*V6);
+   fifr->SetBinError(3,100.*SdQ[2]);
+  }
+  else
+  {
+   cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;  
+  }
   
-  Double_t fV2=0.,fV4=0.,fV6=0.,fV8=0.;
-  Double_t fSdQ[4]={0.};
-  Double_t fChiQ[4]={0.};
-   if (fCumulant[0]>=0.){ 
-    fV2=sqrt(fCumulant[0]);    
-    //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.)>0.){
-     //fChiQ[0]=fAvM*fV2/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.),0.5);
-     //fSdQ[0]=pow(((1./(2.*fAvM*nEvents))*((1.+1.*pow(fChiQ[0],2))/(1.*pow(fChiQ[0],2)))),0.5);
-     cout<<" v_"<<fgkFlow<<"{2} = "<<100.*fV2<<"%, chi{2} = "<<fChiQ[0]<<", sd{2} = "<<100.*fSdQ[0]<<"%"<<endl;//to be improved (2->fgkFlow)
-     //fCommonHistsRes2->FillIntegratedFlow(100.*fV2,100.*fSdQ[0]);
-     //fCommonHistsRes2->FillChi(fChiQ[0]);
-     fifr->SetBinContent(1,100.*fV2);
-    //}
-   //} else {
-    //cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;  
-   } else {
-    //cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl; 
-   }
-  if (fCumulant[1]<=0.){
-    fV4=pow(-fCumulant[1],(1./4.));
-    //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.)>0.){
-     //fChiQ[1]=fAvM*fV4/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.),0.5);
-     //fSdQ[1]=(1./(pow(2.*fAvM*nEvents,.5)))*pow((1.+2.*pow(fChiQ[1],2)+(1./4.)*pow(fChiQ[1],4.)+(1./4.)*pow(fChiQ[1],6.))/((1./4.)*pow(fChiQ[1],6.)),.5);
-     cout<<" v_"<<fgkFlow<<"{4} = "<<100.*fV4<<"%, chi{4} = "<<fChiQ[1]<<", sd{4} = "<<100.*fSdQ[1]<<"%"<<endl;//to be improved (2->fgkFlow)
-     //fCommonHistsRes4->FillInteg8ratedFlow(100.*fV4,100.*fSdQ[1]);
-     //fCommonHistsRes4->FillChi(fChiQ[1]);
-     fifr->SetBinContent(2,100.*fV4);
-    //} else {
-      //cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
-    //}
- // } else {
-   // cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;  
+  //v_2{8}
+  if((cumulant[3]<=0.) && (AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(pow(-(1./33.)*cumulant[3],(1./8.))*AvM,2.)>0.))
+  {  
+   ChiQ[3]=AvM*V8/pow(AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(V8*AvM,2.),0.5);
+   SdQ[3]=(1./(pow(2.*AvM*nEvents,.5)))*pow((12.+96.*pow(ChiQ[3],2)+72.*pow(ChiQ[3],4.)+304.*pow(ChiQ[3],6.)+257.*pow(ChiQ[3],8.)+804.*pow(ChiQ[3],10.)+363.*pow(ChiQ[3],12.)+726.*pow(ChiQ[3],14.))/(726.*pow(ChiQ[3],14.)),.5);
+   cout<<" v_"<<fgkFlow<<"{8} = "<<100.*V8<<"%, chi{8} = "<<ChiQ[3]<<", sd{8} = "<<100.*SdQ[3]<<"%"<<endl;
+   fifr->SetBinContent(4,100.*V8);
+   fifr->SetBinError(4,100.*SdQ[3]);
   } 
-  if (fCumulant[2]>=0.){
-    fV6=pow((1./4.)*fCumulant[2],(1./6.));
-    //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV6*fAvM,2.)>0.){
-     //fChiQ[2]=fAvM*fV6/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV6*fAvM,2.),0.5);
-     //fSdQ[2]=(1./(pow(2.*fAvM*nEvents,.5)))*pow((3.+18.*pow(fChiQ[2],2)+9.*pow(fChiQ[2],4.)+28.*pow(fChiQ[2],6.)+12.*pow(fChiQ[2],8.)+24.*pow(fChiQ[2],10.))/(24.*pow(fChiQ[2],10.)),.5);
-     cout<<" v_"<<fgkFlow<<"{6} = "<<100.*fV6<<"%, chi{6} = "<<fChiQ[2]<<", sd{6} = "<<100.*fSdQ[2]<<"%"<<endl;//to be improved (2->fgkFlow)
-     //fCommonHistsRes6->FillIntegratedFlow(100.*fV6,100.*fSdQ[2]);
-     //fCommonHistsRes6->FillChi(fChiQ[2]);
-     fifr->SetBinContent(3,100.*fV6);
-    //} else {
-     // cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl; 
-   // }
-  } else {
-    //cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;  
+  else 
+  {
+   cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;     
   }
-  if (fCumulant[3]<=0.){
-    fV8=pow(-(1./33.)*fCumulant[3],(1./8.));
-    //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV8*fAvM,2.)>0.){
-     //fChiQ[3]=fAvM*fV8/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV8*fAvM,2.),0.5);
-     //fSdQ[3]=(1./(pow(2.*fAvM*nEvents,.5)))*pow((12.+96.*pow(fChiQ[3],2)+72.*pow(fChiQ[3],4.)+304.*pow(fChiQ[3],6.)+257.*pow(fChiQ[3],8.)+804.*pow(fChiQ[3],10.)+363.*pow(fChiQ[3],12.)+726.*pow(fChiQ[3],14.))/(726.*pow(fChiQ[3],14.)),.5);
-      cout<<" v_"<<fgkFlow<<"{8} = "<<100.*fV8<<"%, chi{8} = "<<fChiQ[3]<<", sd{8} = "<<100.*fSdQ[3]<<"%"<<endl;//to be improved (2->fgkFlow)
-      //fCommonHistsRes8->FillIntegratedFlow(100.*fV8,100.*fSdQ[3]);
-      //fCommonHistsRes8->FillChi(fChiQ[3]);
-      fifr->SetBinContent(4,100.*fV8);
-     //} else {
-       //cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;a
-     //}
-  } else {
-    //cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;     
-  }
-  if (fCumulant[4]>=0.){//fHistProIntFlow
-    cout<<"v_"<<fgkFlow<<"{10} = "<<100.*pow((1./456.)*fCumulant[4],(1./10.))<<"%"<<endl;//to be improved (2->fgkFlow)
-    fifr->SetBinContent(5,100.*pow((1./456.)*fCumulant[4],(1./10.)));
+  
+  //v_2{10}
+  if (cumulant[4]>=0.){
+    cout<<"v_"<<fgkFlow<<"{10} = "<<100.*V10<<"%"<<endl;
+    fifr->SetBinContent(5,100.*pow((1./456.)*cumulant[4],(1./10.)));
   } else {
-      cout<<"v_"<<fgkFlow<<"{10} = Im"<<endl;  //to be improved (2->fgkFlow)
+      cout<<"v_"<<fgkFlow<<"{10} = Im"<<endl; 
   }
-  if (fCumulant[5]<=0.){
-    cout<<"v_"<<fgkFlow<<"{12} = "<<100.*pow(-(1./9460.)*fCumulant[5],(1./12.))<<"%"<<endl;//to be improved (2->fgkFlow)
-    fifr->SetBinContent(6,100.*pow(-(1./9460.)*fCumulant[5],(1./12.)));
+  
+  //v_2{12}
+  if (cumulant[5]<=0.){
+    cout<<"v_"<<fgkFlow<<"{12} = "<<100.*V12<<"%"<<endl;
+    fifr->SetBinContent(6,100.*pow(-(1./9460.)*cumulant[5],(1./12.)));
   } else {
-    cout<<"v_"<<fgkFlow<<"{12} = Im"<<endl;  //to be improved (2->fgkFlow)
+    cout<<"v_"<<fgkFlow<<"{12} = Im"<<endl; 
   }
-  if (fCumulant[6]>=0.){
-    cout<<"v_"<<fgkFlow<<"{14} = "<<100.*pow((1./274800.)*fCumulant[6],(1./14.))<<"%"<<endl;//to be improved (2->fgkFlow)
-    fifr->SetBinContent(7,100.*pow((1./274800.)*fCumulant[6],(1./14.)));
+  
+  //v_2{14}
+  if (cumulant[6]>=0.){
+    cout<<"v_"<<fgkFlow<<"{14} = "<<100.*V14<<"%"<<endl;
+    fifr->SetBinContent(7,100.*pow((1./274800.)*cumulant[6],(1./14.)));
   } else {
-    cout<<"v_"<<fgkFlow<<"{14} = Im"<<endl;  //to be improved (2->fgkFlow)
+    cout<<"v_"<<fgkFlow<<"{14} = Im"<<endl;  
   }
-  if (fCumulant[7]<=0.){
-    cout<<"v_"<<fgkFlow<<"{16} = "<<100.*pow(-(1./10643745.)*fCumulant[7],(1./16.))<<"%"<<endl;//to be improved (2->fgkFlow)
-    fifr->SetBinContent(8,100.*pow(-(1./10643745.)*fCumulant[7],(1./16.)));
+  
+  //v_2{16}
+  if (cumulant[7]<=0.){
+    cout<<"v_"<<fgkFlow<<"{16} = "<<100.*V16<<"%"<<endl;
+    fifr->SetBinContent(8,100.*pow(-(1./10643745.)*cumulant[7],(1./16.)));
   } else {
-    cout<<"v_"<<fgkFlow<<"{16} = Im"<<endl;  //to be improved (2->fgkFlow)
+    cout<<"v_"<<fgkFlow<<"{16} = Im"<<endl;  
   }
-  //cout<<"***************************"<<endl;
-      
-  //DIFFERENTIAL FLOW CALCULATIONS STARTS HERE!!!
+  cout<<"***************************"<<endl;    
+  
   
-  Double_t fX[fgknBins][fgkPmax][fgkQmax]={0.};//see the text bellow relation (11) in PG
-  Double_t fY[fgknBins][fgkPmax][fgkQmax]={0.};
+ /////////////////////////////////////////////////////////////////////////////
+ ///////////////////////DIFFERENTIAL FLOW CALCULATIONS////////////////////////
+ /////////////////////////////////////////////////////////////////////////////
   
+  Double_t X[fgknBins][fgkPmax][fgkQmax]={0.};
+  Double_t Y[fgknBins][fgkPmax][fgkQmax]={0.};
+  
+  /*
+  //3D profiles
   for(Int_t b=0;b<fgknBins;b++){
     for(Int_t p=0;p<fgkPmax;p++){
       for(Int_t q=0;q<fgkQmax;q++){
-       fX[b][p][q]=fDGRe->GetBinContent(b+1,p+1,q+1)/fBvG[p][q];
-       fY[b][p][q]=fDGIm->GetBinContent(b+1,p+1,q+1)/fBvG[p][q];
+       X[b][p][q]=fDiffGenFunRe->GetBinContent(b+1,p+1,q+1)/AvG[p][q];
+       Y[b][p][q]=fDiffGenFunIm->GetBinContent(b+1,p+1,q+1)/AvG[p][q];
       }
     }   
   } 
+  */
+
+  for(Int_t b=0;b<fgknBins;b++){
+    //for(Int_t p=0;p<fgkPmax;p++){
+      for(Int_t q=0;q<fgkQmax;q++){
+       X[b][0][q]=fdRe0->GetBinContent(b+1,q+1)/AvG[0][q];
+       Y[b][0][q]=fdIm0->GetBinContent(b+1,q+1)/AvG[0][q];
+       //--------------------------------------------------
+       X[b][1][q]=fdRe1->GetBinContent(b+1,q+1)/AvG[1][q];
+       Y[b][1][q]=fdIm1->GetBinContent(b+1,q+1)/AvG[1][q];
+       //--------------------------------------------------
+       X[b][2][q]=fdRe2->GetBinContent(b+1,q+1)/AvG[2][q];
+       Y[b][2][q]=fdIm2->GetBinContent(b+1,q+1)/AvG[2][q];
+       //--------------------------------------------------
+       X[b][3][q]=fdRe3->GetBinContent(b+1,q+1)/AvG[3][q];
+       Y[b][3][q]=fdIm3->GetBinContent(b+1,q+1)/AvG[3][q];
+       //--------------------------------------------------
+       X[b][4][q]=fdRe4->GetBinContent(b+1,q+1)/AvG[4][q];
+       Y[b][4][q]=fdIm4->GetBinContent(b+1,q+1)/AvG[4][q];
+       //--------------------------------------------------
+       X[b][5][q]=fdRe5->GetBinContent(b+1,q+1)/AvG[5][q];
+       Y[b][5][q]=fdIm5->GetBinContent(b+1,q+1)/AvG[5][q];
+       //--------------------------------------------------
+       X[b][6][q]=fdRe6->GetBinContent(b+1,q+1)/AvG[6][q];
+       Y[b][6][q]=fdIm6->GetBinContent(b+1,q+1)/AvG[6][q];
+       //--------------------------------------------------
+       X[b][7][q]=fdRe7->GetBinContent(b+1,q+1)/AvG[7][q];
+       Y[b][7][q]=fdIm7->GetBinContent(b+1,q+1)/AvG[7][q];
+      }
+    //}   
+  }
   
-  Double_t fD[fgknBins][fgkPmax]={0.};//implementing relation (11) from PG
+  Double_t D[fgknBins][fgkPmax]={0.};
   
   for (Int_t b=0;b<fgknBins;b++){
     for (Int_t p=0;p<fgkPmax;p++){
       Double_t fTempHere3=0.; 
       for (Int_t q=0;q<fgkQmax;q++){
-       fTempHere3+=cos(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*fX[b][p][q] + sin(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*fY[b][p][q];
+       fTempHere3+=cos(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*X[b][p][q] + sin(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*Y[b][p][q];
       } 
-      fD[b][p]=1.*(pow(fR0*pow(p+1.,.5),fgkMltpl)/fgkQmax)*fTempHere3;
-      if(fD[b][p]){
-       //cout<<"And this "<<b<<" "<<fD[b][p]<<endl;
-      }  
+      D[b][p]=1.*(pow(fR0*pow(p+1.,.5),fgkMltpl)/fgkQmax)*fTempHere3;
     }
   } 
   
-  Double_t fDiffCumulant2[fgknBins]={0.};//implementing relation (12) from PG
-  Double_t fDiffCumulant4[fgknBins]={0.};
-  Double_t fDiffCumulant6[fgknBins]={0.};
-  Double_t fDiffCumulant8[fgknBins]={0.};
+  Double_t DiffCumulant2[fgknBins]={0.};
+  Double_t DiffCumulant4[fgknBins]={0.};
+  Double_t DiffCumulant6[fgknBins]={0.};
+  Double_t DiffCumulant8[fgknBins]={0.};
   
   for (Int_t b=0;b<fgknBins;b++){
-    fDiffCumulant2[b]=(1./(fR0*fR0))*(4.*fD[b][0]-3.*fD[b][1]+(4./3.)*fD[b][2]-(1./4.)*fD[b][3]);
-    fDiffCumulant4[b]=(1./pow(fR0,4.))*((-26./3.)*fD[b][0]+(19./2.)*fD[b][1]-(14./3.)*fD[b][2]+(11./12.)*fD[b][3]);
-    fDiffCumulant6[b]=(1./pow(fR0,6.))*(18.*fD[b][0]-24.*fD[b][1]+14.*fD[b][2]-3.*fD[b][3]);
-    fDiffCumulant8[b]=(1./pow(fR0,8.))*(-24.*fD[b][0]+36.*fD[b][1]-24.*fD[b][2]+6.*fD[b][3]);
+    DiffCumulant2[b]=(1./(fR0*fR0))*(4.*D[b][0]-3.*D[b][1]+(4./3.)*D[b][2]-(1./4.)*D[b][3]);
+    DiffCumulant4[b]=(1./pow(fR0,4.))*((-26./3.)*D[b][0]+(19./2.)*D[b][1]-(14./3.)*D[b][2]+(11./12.)*D[b][3]);
+    DiffCumulant6[b]=(1./pow(fR0,6.))*(18.*D[b][0]-24.*D[b][1]+14.*D[b][2]-3.*D[b][3]);
+    DiffCumulant8[b]=(1./pow(fR0,8.))*(-24.*D[b][0]+36.*D[b][1]-24.*D[b][2]+6.*D[b][3]);
   }
   
-  Double_t fv2[fgknBins],fv4[fgknBins],fv6[fgknBins],fv8[fgknBins];
-  //Double_t fAvPt[fgknBins];
-  Double_t fSddiff2[fgknBins],fSddiff4[fgknBins];
+  Double_t v2[fgknBins],v4[fgknBins],v6[fgknBins],v8[fgknBins];
+  //Double_t Sddiff2[fgknBins],Sddiff4[fgknBins];
 
-  cout<<"number of pt bins: "<<fgknBins<<endl;
-  cout<<"****************************************"<<endl;
+  //cout<<"number of pt bins: "<<fgknBins<<endl;
+  //cout<<"****************************************"<<endl;
   for (Int_t b=0;b<fgknBins;b++){ 
-    //if(fBinNoOfParticles[b]==0)continue;
-    //fAvPt[b]=fBinMeanPt[b]/fBinNoOfParticles[b];
-    cout<<"pt bin: "<<b*fBinWidth<<"-"<<(b+1)*fBinWidth<<" GeV"<<endl;
-    //cout<<"number of particles in this pt bin: "<<tempNo->GetBinContent(b+1,1,1)<<endl;
-    //cout<<"mean pt in this bin: "<<fAvPt[b]<<" GeV"<<endl;
-    if(fCumulant[0]>=0){
-      fv2[b]=100.*fDiffCumulant2[b]/pow(fCumulant[0],.5);
-      //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.)>0.){
-       //fSddiff2[b]=pow((1./(2.*fBinNoOfParticles[b]))*((1.+pow(fChiQ[0],2.))/pow(fChiQ[0],2.)),0.5);
-       cout<<"v'_2/2{2} = "<<fv2[b]<<"%, "<<" "<<"sd{2} = "<<100.*fSddiff2[b]<<"%"<<endl;
-       //fDiffFlowResults2->SetBinContent(b+1,fv2[b],100.*fSddiff2[b]);
-       fdfr2->SetBinContent(b+1,fv2[b]);
-      //} else {
+    //cout<<"pt bin: "<<b*fBinWidth<<"-"<<(b+1)*fBinWidth<<" GeV"<<endl;
+    
+    //v'_{2/2}{2}
+    if(cumulant[0]>=0)
+    {
+      v2[b]=100.*DiffCumulant2[b]/pow(cumulant[0],.5);
+      if (AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(V2*AvM,2.)>0.)
+      {
+       //Sddiff2[b]=pow((1./(2.*fBinNoOfParticles[b]))*((1.+pow(ChiQ[0],2.))/pow(ChiQ[0],2.)),0.5);
+       //cout<<"v'_2/2{2} = "<<v2[b]<<"%, "<<" "<<"sd{2} = "<<100.*Sddiff2[b]<<"%"<<endl;
+       //DiffFlowResults2->SetBinContent(b+1,v2[b],100.*Sddiff2[b]);
+       fdfr2->SetBinContent(b+1,v2[b]);
+      } else {
          //cout<<"v'_2/2{2} = Im"<<endl;
-      //}
+      }
     }else{
-      cout<<"v'_2/2{2} = Im"<<endl;
+      //cout<<"v'_2/2{2} = Im"<<endl;
     } 
-    if(fCumulant[1]<=0){
-      fv4[b]=-100.*fDiffCumulant4[b]/pow(-fCumulant[1],.75);
-      //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.)>0.){
-       //fSddiff4[b]=pow((1./(2.*fBinNoOfParticles[b]))*((2.+6.*pow(fChiQ[1],2.)+pow(fChiQ[1],4.)+pow(fChiQ[1],6.))/pow(fChiQ[1],6.)),0.5);
-       cout<<"v'_2/2{4} = "<<fv4[b]<<"%, "<<" "<<"sd{4} = "<<100.*fSddiff4[b]<<"%"<<endl;
-       //fCommonHistsRes4->FillDifferentialFlow(b+1,fv4[b],100.*fSddiff4[b]);
-       fdfr4->SetBinContent(b+1,fv4[b]);
-      //} else {
-       // cout<<"v'_2/2{4} = Im"<<endl;
-      //} 
+    
+    //v'_{2/2}{4}
+    if(cumulant[1]<=0)
+    {
+      v4[b]=-100.*DiffCumulant4[b]/pow(-cumulant[1],.75);
+      if (AvQ2x+AvQ2y-pow(AvQx,2.)-pow(AvQy,2.)-pow(V4*AvM,2.)>0.)
+      {
+       //Sddiff4[b]=pow((1./(2.*fBinNoOfParticles[b]))*((2.+6.*pow(ChiQ[1],2.)+pow(ChiQ[1],4.)+pow(ChiQ[1],6.))/pow(ChiQ[1],6.)),0.5);
+       //cout<<"v'_2/2{4} = "<<v4[b]<<"%, "<<" "<<"sd{4} = "<<100.*Sddiff4[b]<<"%"<<endl;
+       //fCommonHistsRes4->FillDifferentialFlow(b+1,v4[b],100.*Sddiff4[b]);
+       fdfr4->SetBinContent(b+1,v4[b]);
+      } else {
+         //cout<<"v'_2/2{4} = Im"<<endl;
+      } 
     }else{
-      cout<<"v'_2/2{4} = Im"<<endl;
+      //cout<<"v'_2/2{4} = Im"<<endl;
     }  
-    if(fCumulant[2]>=0){
-      cout<<"v'_2/2{6} = "<<100.*fDiffCumulant6[b]/(4.*pow((1./4.)*fCumulant[2],(5./6.)))<<"%"<<endl;
-      fv6[b]=100.*fDiffCumulant6[b]/(4.*pow((1./4.)*fCumulant[2],(5./6.)));
-      //fCommonHistsRes6->FillDifferentialFlow(b+1,fv6[b],0.);
-      fdfr6->SetBinContent(b+1,fv6[b]);
+    
+    //v'_{2/2}{6}
+    if(cumulant[2]>=0){
+      //cout<<"v'_2/2{6} = "<<100.*DiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)))<<"%"<<endl;
+      v6[b]=100.*DiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)));
+      //fCommonHistsRes6->FillDifferentialFlow(b+1,v6[b],0.);
+      fdfr6->SetBinContent(b+1,v6[b]);
     }else{
-      cout<<"v'_2/2{6} = Im"<<endl;
+      //cout<<"v'_2/2{6} = Im"<<endl;
     }     
-    if(fCumulant[3]<=0){
-      cout<<"v'_2/2{8} = "<<-100.*fDiffCumulant8[b]/(33.*pow(-(1./33.)*fCumulant[3],(7./8.)))<<"%"<<endl;
-      fv8[b]=-100.*fDiffCumulant8[b]/(33.*pow(-(1./33.)*fCumulant[3],(7./8.))); 
-      //fCommonHistsRes8->FillDifferentialFlow(b+1,fv8[b],0.);
-      fdfr8->SetBinContent(b+1,fv8[b]);
+    
+    //v'_{2/2}{8}
+    if(cumulant[3]<=0){
+      //cout<<"v'_2/2{8} = "<<-100.*DiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.)))<<"%"<<endl;
+      v8[b]=-100.*DiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.))); 
+      //fCommonHistsRes8->FillDifferentialFlow(b+1,v8[b],0.);
+      fdfr8->SetBinContent(b+1,v8[b]);
     }else{
-      cout<<"v'_2/2{8} = Im"<<endl;
+      //cout<<"v'_2/2{8} = Im"<<endl;
     }       
-    cout<<"****************************************"<<endl;
+    //cout<<"****************************************"<<endl;
   }  
+  
 }
 
   
index 5ebd4a4..96e579f 100644 (file)
@@ -1,14 +1,21 @@
+/* 
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. 
+ * See cxx source for full Copyright notice 
+ * $Id$ 
+ */
+
+/********************************** 
+ * functions and equations needed * 
+ * for calculation of cumulants   *
+ * and final flow estimates       *
+ *                                *   
+ * author: Ante Bilandzic         * 
+ *          (anteb@nikhef.nl)     *
+ *********************************/ 
+
 #ifndef AliCumulantsFunctions_H
 #define AliCumulantsFunctions_H
 
-//********************************* 
-// functions and equations needed * 
-// for calcualation of cumulants  *
-// and final flow estimates       *   
-// author: Ante Bilandzic         * 
-// email: anteb@nikhef.nl         *
-//********************************* 
-
 #include "AliFlowCommonConstants.h"
 #include "AliFlowCumuConstants.h"
 
@@ -21,29 +28,42 @@ class TObjArray;
 class TList;
 class TFile;
 
+//================================================================================================================
+
 class AliCumulantsFunctions{
  public:
   AliCumulantsFunctions();
   virtual ~AliCumulantsFunctions();
-  AliCumulantsFunctions(TProfile2D *IntGen, TProfile3D *DiffGenRe, TProfile3D *DiffGenIm, TH1D *ifr, TH1D *dfr2, TH1D *dfr4, TH1D *dfr6, TH1D *dfr8, Double_t CvM);
+  AliCumulantsFunctions(TProfile2D *IntGenFun, TProfile3D *DiffGenFunRe, TProfile3D *DiffGenFunIm, TH1D *ifr, TH1D *dfr2, TH1D *dfr4, TH1D *dfr6, TH1D *dfr8, TProfile *AvMult, TProfile *fQVector, TH1D *fQDistrib, TProfile2D *fdRe0, TProfile2D *fdRe1, TProfile2D *fdRe2, TProfile2D *fdRe3, TProfile2D *fdRe4, TProfile2D *fdRe5, TProfile2D *fdRe6, TProfile2D *fdRe7, TProfile2D *fdIm0, TProfile2D *fdIm1, TProfile2D *fdIm2, TProfile2D *fdIm3, TProfile2D *fdIm4, TProfile2D *fdIm5, TProfile2D *fdIm6, TProfile2D *fdIm7);
   void Calculate();
 
  private:
   AliCumulantsFunctions(const AliCumulantsFunctions& fun);
   AliCumulantsFunctions& operator=(const AliCumulantsFunctions& fun);
   
-  TProfile2D *fIG;   //
-  TProfile3D *fDGRe; //
-  TProfile3D *fDGIm; //
-  TH1D *fifr;       //integrated flow final results 
-  TH1D *fdfr2;      //differential flow final results //to be improved
-  TH1D *fdfr4;      //differential flow final results //to be improved
-  TH1D *fdfr6;      //differential flow final results //to be improved
-  TH1D *fdfr8;      //differential flow final results //to be improved
-  Double_t fAvMult; //avarage multiplicity
+  TProfile2D *fIntGenFun;    //avarage value of generating function for int. flow
+  TProfile3D *fDiffGenFunRe; //avarage value of generating function for diff. flow (real part)
+  TProfile3D *fDiffGenFunIm; //avarage value of generating function for diff. flow (imaginary part)
+  TH1D *fifr;                //integrated flow final results 
+  TH1D *fdfr2;               //differential flow final results
+  TH1D *fdfr4;               //differential flow final results
+  TH1D *fdfr6;               //differential flow final results
+  TH1D *fdfr8;               //differential flow final results 
+   
+  TProfile *fAvMult;         //avarage selected multiplicity for int. flow
+  TProfile *fQVector;        //avarage values of Q-vector components (1st bin: <Q_x>, 2nd bin: <Q_y>, 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>) 
+  
+  TH1D     *fQDistrib;       //q-distribution
+  
+  TProfile2D *fdRe0,*fdRe1,*fdRe2,*fdRe3,*fdRe4,*fdRe5,*fdRe6,*fdRe7;//differential flow 
+  TProfile2D *fdIm0,*fdIm1,*fdIm2,*fdIm3,*fdIm4,*fdIm5,*fdIm6,*fdIm7;//differential flow
   
   ClassDef(AliCumulantsFunctions, 0);
 };
+
+//================================================================================================================
+
 #endif
 
 
index a70f428..f95b83d 100644 (file)
@@ -1,3 +1,25 @@
+/*************************************************************************
+* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  * 
+**************************************************************************/
+
+/******************************** 
+ * flow analysis with cumulants * 
+ *                              *
+ * author: Ante Bilandzic       * 
+ *          (anteb@nikhef.nl)   *
+ *******************************/ 
+
 #define AliFlowAnalysisWithCumulants_cxx
 
 #include "Riostream.h"
 #include "TFile.h"
 #include "TList.h" //NEW
 #include "TParticle.h"
-
-
-
-
-
-
-
+#include "TRandom3.h"
 #include "TProfile.h"
 #include "TProfile2D.h" 
 #include "TProfile3D.h"
-
-
-
-
-
-
 #include "AliFlowEventSimple.h"
 #include "AliFlowTrackSimple.h"
 #include "AliFlowAnalysisWithCumulants.h"
 #include "AliFlowCumuConstants.h"
 
 class TH1;
-class TH3;
 class TGraph;
 class TPave;
 class TLatex;
@@ -41,256 +50,369 @@ class TList;
 class TCanvas;
 class TSystem;
 class TROOT;
-
 class AliFlowVector;
 class TVector;
 
-//******************************* 
-// flow analysis with cumulants *   
-// author: Ante Bilandzic       * 
-// email: anteb@nikhef.nl       *
-//******************************* 
+//================================================================================================================
 
 ClassImp(AliFlowAnalysisWithCumulants)
 
-//________________________________________________________________________
-
 AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants():  
-  fTrack(NULL),
-  //fHistFileName("cumulants.root"),//NEW: In the old version theis was not commented out
-  //fHistFilefHistFile(NULL),//NEW: In the old version theis was not commented out
-  fHistList(NULL), //NEW
-  fAvM(0),
-  fR0(0),
-  fPtMax(0),
-  fPtMin(0),
-  fBinWidth(0),
-  fAvQx(0),
-  fAvQy(0),
-  fAvQ2x(0),
-  fAvQ2y(0),
-  fNumberOfEvents(0),
-  
-  
-  
-  fHistProAvM(NULL),
-  fIntFlowResults(NULL),
-  fDiffFlowResults2(NULL),//to be improved
-  fDiffFlowResults4(NULL),
-  fDiffFlowResults6(NULL),
-  fDiffFlowResults8(NULL),
-  fIntFlowGenFun(NULL),
-  fDiffFlowGenFunRe(NULL),
-  fDiffFlowGenFunIm(NULL),
-  
-  
-  
-  fCommonHists(NULL)
-  //fCommonHistsRes2(NULL),//to be improved
-  //fCommonHistsRes4(NULL),//to be improved
-  //fCommonHistsRes6(NULL),//to be improved
-  //fCommonHistsRes8(NULL)//to be improved
- {
-   //constructor 
-   
-   fHistList = new TList(); //NEW
-   
-   fR0=AliFlowCumuConstants::fgR0;
-   fPtMax=AliFlowCommonConstants::GetPtMax(); 
-   fPtMin=AliFlowCommonConstants::GetPtMin();
-   fBinWidth=(fPtMax-fPtMin)/fgknBins;
-   
- }
+ fTrack(NULL),
+ fHistList(NULL),
+ fR0(0),
+ fPtMax(0),
+ fPtMin(0),
+ fBinWidth(0),
+ fAvQx(0),
+ fAvQy(0),
+ fAvQ2x(0),
+ fAvQ2y(0),
+ fAvMultIntFlow(NULL),
+ fQVectorComponents(NULL),
+ fIntFlowResults(NULL),
+ fDiffFlowResults2(NULL),
+ fDiffFlowResults4(NULL),
+ fDiffFlowResults6(NULL),
+ fDiffFlowResults8(NULL),
+ fIntFlowGenFun(NULL),
+ fDiffFlowGenFunRe(NULL),
+ fDiffFlowGenFunIm(NULL),
+ fDiffFlowGenFunRe0(NULL),
+ fDiffFlowGenFunRe1(NULL),
+ fDiffFlowGenFunRe2(NULL),
+ fDiffFlowGenFunRe3(NULL),
+ fDiffFlowGenFunRe4(NULL),
+ fDiffFlowGenFunRe5(NULL),
+ fDiffFlowGenFunRe6(NULL),
+ fDiffFlowGenFunRe7(NULL),
+ fDiffFlowGenFunIm0(NULL),
+ fDiffFlowGenFunIm1(NULL),
+ fDiffFlowGenFunIm2(NULL),
+ fDiffFlowGenFunIm3(NULL),
+ fDiffFlowGenFunIm4(NULL),
+ fDiffFlowGenFunIm5(NULL),
+ fDiffFlowGenFunIm6(NULL),
+ fDiffFlowGenFunIm7(NULL),
+ fCommonHists(NULL),
+ fQDist(NULL)//q-distribution
+{
+ //constructor 
+ fHistList = new TList(); 
+ fR0=AliFlowCumuConstants::fgR0;
+ fPtMax=AliFlowCommonConstants::GetPtMax(); 
+ fPtMin=AliFlowCommonConstants::GetPtMin();
+ fBinWidth=(fPtMax-fPtMin)/fgknBins;  
+}
 
-AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants(){
-  //desctructor
-   delete fHistList; //NEW
+AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants()
+{
+ //desctructor
+ delete fHistList; 
 }
 
-  
-//___________________________________________________________________________
-void AliFlowAnalysisWithCumulants::CreateOutputObjects(){
- //output histograms
+//================================================================================================================
+
+void AliFlowAnalysisWithCumulants::CreateOutputObjects()
+{
+ //various output histograms
  
-  fHistProAvM = new TProfile("fHistProAvM","Avarage Multiplicity",1,0,1,0,100000);
-  fHistProAvM->SetXTitle("");
-  fHistProAvM->SetYTitle("Avarage Multiplicity");
-  fHistList->Add(fHistProAvM);
-  
-  fIntFlowResults = new TH1D("fIntFlowResults","Integrated Flow From Cumulants",8,0,8);
-  fIntFlowResults->SetXTitle("");
-  fIntFlowResults->SetYTitle("Integrated Flow [%]");
-  fHistList->Add(fIntFlowResults);
+ //avarage multiplicity 
+ fAvMultIntFlow = new TProfile("fAvMultIntFlow","Avarage multiplicity of selected particles used for int. flow",1,0,1,0,100000);
+ fAvMultIntFlow->SetXTitle("");
+ fAvMultIntFlow->SetYTitle("<multiplicity>");
+ fHistList->Add(fAvMultIntFlow);
  
-  //to be improved, I should store the results as CommonHistResults
-  fDiffFlowResults2 = new TH1D("fDiffFlowResults2","v'_2/2{2}",fgknBins,fPtMin,fPtMax);
-  fDiffFlowResults2->SetXTitle("pt [GeV]");
-  fDiffFlowResults2->SetYTitle("Differential Flow [%]");
-  fHistList->Add(fDiffFlowResults2);
+ //Q-vector stuff
+ fQVectorComponents = new TProfile("fQVectorComponents","Avarage of Q-vector components",4,0.,4.);
+ fQVectorComponents->SetXTitle("");
+ fQVectorComponents->SetYTitle("< >");
+ fHistList->Add(fQVectorComponents);
  
-  //to be improved, I should store the results as CommonHistResults
-  fDiffFlowResults4 = new TH1D("fDiffFlowResults4","v'_2/2{4}",fgknBins,fPtMin,fPtMax);
-  fDiffFlowResults4->SetXTitle("pt [GeV]");
-  fDiffFlowResults4->SetYTitle("Differential Flow [%]");
-  fHistList->Add(fDiffFlowResults4);
+ //final results for integrated flow (v_n{2}, v_n{4},..., v_n{16}) from cumulants (by default n=2) 
+ fIntFlowResults = new TH1D("fIntFlowResults","Integrated Flow From Cumulants",8,0,8);
+ fIntFlowResults->SetXTitle("");
+ fIntFlowResults->SetYTitle("Integrated Flow [%]");
+ fHistList->Add(fIntFlowResults);
   
-  //to be improved, I should store the results as CommonHistResults
-  fDiffFlowResults6 = new TH1D("fDiffFlowResults6","v'_2/2{6}",fgknBins,fPtMin,fPtMax);
-  fDiffFlowResults6->SetXTitle("pt [GeV]");
-  fDiffFlowResults6->SetYTitle("Differential Flow [%]");
-  fHistList->Add(fDiffFlowResults6);
-  
-  //to be improved, I should store the results as CommonHistResults
-  fDiffFlowResults8 = new TH1D("fDiffFlowResults8","v'_2/2{8}",fgknBins,fPtMin,fPtMax);
-  fDiffFlowResults8->SetXTitle("pt [GeV]");
-  fDiffFlowResults8->SetYTitle("Differential Flow [%]");
-  fHistList->Add(fDiffFlowResults8);
+ //final results for differential flow v_p/n{2} (by default p=n=2)
+ fDiffFlowResults2 = new TH1D("fDiffFlowResults2","v'_2/2{2}",fgknBins,fPtMin,fPtMax);
+ fDiffFlowResults2->SetXTitle("pt [GeV]");
+ fDiffFlowResults2->SetYTitle("Differential Flow [%]");
+ fHistList->Add(fDiffFlowResults2);
+
+ //final results for differential flow v_p/n{4} (by default p=n=2) 
+ fDiffFlowResults4 = new TH1D("fDiffFlowResults4","v'_2/2{4}",fgknBins,fPtMin,fPtMax);
+ fDiffFlowResults4->SetXTitle("pt [GeV]");
+ fDiffFlowResults4->SetYTitle("Differential Flow [%]");
+ fHistList->Add(fDiffFlowResults4);
  
-  fIntFlowGenFun = new TProfile2D("fIntFlowGenFun","G[p][q]",8,0.,8.,17,0.,17.);//to be improved (z_down and z_up)
-  fIntFlowGenFun->SetXTitle("p");
-  fIntFlowGenFun->SetYTitle("q");
-  fHistList->Add(fIntFlowGenFun);
-     
-  fDiffFlowGenFunRe = new TProfile3D("fDiffFlowGenFunRe","Re(D[p][q])",fgknBins,fPtMin/fBinWidth,fPtMax/fBinWidth,8,0.,8.,17,0.,17.);
-  fHistList->Add(fDiffFlowGenFunRe);
+ //final results for differential flow v_p/n{6} (by default p=n=2)  
+ fDiffFlowResults6 = new TH1D("fDiffFlowResults6","v'_2/2{6}",fgknBins,fPtMin,fPtMax);
+ fDiffFlowResults6->SetXTitle("pt [GeV]");
+ fDiffFlowResults6->SetYTitle("Differential Flow [%]");
+ fHistList->Add(fDiffFlowResults6);
  
-  fDiffFlowGenFunIm = new TProfile3D("fDiffFlowGenFunIm","Im(D[p][q])",fgknBins,fPtMin/fBinWidth,fPtMax/fBinWidth,8,0.,8.,17,0.,17.);
-  fHistList->Add(fDiffFlowGenFunIm);
+ //final results for differential flow v_p/n{8} (by default p=n=2)
+ fDiffFlowResults8 = new TH1D("fDiffFlowResults8","v'_2/2{8}",fgknBins,fPtMin,fPtMax);
+ fDiffFlowResults8->SetXTitle("pt [GeV]");
+ fDiffFlowResults8->SetYTitle("Differential Flow [%]");
+ fHistList->Add(fDiffFlowResults8);
   
-
- /* 
- fCommonHistsRes2 = new AliFlowCommonHistResults("Cumulants2");//to be improved
- fCommonHistsRes4 = new AliFlowCommonHistResults("Cumulants4");//to be improved
- fCommonHistsRes6 = new AliFlowCommonHistResults("Cumulants6");//to be improved
- fCommonHistsRes8 = new AliFlowCommonHistResults("Cumulants8");//to be improved
+ //avarage of the generating function for integrated flow <G[p][q]>
+ fIntFlowGenFun = new TProfile2D("fIntFlowGenFun","<G[p][q]>",fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
+ fIntFlowGenFun->SetXTitle("p");
+ fIntFlowGenFun->SetYTitle("q");
+ fHistList->Add(fIntFlowGenFun);
  
+ /*
+ //avarage of the real part of generating function for differential flow <Re(D[b][p][q])>
+ fDiffFlowGenFunRe = new TProfile3D("fDiffFlowGenFunRe","<Re(D[b][p][q])>",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunRe->SetXTitle("b");
+ fDiffFlowGenFunRe->SetYTitle("p");
+ fDiffFlowGenFunRe->SetZTitle("q");
+ fHistList->Add(fDiffFlowGenFunRe);
  
- fHistList->Add(fCommonHistsRes2->GetHistList()); //NEW //to be improved
- fHistList->Add(fCommonHistsRes4->GetHistList()); //NEW //to be improved
- fHistList->Add(fCommonHistsRes6->GetHistList()); //NEW //to be improved 
- fHistList->Add(fCommonHistsRes8->GetHistList()); //NEW //to be improved
+ //avarage of the imaginary part of generating function for differential flow <Im(D[b][p][q])>
+ fDiffFlowGenFunIm = new TProfile3D("fDiffFlowGenFunIm","<Im(D[b][p][q])>",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunIm->SetXTitle("b");
+ fDiffFlowGenFunIm->SetYTitle("p");
+ fDiffFlowGenFunIm->SetZTitle("q");
+ fHistList->Add(fDiffFlowGenFunIm);
  */
-
- //control histograms
+ fDiffFlowGenFunRe0 = new TProfile2D("fDiffFlowGenFunRe0","Re(<D[b][0][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunRe0->SetXTitle("b");
+ fDiffFlowGenFunRe0->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunRe0);
+
+ fDiffFlowGenFunIm0 = new TProfile2D("fDiffFlowGenFunIm0","Im(<D[b][0][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunIm0->SetXTitle("b");
+ fDiffFlowGenFunIm0->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunIm0);
+ fDiffFlowGenFunRe1 = new TProfile2D("fDiffFlowGenFunRe1","Re(<D[b][1][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunRe1->SetXTitle("b");
+ fDiffFlowGenFunRe1->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunRe1);
+ fDiffFlowGenFunIm1 = new TProfile2D("fDiffFlowGenFunIm1","Im(<D[b][1][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunIm1->SetXTitle("b");
+ fDiffFlowGenFunIm1->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunIm1);
+ fDiffFlowGenFunRe2 = new TProfile2D("fDiffFlowGenFunRe2","Re(<D[b][2][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunRe2->SetXTitle("b");
+ fDiffFlowGenFunRe2->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunRe2);
+ fDiffFlowGenFunIm2 = new TProfile2D("fDiffFlowGenFunIm2","Im(<D[b][2][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunIm2->SetXTitle("b");
+ fDiffFlowGenFunIm2->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunIm2);
+ fDiffFlowGenFunRe3 = new TProfile2D("fDiffFlowGenFunRe3","Re(<D[b][3][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunRe3->SetXTitle("b");
+ fDiffFlowGenFunRe3->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunRe3);
+ fDiffFlowGenFunIm3 = new TProfile2D("fDiffFlowGenFunIm3","Im(<D[b][3][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunIm3->SetXTitle("b");
+ fDiffFlowGenFunIm3->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunIm3);
+ fDiffFlowGenFunRe4 = new TProfile2D("fDiffFlowGenFunRe4","Re(<D[b][4][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunRe4->SetXTitle("b");
+ fDiffFlowGenFunRe4->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunRe4);
+ fDiffFlowGenFunIm4 = new TProfile2D("fDiffFlowGenFunIm4","Im(<D[b][4][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunIm4->SetXTitle("b");
+ fDiffFlowGenFunIm4->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunIm4);
+ fDiffFlowGenFunRe5 = new TProfile2D("fDiffFlowGenFunRe5","Re(<D[b][5][q]>)",fgkQmax,0.,(Double_t)fgkQmax,fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth));
+ fDiffFlowGenFunRe5->SetXTitle("b");
+ fDiffFlowGenFunRe5->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunRe5);
+ fDiffFlowGenFunIm5 = new TProfile2D("fDiffFlowGenFunIm5","Im(<D[b][5][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunIm5->SetXTitle("b");
+ fDiffFlowGenFunIm5->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunIm5);
+ fDiffFlowGenFunRe6 = new TProfile2D("fDiffFlowGenFunRe6","Re(<D[b][6][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunRe6->SetXTitle("b");
+ fDiffFlowGenFunRe6->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunRe6);
+ fDiffFlowGenFunIm6 = new TProfile2D("fDiffFlowGenFunIm6","Im(<D[b][6][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunIm6->SetXTitle("b");
+ fDiffFlowGenFunIm6->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunIm6);
+ fDiffFlowGenFunRe7 = new TProfile2D("fDiffFlowGenFunRe7","Re(<D[b][7][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunRe7->SetXTitle("b");
+ fDiffFlowGenFunRe7->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunRe7);
+ fDiffFlowGenFunIm7 = new TProfile2D("fDiffFlowGenFunIm7","Im(<D[b][7][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
+ fDiffFlowGenFunIm7->SetXTitle("b");
+ fDiffFlowGenFunIm7->SetYTitle("q");
+ fHistList->Add(fDiffFlowGenFunIm7);
+ //common control histograms
  fCommonHists = new AliFlowCommonHist("Cumulants");
- fHistList->Add(fCommonHists->GetHistList()); 
+ fHistList->Add(fCommonHists->GetHistList());  
  
-}
-
-//________________________________________________________________________
-void AliFlowAnalysisWithCumulants::Make(AliFlowEventSimple* anEvent) {
-  //running over data
+ //q-distribution 
+ fQDist = new TH1D("fQDist","q-distribution",100,0,10);
+ fQDist->SetXTitle("q=Q/#sqrt{M}");
+ fQDist->SetYTitle("dN_{event}/dq");
+ fHistList->Add(fQDist);
   
-  //---------------------------------------------------------  
-  //fill the common control histograms
-  fCommonHists->FillControlHistograms(anEvent);   
-  //---------------------------------------------------------
+}//end of CreateOutputObjects()
+
+//================================================================================================================
+
+void AliFlowAnalysisWithCumulants::Make(AliFlowEventSimple* anEvent)
+{
+ //running over data
   
+ //fill the common control histograms
+ fCommonHists->FillControlHistograms(anEvent);   
   
-  //---------------------------------------------------------
-  //initializing the generating function for integrated flow
-  Double_t fG[fgkPmax][fgkQmax];
-  for(Int_t p=0;p<fgkPmax;p++){
-   for(Int_t q=0;q<fgkQmax;q++){
-    fG[p][q]=1.; 
-   }
+ //initializing the generating function G[p][q] for integrated flow 
+ Double_t G[fgkPmax][fgkQmax];
+ for(Int_t p=0;p<fgkPmax;p++){
+  for(Int_t q=0;q<fgkQmax;q++){
+   G[p][q]=1.; 
   }
-  //---------------------------------------------------------
-  
-  
-  //---------------------------------------------------------
-  //Q vector stuff 
-  AliFlowVector fQVector;
-  fQVector.Set(0.,0.);
-  fQVector.SetMult(0);
-  
-  fQVector=anEvent->GetQ();//get the Q vector for this event
+ }
   
-  fAvQx+=fQVector.X();
-  fAvQy+=fQVector.Y();
-  fAvQ2x+=pow(fQVector.X(),2.);
-  fAvQ2y+=pow(fQVector.Y(),2.);
-  //----------------------------------------------------------
-    
-  Int_t nPrim = anEvent->NumberOfTracks();
-  Int_t fEventNSelTracksIntFlow = anEvent->GetEventNSelTracksIntFlow();
-  Int_t fSelTracksIntFlow = 0;
-    
-  cout<<"Number of input tracks for cumulant analysis: "<<nPrim<<endl;
-  cout<<"Number of selected tracks for cumulant analysis: "<<fEventNSelTracksIntFlow<<endl;
+ Int_t nPrim = anEvent->NumberOfTracks();//total multiplicity
+ Int_t nEventNSelTracksIntFlow = anEvent->GetEventNSelTracksIntFlow();//selected multiplicity (parrticles used for int. flow)
+ Int_t nSelTracksIntFlow = 0;//cross-checking the selected multiplicity
   
-  //------------------------------------------------------------------------------------
-  //STARTING THE FIRST LOOP (CALCULATING THE GENERATING FUNCTION FOR INTEGRATED FLOW)
-  for(Int_t i=0;i<nPrim;i++){
-    fTrack=anEvent->GetTrack(i);
-    if(fTrack&&fTrack->UseForIntegratedFlow()){
-      fSelTracksIntFlow++;
-      for(Int_t p=0;p<fgkPmax;p++){
-       for(Int_t q=0;q<fgkQmax;q++){
-         fG[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/fEventNSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)); 
-       }
-      }
+ //first loop over data: evaluating the generating function G[p][q] for integrated flow 
+ for(Int_t i=0;i<nPrim;i++){
+  fTrack=anEvent->GetTrack(i);
+  if(fTrack&&fTrack->UseForIntegratedFlow()){
+   nSelTracksIntFlow++;
+   for(Int_t p=0;p<fgkPmax;p++){
+    for(Int_t q=0;q<fgkQmax;q++){
+     G[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/nEventNSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)); 
     }
-  }// ENDING THE FIRST LOOP OVER TRACKS
-  //------------------------------------------------------------------------------------ 
-  
-  
-  fHistProAvM->Fill(0.,fSelTracksIntFlow,1);
-  fAvM=fHistProAvM->GetBinContent(1);
-
-  
-  //------------------------------------------------------------------------------------
-  //STORING THE VALUE OF GENERATING FUNCTION FOR INTEGRATED FLOW INTO 2D PROFILE
-  for(Int_t p=0;p<fgkPmax;p++){
-   for(Int_t q=0;q<fgkQmax;q++){
-    fIntFlowGenFun->Fill((Double_t)p,(Double_t)q,fG[p][q],1);
-   } 
+   }
   }
-  //------------------------------------------------------------------------------------
-   
-                
-  //------------------------------------------------------------------------------------
-  //STARTING THE SECOND LOOP (CALCULATING THE GENERATING FUNCTION FOR DIFFERENTIAL FLOW)
-  // Remark 0: generating function for diff. flow is complex number, I need to calcuate separately real and imaginary part
-  // Remark 1: note that I need here fG[p][q], the value of generating function for integrated flow for the CURRENT event
-  // Remark 2: results are immediately stored in two 3D profiles, one for real and one for imaginary part
-  for(Int_t i=0;i<nPrim;i++){
-   fTrack=anEvent->GetTrack(i);
-   if (fTrack && fTrack->UseForDifferentialFlow()){
-    for(Int_t p=0;p<fgkPmax;p++){
-     for(Int_t q=0;q<fgkQmax;q++){
-      fDiffFlowGenFunRe->Fill(fTrack->Pt()/fBinWidth,(Double_t)p,(Double_t)q,fG[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/
-           (1.+(2.*fR0*sqrt(p+1.)/fSelTracksIntFlow) *
-            cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
-      fDiffFlowGenFunIm->Fill(fTrack->Pt()/fBinWidth,(Double_t)p,(Double_t)q,fG[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/
-           (1.+(2.*fR0*sqrt(p+1.)/fSelTracksIntFlow) * 
-            cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);        
-     }
+ }//ending the first loop over data
+ //storing the value of G[p][q] in 2D profile in order to get automatically the avarage <G[p][q]>
+ for(Int_t p=0;p<fgkPmax;p++){
+  for(Int_t q=0;q<fgkQmax;q++){
+   fIntFlowGenFun->Fill((Double_t)p,(Double_t)q,G[p][q],1);
+  } 
+ }
+ //storing the selected multiplicity (if fAvMultIntFlow is not filled here then you had wrongly selected the particles used for integrated flow)
+ if(nSelTracksIntFlow==nEventNSelTracksIntFlow)
+ {
+  fAvMultIntFlow->Fill(0.,nSelTracksIntFlow,1);
+ }
+ //calculating Q-vector of event (needed for errors)
+ AliFlowVector fQVector;
+ fQVector.Set(0.,0.);
+ fQVector.SetMult(0);
+ fQVector=anEvent->GetQ();                            //get the Q vector for this event
+ fQVectorComponents->Fill(0.,fQVector.X(),1);         //in the 1st bin fill Q_x
+ fQVectorComponents->Fill(1.,fQVector.Y(),1);         //in the 2nd bin fill Q_y
+ fQVectorComponents->Fill(2.,pow(fQVector.X(),2.),1); //in the 3rd bin fill (Q_x)^2
+ fQVectorComponents->Fill(3.,pow(fQVector.Y(),2.),1); //in the 4th bin fill (Q_y)^2
+ //q-distribution
+ Double_t qDist = fQVector.Mod()/sqrt(fQVector.GetMult());
+ fQDist->Fill(qDist,1); 
+
+ /*
+ //two 3D profiles for differential flow
+ //second loop over data: evaluating the generating function D[b][p][q] for differential flow 
+ //remark 0: D[b][p][q] is a complex number => real and imaginary part are calculated separately
+ //remark 1: note that below G[p][q] is needed, the value of generating function for integrated flow for the CURRENT event
+ //remark 2: results are stored in two 3D profiles in order to automatically get <Re(D[b][p][q])> and <Im(D[b][p][q])>
+ for(Int_t i=0;i<nPrim;i++){
+  fTrack=anEvent->GetTrack(i);
+  if (fTrack && fTrack->UseForDifferentialFlow()){
+   for(Int_t p=0;p<fgkPmax;p++){
+    for(Int_t q=0;q<fgkQmax;q++){
+     //real part
+     fDiffFlowGenFunRe->Fill(fTrack->Pt()/fBinWidth,(Double_t)p,(Double_t)q,G[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //imaginary part
+     fDiffFlowGenFunIm->Fill(fTrack->Pt()/fBinWidth,(Double_t)p,(Double_t)q,G[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);        
     }
-   }  
-  }// ENDING THE SECOND LOOP OVER TRACKS
-  //------------------------------------------------------------------------------------
-  
+   }
+  }  
+ }//ending the second loop over data    
+ */
+ //sixteen 2D profiles for differential flow          
+ for(Int_t i=0;i<nPrim;i++){
+  fTrack=anEvent->GetTrack(i);
+  if (fTrack && fTrack->UseForDifferentialFlow()){
+   //for(Int_t p=0;p<fgkPmax;p++){
+    for(Int_t q=0;q<fgkQmax;q++){
+     //real part
+     fDiffFlowGenFunRe0->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[0][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(0.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //imaginary part
+     fDiffFlowGenFunIm0->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[0][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(0.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); 
+     //-----------------------------------------------------------------------
+     //real part
+     fDiffFlowGenFunRe1->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[1][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(1.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //imaginary part
+     fDiffFlowGenFunIm1->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[1][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(1.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);             
+    //-----------------------------------------------------------------------
+     //real part
+     fDiffFlowGenFunRe2->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[2][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(2.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //imaginary part
+     fDiffFlowGenFunIm2->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[2][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(2.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);  
+     //-----------------------------------------------------------------------
+     //real part
+     fDiffFlowGenFunRe3->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[3][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(3.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //imaginary part
+     fDiffFlowGenFunIm3->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[3][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(3.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);  
+     //-----------------------------------------------------------------------
+     //real part
+     fDiffFlowGenFunRe4->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[4][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(4.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //imaginary part
+     fDiffFlowGenFunIm4->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[4][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(4.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);  
+     //-----------------------------------------------------------------------
+     //real part
+     fDiffFlowGenFunRe5->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[5][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(5.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //imaginary part
+     fDiffFlowGenFunIm5->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[5][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(5.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);  
+     //-----------------------------------------------------------------------
+     //real part
+     fDiffFlowGenFunRe6->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[6][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(6.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //imaginary part
+     fDiffFlowGenFunIm6->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[6][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(6.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //-----------------------------------------------------------------------
+     //real part
+     fDiffFlowGenFunRe7->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[7][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(7.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
+     //imaginary part
+     fDiffFlowGenFunIm7->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,G[7][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(7.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);    
+   }
+   //}
+  }  
+ }//ending the second loop over data            
+      
 }//end of Make()
 
-//________________________________________________________________________
-void AliFlowAnalysisWithCumulants::Finish(){
-  //final results
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
+//================================================================================================================
 
+void AliFlowAnalysisWithCumulants::Finish()
+{
+ //not needed for the time being...
+}
 
 
 
index ccce05e..e63bd00 100644 (file)
+/* 
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. 
+ * See cxx source for full Copyright notice 
+ * $Id$ 
+ */
+
+/******************************** 
+ * flow analysis with cumulants * 
+ *                              * 
+ * author: Ante Bilandzic       * 
+ *          (anteb@nikhef.nl)   *
+ *******************************/ 
+
 #ifndef AliFlowAnalysisWithCumulants_H
 #define AliFlowAnalysisWithCumulants_H
 
-//******************************* 
-// flow analysis with cumulants *   
-// author: Ante Bilandzic       * 
-// email: anteb@nikhef.nl       *
-//******************************* 
-
 #include "AliFlowCommonConstants.h"
 #include "AliFlowCumuConstants.h"
 
-class TH1;
-class TH3D;
-
 class TObjArray;
+class TList;
+class TFile;
 
-
+class TH1;
 class TProfile;
 class TProfile2D;
 class TProfile3D;
 
-
-class TList;
-class TFile;
-
 class AliFlowEventSimple;
 class AliFlowTrackSimple;
 class AliFlowCommonHist;
 class AliFlowCommonHistResults;
 class AliFlowVector;
 
-class AliFlowAnalysisWithCumulants {
+//================================================================================================================
+
+class AliFlowAnalysisWithCumulants{
  public:
   AliFlowAnalysisWithCumulants();
-  virtual ~AliFlowAnalysisWithCumulants();
+  virtual ~AliFlowAnalysisWithCumulants(); 
   
   virtual void CreateOutputObjects();
   virtual void Make(AliFlowEventSimple* anEvent);
   virtual void Finish();
-
-  // Output 
   
-  // Output
-  TList*   GetHistList() const             { return this->fHistList ; }     // Gets output histogram list //NEW
-
- private:
-  AliFlowAnalysisWithCumulants(const AliFlowAnalysisWithCumulants& aAnalysis);
-  AliFlowAnalysisWithCumulants& operator=(const AliFlowAnalysisWithCumulants& aAnalysis);
-
-  AliFlowTrackSimple* fTrack;//track
-  static const Int_t fgkQmax=AliFlowCumuConstants::kQmax;//needed for numerics
-  static const Int_t fgkPmax=AliFlowCumuConstants::kPmax;//needed for numerics  
-  static const Int_t fgkFlow=AliFlowCumuConstants::kFlow;//integrated flow coefficient to be calculated
-  static const Int_t fgkMltpl=AliFlowCumuConstants::kMltpl;//the multiple in p=m*n (diff. flow) 
-  static const Int_t fgknBins=100;//number of pt bins
+  TList* GetHistList() const {return this->fHistList;}      //output histogram list
  
-  //TString  fHistFileName;      //! The output file name  // NEW: In the old version this was not commented out     
-  //TFile* fHistFile; // histogram file for Cumulants // NEW: In the old version this was not commented out
-  
-  TList* fHistList; //list to hold all output histograms //NEW
-  
-  Double_t fAvM;//avarage SELECTED multiplicity
-
-  Double_t fR0;//needed for numerics
-  Double_t fPtMax;//maximum pt
-  Double_t fPtMin;//minimum pt
-  Double_t fBinWidth;//width of pt bin (in GeV)
+ private:
+  AliFlowAnalysisWithCumulants(const AliFlowAnalysisWithCumulants& afawc);
+  AliFlowAnalysisWithCumulants& operator=(const AliFlowAnalysisWithCumulants& afawc);
+  AliFlowTrackSimple* fTrack;                               //track
+  static const Int_t fgkQmax=AliFlowCumuConstants::kQmax;   //needed for numerics
+  static const Int_t fgkPmax=AliFlowCumuConstants::kPmax;   //needed for numerics  
+  static const Int_t fgkFlow=AliFlowCumuConstants::kFlow;   //integrated flow coefficient to be calculated
+  static const Int_t fgkMltpl=AliFlowCumuConstants::kMltpl; //the multiple in p=m*n (diff. flow) 
+  static const Int_t fgknBins=100;                          //number of pt bins
+  TList* fHistList;                                         //list to hold all output histograms
+
+  Double_t fR0;       //needed for numerics
+  Double_t fPtMax;    //maximum pt
+  Double_t fPtMin;    //minimum pt
+  Double_t fBinWidth; //width of pt bin (in GeV)
       
-  Double_t fAvQx;//<Q_x>
-  Double_t fAvQy;//<Q_y>
-  Double_t fAvQ2x;//<(Q_x)^2>
-  Double_t fAvQ2y;//<(Q_y)^2>
-   
-  Int_t fNumberOfEvents;//number of events 
-   
-  
-  TProfile*          fHistProAvM;       //avarage multiplicity
-  
+  Double_t fAvQx;  //<Q_x>
+  Double_t fAvQy;  //<Q_y>
+  Double_t fAvQ2x; //<(Q_x)^2>
+  Double_t fAvQ2y; //<(Q_y)^2>
+
+  TProfile*          fAvMultIntFlow;     //avarage selected multiplicity
+  TProfile*          fQVectorComponents; //avarages of Q-vector components (1st bin: <Q_x>, 2nd bin: <Q_y>, 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>)
+    
   TH1D*              fIntFlowResults;   //integrated flow final results
   
-  TH1D*              fDiffFlowResults2; //differential flow final results //to be improved
-  TH1D*              fDiffFlowResults4; //differential flow final results //to be improved
-  TH1D*              fDiffFlowResults6; //differential flow final results //to be improved
-  TH1D*              fDiffFlowResults8; //differential flow final results //to be improved
+  TH1D*              fDiffFlowResults2; //differential flow final results (2nd order estimate) 
+  TH1D*              fDiffFlowResults4; //differential flow final results (4th order estimate)
+  TH1D*              fDiffFlowResults6; //differential flow final results (6th order estimate)
+  TH1D*              fDiffFlowResults8; //differential flow final results (8th order estimate)
   
-  TProfile2D*        fIntFlowGenFun;    // avarage of the generating function for integrated flow 
-  TProfile3D*        fDiffFlowGenFunRe; // avarage of the generating function for differential flow (real part)
-  TProfile3D*        fDiffFlowGenFunIm; // avarage of the generating function for differential flow (imaginary part)
+  TProfile2D*        fIntFlowGenFun;    //avarage of the generating function for integrated flow 
+  TProfile3D*        fDiffFlowGenFunRe; //avarage of the generating function for differential flow (real part)
+  TProfile3D*        fDiffFlowGenFunIm; //avarage of the generating function for differential flow (imaginary part)
   
-  AliFlowCommonHist* fCommonHists;//control histograms
-  //AliFlowCommonHistResults *fCommonHistsRes2, *fCommonHistsRes4, *fCommonHistsRes6, *fCommonHistsRes8;//final results//to be improved 
+  TProfile2D *fDiffFlowGenFunRe0,*fDiffFlowGenFunRe1,*fDiffFlowGenFunRe2,*fDiffFlowGenFunRe3;//differential flow 
+  TProfile2D *fDiffFlowGenFunRe4,*fDiffFlowGenFunRe5,*fDiffFlowGenFunRe6,*fDiffFlowGenFunRe7;//differential flow
+  TProfile2D *fDiffFlowGenFunIm0,*fDiffFlowGenFunIm1,*fDiffFlowGenFunIm2,*fDiffFlowGenFunIm3;//differential flow
+  TProfile2D *fDiffFlowGenFunIm4,*fDiffFlowGenFunIm5,*fDiffFlowGenFunIm6,*fDiffFlowGenFunIm7;//differential flow
+  AliFlowCommonHist* fCommonHists;      //common control histograms
   
-  Double_t fAvG[fgkPmax][fgkQmax];//avarage of the generating function used for integrated flow
-  Int_t fBinEventEntries[fgknBins];//counts how many events have at least 1 particle in particular bin
-  Int_t fBinNoOfParticles[fgknBins];//number of particles per bin
-  Double_t fBinMeanPt[fgknBins];//mean pt per bin
-  Double_t fBinEventDRe[fgknBins][fgkPmax][fgkQmax];//real part of the generating function used for differential flow
-  Double_t fBinEventDIm[fgknBins][fgkPmax][fgkQmax];//imaginary part of the generating function used for differential flow
-    
+  TH1D*              fQDist;            //q-distribution
+      
   ClassDef(AliFlowAnalysisWithCumulants, 0);
 };
+
+//================================================================================================================
+
 #endif
 
 
index 669cf9a..0e5d0d0 100644 (file)
@@ -1,43 +1,11 @@
-//RUN SETTINGS
-//analysis type can be ESD, AOD, MC, ESDMC0, ESDMC1
-const TString type = "ESD";
-
-
-//SETTING THE CUTS
-
-//for integrated flow
-const Double_t ptmin1 = 0.0;
-const Double_t ptmax1 = 1000.0;
-const Double_t ymin1  = -2.;
-const Double_t ymax1  = 2.;
-const Int_t mintrackrefsTPC1 = 2;
-const Int_t mintrackrefsITS1 = 3;
-const Int_t charge1 = 1;
-const Int_t PDG1 = 211;
-const Int_t minclustersTPC1 = 50;
-const Int_t maxnsigmatovertex1 = 3;
-
-//for differential flow
-const Double_t ptmin2 = 0.0;
-const Double_t ptmax2 = 1000.0;
-const Double_t ymin2  = -2.;
-const Double_t ymax2  = 2.;
-const Int_t mintrackrefsTPC2 = 2;
-const Int_t mintrackrefsITS2 = 3;
-const Int_t charge2 = 1;
-const Int_t PDG2 = 321;
-const Int_t minclustersTPC2 = 50;
-const Int_t maxnsigmatovertex2 = 3;
-
-
 // from CreateESDChain.C - instead of  #include "CreateESDChain.C"
 TChain* CreateESDChain(const char* aDataDir = "ESDfiles.txt", Int_t aRuns = 20, Int_t offset = 0) ;
 void LookupWrite(TChain* chain, const char* target) ;
 
-//void runAliAnalysisTaskCumulants(Int_t nRuns = -1, TString type = "ESD", const Char_t* dataDir="/data/alice2/ante/ab", Int_t offset = 0) 
+void runAliAnalysisTaskCumulants(Int_t nRuns =-1, TString type = "ESD", const Char_t* dataDir="/data/alice2/ante/ab2", Int_t offset = 0) 
 //void runAliAnalysisTaskCumulants(Int_t nRuns = -1, TString type = "ESD", const Char_t* dataDir="/data/alice2/ab2", Int_t offset = 0) 
 //void runAliAnalysisTaskCumulants(Int_t nRuns = -1, TString type = "ESD", const Char_t* dataDir="/data/alice2/LHyquid3_rot", Int_t offset = 0)
-void runAliAnalysisTaskCumulants(Int_t nRuns = -1, TString type = "ESD", const Char_t* dataDir="/Users/snelling/alice_data/TherminatorFix", Int_t offset = 0) 
+//void runAliAnalysisTaskCumulants(Int_t nRuns = 2, TString type = "MC", const Char_t* dataDir="/Users/snelling/alice_data/TherminatorFix", Int_t offset = 0) 
 {
   TStopwatch timer;
   timer.Start();
@@ -52,139 +20,8 @@ void runAliAnalysisTaskCumulants(Int_t nRuns = -1, TString type = "ESD", const C
   cerr<<"libESD loaded..."<<endl;
   gSystem->Load("libANALYSIS.so");
   cerr<<"libANALYSIS.so loaded..."<<endl;
-  gSystem->Load("libCORRFW.so");
-  cerr<<"libCORRFW.so loaded..."<<endl;
   gSystem->Load("libPWG2flow.so");
   cerr<<"libPWG2flow.so loaded..."<<endl;
-
-
-
-
-//____________________________________________//
- //Create cuts using correction framework
-
- //############# cuts on MC
- AliCFTrackKineCuts* mcKineCuts1 = new AliCFTrackKineCuts("mcKineCuts1","MC-level kinematic cuts");
- mcKineCuts1->SetPtRange(ptmin1,ptmax1);
- mcKineCuts1->SetRapidityRange(ymin1,ymax1);
- mcKineCuts1->SetChargeMC(charge1);
-
- AliCFTrackKineCuts* mcKineCuts2 = new AliCFTrackKineCuts("mcKineCuts2","MC-level kinematic cuts");
- mcKineCuts2->SetPtRange(ptmin2,ptmax2);
- mcKineCuts2->SetRapidityRange(ymin2,ymax2);
- mcKineCuts2->SetChargeMC(charge2);
-
- AliCFParticleGenCuts* mcGenCuts1 = new AliCFParticleGenCuts("mcGenCuts1","MC particle generation cuts");
- mcGenCuts1->SetRequireIsPrimary();
- mcGenCuts1->SetRequirePdgCode(PDG1);
-
- AliCFParticleGenCuts* mcGenCuts2 = new AliCFParticleGenCuts("mcGenCuts2","MC particle generation cuts");
- mcGenCuts2->SetRequireIsPrimary();
- mcGenCuts2->SetRequirePdgCode(PDG2);
-
- //############# Acceptance Cuts  ????????
- AliCFAcceptanceCuts *mcAccCuts = new AliCFAcceptanceCuts("mcAccCuts","MC acceptance cuts");
- mcAccCuts->SetMinNHitITS(mintrackrefsITS1);
- mcAccCuts->SetMinNHitTPC(mintrackrefsTPC1);
-
- //############# Rec-Level kinematic cuts
- AliCFTrackKineCuts *recKineCuts1 = new AliCFTrackKineCuts("recKineCuts1","rec-level kine cuts");
- recKineCuts1->SetPtRange(ptmin1,ptmax1);
- recKineCuts1->SetRapidityRange(ymin1,ymax1);
- recKineCuts1->SetChargeRec(charge1);
-
- AliCFTrackKineCuts *recKineCuts2 = new AliCFTrackKineCuts("recKineCuts2","rec-level kine cuts");
- recKineCuts2->SetPtRange(ptmin2,ptmax2);
- recKineCuts2->SetRapidityRange(ymin2,ymax2);
- recKineCuts2->SetChargeRec(charge2);
-
- AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
- recQualityCuts->SetMinNClusterTPC(minclustersTPC1);
- recQualityCuts->SetRequireITSRefit(kTRUE);
-
- AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
- recIsPrimaryCuts->SetMaxNSigmaToVertex(maxnsigmatovertex1);
-
- AliCFTrackCutPid* cutPID1 = new AliCFTrackCutPid("cutPID1","ESD_PID") ;
- AliCFTrackCutPid* cutPID2 = new AliCFTrackCutPid("cutPID2","ESD_PID") ;
- int n_species = AliPID::kSPECIES ;
- Double_t* prior = new Double_t[n_species];
-
- prior[0] = 0.0244519 ;
- prior[1] = 0.0143988 ;
- prior[2] = 0.805747  ;
- prior[3] = 0.0928785 ;
- prior[4] = 0.0625243 ;
-
- cutPID1->SetPriors(prior);
- cutPID1->SetProbabilityCut(0.0);
- cutPID1->SetDetectors("TPC TOF");
- switch(TMath::Abs(PDG1)) {
- case 11   : cutPID1->SetParticleType(AliPID::kElectron, kTRUE); break;
- case 13   : cutPID1->SetParticleType(AliPID::kMuon    , kTRUE); break;
- case 211  : cutPID1->SetParticleType(AliPID::kPion    , kTRUE); break;
- case 321  : cutPID1->SetParticleType(AliPID::kKaon    , kTRUE); break;
- case 2212 : cutPID1->SetParticleType(AliPID::kProton  , kTRUE); break;
- default   : printf("UNDEFINED PID\n"); break;
- }
-
- cutPID2->SetPriors(prior);
- cutPID2->SetProbabilityCut(0.0);
- cutPID2->SetDetectors("TPC TOF");
- switch(TMath::Abs(PDG2)) {
- case 11   : cutPID2->SetParticleType(AliPID::kElectron, kTRUE); break;
- case 13   : cutPID2->SetParticleType(AliPID::kMuon    , kTRUE); break;
- case 211  : cutPID2->SetParticleType(AliPID::kPion    , kTRUE); break;
- case 321  : cutPID2->SetParticleType(AliPID::kKaon    , kTRUE); break;
- case 2212 : cutPID2->SetParticleType(AliPID::kProton  , kTRUE); break;
- default   : printf("UNDEFINED PID\n"); break;
- }
-
- printf("CREATE MC KINE CUTS\n");
- TObjArray* mcList1 = new TObjArray(0);
- mcList1->AddLast(mcKineCuts1);
- mcList1->AddLast(mcGenCuts1);
-
- TObjArray* mcList2 = new TObjArray(0);
- mcList2->AddLast(mcKineCuts2);
- mcList2->AddLast(mcGenCuts2);
-
- printf("CREATE ACCEPTANCE CUTS\n");
- TObjArray* accList = new TObjArray(0) ;
- accList->AddLast(mcAccCuts);
-
- printf("CREATE RECONSTRUCTION CUTS\n");
- TObjArray* recList1 = new TObjArray(0) ;
- recList1->AddLast(recKineCuts1);
- recList1->AddLast(recQualityCuts);
- recList1->AddLast(recIsPrimaryCuts);
-
- TObjArray* recList2 = new TObjArray(0) ;
- recList2->AddLast(recKineCuts2);
- recList2->AddLast(recQualityCuts);
- recList2->AddLast(recIsPrimaryCuts);
-
- printf("CREATE PID CUTS\n");
- TObjArray* fPIDCutList1 = new TObjArray(0) ;
- fPIDCutList1->AddLast(cutPID1);
-
- TObjArray* fPIDCutList2 = new TObjArray(0) ;
- fPIDCutList2->AddLast(cutPID2);
-
- printf("CREATE INTERFACE AND CUTS\n");
- AliCFManager* cfmgr1 = new AliCFManager();
- cfmgr1->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList1);
- //cfmgr1->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
- cfmgr1->SetParticleCutsList(AliCFManager::kPartRecCuts,recList1);
- cfmgr1->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList1);
-
- AliCFManager* cfmgr2 = new AliCFManager();
- cfmgr2->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList2);
- //cfmgr2->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
- cfmgr2->SetParticleCutsList(AliCFManager::kPartRecCuts,recList2);
- cfmgr2->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList2);
-
-
   
   // create the TChain. CreateESDChain() is defined in CreateESDChain.C
   TChain* chain = CreateESDChain(dataDir, nRuns, offset);
@@ -214,9 +51,7 @@ void runAliAnalysisTaskCumulants(Int_t nRuns = -1, TString type = "ESD", const C
   //____________________________________________//
   // 1st cumulant task
   AliAnalysisTaskCumulants *task1 = new AliAnalysisTaskCumulants("TaskCumulants");
-  task1->SetAnalysisType("ESD");
-  task1->SetCFManager1(cfmgr1);
-  task1->SetCFManager2(cfmgr2);
+  task1->SetAnalysisType(type);
   mgr->AddTask(task1);
 
   // Create containers for input/output
@@ -348,3 +183,4 @@ void LookupWrite(TChain* chain, const char* target)
   
   delete iter;
 }
+
index c8b01da..f0df43a 100644 (file)
@@ -30,31 +30,47 @@ const Int_t minclustersTPC2 = 50;
 const Int_t maxnsigmatovertex2 = 3;
 
 
-  void runProofCumulants(const Char_t* data="/PWG2/akisiel/LHC500C2030", Int_t nRuns=100, Int_t offset=0) {
+  void runProofCumulants(const Char_t* data="/PWG2/akisiel/LHC500C2030", Int_t nRuns=4, Int_t offset=0) {
+  //void runProofCumulants(const Char_t* data="/PWG2/akisiel/LHC500C0005", Int_t nRuns=-1, Int_t offset=0){ 
+
+  //void runProofCumulants(const Char_t* data="/PWG2/pganoti/Pythia6At10TeV_05T_b", Int_t nRuns=1000, Int_t offset=0){ 
+  //void runProofCumulants(const Char_t* data="/PWG2/jgrosseo/sim_1600XX_esd", Int_t nRuns=-1, Int_t offset=0){ 
+  
+  //void runProofCumulants(const Char_t* data="/COMMON/COMMON/run15035_PbPb", Int_t nRuns=-1, Int_t offset=0){
+
+  //void runProofCumulants(const Char_t* data="/PWG2/hricaud/LHC07f_160033DataSet", Int_t nRuns=-1, Int_t offset=0){
+
+  //void runProofCumulants(const Char_t* data="/PWG2/hricaud/LHC07f_160038_root_archiveDataSet", Int_t nRuns=-1, Int_t offset=0){
+
+ //void runProofCumulants(const Char_t* data="/PWG2/belikov/40825", Int_t nRuns=-1, Int_t offset=0){
+
+
   TStopwatch timer;
   timer.Start();
 
   printf("*** Connect to PROOF ***\n");
-  //  TProof::Open("snelling@lxb6046.cern.ch");
-  TProof::Open("snelling@localhost");
-
-  gProof->UploadPackage("STEERBase.par");
-  gProof->EnablePackage("STEERBase");
-  gProof->UploadPackage("ESD.par");
-  gProof->EnablePackage("ESD");
-  gProof->UploadPackage("AOD.par");
-  gProof->EnablePackage("AOD");
-  gProof->UploadPackage("ANALYSIS.par");
-  gProof->EnablePackage("ANALYSIS");
-  gProof->UploadPackage("ANALYSISalice.par");
-  gProof->EnablePackage("ANALYSISalice");
-  gProof->UploadPackage("PWG2AOD.par");
-  gProof->EnablePackage("PWG2AOD");
-  gProof->UploadPackage("CORRFW.par");
-  gProof->EnablePackage("CORRFW");
-  //  gProof->ClearPackage("PWG2flow");
-  gProof->UploadPackage("PWG2flow.par");
-  gProof->EnablePackage("PWG2flow");
+  TProof::Open("abilandz@lxb6046.cern.ch");
+  //TProof::Open("snelling@localhost");
+
+ gProof->UploadPackage("STEERBase.par");
+ gProof->EnablePackage("STEERBase");
+ gProof->UploadPackage("ESD.par");
+ gProof->EnablePackage("ESD");
+ gProof->UploadPackage("AOD.par");
+ gProof->EnablePackage("AOD");
+ gProof->UploadPackage("ANALYSIS.par");
+ gProof->EnablePackage("ANALYSIS");
+ gProof->UploadPackage("ANALYSISalice.par");
+ gProof->EnablePackage("ANALYSISalice");
+ gProof->UploadPackage("PWG2AOD.par");
+ gProof->EnablePackage("PWG2AOD");
+ gProof->UploadPackage("CORRFW.par");
+ gProof->EnablePackage("CORRFW");
+ gProof->ClearPackage("PWG2flow");
+ gProof->UploadPackage("PWG2flow.par");
+ gProof->EnablePackage("PWG2flow");
+
 
 
 //____________________________________________//
@@ -106,7 +122,7 @@ const Int_t maxnsigmatovertex2 = 3;
  AliCFTrackCutPid* cutPID2 = new AliCFTrackCutPid("cutPID2","ESD_PID") ;
  int n_species = AliPID::kSPECIES ;
  Double_t* prior = new Double_t[n_species];
-
  prior[0] = 0.0244519 ;
  prior[1] = 0.0143988 ;
  prior[2] = 0.805747  ;
@@ -168,7 +184,7 @@ const Int_t maxnsigmatovertex2 = 3;
  TObjArray* fPIDCutList2 = new TObjArray(0) ;
  fPIDCutList2->AddLast(cutPID2);
 
- printf("CREATE INTERFACE AND CUTS\n");
+printf("CREATE INTERFACE AND CUTS\n");
  AliCFManager* cfmgr1 = new AliCFManager();
  cfmgr1->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList1);
  //cfmgr1->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
@@ -180,27 +196,57 @@ const Int_t maxnsigmatovertex2 = 3;
  //cfmgr2->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
  cfmgr2->SetParticleCutsList(AliCFManager::kPartRecCuts,recList2);
  cfmgr2->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList2);
+ //____________________________________________//
+ // Make the analysis manager
+ AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+
+ if (type == "ESD"){
+   AliVEventHandler* esdH = new AliESDInputHandler;
+   //esdH->SetInactiveBranches("FMD CaloCluster");   //needed?
+   mgr->SetInputEventHandler(esdH); }
+
+ if (type == "AOD"){
+   AliVEventHandler* aodH = new AliAODInputHandler;
+   mgr->SetInputEventHandler(aodH); }
+
+ if (type == "MC" || type == "ESDMC0" || type == "ESDMC1"){
+   AliVEventHandler* esdH = new AliESDInputHandler;
+   mgr->SetInputEventHandler(esdH);
+
+   AliMCEventHandler *mc = new AliMCEventHandler();
+   mgr->SetMCtruthEventHandler(mc); }
 
-
-  //____________________________________________//
-  // Make the analysis manager
-  AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
-  AliESDInputHandler* esdH = new AliESDInputHandler;
-  esdH->SetInactiveBranches("FMD CaloCluster");
-  mgr->SetInputEventHandler(esdH);  
+  
+  
+  
+  
+  
+  
+  
+  
   //____________________________________________//
   // 1st Pt task
   AliAnalysisTaskCumulants *task1 = new AliAnalysisTaskCumulants("TaskCumulants");
-  task1->SetAnalysisType("ESD");
+  task1->SetAnalysisType(type);
   task1->SetCFManager1(cfmgr1);
   task1->SetCFManager2(cfmgr2);
   mgr->AddTask(task1);
 
-
   // Create containers for input/output
-  AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
-  //  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("clist1", TList::Class(),AliAnalysisManager::kOutputContainer,"OutputFromCumulantAnlysisESD.root");
-  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("clist1", TList::Class(),AliAnalysisManager::kOutputContainer,"outputFromCumulantAnalysisESD.root");
+ AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
+ TString outputName = "outputFromCumulantAnalysis";
+ outputName+=type;
+ outputName+=".root";
+ AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("clist1", TList::Class(),AliAnalysisManager::kOutputContainer,outputName);
+  
+ //added according to Andrei:
+ //coutput1->SetSpecialOutput();
+  
   
   //____________________________________________//
   mgr->ConnectInput(task1,0,cinput1);