PWG2/FLOW/AliAnalysisTaskLeeYangZeros.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 #include "Riostream.h" //needed as include
  17 #include "TChain.h"
  18 #include "TTree.h"
  19 #include "TFile.h"
  20 #include "TList.h"
  21 #include "TProfile.h"
  22
  23 class AliAnalysisTask;
  24 #include "AliAnalysisManager.h"
  25
  26 #include "AliESDEvent.h"
  27 #include "AliESDInputHandler.h"
  28
  29 #include "AliAODEvent.h"
  30 #include "AliAODInputHandler.h"
  31
  32 #include "AliMCEventHandler.h"
  33 #include "AliMCEvent.h"
  34
  35 #include "../../CORRFW/AliCFManager.h"
  36
  37 #include "AliFlowLYZConstants.h"
  38 #include "AliAnalysisTaskLeeYangZeros.h"
  39 #include "AliFlowEventSimpleMaker.h"
  40 #include "AliFlowCommonHist.h"
  41 #include "AliFlowCommonHistResults.h"
  42 #include "AliFlowLYZHist1.h"
  43 #include "AliFlowLYZHist2.h"
  44 #include "AliFlowAnalysisWithLeeYangZeros.h"
  45
  46 // AliAnalysisTaskLeeYangZeros:
  47 // analysis task for Lee Yang Zeros method
  48 // Author: Naomi van der Kolk (kolk@nikhef.nl)
  49
  50 ClassImp(AliAnalysisTaskLeeYangZeros)
  51
  52 //________________________________________________________________________
  53 AliAnalysisTaskLeeYangZeros::AliAnalysisTaskLeeYangZeros(const char *name, Bool_t firstrun, Bool_t on) :
  54   AliAnalysisTask(name, ""),
  55   fESD(0),
  56   fAOD(0),
  57   fAnalysisType("ESD"),
  58   fCFManager1(NULL),
  59   fCFManager2(NULL),
  60   fLyz(0),
  61   fEventMaker(0),
  62   fFirstRunFile(0),
  63   fListHistos(NULL),
  64   fQAInt(NULL),
  65   fQADiff(NULL),
  66   fFirstRunLYZ(firstrun),  //set boolean for firstrun to initial value
  67   fUseSumLYZ(kTRUE),       //set boolean for use sum to initial value
  68   fQA(on)
  69 {
  70   // Constructor
  71   cout<<"AliAnalysisTaskLeeYangZeros::AliAnalysisTaskLeeYangZeros(const char *name)"<<endl;
  72
  73   // Define input and output slots here
  74   // Input slot #0 works with a TChain
  75   DefineInput(0, TChain::Class());
  76   if (!firstrun) DefineInput(1, TList::Class()); //for second loop
  77   // Output slot #0 writes into a TList container
  78   DefineOutput(0, TList::Class());
  79   if(on) {
  80     DefineOutput(1, TList::Class());
  81     DefineOutput(2, TList::Class()); }
  82 }
  83
  84
  85
  86 //________________________________________________________________________
  87 AliAnalysisTaskLeeYangZeros::AliAnalysisTaskLeeYangZeros() :
  88   fESD(0),
  89   fAOD(0),
  90   fAnalysisType("ESD"),
  91   fCFManager1(NULL),
  92   fCFManager2(NULL),
  93   fLyz(0),
  94   fEventMaker(0),
  95   fFirstRunFile(0),
  96   fListHistos(NULL),
  97   fQAInt(NULL),
  98   fQADiff(NULL),
  99   fFirstRunLYZ(kTRUE), //set boolean for firstrun to initial value
 100   fUseSumLYZ(kTRUE),    //set boolean for use sum to initial value
 101   fQA(kFALSE)
 102 {
 103   // Constructor
 104   cout<<"AliAnalysisTaskLeeYangZeros::AliAnalysisTaskLeeYangZeros()"<<endl;
 105
 106 }
 107
 108 //________________________________________________________________________
 109 AliAnalysisTaskLeeYangZeros::~AliAnalysisTaskLeeYangZeros()
 110 {
 111
 112   //destructor
 113
 114 }
 115
 116 //________________________________________________________________________
 117 void AliAnalysisTaskLeeYangZeros::ConnectInputData(Option_t *)
 118 {
 119   // Connect ESD or AOD here
 120   // Called once
 121   cout<<"AliAnalysisTaskLeeYangZeros::ConnectInputData(Option_t *)"<<endl;
 122
 123   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 124   if (!tree) {
 125     Printf("ERROR: Could not read chain from input slot 0");
 126   }
 127   else {
 128     // Disable all branches and enable only the needed ones
 129     if (fAnalysisType == "MC") {
 130       // we want to process only MC
 131       tree->SetBranchStatus("*", kFALSE);
 132
 133       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 134
 135       if (!esdH) {
 136         Printf("ERROR: Could not get ESDInputHandler");
 137       } else {
 138         fESD = esdH->GetEvent();
 139       }
 140     }
 141     else if (fAnalysisType == "ESD" || fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1"  ) {
 142       tree->SetBranchStatus("*", kFALSE);
 143       tree->SetBranchStatus("Tracks.*", kTRUE);
 144
 145       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 146
 147       if (!esdH) {
 148         Printf("ERROR: Could not get ESDInputHandler");
 149       } else
 150         fESD = esdH->GetEvent();
 151     }
 152     else if (fAnalysisType == "AOD") {
 153       AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 154
 155       if (!aodH) {
 156         Printf("ERROR: Could not get AODInputHandler");
 157       }
 158       else {
 159         fAOD = aodH->GetEvent();
 160       }
 161     }
 162     else {
 163       Printf("Wrong analysis type: Only ESD, ESDMC0, ESDMC1, AOD and MC types are allowed!");
 164
 165     }
 166   }
 167 }
 168
 169 //________________________________________________________________________
 170 void AliAnalysisTaskLeeYangZeros::CreateOutputObjects()
 171 {
 172   // Called once
 173   cout<<"AliAnalysisTaskLeeYangZeros::CreateOutputObjects()"<<endl;
 174
 175   if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMC0"  || fAnalysisType == "ESDMC1" || fAnalysisType == "MC")) {
 176     cout<<"WRONG ANALYSIS TYPE! only ESD, ESDMC0, ESDMC1, AOD and MC are allowed."<<endl;
 177     exit(1);
 178   }
 179
 180   //event maker
 181   fEventMaker = new AliFlowEventSimpleMaker();
 182   //Analyser
 183   fLyz = new AliFlowAnalysisWithLeeYangZeros() ;
 184
 185   fLyz -> SetFirstRun(GetFirstRunLYZ());   //set first run true or false
 186   fLyz -> SetUseSum(GetUseSumLYZ());       //set use sum true or false
 187
 188   // Get data from input slot 1
 189   if (GetNinputs() == 2) {                   //if there are two input slots
 190     TList* pFirstRunList = (TList*)GetInputData(1);
 191     if (pFirstRunList) {
 192       fLyz -> SetFirstRunList(pFirstRunList);
 193     } else { cout<<"No first run List!"<<endl; exit(0); }
 194   }
 195
 196   fLyz -> Init();
 197
 198   if (fLyz->GetHistList()) {
 199     fListHistos = fLyz->GetHistList();
 200     //    fListHistos->Print();
 201   }
 202   else {Printf("ERROR: Could not retrieve histogram list"); }
 203
 204 }
 205
 206 //________________________________________________________________________
 207 void AliAnalysisTaskLeeYangZeros::Exec(Option_t *)
 208 {
 209   // Main loop
 210   // Called for each event
 211   if (fAnalysisType == "MC") {
 212     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 213     // This handler can return the current MC event
 214
 215     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 216     if (!eventHandler) {
 217       Printf("ERROR: Could not retrieve MC event handler");
 218       return;
 219     }
 220
 221     AliMCEvent* mcEvent = eventHandler->MCEvent();
 222     if (!mcEvent) {
 223       Printf("ERROR: Could not retrieve MC event");
 224       return;
 225     }
 226
 227     fCFManager1->SetEventInfo(mcEvent);
 228     fCFManager2->SetEventInfo(mcEvent);
 229
 230     Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 231
 232     //lee yang zeros analysis
 233     //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(mcEvent);
 234     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(mcEvent,fCFManager1,fCFManager2);
 235     fLyz->Make(fEvent);
 236     delete fEvent;
 237   }
 238   else if (fAnalysisType == "ESD") {
 239     if (!fESD) {
 240       Printf("ERROR: fESD not available");
 241       return;
 242     }
 243     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 244
 245     //lee yang zeros analysis
 246     //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD);
 247     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD,fCFManager1,fCFManager2);
 248     fLyz->Make(fEvent);
 249     delete fEvent;
 250   }
 251   else if (fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1") {
 252     if (!fESD) {
 253       Printf("ERROR: fESD not available");
 254       return;
 255     }
 256     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 257
 258     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 259     if (!eventHandler) {
 260       Printf("ERROR: Could not retrieve MC event handler");
 261       return;
 262     }
 263
 264     AliMCEvent* mcEvent = eventHandler->MCEvent();
 265     if (!mcEvent) {
 266       Printf("ERROR: Could not retrieve MC event");
 267       return;
 268     }
 269
 270     fCFManager1->SetEventInfo(mcEvent);
 271     fCFManager2->SetEventInfo(mcEvent);
 272
 273     //lee yang zeros analysis
 274     //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fESD,mcEvent,0); //0 = kine from ESD, 1 = kine from MC
 275     AliFlowEventSimple* fEvent=NULL;
 276     if (fAnalysisType == "ESDMC0") {
 277       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 0); //0 = kine from ESD, 1 = kine from MC
 278     } else if (fAnalysisType == "ESDMC1") {
 279       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 1); //0 = kine from ESD, 1 = kine from MC
 280     }
 281     fLyz->Make(fEvent);
 282     delete fEvent;
 283     //delete mcEvent;
 284   }
 285
 286   else if (fAnalysisType == "AOD") {
 287     if (!fAOD) {
 288       Printf("ERROR: fAOD not available");
 289       return;
 290     }
 291     Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
 292
 293     //lee yang zeros analysis
 294     AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD); //no CF yet!
 295     fLyz->Make(fEvent);
 296     delete fEvent;
 297   }
 298
 299   PostData(0,fListHistos); //here for CAF
 300   if (fQA) {
 301     PostData(1,fQAInt);
 302     PostData(2,fQADiff); }
 303 }
 304
 305 //________________________________________________________________________
 306 void AliAnalysisTaskLeeYangZeros::Terminate(Option_t *)
 307 {
 308   // Called once at the end of the query
 309
 310   const Int_t iNtheta = AliFlowLYZConstants::kTheta;
 311
 312   AliFlowAnalysisWithLeeYangZeros* fLyzTerm = new AliFlowAnalysisWithLeeYangZeros() ;
 313   fLyzTerm -> SetFirstRun(GetFirstRunLYZ());   //set first run true or false
 314   fLyzTerm -> SetUseSum(GetUseSumLYZ());       //set use sum true or false
 315
 316   fListHistos = (TList*)GetOutputData(0);
 317   //cout << "histogram list in Terminate" << endl;
 318
 319   if (fListHistos) {
 320
 321     //define histograms for first and second run
 322     AliFlowCommonHist *pCommonHist = NULL;
 323     AliFlowCommonHistResults *pCommonHistResults = NULL;
 324     TProfile* pHistProVtheta = NULL;
 325     TProfile* pHistProReDenom = NULL;
 326     TProfile* pHistProImDenom = NULL;
 327     TProfile* pHistProReDtheta = NULL;
 328     TProfile* pHistProImDtheta = NULL;
 329     TProfile* pHistProVeta = NULL;
 330     TProfile* pHistProVPt  = NULL;
 331     AliFlowLYZHist1 *pLYZHist1[iNtheta] = {NULL};      //array of pointers to AliFlowLYZHist1
 332     AliFlowLYZHist2 *pLYZHist2[iNtheta] = {NULL};      //array of pointers to AliFlowLYZHist2
 333
 334     if (GetFirstRunLYZ()) { //first run
 335       //Get the common histograms from the output list
 336       pCommonHist = dynamic_cast<AliFlowCommonHist*>
 337         (fListHistos->FindObject("AliFlowCommonHistLYZ1"));
 338       pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 339         (fListHistos->FindObject("AliFlowCommonHistResultsLYZ1"));
 340     }
 341     else { //second run
 342       //Get the common histograms from the output list
 343       pCommonHist = dynamic_cast<AliFlowCommonHist*>
 344         (fListHistos->FindObject("AliFlowCommonHistLYZ2"));
 345       pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 346         (fListHistos->FindObject("AliFlowCommonHistResultsLYZ2"));
 347     }
 348
 349     TProfile* pHistProR0theta = dynamic_cast<TProfile*>
 350       (fListHistos->FindObject("First_FlowPro_r0theta_LYZ"));
 351
 352     TH1F* pHistQsumforChi = dynamic_cast<TH1F*>
 353       (fListHistos->FindObject("Flow_QsumforChi_LYZ"));
 354
 355
 356     if (GetFirstRunLYZ()) { //for firstrun
 357       //Get the histograms from the output list
 358       for(Int_t theta = 0;theta<iNtheta;theta++){
 359         TString name = "AliFlowLYZHist1_";
 360         name += theta;
 361         pLYZHist1[theta] = dynamic_cast<AliFlowLYZHist1*>
 362           (fListHistos->FindObject(name));
 363       }
 364       pHistProVtheta = dynamic_cast<TProfile*>
 365           (fListHistos->FindObject("First_FlowPro_Vtheta_LYZ"));
 366
 367       //Set the histogram pointers and call Finish()
 368       if (pCommonHist && pCommonHistResults && pLYZHist1[0] &&
 369           pHistProVtheta && pHistProR0theta && pHistQsumforChi ) {
 370         fLyzTerm->SetCommonHists(pCommonHist);
 371         fLyzTerm->SetCommonHistsRes(pCommonHistResults);
 372         fLyzTerm->SetHist1(pLYZHist1);
 373         fLyzTerm->SetHistProVtheta(pHistProVtheta);
 374         fLyzTerm->SetHistProR0theta(pHistProR0theta);
 375         fLyzTerm->SetHistQsumforChi(pHistQsumforChi);
 376         fLyzTerm->Finish();
 377         PostData(0,fListHistos);
 378       } else {
 379         cout<<"WARNING: Histograms needed to run Finish() firstrun are not accessable!"<<endl;
 380       }
 381     } else { //for second run
 382       //Get the histograms from the output list
 383       for(Int_t theta = 0;theta<iNtheta;theta++){
 384         TString name = "AliFlowLYZHist2_";
 385         name += theta;
 386         pLYZHist2[theta] = dynamic_cast<AliFlowLYZHist2*>
 387           (fListHistos->FindObject(name));
 388       }
 389
 390       pHistProReDenom = dynamic_cast<TProfile*>
 391         (fListHistos->FindObject("Second_FlowPro_ReDenom_LYZ"));
 392       pHistProImDenom = dynamic_cast<TProfile*>
 393         (fListHistos->FindObject("Second_FlowPro_ImDenom_LYZ"));
 394
 395       pHistProReDtheta = dynamic_cast<TProfile*>
 396         (fListHistos->FindObject("Second_FlowPro_ReDtheta_LYZ"));
 397       pHistProImDtheta = dynamic_cast<TProfile*>
 398         (fListHistos->FindObject("Second_FlowPro_ImDtheta_LYZ"));
 399
 400       pHistProVeta = dynamic_cast<TProfile*>
 401         (fListHistos->FindObject("Second_FlowPro_Veta_LYZ"));
 402       pHistProVPt = dynamic_cast<TProfile*>
 403         (fListHistos->FindObject("Second_FlowPro_VPt_LYZ"));
 404
 405       //Set the histogram pointers and call Finish()
 406       if (pCommonHist && pCommonHistResults && pLYZHist2[0] && pHistProR0theta &&
 407           pHistProReDenom && pHistProImDenom && pHistProVeta && pHistProVPt) {
 408         fLyzTerm->SetCommonHists(pCommonHist);
 409         fLyzTerm->SetCommonHistsRes(pCommonHistResults);
 410         fLyzTerm->SetHist2(pLYZHist2);
 411         fLyzTerm->SetHistProR0theta(pHistProR0theta);
 412         fLyzTerm->SetHistProReDenom(pHistProReDenom);
 413         fLyzTerm->SetHistProImDenom(pHistProImDenom);
 414         fLyzTerm->SetHistProReDtheta(pHistProReDtheta);
 415         fLyzTerm->SetHistProImDtheta(pHistProImDtheta);
 416         fLyzTerm->SetHistProVeta(pHistProVeta);
 417         fLyzTerm->SetHistProVPt(pHistProVPt);
 418         fLyzTerm->SetHistQsumforChi(pHistQsumforChi);
 419         fLyzTerm->Finish();
 420         PostData(0,fListHistos);
 421       } else {
 422         cout<<"WARNING: Histograms needed to run Finish() secondrun are not accessable!"<<endl;
 423       }
 424     }
 425
 426     //    fListHistos->Print();
 427   }
 428   else { cout << "histogram list pointer is empty" << endl;}
 429
 430   cout<<".....finished LYZ"<<endl;
 431 }