PWG2/FLOW/AliAnalysisTaskFittingQDistribution.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 fitting       *
  18  *         q-distribution             *
  19  *                                    *
  20  * authors: Naomi van der Kolk        *
  21  *           (kolk@nikhef.nl)         *
  22  *          Raimond Snellings         *
  23  *           (snelling@nikhef.nl)     *
  24  *          Ante Bilandzic            *
  25  *           (anteb@nikhef.nl)        *
  26  * ***********************************/
  27
  28 #include "Riostream.h"
  29 #include "TChain.h"
  30 #include "TTree.h"
  31 #include "TFile.h"
  32 #include "TList.h"
  33 #include "TH1.h"
  34 #include "TProfile.h"
  35
  36 #include "AliAnalysisTask.h"
  37 #include "AliAnalysisDataSlot.h"
  38 #include "AliAnalysisDataContainer.h"
  39 #include "AliAnalysisManager.h"
  40
  41 #include "AliESDEvent.h"
  42 #include "AliESDInputHandler.h"
  43
  44 #include "AliAODEvent.h"
  45 #include "AliAODInputHandler.h"
  46
  47 #include "AliMCEventHandler.h"
  48 #include "AliMCEvent.h"
  49
  50 #include "../../CORRFW/AliCFManager.h"
  51
  52 #include "AliAnalysisTaskFittingQDistribution.h"
  53 #include "AliFlowEventSimpleMaker.h"
  54 #include "AliFittingQDistribution.h"
  55 #include "AliFlowCommonConstants.h"
  56 #include "AliFlowCommonHistResults.h"
  57 #include "AliFittingFunctionsForQDistribution.h"
  58
  59 ClassImp(AliAnalysisTaskFittingQDistribution)
  60
  61 //================================================================================================================
  62
  63 AliAnalysisTaskFittingQDistribution::AliAnalysisTaskFittingQDistribution(const char *name, Bool_t on):
  64  AliAnalysisTask(name,""),
  65  fESD(NULL),
  66  fAOD(NULL),
  67  fFQDA(NULL),//Fitting Q_Distribution Analysis (FQDA) object
  68  fEventMaker(NULL),
  69  fAnalysisType("ESD"),
  70  fCFManager1(NULL),
  71  fCFManager2(NULL),
  72  fListHistos(NULL),
  73  fQAInt(NULL),
  74  fQADiff(NULL),
  75  fQA(on)
  76 {
  77  //constructor
  78  cout<<"AliAnalysisTaskFittingQDistribution::AliAnalysisTaskFittingQDistribution(const char *name)"<<endl;
  79
  80  // Define input and output slots here
  81  // Input slot #0 works with a TChain
  82  DefineInput(0, TChain::Class());
  83
  84  // Output slot #0 writes into a TList container
  85  DefineOutput(0, TList::Class());
  86  if(on)
  87  {
  88   DefineOutput(1, TList::Class());
  89   DefineOutput(2, TList::Class());
  90  }
  91 }
  92
  93 AliAnalysisTaskFittingQDistribution::AliAnalysisTaskFittingQDistribution():
  94  fESD(NULL),
  95  fAOD(NULL),
  96  fFQDA(NULL),//Fitting q-distribution Analysis (FQDA) object
  97  fEventMaker(NULL),
  98  fAnalysisType("ESD"),
  99  fCFManager1(NULL),
 100  fCFManager2(NULL),
 101  fListHistos(NULL),
 102  fQAInt(NULL),
 103  fQADiff(NULL),
 104  fQA(kFALSE)
 105 {
 106  //dummy constructor
 107  cout<<"AliAnalysisTaskFittingQDistribution::AliAnalysisTaskFittingQDistribution()"<<endl;
 108 }
 109
 110 //================================================================================================================
 111
 112 void AliAnalysisTaskFittingQDistribution::ConnectInputData(Option_t *)
 113 {
 114  //connect ESD or AOD (called once)
 115  cout<<"AliAnalysisTaskFittingQDistribution::ConnectInputData(Option_t *)"<<endl;
 116
 117  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 118  if (!tree)
 119  {
 120   Printf("ERROR: Could not read chain from input slot 0");
 121  }
 122  else
 123  {
 124  //disable all branches and enable only the needed ones
 125   if (fAnalysisType == "MC") {
 126      // we want to process only MC
 127       tree->SetBranchStatus("*", kFALSE);
 128
 129       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 130
 131       if (!esdH) {
 132         Printf("ERROR: Could not get ESDInputHandler");
 133       } else {
 134         fESD = esdH->GetEvent();
 135       }
 136     }
 137     else if (fAnalysisType == "ESD" || fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1"  ) {
 138       tree->SetBranchStatus("*", kFALSE);
 139       tree->SetBranchStatus("Tracks.*", kTRUE);
 140
 141       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 142
 143       if (!esdH) {
 144         Printf("ERROR: Could not get ESDInputHandler");
 145       } else
 146         fESD = esdH->GetEvent();
 147     }
 148     else if (fAnalysisType == "AOD") {
 149       AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 150
 151       if (!aodH) {
 152         Printf("ERROR: Could not get AODInputHandler");
 153       }
 154       else {
 155         fAOD = aodH->GetEvent();
 156       }
 157     }
 158     else {
 159       Printf("Wrong analysis type: Only ESD, ESDMC0, ESDMC1, AOD and MC types are allowed!");
 160
 161     }
 162   }
 163 }
 164
 165 //================================================================================================================
 166
 167 void AliAnalysisTaskFittingQDistribution::CreateOutputObjects()
 168 {
 169  //called at every worker node to initialize
 170  cout<<"AliAnalysisTaskFittingQDistribution::CreateOutputObjects()"<<endl;
 171
 172
 173  //OpenFile(0);
 174
 175
 176  if(!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1" || fAnalysisType == "MC"))
 177  {
 178   cout<<"WRONG ANALYSIS TYPE! only ESD, ESDMC0, ESDMC1, AOD and MC are allowed."<<endl;
 179   exit(1);
 180  }
 181
 182  //event maker
 183  fEventMaker = new AliFlowEventSimpleMaker();
 184
 185  //analyser
 186  fFQDA = new AliFittingQDistribution();
 187  fFQDA->CreateOutputObjects();
 188
 189  if(fFQDA->GetHistList())
 190  {
 191   fListHistos = fFQDA->GetHistList();
 192   //fListHistos->Print();
 193  }
 194  else
 195  {
 196   Printf("ERROR: Could not retrieve histogram list");
 197  }
 198
 199  //PostData(0,fListHistos);
 200
 201 }
 202
 203 //================================================================================================================
 204
 205 void AliAnalysisTaskFittingQDistribution::Exec(Option_t *)
 206 {
 207  //main loop (called for each event)
 208  if (fAnalysisType == "MC") {
 209     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 210     // This handler can return the current MC event
 211
 212     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 213     if (!eventHandler) {
 214       Printf("ERROR: Could not retrieve MC event handler");
 215       return;
 216     }
 217
 218     AliMCEvent* mcEvent = eventHandler->MCEvent();
 219     if (!mcEvent) {
 220       Printf("ERROR: Could not retrieve MC event");
 221       return;
 222     }
 223
 224     Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 225     fCFManager1->SetEventInfo(mcEvent);
 226     fCFManager2->SetEventInfo(mcEvent);
 227
 228     //fitting q-distribution
 229     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(mcEvent,fCFManager1,fCFManager2);
 230     fFQDA->Make(fEvent);
 231     delete fEvent;
 232   }
 233   else if (fAnalysisType == "ESD") {
 234     if (!fESD) {
 235       Printf("ERROR: fESD not available");
 236       return;
 237     }
 238     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 239
 240     //fitting q-distribution
 241     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD,fCFManager1,fCFManager2);
 242     //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD);
 243     fFQDA->Make(fEvent);
 244     delete fEvent;
 245   }
 246   else if (fAnalysisType == "ESDMC0") {
 247     if (!fESD) {
 248       Printf("ERROR: fESD not available");
 249       return;
 250     }
 251     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 252
 253     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 254     if (!eventHandler) {
 255       Printf("ERROR: Could not retrieve MC event handler");
 256       return;
 257     }
 258
 259     AliMCEvent* mcEvent = eventHandler->MCEvent();
 260     if (!mcEvent) {
 261       Printf("ERROR: Could not retrieve MC event");
 262       return;
 263     }
 264
 265     fCFManager1->SetEventInfo(mcEvent);
 266     fCFManager2->SetEventInfo(mcEvent);
 267
 268     //fitting q-distribution
 269     AliFlowEventSimple* fEvent=NULL;
 270     if (fAnalysisType == "ESDMC0") {
 271       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 0); //0 = kine from ESD, 1 = kine from MC
 272     } else if (fAnalysisType == "ESDMC1") {
 273       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 1); //0 = kine from ESD, 1 = kine from MC
 274     }
 275     fFQDA->Make(fEvent);
 276     delete fEvent;
 277     //delete mcEvent;
 278   }
 279
 280   else if (fAnalysisType == "AOD") {
 281     if (!fAOD) {
 282       Printf("ERROR: fAOD not available");
 283       return;
 284     }
 285     Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
 286
 287     // analysis
 288     //For the moment don't use CF //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD,fCFManager1,fCFManager2);
 289     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD);
 290     fFQDA->Make(fEvent);
 291     delete fEvent;
 292   }
 293
 294   PostData(0,fListHistos);
 295   if(fQA)
 296   {
 297    PostData(1,fQAInt);
 298    PostData(2,fQADiff);
 299   }
 300 }
 301
 302 //================================================================================================================
 303
 304 void AliAnalysisTaskFittingQDistribution::Terminate(Option_t *)
 305 {
 306  //accessing the output list
 307  fListHistos = (TList*)GetOutputData(0);
 308  //fListHistos->Print();
 309
 310  if(fListHistos)
 311  {
 312   //final results (integrated flow)
 313   TH1D *intFlowResults = dynamic_cast<TH1D*>(fListHistos->FindObject("fIntFlowResultsFQD"));
 314
 315   //sigma^2
 316   TH1D *sigma2 = dynamic_cast<TH1D*>(fListHistos->FindObject("fSigma2"));
 317
 318   //common histograms to store the final results for the integrated flow
 319   AliFlowCommonHistResults *commonHistRes = dynamic_cast<AliFlowCommonHistResults*>(fListHistos->FindObject("AliFlowCommonHistResultsFQD"));
 320
 321   //average selected multiplicity (for int. flow)
 322   TProfile *AvMult = dynamic_cast<TProfile*>(fListHistos->FindObject("fAvMultIntFlowFQD"));
 323
 324   //q-distribution
 325   TH1D *qDist = dynamic_cast<TH1D*>(fListHistos->FindObject("fQDistributionFQD"));
 326
 327   //----------------------------------------------------
 328
 329   fFQDA = new AliFittingQDistribution();
 330
 331   fFQDA->SetIntFlowResults(intFlowResults);
 332   fFQDA->SetSigma2(sigma2);
 333   fFQDA->SetCommonHistsResults(commonHistRes);
 334
 335   fFQDA->SetAverageMultiplicity(AvMult);
 336   fFQDA->SetQDistribution(qDist);
 337
 338   fFQDA->Finish();
 339
 340   //----------------------------------------------------
 341  }
 342  else
 343  {
 344   cout<<"histogram list pointer is empty"<<endl;
 345  }
 346 }
 347
 348 //================================================================================================================
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367