PWG2/FLOW/AliFlowTasks/AliAnalysisTaskFlowEvent.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 // AliAnalysisTaskFlowEvent:
  18 //
  19 // analysis task for filling the flow event
  20 // from MCEvent, ESD, AOD ....
  21 // and put it in an output stream so it can
  22 // be used by the various flow analysis methods
  23 // for cuts the correction framework is used
  24 // which also outputs QA histograms to view
  25 // the effects of the cuts
  26 ////////////////////////////////////////////////////
  27
  28 #include "Riostream.h" //needed as include
  29 #include "TChain.h"
  30 #include "TTree.h"
  31 #include "TFile.h" //needed as include
  32 #include "TList.h"
  33 #include "TRandom3.h"
  34 #include "TTimeStamp.h"
  35
  36 // ALICE Analysis Framework
  37 class AliAnalysisTask;
  38 #include "AliAnalysisManager.h"
  39
  40 // ESD interface
  41 #include "AliESDEvent.h"
  42 #include "AliESDInputHandler.h"
  43
  44 // AOD interface
  45 #include "AliAODEvent.h"
  46 #include "AliAODInputHandler.h"
  47
  48 // Monte Carlo Event
  49 #include "AliMCEventHandler.h"
  50 #include "AliMCEvent.h"
  51
  52 // ALICE Correction Framework
  53 #include "AliCFManager.h"
  54
  55 // Interface to Event generators to get Reaction Plane Angle
  56 #include "AliGenCocktailEventHeader.h"
  57 #include "AliGenHijingEventHeader.h"
  58 #include "AliGenGeVSimEventHeader.h"
  59 #include "AliGenEposEventHeader.h"
  60
  61 // Interface to make the Flow Event Simple used in the flow analysis methods
  62 #include "AliFlowEventSimpleMaker.h"
  63
  64 #include "AliAnalysisTaskFlowEvent.h"
  65
  66
  67 ClassImp(AliAnalysisTaskFlowEvent)
  68
  69 //________________________________________________________________________
  70 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on) :
  71   AliAnalysisTask(name, ""),
  72 //  fOutputFile(NULL),
  73   fESD(NULL),
  74   fAOD(NULL),
  75   fEventMaker(NULL),
  76   fAnalysisType("ESD"),
  77   fCFManager1(NULL),
  78   fCFManager2(NULL),
  79   fQAInt(NULL),
  80   fQADiff(NULL),
  81   fMinMult(0),
  82   fMaxMult(10000000),
  83   fQA(on),
  84   fMCReactionPlaneAngle(0.),
  85   fCount(0),
  86   fNoOfLoops(1),
  87   fEllipticFlowValue(0.),
  88   fSigmaEllipticFlowValue(0.),
  89   fMultiplicityOfEvent(1000000000),
  90   fSigmaMultiplicityOfEvent(0),
  91   fMyTRandom3(NULL)
  92 {
  93   // Constructor
  94   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name)"<<endl;
  95
  96   // Define input and output slots here
  97   // Input slot #0 works with a TChain
  98   DefineInput(0, TChain::Class());
  99   // Define here the flow event output
 100   DefineOutput(0, AliFlowEventSimple::Class());
 101   if(on) {
 102     DefineOutput(1, TList::Class());
 103     DefineOutput(2, TList::Class()); }
 104   // and for testing open an output file
 105   //  fOutputFile = new TFile("FlowEvents.root","RECREATE");
 106
 107 }
 108
 109 //________________________________________________________________________
 110 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent() :
 111   //  fOutputFile(NULL),
 112   fESD(NULL),
 113   fAOD(NULL),
 114   fEventMaker(NULL),
 115   fAnalysisType("ESD"),
 116   fCFManager1(NULL),
 117   fCFManager2(NULL),
 118   fQAInt(NULL),
 119   fQADiff(NULL),
 120   fMinMult(0),
 121   fMaxMult(10000000),
 122   fQA(kFALSE),
 123   fMCReactionPlaneAngle(0.),
 124   fCount(0),
 125   fNoOfLoops(1),
 126   fEllipticFlowValue(0.),
 127   fSigmaEllipticFlowValue(0.),
 128   fMultiplicityOfEvent(1000000000),
 129   fSigmaMultiplicityOfEvent(0),
 130   fMyTRandom3(NULL)
 131 {
 132   // Constructor
 133   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()"<<endl;
 134 }
 135
 136 //________________________________________________________________________
 137 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed) :
 138   AliAnalysisTask(name, ""),
 139 //  fOutputFile(NULL),
 140   fESD(NULL),
 141   fAOD(NULL),
 142   fEventMaker(NULL),
 143   fAnalysisType("ESD"),
 144   fCFManager1(NULL),
 145   fCFManager2(NULL),
 146   fQAInt(NULL),
 147   fQADiff(NULL),
 148   fMinMult(0),
 149   fMaxMult(10000000),
 150   fQA(on),
 151   fMCReactionPlaneAngle(0.),
 152   fCount(0),
 153   fNoOfLoops(1),
 154   fEllipticFlowValue(0.),
 155   fSigmaEllipticFlowValue(0.),
 156   fMultiplicityOfEvent(1000000000),
 157   fSigmaMultiplicityOfEvent(0),
 158   fMyTRandom3(NULL)
 159 {
 160   // Constructor
 161   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)"<<endl;
 162
 163   fMyTRandom3 = new TRandom3(iseed);
 164   gRandom->SetSeed(fMyTRandom3->Integer(65539));
 165
 166   // Define input and output slots here
 167   // Input slot #0 works with a TChain
 168   DefineInput(0, TChain::Class());
 169   // Define here the flow event output
 170   DefineOutput(0, AliFlowEventSimple::Class());
 171   if(on) {
 172     DefineOutput(1, TList::Class());
 173     DefineOutput(2, TList::Class()); }
 174   // and for testing open an output file
 175   //  fOutputFile = new TFile("FlowEvents.root","RECREATE");
 176
 177 }
 178
 179 //________________________________________________________________________
 180 AliAnalysisTaskFlowEvent::~AliAnalysisTaskFlowEvent()
 181 {
 182   //
 183   // Destructor
 184   //
 185   if (fMyTRandom3) delete fMyTRandom3;
 186   // objects in the output list are deleted
 187   // by the TSelector dtor (I hope)
 188
 189 }
 190
 191 //________________________________________________________________________
 192 void AliAnalysisTaskFlowEvent::ConnectInputData(Option_t *)
 193 {
 194   // Connect ESD or AOD here
 195   // Called once
 196   cout<<"AliAnalysisTaskFlowEvent::ConnectInputData(Option_t *)"<<endl;
 197
 198   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 199   if (!tree) {
 200     Printf("ERROR: Could not read chain from input slot 0");
 201   }
 202   else {
 203     if (fAnalysisType == "ESD" || fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1" || fAnalysisType == "MC") {
 204       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 205       if (!esdH) {
 206         Printf("ERROR: Could not get ESDInputHandler");
 207       }
 208       else {fESD = esdH->GetEvent();}
 209     }
 210     else if (fAnalysisType == "AOD") {
 211       AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 212       if (!aodH) {
 213         Printf("ERROR: Could not get AODInputHandler");
 214       }
 215       else { fAOD = aodH->GetEvent();}
 216     }
 217     else {
 218       Printf("!!!!!Wrong analysis type: Only ESD, ESDMC0, ESDMC1, AOD and MC types are allowed!");
 219       exit(1);
 220     }
 221   }
 222 }
 223
 224 //________________________________________________________________________
 225 void AliAnalysisTaskFlowEvent::CreateOutputObjects()
 226 {
 227   // Called at every worker node to initialize
 228   cout<<"AliAnalysisTaskFlowEvent::CreateOutputObjects()"<<endl;
 229
 230   if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMC0"  || fAnalysisType == "ESDMC1" || fAnalysisType == "MC")) {
 231     cout<<"WRONG ANALYSIS TYPE! only ESD, ESDMC0, ESDMC1, AOD and MC are allowed."<<endl;
 232     exit(1);
 233   }
 234
 235   // Flow Event maker
 236   fEventMaker = new AliFlowEventSimpleMaker();
 237 }
 238
 239 //________________________________________________________________________
 240 void AliAnalysisTaskFlowEvent::Exec(Option_t *)
 241 {
 242   // Main loop
 243   // Called for each event
 244   AliFlowEventSimple* fEvent = NULL;
 245   Double_t fRP = 0.; // the monte carlo reaction plane angle
 246   AliMCEvent* mcEvent = NULL;
 247   // See if we can get Monte Carlo Information and if so get the reaction plane
 248
 249   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 250   if (eventHandler) {
 251     mcEvent = eventHandler->MCEvent();
 252     if (mcEvent) {
 253
 254       //COCKTAIL with HIJING
 255       if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Cocktail Header")) { //returns 0 if matches
 256         AliGenCocktailEventHeader *headerC = dynamic_cast<AliGenCocktailEventHeader *> (mcEvent-> GenEventHeader());
 257         if (headerC) {
 258           TList *lhd = headerC->GetHeaders();
 259           if (lhd) {
 260             AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(0));
 261             if (hdh) {
 262               fRP = hdh->ReactionPlaneAngle();
 263               //cout<<"The reactionPlane from Hijing is: "<< fRP <<endl;
 264             }
 265           }
 266         }
 267         //else { cout<<"headerC is NULL"<<endl; }
 268       }
 269
 270       //GEVSIM
 271       else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"GeVSim header")) { //returns 0 if matches
 272         AliGenGeVSimEventHeader* headerG = (AliGenGeVSimEventHeader*)(mcEvent->GenEventHeader());
 273         if (headerG) {
 274           fRP = headerG->GetEventPlane();
 275           //cout<<"The reactionPlane from GeVSim is: "<< fRP <<endl;
 276         }
 277         //else { cout<<"headerG is NULL"<<endl; }
 278       }
 279
 280       //HIJING
 281       else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Hijing")) { //returns 0 if matches
 282         AliGenHijingEventHeader* headerH = (AliGenHijingEventHeader*)(mcEvent->GenEventHeader());
 283         if (headerH) {
 284           fRP = headerH->ReactionPlaneAngle();
 285           //cout<<"The reactionPlane from Hijing is: "<< fRP <<endl;
 286         }
 287         //else { cout<<"headerH is NULL"<<endl; }
 288       }
 289
 290       //EPOS
 291       else if (!strcmp(mcEvent->GenEventHeader()->GetName(),"EPOS")) {
 292         AliGenEposEventHeader* headerE = (AliGenEposEventHeader*)(mcEvent->GenEventHeader());
 293         if (headerE) {
 294           fRP = headerE->ReactionPlaneAngle();
 295           //cout<<"The reactionPlane from EPOS is: "<< fR <<endl;
 296         }
 297         //else { cout<<"headerE is NULL"<<endl; }
 298       }
 299
 300     }
 301     //else {cout<<"No MC event!"<<endl; }
 302
 303   }
 304   //else {cout<<"No eventHandler!"<<endl; }
 305
 306
 307   fEventMaker->SetMCReactionPlaneAngle(fRP);
 308   //setting event cuts
 309   fEventMaker->SetMinMult(fMinMult);
 310   fEventMaker->SetMaxMult(fMaxMult);
 311
 312   if (fEllipticFlowValue != 0.) {
 313     // set the value of the monte carlo event plane for the flow event
 314     cout << "settings for afterburner in TaskFlowEvent.cxx:" << endl;
 315     cout << "fCount = " << fCount << endl;
 316     cout << "fNoOfLoops = " << fNoOfLoops << endl;
 317     cout << "fEllipticFlowValue = " << fEllipticFlowValue << endl;
 318     cout << "fSigmaEllipticFlowValue = " << fSigmaEllipticFlowValue << endl;
 319     cout << "fMultiplicityOfEvent = " << fMultiplicityOfEvent << endl;
 320     cout << "fSigmaMultiplicityOfEvent = " << fSigmaMultiplicityOfEvent << endl;
 321
 322     Double_t xRPangle=0.;
 323     Double_t xNewFlowValue = 0.;
 324     Int_t nNewMultOfEvent = 100000000;
 325
 326     if (fMyTRandom3) {
 327       xRPangle = TMath::TwoPi()*(fMyTRandom3->Rndm());
 328       xNewFlowValue = fMyTRandom3->Gaus(fEllipticFlowValue,fSigmaEllipticFlowValue);
 329       nNewMultOfEvent =  TMath::Nint(fMyTRandom3->Gaus(fMultiplicityOfEvent,fSigmaMultiplicityOfEvent));
 330     }
 331     else {
 332       cout << "no random generator pointer initialized " << endl;
 333     }
 334     cout << "xRPangle = " << xRPangle << endl;
 335     cout << "xNewFlowValue = " << xNewFlowValue << endl;
 336     cout << "nNewMultOfEvent = " << nNewMultOfEvent << endl;
 337     cout << "settings for after burner" << endl;
 338
 339     fEventMaker->SetMCReactionPlaneAngle(xRPangle);
 340     fEventMaker->SetNoOfLoops(fNoOfLoops);
 341     fEventMaker->SetEllipticFlowValue(xNewFlowValue);
 342     fEventMaker->SetMultiplicityOfEvent(nNewMultOfEvent);
 343     //end settings afterburner
 344   }
 345
 346
 347   // Fill the FlowEventSimple for MC input
 348   if (fAnalysisType == "MC") {
 349     if (!fCFManager1) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 350     if (!fCFManager2) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 351
 352     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 353     // This handler can return the current MC event
 354     if (!mcEvent) { Printf("ERROR: Could not retrieve MC event"); return;}
 355
 356     fCFManager1->SetMCEventInfo(mcEvent);
 357     fCFManager2->SetMCEventInfo(mcEvent);
 358
 359     // analysis
 360     Printf("Number of MC particles: %d", mcEvent->GetNumberOfTracks());
 361     fEvent = fEventMaker->FillTracks(mcEvent,fCFManager1,fCFManager2);
 362     // here we have the fEvent and want to make it available as an output stream
 363     // so no delete fEvent;
 364   }
 365   // Fill the FlowEventSimple for ESD input
 366   else if (fAnalysisType == "ESD") {
 367     if (!fCFManager1) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 368     if (!fCFManager2) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 369
 370     if (!fESD) { Printf("ERROR: fESD not available"); return;}
 371     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 372
 373     // analysis
 374     fEvent = fEventMaker->FillTracks(fESD,fCFManager1,fCFManager2);
 375   }
 376   // Fill the FlowEventSimple for ESD input combined with MC info
 377   else if (fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1" ) {
 378     if (!fCFManager1) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 379     if (!fCFManager2) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 380     if (!fESD) { Printf("ERROR: fESD not available"); return;}
 381     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 382
 383     if (!mcEvent) {Printf("ERROR: Could not retrieve MC event"); return;}
 384
 385     fCFManager1->SetMCEventInfo(mcEvent);
 386     fCFManager2->SetMCEventInfo(mcEvent);
 387
 388
 389     if (fAnalysisType == "ESDMC0") {
 390       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 0); //0 = kine from ESD, 1 = kine from MC
 391     } else if (fAnalysisType == "ESDMC1") {
 392       fEvent = fEventMaker->FillTracks(fESD, mcEvent, fCFManager1, fCFManager2, 1); //0 = kine from ESD, 1 = kine from MC
 393     }
 394   }
 395   // Fill the FlowEventSimple for AOD input
 396   else if (fAnalysisType == "AOD") {
 397     if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
 398     Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
 399
 400     // analysis
 401     //For the moment don't use CF //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fAOD,fCFManager1,fCFManager2);
 402     fEvent = fEventMaker->FillTracks(fAOD);
 403   }
 404
 405   //fListHistos->Print();
 406   //  fOutputFile->WriteObject(fEvent,"myFlowEventSimple");
 407   PostData(0,fEvent);
 408   if (fQA) {
 409     PostData(1,fQAInt);
 410     PostData(2,fQADiff); }
 411 }
 412
 413 //________________________________________________________________________
 414 void AliAnalysisTaskFlowEvent::Terminate(Option_t *)
 415 {
 416   // Called once at the end of the query -- do not call in case of CAF
 417
 418 }
 419
 420