]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/FLOW/AliAnalysisTaskCumulants.cxx
make cumulants work on caf (not using tprof3d for the time being)
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliAnalysisTaskCumulants.cxx
index 8ac3914875168dd2d4eff80caa3d9121b38df073..5808406afab3dad6b30fce21343bce69a86a6514 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;
+ }
 }