PWG2/FLOW/AliAnalysisTaskCumulants.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /**************************************
  17  * analysis task for cumulant method  *
  18  *                                    *
  19  * authors: Naomi van der Kolk        *
  20  *           (kolk@nikhef.nl)         *
  21  *          Raimond Snellings         *
  22  *           (snelling@nikhef.nl)     *
  23  *          Ante Bilandzic            *
  24  *           (anteb@nikhef.nl)        *
  25  * ***********************************/
  26
  27 #include "Riostream.h"
  28 #include "TChain.h"
  29 #include "TTree.h"
  30 #include "TFile.h"
  31 #include "TList.h"
  32 #include "TH1.h"
  33 #include "TProfile.h"
  34 #include "TProfile2D.h"
  35 #include "TProfile3D.h"
  36
  37 #include "AliAnalysisTask.h"
  38 #include "AliAnalysisDataSlot.h"
  39 #include "AliAnalysisDataContainer.h"
  40 #include "AliAnalysisManager.h"
  41
  42 #include "AliESDEvent.h"
  43 #include "AliESDInputHandler.h"
  44
  45 #include "AliAODEvent.h"
  46 #include "AliAODInputHandler.h"
  47
  48 #include "AliMCEventHandler.h"
  49 #include "AliMCEvent.h"
  50
  51 #include "../../CORRFW/AliCFManager.h"
  52
  53 #include "AliAnalysisTaskCumulants.h"
  54 #include "AliFlowEventSimpleMaker.h"
  55 #include "AliFlowAnalysisWithCumulants.h"
  56 #include "AliFlowCumuConstants.h"
  57 #include "AliFlowCommonConstants.h"
  58 #include "AliFlowCommonHistResults.h"
  59 #include "AliCumulantsFunctions.h"
  60
  61 ClassImp(AliAnalysisTaskCumulants)
  62
  63 //================================================================================================================
  64
  65 AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name):
  66  AliAnalysisTask(name,""),
  67  fESD(NULL),
  68  fAOD(NULL),
  69  fCA(NULL),//Cumulant Analysis (CA) object
  70  fEventMaker(NULL),
  71  fAnalysisType("ESD"),
  72  fCFManager1(NULL),
  73  fCFManager2(NULL),
  74  fListHistos(NULL)
  75 {
  76  //constructor
  77  cout<<"AliAnalysisTaskCumulants::AliAnalysisTaskCumulants(const char *name)"<<endl;
  78  //input and output slots:
  79  //input slot #0 works with a TChain
  80  DefineInput(0, TChain::Class());
  81  //output slot #0 writes into a TList container
  82  DefineOutput(0, TList::Class());
  83 }
  84
  85 AliAnalysisTaskCumulants::AliAnalysisTaskCumulants():
  86  fESD(NULL),
  87  fAOD(NULL),
  88  fCA(NULL),//Cumulant Analysis (CA) object
  89  fEventMaker(NULL),
  90  fAnalysisType("ESD"),
  91  fCFManager1(NULL),
  92  fCFManager2(NULL),
  93  fListHistos(NULL)
  94 {
  95  //dummy constructor
  96  cout<<"AliAnalysisTaskCumulants::AliAnalysisTaskCumulants()"<<endl;
  97 }
  98
  99 //================================================================================================================
 100
 101 void AliAnalysisTaskCumulants::ConnectInputData(Option_t *)
 102 {
 103  //connect ESD or AOD (called once)
 104  cout<<"AliAnalysisTaskCumulants::ConnectInputData(Option_t *)"<<endl;
 105
 106  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 107  if (!tree)
 108  {
 109   Printf("ERROR: Could not read chain from input slot 0");
 110  }
 111  else
 112  {
 113  //disable all branches and enable only the needed ones
 114   if (fAnalysisType == "MC") {
 115      // we want to process only MC
 116       tree->SetBranchStatus("*", kFALSE);
 117
 118       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 119
 120       if (!esdH) {
 121         Printf("ERROR: Could not get ESDInputHandler");
 122       } else {
 123         fESD = esdH->GetEvent();
 124       }
 125     }
 126     else if (fAnalysisType == "ESD" || fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1"  ) {
 127       tree->SetBranchStatus("*", kFALSE);
 128       tree->SetBranchStatus("Tracks.*", kTRUE);
 129
 130       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 131
 132       if (!esdH) {
 133         Printf("ERROR: Could not get ESDInputHandler");
 134       } else
 135         fESD = esdH->GetEvent();
 136     }
 137     else if (fAnalysisType == "AOD") {
 138       AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 139
 140       if (!aodH) {
 141         Printf("ERROR: Could not get AODInputHandler");
 142       }
 143       else {
 144         fAOD = aodH->GetEvent();
 145       }
 146     }
 147     else {
 148       Printf("Wrong analysis type: Only ESD, ESDMC0, ESDMC1, AOD and MC types are allowed!");
 149
 150     }
 151   }
 152 }
 153
 154 //================================================================================================================
 155
 156 void AliAnalysisTaskCumulants::CreateOutputObjects()
 157 {
 158  //called at every worker node to initialize
 159  cout<<"AliAnalysisTaskCumulants::CreateOutputObjects()"<<endl;
 160
 161
 162  //OpenFile(0);
 163
 164
 165  if(!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1" || fAnalysisType == "MC"))
 166  {
 167   cout<<"WRONG ANALYSIS TYPE! only ESD, ESDMC0, ESDMC1, AOD and MC are allowed."<<endl;
 168   exit(1);
 169  }
 170
 171  //event maker
 172  fEventMaker = new AliFlowEventSimpleMaker();
 173
 174  //analyser
 175  fCA = new AliFlowAnalysisWithCumulants();
 176  fCA->CreateOutputObjects();
 177
 178  if(fCA->GetHistList())
 179  {
 180   fListHistos = fCA->GetHistList();
 181   //fListHistos->Print();
 182  }
 183  else
 184  {
 185   Printf("ERROR: Could not retrieve histogram list");
 186  }
 187
 188  //PostData(0,fListHistos);
 189
 190 }
 191
 192 //================================================================================================================
 193
 194 void AliAnalysisTaskCumulants::Exec(Option_t *)
 195 {
 196  //main loop (called for each event)
 197  if (fAnalysisType == "MC") {
 198     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 199     // This handler can return the current MC event
 200
 201     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 202     if (!eventHandler) {
 203       Printf("ERROR: Could not retrieve MC event handler");
 204       return;
 205     }
 206
 207     AliMCEvent* mcEvent = eventHandler->MCEvent();
 208     if (!mcEvent) {
 209       Printf("ERROR: Could not retrieve MC event");
 210       return;
 211     }
 212
 213     Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 214     fCFManager1->SetEventInfo(mcEvent);
 215     fCFManager2->SetEventInfo(mcEvent);
 216
 217     //cumulant analysis
 218     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(mcEvent,fCFManager1,fCFManager2);
 219     fCA->Make(fEvent);
 220     delete fEvent;
 221   }
 222   else if (fAnalysisType == "ESD") {
 223     if (!fESD) {
 224       Printf("ERROR: fESD not available");
 225       return;
 226     }
 227     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 228
 229     //cumulant analysis
 230     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD);
 231     fCA->Make(fEvent);
 232     delete fEvent;
 233   }
 234   else if (fAnalysisType == "ESDMC0") {
 235     if (!fESD) {
 236       Printf("ERROR: fESD not available");
 237       return;
 238     }
 239     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 240
 241     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 242     if (!eventHandler) {
 243       Printf("ERROR: Could not retrieve MC event handler");
 244       return;
 245     }
 246
 247     AliMCEvent* mcEvent = eventHandler->MCEvent();
 248     if (!mcEvent) {
 249       Printf("ERROR: Could not retrieve MC event");
 250       return;
 251     }
 252
 253     fCFManager1->SetEventInfo(mcEvent);
 254     fCFManager2->SetEventInfo(mcEvent);
 255
 256     //cumulant analysis
 257     AliFlowEventSimple* fEvent=NULL;
 258     if (fAnalysisType == "ESDMC0") {
 259       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 0); //0 = kine from ESD, 1 = kine from MC
 260     } else if (fAnalysisType == "ESDMC1") {
 261       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 1); //0 = kine from ESD, 1 = kine from MC
 262     }
 263     fCA->Make(fEvent);
 264     delete fEvent;
 265     //delete mcEvent;
 266   }
 267
 268   else if (fAnalysisType == "AOD") {
 269     if (!fAOD) {
 270       Printf("ERROR: fAOD not available");
 271       return;
 272     }
 273     Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
 274
 275     // analysis
 276     //For the moment don't use CF //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD,fCFManager1,fCFManager2);
 277     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD);
 278     fCA->Make(fEvent);
 279     delete fEvent;
 280   }
 281
 282
 283   PostData(0,fListHistos);
 284 }
 285
 286 //================================================================================================================
 287
 288 void AliAnalysisTaskCumulants::Terminate(Option_t *)
 289 {
 290  //accessing the output list which contains the merged 2D and 3D profiles from all worker nodes
 291  fListHistos = (TList*)GetOutputData(0);
 292  //fListHistos->Print();
 293
 294  if(fListHistos)
 295  {
 296   //profiles with avarage values of generating functions for int. and diff. flow
 297   TProfile2D *intFlowGenFun = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fIntFlowGenFun"));
 298   //TProfile3D *diffFlowGenFunRe=dynamic_cast<TProfile3D*>(fListHistos->FindObject("fDiffFlowGenFunRe"));
 299   //TProfile3D *diffFlowGenFunIm=dynamic_cast<TProfile3D*>(fListHistos->FindObject("fDiffFlowGenFunIm"));
 300
 301   TProfile2D *diffFlowGenFunRe0 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe0"));
 302   TProfile2D *diffFlowGenFunRe1 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe1"));
 303   TProfile2D *diffFlowGenFunRe2 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe2"));
 304   TProfile2D *diffFlowGenFunRe3 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe3"));
 305   TProfile2D *diffFlowGenFunRe4 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe4"));
 306   TProfile2D *diffFlowGenFunRe5 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe5"));
 307   TProfile2D *diffFlowGenFunRe6 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe6"));
 308   TProfile2D *diffFlowGenFunRe7 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunRe7"));
 309   TProfile2D *diffFlowGenFunIm0 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm0"));
 310   TProfile2D *diffFlowGenFunIm1 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm1"));
 311   TProfile2D *diffFlowGenFunIm2 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm2"));
 312   TProfile2D *diffFlowGenFunIm3 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm3"));
 313   TProfile2D *diffFlowGenFunIm4 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm4"));
 314   TProfile2D *diffFlowGenFunIm5 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm5"));
 315   TProfile2D *diffFlowGenFunIm6 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm6"));
 316   TProfile2D *diffFlowGenFunIm7 = dynamic_cast<TProfile2D*>(fListHistos->FindObject("fDiffFlowGenFunIm7"));
 317
 318   //histograms to store the final results
 319   TH1D *intFlowResults   = dynamic_cast<TH1D*>(fListHistos->FindObject("fIntFlowResults"));
 320   TH1D *diffFlowResults2 = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults2"));
 321   TH1D *diffFlowResults4 = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults4"));
 322   TH1D *diffFlowResults6 = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults6"));
 323   TH1D *diffFlowResults8 = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults8"));
 324
 325   //profile with avarage selected multiplicity for int. flow
 326   TProfile *avMult = dynamic_cast<TProfile*>(fListHistos->FindObject("fAvMultIntFlow"));
 327
 328   //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>)
 329   TProfile *QVectorComponents = dynamic_cast<TProfile*>(fListHistos->FindObject("fQVectorComponents"));
 330
 331   //q-distribution
 332   TH1D *qDist = dynamic_cast<TH1D*>(fListHistos->FindObject("fQDist"));
 333
 334   AliCumulantsFunctions finalResults(intFlowGenFun,NULL,NULL, intFlowResults,diffFlowResults2,diffFlowResults4,diffFlowResults6,diffFlowResults8,avMult,QVectorComponents,qDist,diffFlowGenFunRe0,diffFlowGenFunRe1,diffFlowGenFunRe2,
 335 diffFlowGenFunRe3,diffFlowGenFunRe4,diffFlowGenFunRe5,diffFlowGenFunRe6,diffFlowGenFunRe7,diffFlowGenFunIm0,diffFlowGenFunIm1,
 336 diffFlowGenFunIm2,diffFlowGenFunIm3,diffFlowGenFunIm4,diffFlowGenFunIm5,diffFlowGenFunIm6,diffFlowGenFunIm7);
 337
 338   finalResults.Calculate();
 339
 340  }
 341  else
 342  {
 343   cout<<"histogram list pointer is empty"<<endl;
 344  }
 345 }
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366