PWG2/FLOW/AliAnalysisTaskQCumulants.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   *f
   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 Q-cumulants      *
  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 "AliAnalysisTaskQCumulants.h"
  54 #include "AliFlowEventSimpleMaker.h"
  55 #include "AliFlowAnalysisWithQCumulants.h"
  56 #include "AliFlowCumuConstants.h"
  57 #include "AliFlowCommonConstants.h"
  58 #include "AliFlowCommonHistResults.h"
  59 #include "AliQCumulantsFunctions.h"
  60
  61 ClassImp(AliAnalysisTaskQCumulants)
  62
  63 //================================================================================================================
  64
  65 AliAnalysisTaskQCumulants::AliAnalysisTaskQCumulants(const char *name, Bool_t on):
  66  AliAnalysisTask(name,""),
  67  fESD(NULL),
  68  fAOD(NULL),
  69  fQCA(NULL),//Q-cumulant Analysis (QCA) object
  70  fEventMaker(NULL),
  71  fAnalysisType("ESD"),
  72  fCFManager1(NULL),
  73  fCFManager2(NULL),
  74  fListHistos(NULL),
  75  fQAInt(NULL),
  76  fQADiff(NULL),
  77  fQA(on)
  78 {
  79  //constructor
  80  cout<<"AliAnalysisTaskQCumulants::AliAnalysisTaskQCumulants(const char *name)"<<endl;
  81
  82  // Define input and output slots here
  83  // Input slot #0 works with a TChain
  84  DefineInput(0, TChain::Class());
  85
  86  // Output slot #0 writes into a TList container
  87  DefineOutput(0, TList::Class());
  88  if(on)
  89  {
  90   DefineOutput(1, TList::Class());
  91   DefineOutput(2, TList::Class());
  92  }
  93 }
  94
  95 AliAnalysisTaskQCumulants::AliAnalysisTaskQCumulants():
  96  fESD(NULL),
  97  fAOD(NULL),
  98  fQCA(NULL),//Q-cumulant Analysis (QCA) object
  99  fEventMaker(NULL),
 100  fAnalysisType("ESD"),
 101  fCFManager1(NULL),
 102  fCFManager2(NULL),
 103  fListHistos(NULL),
 104  fQAInt(NULL),
 105  fQADiff(NULL),
 106  fQA(kFALSE)
 107 {
 108  //dummy constructor
 109  cout<<"AliAnalysisTaskQCumulants::AliAnalysisTaskQCumulants()"<<endl;
 110 }
 111
 112 //================================================================================================================
 113
 114 void AliAnalysisTaskQCumulants::ConnectInputData(Option_t *)
 115 {
 116  //connect ESD or AOD (called once)
 117  cout<<"AliAnalysisTaskQCumulants::ConnectInputData(Option_t *)"<<endl;
 118
 119  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 120  if (!tree)
 121  {
 122   Printf("ERROR: Could not read chain from input slot 0");
 123  }
 124  else
 125  {
 126  //disable all branches and enable only the needed ones
 127   if (fAnalysisType == "MC") {
 128      // we want to process only MC
 129       tree->SetBranchStatus("*", kFALSE);
 130
 131       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 132
 133       if (!esdH) {
 134         Printf("ERROR: Could not get ESDInputHandler");
 135       } else {
 136         fESD = esdH->GetEvent();
 137       }
 138     }
 139     else if (fAnalysisType == "ESD" || fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1"  ) {
 140       tree->SetBranchStatus("*", kFALSE);
 141       tree->SetBranchStatus("Tracks.*", kTRUE);
 142
 143       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 144
 145       if (!esdH) {
 146         Printf("ERROR: Could not get ESDInputHandler");
 147       } else
 148         fESD = esdH->GetEvent();
 149     }
 150     else if (fAnalysisType == "AOD") {
 151       AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 152
 153       if (!aodH) {
 154         Printf("ERROR: Could not get AODInputHandler");
 155       }
 156       else {
 157         fAOD = aodH->GetEvent();
 158       }
 159     }
 160     else {
 161       Printf("Wrong analysis type: Only ESD, ESDMC0, ESDMC1, AOD and MC types are allowed!");
 162
 163     }
 164   }
 165 }
 166
 167 //================================================================================================================
 168
 169 void AliAnalysisTaskQCumulants::CreateOutputObjects()
 170 {
 171  //called at every worker node to initialize
 172  cout<<"AliAnalysisTaskQCumulants::CreateOutputObjects()"<<endl;
 173
 174
 175  //OpenFile(0);
 176
 177
 178  if(!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1" || fAnalysisType == "MC"))
 179  {
 180   cout<<"WRONG ANALYSIS TYPE! only ESD, ESDMC0, ESDMC1, AOD and MC are allowed."<<endl;
 181   exit(1);
 182  }
 183
 184  //event maker
 185  fEventMaker = new AliFlowEventSimpleMaker();
 186
 187  //analyser
 188  fQCA = new AliFlowAnalysisWithQCumulants();
 189  fQCA->CreateOutputObjects();
 190
 191  if(fQCA->GetHistList())
 192  {
 193   fListHistos = fQCA->GetHistList();
 194   //fListHistos->Print();
 195  }
 196  else
 197  {
 198   Printf("ERROR: Could not retrieve histogram list");
 199  }
 200
 201  //PostData(0,fListHistos);
 202
 203 }
 204
 205 //================================================================================================================
 206
 207 void AliAnalysisTaskQCumulants::Exec(Option_t *)
 208 {
 209  //main loop (called for each event)
 210  if (fAnalysisType == "MC") {
 211     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 212     // This handler can return the current MC event
 213
 214     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 215     if (!eventHandler) {
 216       Printf("ERROR: Could not retrieve MC event handler");
 217       return;
 218     }
 219
 220     AliMCEvent* mcEvent = eventHandler->MCEvent();
 221     if (!mcEvent) {
 222       Printf("ERROR: Could not retrieve MC event");
 223       return;
 224     }
 225
 226     Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 227     fCFManager1->SetEventInfo(mcEvent);
 228     fCFManager2->SetEventInfo(mcEvent);
 229
 230     //Q-cumulant analysis
 231     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(mcEvent,fCFManager1,fCFManager2);
 232     fQCA->Make(fEvent);
 233     delete fEvent;
 234   }
 235   else if (fAnalysisType == "ESD") {
 236     if (!fESD) {
 237       Printf("ERROR: fESD not available");
 238       return;
 239     }
 240     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 241
 242     //Q-cumulant analysis
 243     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD,fCFManager1,fCFManager2);
 244     //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD);
 245
 246     fQCA->Make(fEvent);
 247
 248     delete fEvent;
 249   }
 250   else if (fAnalysisType == "ESDMC0") {
 251     if (!fESD) {
 252       Printf("ERROR: fESD not available");
 253       return;
 254     }
 255     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 256
 257     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 258     if (!eventHandler) {
 259       Printf("ERROR: Could not retrieve MC event handler");
 260       return;
 261     }
 262
 263     AliMCEvent* mcEvent = eventHandler->MCEvent();
 264     if (!mcEvent) {
 265       Printf("ERROR: Could not retrieve MC event");
 266       return;
 267     }
 268
 269     fCFManager1->SetEventInfo(mcEvent);
 270     fCFManager2->SetEventInfo(mcEvent);
 271
 272     //Q-cumulant analysis
 273     AliFlowEventSimple* fEvent=NULL;
 274     if (fAnalysisType == "ESDMC0") {
 275       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 0); //0 = kine from ESD, 1 = kine from MC
 276     } else if (fAnalysisType == "ESDMC1") {
 277       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 1); //0 = kine from ESD, 1 = kine from MC
 278     }
 279     fQCA->Make(fEvent);
 280     delete fEvent;
 281     //delete mcEvent;
 282   }
 283
 284   else if (fAnalysisType == "AOD") {
 285     if (!fAOD) {
 286       Printf("ERROR: fAOD not available");
 287       return;
 288     }
 289     Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
 290
 291     // analysis
 292     //For the moment don't use CF //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD,fCFManager1,fCFManager2);
 293     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD);
 294     fQCA->Make(fEvent);
 295     delete fEvent;
 296   }
 297
 298   PostData(0,fListHistos);
 299   if(fQA)
 300   {
 301    PostData(1,fQAInt);
 302    PostData(2,fQADiff);
 303   }
 304 }
 305
 306 //================================================================================================================
 307
 308 void AliAnalysisTaskQCumulants::Terminate(Option_t *)
 309 {
 310  //accessing the output list which contains the merged 2D and 3D profiles from all worker nodes
 311  fListHistos = (TList*)GetOutputData(0);
 312  //fListHistos->Print();
 313
 314  if(fListHistos)
 315  {
 316   //final results (integrated flow)
 317   TH1D *intFlowResults = dynamic_cast<TH1D*>(fListHistos->FindObject("fIntFlowResultsQC"));
 318
 319   //final results (differential flow)
 320   TH1D *diffFlowResults2ndOrder = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults2ndOrderQC"));
 321   TH1D *diffFlowResults4thOrder = dynamic_cast<TH1D*>(fListHistos->FindObject("fDiffFlowResults4thOrderQC"));
 322
 323   //final results for covariances (1st bin <2*4>-<2>*<4>, 2nd bin <2*6>-<2>*<6>, ...)
 324   TH1D *covariances = dynamic_cast<TH1D*>(fListHistos->FindObject("fCovariances"));
 325
 326   //common histograms to store the final results for the 2nd order integrated and differential flow
 327   AliFlowCommonHistResults *commonHistRes2nd = dynamic_cast<AliFlowCommonHistResults*>(fListHistos->FindObject("AliFlowCommonHistResults2ndOrderQC"));
 328
 329   //common histograms to store the final results for the 4th order integrated and differential flow
 330   AliFlowCommonHistResults *commonHistRes4th = dynamic_cast<AliFlowCommonHistResults*>(fListHistos->FindObject("AliFlowCommonHistResults4thOrderQC"));
 331
 332   //common histograms to store the final results for the 6th order integrated and differential flow
 333   AliFlowCommonHistResults *commonHistRes6th = dynamic_cast<AliFlowCommonHistResults*>(fListHistos->FindObject("AliFlowCommonHistResults6thOrderQC"));
 334
 335   //common histograms to store the final results for the 8th order integrated and differential flow
 336   AliFlowCommonHistResults *commonHistRes8th = dynamic_cast<AliFlowCommonHistResults*>(fListHistos->FindObject("AliFlowCommonHistResults8thOrderQC"));
 337
 338   //average selected multiplicity (for int. flow)
 339   TProfile *AvMult = dynamic_cast<TProfile*>(fListHistos->FindObject("fAvMultIntFlowQC"));
 340
 341   //multi-particle correlations calculated from Q-vectors
 342   TProfile *QCorrelations = dynamic_cast<TProfile*>(fListHistos->FindObject("fQCorrelations"));
 343
 344   //average of products: 1st bin: <2*4>, 2nd bin: <2*6>, ...
 345   TProfile *QProduct = dynamic_cast<TProfile*>(fListHistos->FindObject("fQProduct"));
 346
 347   //average 2-, 3- and 4-particle correlations per bin
 348   TProfile *binned2p_1n1n = dynamic_cast<TProfile*>(fListHistos->FindObject("f2_1n1n"));
 349   TProfile *binned2p_2n2n = dynamic_cast<TProfile*>(fListHistos->FindObject("f2_2n2n"));
 350   TProfile *binned3p_2n1n1n = dynamic_cast<TProfile*>(fListHistos->FindObject("f3_2n1n1n"));
 351   TProfile *binned3p_1n1n2n = dynamic_cast<TProfile*>(fListHistos->FindObject("f3_1n1n2n"));
 352   TProfile *binned4p_1n1n1n1n = dynamic_cast<TProfile*>(fListHistos->FindObject("f4_1n1n1n1n"));
 353
 354   //average values of Q-vector components (1st bin: <Q_x>, 2nd bin: <Q_y>, 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>)
 355   TProfile *QVectorComponents = dynamic_cast<TProfile*>(fListHistos->FindObject("fQvectorComponents"));
 356
 357   //multi-particle correlations calculated with nested loop
 358   TProfile *DirectCorrelations = dynamic_cast<TProfile*>(fListHistos->FindObject("fDirectCorrelations"));
 359
 360
 361   fQCA = new AliFlowAnalysisWithQCumulants();
 362
 363   fQCA->SetIntFlowResults(intFlowResults);
 364   fQCA->SetDiffFlowResults2nd(diffFlowResults2ndOrder);
 365   fQCA->SetDiffFlowResults4th(diffFlowResults4thOrder);
 366   fQCA->SetCovariances(covariances);
 367
 368   fQCA->SetCommonHistsResults2nd(commonHistRes2nd);
 369   fQCA->SetCommonHistsResults4th(commonHistRes4th);
 370   fQCA->SetCommonHistsResults6th(commonHistRes6th);
 371   fQCA->SetCommonHistsResults8th(commonHistRes8th);
 372
 373   fQCA->SetAverageMultiplicity(AvMult);
 374   fQCA->SetQCorrelations(QCorrelations);
 375   fQCA->SetQProduct(QProduct);
 376   fQCA->SetQVectorComponents(QVectorComponents);
 377
 378   fQCA->SetTwo_1n1nPerBin(binned2p_1n1n);
 379   fQCA->SetTwo_2n2nPerBin(binned2p_2n2n);
 380   fQCA->SetThree_2n1n1nPerBin(binned3p_2n1n1n);
 381   fQCA->SetThree_1n1n2nPerBin(binned3p_1n1n2n);
 382   fQCA->SetFour_1n1n1n1nPerBin(binned4p_1n1n1n1n);
 383
 384   fQCA->SetDirectCorrelations(DirectCorrelations);
 385
 386   fQCA->Finish();
 387  }
 388  else
 389  {
 390   cout<<"histogram list pointer is empty"<<endl;
 391  }
 392 }
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413