]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
added support for correction framework in cumulants
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Jul 2008 15:22:29 +0000 (15:22 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Jul 2008 15:22:29 +0000 (15:22 +0000)
PWG2/FLOW/AliAnalysisTaskCumulants.cxx
PWG2/FLOW/AliAnalysisTaskCumulants.h
PWG2/FLOW/AliCumulantsFunctions.cxx [new file with mode: 0644]
PWG2/FLOW/AliCumulantsFunctions.h [new file with mode: 0644]
PWG2/FLOW/AliFlowAnalysisWithCumulants.cxx
PWG2/FLOW/AliFlowAnalysisWithCumulants.h
PWG2/FLOW/macros/runAliAnalysisTaskCumulants.C
PWG2/FLOW/macros/runProofCumulants.C [new file with mode: 0644]
PWG2/PWG2flowLinkDef.h
PWG2/libPWG2flow.pkg

index e6c92ea30fb9dc9403718fadfe8e24351a9b779f..8ac3914875168dd2d4eff80caa3d9121b38df073 100644 (file)
 #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"
 
-class AliAnalysisTask;
+#include "AliAnalysisTask.h"
+#include "AliAnalysisDataSlot.h"
+#include "AliAnalysisDataContainer.h"
+
+//class AliAnalysisTask;
 #include "AliAnalysisManager.h"
 
 #include "AliESDEvent.h"
@@ -31,13 +42,22 @@ class AliAnalysisTask;
 #include "AliMCEventHandler.h"
 #include "AliMCEvent.h"
 
+#include "../../CORRFW/AliCFManager.h"
+
 #include "AliAnalysisTaskCumulants.h"
 #include "AliFlowEventSimpleMaker.h"
 #include "AliFlowAnalysisWithCumulants.h"
+#include "AliFlowCumuConstants.h"
+#include "AliFlowCommonConstants.h"
+#include "AliFlowCommonHistResults.h"
+
+#include "AliCumulantsFunctions.h"
 
 // AliAnalysisTaskCumulants:
-// analysis task for Lee Yang Zeros method
-// Author: Naomi van der Kolk (kolk@nikhef.nl)
+// analysis task for 
+// Cumulant method
+// with many authors (N.K. R.S. A.B.)
+// who do something
 
 ClassImp(AliAnalysisTaskCumulants)
 
@@ -46,9 +66,12 @@ AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name) :
   AliAnalysisTask(name, ""), 
   fESD(NULL),
   fAOD(NULL),
-  fAnalysisType("ESD"), 
   fMyCumuAnalysis(NULL),
-  fEventMaker(NULL)
+  fEventMaker(NULL),
+  fAnalysisType("ESD"), 
+  fCFManager1(NULL),
+  fCFManager2(NULL),
+  fListHistos(NULL)
 {
   // Constructor
   cout<<"AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name)"<<endl;
@@ -60,6 +83,20 @@ AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name) :
   // 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)
+{
+  // Constructor
+  cout<<"AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name)"<<endl;
+  }
 
 //________________________________________________________________________
 void AliAnalysisTaskCumulants::ConnectInputData(Option_t *) 
@@ -117,30 +154,28 @@ void AliAnalysisTaskCumulants::ConnectInputData(Option_t *)
 //________________________________________________________________________
 void AliAnalysisTaskCumulants::CreateOutputObjects() 
 {
-  // Called once
+  // 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
+  //analyser
   fMyCumuAnalysis = new AliFlowAnalysisWithCumulants() ;
    
-
-  //output file
-  TString outputName = "outputFromCumulantsAnalysis" ;
-  outputName += fAnalysisType.Data();
-  outputName += ".root";
-  fMyCumuAnalysis->SetHistFileName( outputName.Data() ); //Ante please implement
-  
-  
   fMyCumuAnalysis->CreateOutputObjects();
 
-}
+  if (fMyCumuAnalysis->GetHistList()) {
+    // fSP->GetHistList()->Print();
+    fListHistos = fMyCumuAnalysis->GetHistList();
+    // fListHistos->Print();
+  }
+  else {Printf("ERROR: Could not retrieve histogram list"); }
+ }
 
 //________________________________________________________________________
 void AliAnalysisTaskCumulants::Exec(Option_t *) 
@@ -164,9 +199,11 @@ void AliAnalysisTaskCumulants::Exec(Option_t *)
     }
 
     Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
+    fCFManager1->SetEventInfo(mcEvent);
+    fCFManager2->SetEventInfo(mcEvent);
 
-    //Cumulants analysis 
-    AliFlowEventSimple* fEvent = fEventMaker->FillTracks(mcEvent);
+    //cumulant analysis 
+    AliFlowEventSimple* fEvent = fEventMaker->FillTracks(mcEvent,fCFManager1,fCFManager2);
     fMyCumuAnalysis->Make(fEvent);
     delete fEvent;
   }
@@ -201,37 +238,21 @@ void AliAnalysisTaskCumulants::Exec(Option_t *)
       return;
     }
 
-    //Cumulant analysis 
-    AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD,mcEvent,0); //0 = kine from ESD, 1 = kine from MC
-    fMyCumuAnalysis->Make(fEvent);
-    delete fEvent;
-    //delete mcEvent;
-  }
-  else if (fAnalysisType == "ESDMC1") {
-    if (!fESD) {
-      Printf("ERROR: fESD not available");
-      return;
-    }
-    Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
-    
-    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-    if (!eventHandler) {
-      Printf("ERROR: Could not retrieve MC event handler");
-      return;
-    }
-
-    AliMCEvent* mcEvent = eventHandler->MCEvent();
-    if (!mcEvent) {
-      Printf("ERROR: Could not retrieve MC event");
-      return;
-    }
+    fCFManager1->SetEventInfo(mcEvent);
+    fCFManager2->SetEventInfo(mcEvent);
 
     //Cumulant analysis 
-    AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD,mcEvent,1); //0 = kine from ESD, 1 = kine from MC
+    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);
     delete fEvent;
     //delete mcEvent;
   }
+  
   else if (fAnalysisType == "AOD") {
     if (!fAOD) {
       Printf("ERROR: fAOD not available");
@@ -239,27 +260,71 @@ void AliAnalysisTaskCumulants::Exec(Option_t *)
     }
     Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
 
-    //Cumulant analysis 
+    // analysis 
+    //For the moment don't use CF //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD,fCFManager1,fCFManager2);
     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD);
     fMyCumuAnalysis->Make(fEvent);
     delete fEvent;
   }
-  
+
+
+  PostData(0,fListHistos); 
 }
 
 //________________________________________________________________________
 void AliAnalysisTaskCumulants::Terminate(Option_t *) 
-{
-  // Called once at the end of the query
-  cerr<<"fMyCumuAnalysis->GetHistFile() -> IsOpen() = "<<fMyCumuAnalysis->GetHistFile() -> IsOpen()<<endl;
+{  
+  //=====================================================================================================
+  // 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"));
+      
+      // 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;
+    }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
-  fMyCumuAnalysis->Finish(); 
 
-  PostData(0,fMyCumuAnalysis->GetHistFile());
 
-  delete fMyCumuAnalysis;
-  delete fEventMaker;
 
-  cout<<".....finished"<<endl;
-}
 
index 0cf598ed7da4636d7637b34653a072ec315a9096..2f51e036f6e16e5293caf714772d3c8ff11fe4c9 100644 (file)
@@ -9,11 +9,12 @@
 // AliAnalysisTaskCumulants:
 // analysis task for 
 // Cumulant method
-// with many authors
+// with many authors (N.K. R.S. A.B.)
 // who do something
 
 class AliESDEvent;
 class AliAODEvent;
+class AliCFManager;
 class AliFlowAnalysisWithCumulants;
 class AliFlowEventSimpleMaker;
 class TFile;
@@ -22,7 +23,9 @@ class TFile;
 
 class AliAnalysisTaskCumulants : public AliAnalysisTask {
  public:
-  AliAnalysisTaskCumulants(const char *name = "AliAnalysisTaskCumulants");
+  //AliAnalysisTaskCumulants(const char *name = "AliAnalysisTaskCumulants");
+  AliAnalysisTaskCumulants();
+  AliAnalysisTaskCumulants(const char *name);
   virtual ~AliAnalysisTaskCumulants() {}
   
   virtual void   ConnectInputData(Option_t *);
@@ -31,7 +34,12 @@ class AliAnalysisTaskCumulants : public AliAnalysisTask {
   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; } 
  
  private:
  
@@ -40,11 +48,25 @@ class AliAnalysisTaskCumulants : public AliAnalysisTask {
 
   AliESDEvent *fESD;                      //ESD object
   AliAODEvent* fAOD;                      //AOD object
-  TString fAnalysisType;                  //string to select which kind of input to analyse: ESD, AOD or MC
   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
 };
 
 #endif
+
+
+
+
+
+
+
+
+
+
+
diff --git a/PWG2/FLOW/AliCumulantsFunctions.cxx b/PWG2/FLOW/AliCumulantsFunctions.cxx
new file mode 100644 (file)
index 0000000..bf9f450
--- /dev/null
@@ -0,0 +1,364 @@
+//******************************* 
+// flow analysis with cumulants *   
+// author: Ante Bilandzic       * 
+// email: anteb@nikhef.nl       *
+//******************************* 
+
+#define AliCumulantsFunctions_cxx
+
+#include "Riostream.h"
+#include "AliFlowCommonConstants.h"
+#include "TChain.h"
+#include "TFile.h"
+#include "TList.h" //NEW
+#include "TParticle.h"
+
+#include "TProfile.h"
+#include "TProfile2D.h" 
+#include "TProfile3D.h"
+
+#include "AliFlowEventSimple.h"
+#include "AliFlowTrackSimple.h"
+#include "AliFlowAnalysisWithCumulants.h"
+#include "AliFlowCumuConstants.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 
+ }
+
+AliCumulantsFunctions::~AliCumulantsFunctions(){
+  //desctructor
+}
+
+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 
+ };
+  
+//___________________________________________________________________________
+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)   
+     
+  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
+    }
+  }
+  /////////////////////////////////////////////////////////////////////////////
+  ///////avaraging the gen. function for the cumulants over azimuth////////////
+  /////////////////////////////////////////////////////////////////////////////
+  
+  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;
+  }
+  
+  /////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////////final results//////////////////////////////
+  /////////////////////////////////////////////////////////////////////////////
+  
+  Double_t fCumulant[fgkPmax];//c_{iFlow}\{2(p+1)\}
+  
+  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]);
+
+  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]);
+
+  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]);
+
+  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]);
+
+  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]);
+
+  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]);
+
+  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]);
+
+  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]);
+  
+  
+  cout<<""<<endl;
+  cout<<"***************************"<<endl;
+  cout<<"cumulants:"<<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;
+  
+  
+  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;  
+  } 
+  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;  
+  }
+  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.)));
+  } else {
+      cout<<"v_"<<fgkFlow<<"{10} = Im"<<endl;  //to be improved (2->fgkFlow)
+  }
+  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.)));
+  } else {
+    cout<<"v_"<<fgkFlow<<"{12} = Im"<<endl;  //to be improved (2->fgkFlow)
+  }
+  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.)));
+  } else {
+    cout<<"v_"<<fgkFlow<<"{14} = Im"<<endl;  //to be improved (2->fgkFlow)
+  }
+  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.)));
+  } else {
+    cout<<"v_"<<fgkFlow<<"{16} = Im"<<endl;  //to be improved (2->fgkFlow)
+  }
+  //cout<<"***************************"<<endl;
+      
+  //DIFFERENTIAL FLOW CALCULATIONS STARTS HERE!!!
+  
+  Double_t fX[fgknBins][fgkPmax][fgkQmax]={0.};//see the text bellow relation (11) in PG
+  Double_t fY[fgknBins][fgkPmax][fgkQmax]={0.};
+  
+  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];
+      }
+    }   
+  } 
+  
+  Double_t fD[fgknBins][fgkPmax]={0.};//implementing relation (11) from PG
+  
+  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];
+      } 
+      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;
+      }  
+    }
+  } 
+  
+  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.};
+  
+  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]);
+  }
+  
+  Double_t fv2[fgknBins],fv4[fgknBins],fv6[fgknBins],fv8[fgknBins];
+  //Double_t fAvPt[fgknBins];
+  Double_t fSddiff2[fgknBins],fSddiff4[fgknBins];
+
+  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<<"v'_2/2{2} = Im"<<endl;
+      //}
+    }else{
+      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;
+      //} 
+    }else{
+      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]);
+    }else{
+      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]);
+    }else{
+      cout<<"v'_2/2{8} = Im"<<endl;
+    }       
+    cout<<"****************************************"<<endl;
+  }  
+}
+
+  
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/PWG2/FLOW/AliCumulantsFunctions.h b/PWG2/FLOW/AliCumulantsFunctions.h
new file mode 100644 (file)
index 0000000..5ebd4a4
--- /dev/null
@@ -0,0 +1,52 @@
+#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"
+
+class TH1;
+class TProfile;
+class TProfile2D;
+class TProfile3D;
+
+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);
+  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
+  
+  ClassDef(AliCumulantsFunctions, 0);
+};
+#endif
+
+
+
+
+
index dd4ecc795a987b74b9cfc47dda767e8dc9f45374..a70f428f4829749a1f6c3dd57f8f73053c2c5fa3 100644 (file)
@@ -6,14 +6,31 @@
 #include "AliFlowCommonHistResults.h"
 #include "TChain.h"
 #include "TFile.h"
+#include "TList.h" //NEW
 #include "TParticle.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;
@@ -40,8 +57,9 @@ ClassImp(AliFlowAnalysisWithCumulants)
 
 AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants():  
   fTrack(NULL),
-  fHistFileName("cumulants.root"),
-  fHistFile(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),
@@ -52,61 +70,133 @@ AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants():
   fAvQ2x(0),
   fAvQ2y(0),
   fNumberOfEvents(0),
-  fCommonHists(NULL),
-  fCommonHistsRes2(NULL),
-  fCommonHistsRes4(NULL),
-  fCommonHistsRes6(NULL),
-  fCommonHistsRes8(NULL)
+  
+  
+  
+  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;
    
-   for(Int_t n=0;n<fgknBins;n++){
-     fBinEventEntries[n]=0;
-     fBinNoOfParticles[n]=0;
-     fBinMeanPt[n]=0;
-     for(Int_t p=0;p<fgkPmax;p++){
-       for(Int_t q=0;q<fgkQmax;q++){
-        fAvG[p][q]=0;
-        fBinEventDRe[n][p][q]=0; 
-        fBinEventDIm[n][p][q]=0;
-       }
-     }
-   }
  }
 
 AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants(){
   //desctructor
+   delete fHistList; //NEW
 }
 
   
 //___________________________________________________________________________
 void AliFlowAnalysisWithCumulants::CreateOutputObjects(){
  //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);
+  //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);
+  //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);
+  
+  //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);
+  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);
+  fDiffFlowGenFunIm = new TProfile3D("fDiffFlowGenFunIm","Im(D[p][q])",fgknBins,fPtMin/fBinWidth,fPtMax/fBinWidth,8,0.,8.,17,0.,17.);
+  fHistList->Add(fDiffFlowGenFunIm);
+  
+
+ /* 
+ 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
+ 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
+ */
 
- fHistFile = new TFile(fHistFileName.Data(),"RECREATE");
- fCommonHists = new AliFlowCommonHist("Cumulants");//control histograms
- fCommonHistsRes2 = new AliFlowCommonHistResults("Cumulants2");
- fCommonHistsRes4 = new AliFlowCommonHistResults("Cumulants4");
- fCommonHistsRes6 = new AliFlowCommonHistResults("Cumulants6");
- fCommonHistsRes8 = new AliFlowCommonHistResults("Cumulants8");
+ //control histograms
+ fCommonHists = new AliFlowCommonHist("Cumulants");
+ fHistList->Add(fCommonHists->GetHistList()); 
 }
 
 //________________________________________________________________________
 void AliFlowAnalysisWithCumulants::Make(AliFlowEventSimple* anEvent) {
   //running over data
-   
+  
+  //---------------------------------------------------------  
+  //fill the common control histograms
   fCommonHists->FillControlHistograms(anEvent);   
+  //---------------------------------------------------------
+  
   
-  Double_t fG[fgkPmax][fgkQmax];//generating function for integrated flow
+  //---------------------------------------------------------
+  //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.; 
    }
   }
+  //---------------------------------------------------------
+  
   
   //---------------------------------------------------------
   //Q vector stuff 
@@ -129,8 +219,8 @@ void AliFlowAnalysisWithCumulants::Make(AliFlowEventSimple* anEvent) {
   cout<<"Number of input tracks for cumulant analysis: "<<nPrim<<endl;
   cout<<"Number of selected tracks for cumulant analysis: "<<fEventNSelTracksIntFlow<<endl;
   
-    //------------------------------------------------------------------------------------
-    //STARTING THE FIRST LOOP (CALCULATING THE GENERATING FUNCTION FOR INTEGRATED FLOW)
+  //------------------------------------------------------------------------------------
+  //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()){
@@ -141,366 +231,67 @@ void AliFlowAnalysisWithCumulants::Make(AliFlowEventSimple* anEvent) {
        }
       }
     }
-  }
-  // ENDING THE FIRST LOOP OVER TRACKS
-  //------------------------------------------------------------------------------------
+  }// ENDING THE FIRST LOOP OVER TRACKS
+  //------------------------------------------------------------------------------------ 
   
   
+  fHistProAvM->Fill(0.,fSelTracksIntFlow,1);
+  fAvM=fHistProAvM->GetBinContent(1);
+
+  
   //------------------------------------------------------------------------------------
-  //avarage multiplicity
-  fAvM+=fSelTracksIntFlow;
-  //avarage of the generating function for integrated flow
+  //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++){
-      fAvG[p][q]+=1.*fG[p][q];
+   for(Int_t q=0;q<fgkQmax;q++){
+    fIntFlowGenFun->Fill((Double_t)p,(Double_t)q,fG[p][q],1);
    } 
-  }  
+  }
   //------------------------------------------------------------------------------------
-  
-  
-  //STARTING WITH DIFFERENTIAL FLOW...  
-  Int_t fBinEntries[fgknBins]={0};//stores the number of particles per bin for the current event
-  Double_t fNumDRe[fgknBins][fgkPmax][fgkQmax]={0.};//real part of the numerator of D (see relation (10) in PG) 
-  Double_t fNumDIm[fgknBins][fgkPmax][fgkQmax]={0.};//imaginary part of the numerator D   
-  
-  //------------------------------------------------------------------------------------------------
-  //STARTING THE SECOND LOOP OVER TRACKS (CALCULATING THE GENERATING FUNCTION FOR DIFFERENTIAL FLOW)
+   
+                
+  //------------------------------------------------------------------------------------
+  //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()) {
-      Int_t fBin=TMath::Nint(floor(fTrack->Pt()/fBinWidth));
-      if(fBin>=fgknBins)continue;//ignoring the particles with pt>ptMax
-      fBinNoOfParticles[fBin]++;
-      fBinEntries[fBin]++;
-      fBinMeanPt[fBin]+=fTrack->Pt();
-      for(Int_t p=0;p<fgkPmax;p++){
-       for(Int_t q=0;q<fgkQmax;q++){
-         fNumDRe[fBin][p][q]+=fG[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/
+   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))
-         fNumDIm[fBin][p][q]+=fG[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/
+            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)); 
-       }
-      }
+            cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);        
+     }
     }
-  } 
-  //ENDING THE SECOND LOOP OVER TRACKS 
-  //----------------------------------------------------------------------------------------------- 
-  
-  
+   }  
+  }// ENDING THE SECOND LOOP OVER TRACKS
+  //------------------------------------------------------------------------------------
   
-  //----------------------------------------------------------
-  //AVARAGING OVER ALL pt BINS WITHIN ONE EVENT 
-  for(Int_t b=0;b<fgknBins;b++){
-    if(fBinEntries[b]==0)continue;
-    fBinEventEntries[b]++;
-    for(Int_t p=0;p<fgkPmax;p++){
-      for(Int_t q=0;q<fgkQmax;q++){
-       fBinEventDRe[b][p][q]+=fNumDRe[b][p][q]/fBinEntries[b];
-       fBinEventDIm[b][p][q]+=fNumDIm[b][p][q]/fBinEntries[b];
-      }
-    }
-  }
-  //----------------------------------------------------------
-
-fNumberOfEvents++;
-
-}
+}//end of Make()
 
 //________________________________________________________________________
 void AliFlowAnalysisWithCumulants::Finish(){
   //final results
-    
-  Int_t nEvents=fNumberOfEvents;
-  
-  cout<<""<<endl;
-  cout<<"***************************************"<<endl;
-  cout<<"**** results of cumulant analysis: ****"<<endl;
-  cout<<"***************************************"<<endl;
-  cout<<""<<endl;
-  cout<<"number of events = "<<nEvents<<endl;
-  
-  //final avarage multiplicity
-  fAvM/=nEvents;
-  
-  //final avarage of generating function for the integrated flow
-  for(Int_t p=0;p<fgkPmax;p++){
-   for(Int_t q=0;q<fgkQmax;q++){
-    fAvG[p][q]/=nEvents;
-   }
-  }    
-  
-  //final avarage of the Q vector stuff
-  fAvQx/=nEvents;
-  fAvQy/=nEvents;
-  fAvQ2x/=nEvents;
-  fAvQ2y/=nEvents;
-  
-  /////////////////////////////////////////////////////////////////////////////      
-  //////////////////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.*fAvM*(pow(fAvG[p][q],(1./fAvM))-1.);
-    }
-  }
-  
-  /////////////////////////////////////////////////////////////////////////////
-  ///////avaraging the gen. function for the cumulants over azimuth////////////
-  /////////////////////////////////////////////////////////////////////////////
-  
-  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;
-  }
-  
-  /////////////////////////////////////////////////////////////////////////////
-  //////////////////////////////////final results//////////////////////////////
-  /////////////////////////////////////////////////////////////////////////////
-  
-  Double_t fCumulant[fgkPmax];//c_{iFlow}\{2(p+1)\}
-  
-  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]);
+}
+
+
+
+
+
+
+
+
 
-  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]);
 
-  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]);
 
-  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]);
 
-  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]);
 
-  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]);
 
-  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]);
 
-  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]);
-  
-  cout<<""<<endl;
-  cout<<"***************************"<<endl;
-  cout<<"cumulants:"<<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;
-  
-  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;
-     fCommonHistsRes2->FillIntegratedFlow(100.*fV2,100.*fSdQ[0]);
-     fCommonHistsRes2->FillChi(fChiQ[0]);
-    }
-   } else {
-    cout<<" v_"<<fgkFlow<<"{2} = 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;
-     fCommonHistsRes4->FillIntegratedFlow(100.*fV4,100.*fSdQ[1]);
-     fCommonHistsRes4->FillChi(fChiQ[1]);
-    } else {
-      cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
-    }
-  } else {
-    cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;  
-  } 
-  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;
-     fCommonHistsRes6->FillIntegratedFlow(100.*fV6,100.*fSdQ[2]);
-     fCommonHistsRes6->FillChi(fChiQ[2]);
-    } else {
-      cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl; 
-    }
-  } else {
-    cout<<" v_"<<fgkFlow<<"{6} = 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;
-      fCommonHistsRes8->FillIntegratedFlow(100.*fV8,100.*fSdQ[3]);
-      fCommonHistsRes8->FillChi(fChiQ[3]);
-     } else {
-       cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;
-     }
-  } else {
-    cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;     
-  }
-  if (fCumulant[4]>=0.){
-    cout<<"v_"<<fgkFlow<<"{10} = "<<100.*pow((1./456.)*fCumulant[4],(1./10.))<<"%"<<endl;
-  } else {
-    cout<<"v_"<<fgkFlow<<"{10} = Im"<<endl;  
-  }
-  if (fCumulant[5]<=0.){
-    cout<<"v_"<<fgkFlow<<"{12} = "<<100.*pow(-(1./9460.)*fCumulant[5],(1./12.))<<"%"<<endl;
-  } else {
-    cout<<"v_"<<fgkFlow<<"{12} = Im"<<endl;  
-  }
-  if (fCumulant[6]>=0.){
-    cout<<"v_"<<fgkFlow<<"{14} = "<<100.*pow((1./274800.)*fCumulant[6],(1./14.))<<"%"<<endl;
-  } else {
-    cout<<"v_"<<fgkFlow<<"{14} = Im"<<endl;  
-  }
-  if (fCumulant[7]<=0.){
-    cout<<"v_"<<fgkFlow<<"{16} = "<<100.*pow(-(1./10643745.)*fCumulant[7],(1./16.))<<"%"<<endl;
-  } else {
-    cout<<"v_"<<fgkFlow<<"{16} = Im"<<endl;  
-  }
-  cout<<"***************************"<<endl;
-  
-  //cout<<""<<endl;
-  //cout<<"continuing with calculations for differential flow..."<<endl;
-  //cout<<""<<endl;
-  Double_t fBinEventDReAv[fgknBins][fgkPmax][fgkQmax]={0.};
-  Double_t fBinEventDImAv[fgknBins][fgkPmax][fgkQmax]={0.};
-  
-  for(Int_t b=0;b<fgknBins;b++){
-    if(fBinEventEntries[b]==0) continue;
-    for(Int_t p=0;p<fgkPmax;p++){
-      for(Int_t q=0;q<fgkQmax;q++){
-       fBinEventDReAv[b][p][q]=fBinEventDRe[b][p][q]/fBinEventEntries[b];//avarage of the real part of numerator in relation (10) in PG 
-       fBinEventDImAv[b][p][q]=fBinEventDIm[b][p][q]/fBinEventEntries[b];//avarage of the imaginary part of numerator in relation (10) in PG
-      }
-    }                                                                                                                                                                     
-  }
 
-  Double_t fX[fgknBins][fgkPmax][fgkQmax]={0.};//see the text bellow relation (11) in PG
-  Double_t fY[fgknBins][fgkPmax][fgkQmax]={0.};
-  
-  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]=fBinEventDReAv[b][p][q]/fAvG[p][q];
-       fY[b][p][q]=fBinEventDImAv[b][p][q]/fAvG[p][q];
-      }
-    }   
-  } 
-  
-  //cout<<""<<endl;
-  //cout<<"I have calculated X and Y."<<endl;
-  //cout<<""<<endl;
-  
-  Double_t fD[fgknBins][fgkPmax]={0.};//implementing relation (11) from PG
-  
-  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];
-      } 
-      fD[b][p]=1.*(pow(fR0*pow(p+1.,.5),fgkMltpl)/fgkQmax)*fTempHere3;
-    }
-  } 
-  
-  //cout<<""<<endl;
-  //cout<<"calculating differential cumulants now..."<<endl;
-  //cout<<""<<endl;
-  
-  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.};
-  
-  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]);
-  }
-  
-  Double_t fv2[fgknBins],fv4[fgknBins],fv6[fgknBins],fv8[fgknBins];
-  Double_t fAvPt[fgknBins];
-  Double_t fSddiff2[fgknBins],fSddiff4[fgknBins];
-
-  //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: "<<fBinNoOfParticles[b]<<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;
-       fCommonHistsRes2->FillDifferentialFlow(b+1,fv2[b],100.*fSddiff2[b]);
-      } else {
-        //cout<<"v'_2/2{2} = Im"<<endl;
-      }
-    }else{
-      //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]);
-      } else {
-        //cout<<"v'_2/2{4} = Im"<<endl;
-      } 
-    }else{
-      //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.);
-    }else{
-      //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.);
-    }else{
-      //cout<<"v'_2/2{8} = Im"<<endl;
-    }       
-    //cout<<"****************************************"<<endl;
-  }  
-}
 
 
index 219b52f9c95ccd79959f1a3bbe689c3a1959eded..ccce05e5377f4cba951786e3da048385ff8a5a01 100644 (file)
 #include "AliFlowCumuConstants.h"
 
 class TH1;
+class TH3D;
+
 class TObjArray;
+
+
+class TProfile;
+class TProfile2D;
+class TProfile3D;
+
+
 class TList;
 class TFile;
 
@@ -31,10 +40,9 @@ class AliFlowAnalysisWithCumulants {
   virtual void Finish();
 
   // Output 
-  void SetHistFileName(TString name)  {this->fHistFileName = name ;} // Sets output file name
-  TString GetHistFileName() const { return this->fHistFileName ; } // Gets output file name
-  TFile* GetHistFile() const { return this->fHistFile ; }     // Gets output file
-
+  
+  // Output
+  TList*   GetHistList() const             { return this->fHistList ; }     // Gets output histogram list //NEW
 
  private:
   AliFlowAnalysisWithCumulants(const AliFlowAnalysisWithCumulants& aAnalysis);
@@ -47,8 +55,11 @@ class AliFlowAnalysisWithCumulants {
   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
  
-  TString  fHistFileName;      //! The output file name     
-  TFile* fHistFile; // histogram file for Cumulants
+  //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
@@ -63,8 +74,23 @@ class AliFlowAnalysisWithCumulants {
    
   Int_t fNumberOfEvents;//number of events 
    
+  
+  TProfile*          fHistProAvM;       //avarage multiplicity
+  
+  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
+  
+  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;//histograms with various order final results 
+  //AliFlowCommonHistResults *fCommonHistsRes2, *fCommonHistsRes4, *fCommonHistsRes6, *fCommonHistsRes8;//final results//to be improved 
   
   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
@@ -78,3 +104,6 @@ class AliFlowAnalysisWithCumulants {
 #endif
 
 
+
+
+
index 8721bf831f4fe339978e361701deed26e6f870d6..669cf9a329ab932b071d705f8832b50388693f8d 100644 (file)
@@ -1,9 +1,43 @@
+//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 = 100, TString type = "ESD", const Char_t* dataDir="/data/alice2/ab2/", 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) 
+//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/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) 
 {
   TStopwatch timer;
   timer.Start();
@@ -18,21 +52,140 @@ void runAliAnalysisTaskCumulants(Int_t nRuns = 100, TString type = "ESD", const
   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");
-  
-  //  gROOT->LoadMacro("AliFlowLYZConstants.cxx+");
-  //  gROOT->LoadMacro("AliFlowCommonConstants.cxx+");
-  //  gROOT->LoadMacro("AliFlowVector.cxx+");
-  //  gROOT->LoadMacro("AliFlowTrackSimple.cxx+");
-  //  gROOT->LoadMacro("AliFlowEventSimple.cxx+");
-  //  gROOT->LoadMacro("AliFlowEventSimpleMaker.cxx+");
-  //  gROOT->LoadMacro("AliFlowCommonHist.cxx+");
-  //  gROOT->LoadMacro("AliFlowCommonHistResults.cxx+");
-  //  gROOT->LoadMacro("AliFlowLYZHist1.cxx+");
-  //  gROOT->LoadMacro("AliFlowLYZHist2.cxx+");
-  //  gROOT->LoadMacro("AliFlowAnalysisWithLeeYangZeros.cxx+"); 
-  //  gROOT->LoadMacro("AliAnalysisTaskLeeYangZeros.cxx+");
+  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);
   cout<<"chain ("<<chain<<")"<<endl;
@@ -40,21 +193,37 @@ void runAliAnalysisTaskCumulants(Int_t nRuns = 100, TString type = "ESD", const
   //____________________________________________//
   // Make the analysis manager
   AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+   
+  //AliVEventHandler* esdH = new AliESDInputHandler;
+  //mgr->SetInputEventHandler(esdH);  
+  //AliMCEventHandler *mc = new AliMCEventHandler();
+  //mgr->SetMCtruthEventHandler(mc);
+  
+  if (type == "ESD"){
   AliVEventHandler* esdH = new AliESDInputHandler;
-  mgr->SetInputEventHandler(esdH);  
+  mgr->SetInputEventHandler(esdH); }
+   
+  if (type == "AOD"){
+  AliVEventHandler* aodH = new AliAODInputHandler;
+  mgr->SetInputEventHandler(aodH); }
+  
+  if (type == "MC"){
   AliMCEventHandler *mc = new AliMCEventHandler();
-  mgr->SetMCtruthEventHandler(mc);
+  mgr->SetMCtruthEventHandler(mc); }
+  
   //____________________________________________//
-  // 1st Cumulant task
+  // 1st cumulant task
   AliAnalysisTaskCumulants *task1 = new AliAnalysisTaskCumulants("TaskCumulants");
-  task1->SetAnalysisType(type);
+  task1->SetAnalysisType("ESD");
+  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("cobj1", TList::Class(),AliAnalysisManager::kOutputContainer);
+    mgr->CreateContainer("cobj1", TList::Class(),AliAnalysisManager::kOutputContainer,"outputFromCumulantAnalysisESD.root");
 
   //____________________________________________//
   mgr->ConnectInput(task1,0,cinput1);
@@ -158,7 +327,7 @@ TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
   
   return chain;
 }
-
+1
 void LookupWrite(TChain* chain, const char* target)
 {
   // looks up the chain and writes the remaining files to the text file target
diff --git a/PWG2/FLOW/macros/runProofCumulants.C b/PWG2/FLOW/macros/runProofCumulants.C
new file mode 100644 (file)
index 0000000..c8b01da
--- /dev/null
@@ -0,0 +1,217 @@
+//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;
+
+
+  void runProofCumulants(const Char_t* data="/PWG2/akisiel/LHC500C2030", Int_t nRuns=100, 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");
+
+
+//____________________________________________//
+ //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);
+
+
+  //____________________________________________//
+  // 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->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");
+  
+  //____________________________________________//
+  mgr->ConnectInput(task1,0,cinput1);
+  mgr->ConnectOutput(task1,0,coutput1);
+
+  if (!mgr->InitAnalysis()) return;
+  mgr->PrintStatus();
+  // old way with a chain
+  //  mgr->StartAnalysis("proof",chain);
+  mgr->StartAnalysis("proof",data,nRuns,offset);
+
+  timer.Stop();
+  timer.Print();
+}
index cdbc748e7382bf9af29c9f382d534c585d612c5c..370551295884c56ac530c99ebbce7a33c488ad29 100644 (file)
@@ -27,5 +27,6 @@
 #pragma link C++ class AliAnalysisTaskLYZEventPlane+;
 #pragma link C++ class AliAnalysisTaskLeeYangZeros+;
 #pragma link C++ class AliAnalysisTaskCumulants+;
+#pragma link C++ class AliCumulantsFunctions+;
 
 #endif
index a731c0e406d2abf3f32c5e1b442c9a606ed0a648..6116cd59441c9925e1a4b0b23cc8c6d57c98ed8c 100644 (file)
@@ -17,6 +17,7 @@ SRCS= FLOW/AliFlowEventSimple.cxx \
       FLOW/AliFlowAnalysisWithLYZEventPlane.cxx \
       FLOW/AliFlowAnalysisWithLeeYangZeros.cxx \
       FLOW/AliFlowAnalysisWithCumulants.cxx \
+      FLOW/AliCumulantsFunctions.cxx \
       FLOW/AliAnalysisTaskScalarProduct.cxx \
       FLOW/AliAnalysisTaskMCEventPlane.cxx \
       FLOW/AliAnalysisTaskLYZEventPlane.cxx \