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() :
  71   AliAnalysisTaskSE(),
  72   //  fOutputFile(NULL),
  73   fEventMaker(NULL),
  74   fAnalysisType("ESD"),
  75   fCFManager1(NULL),
  76   fCFManager2(NULL),
  77   fQAInt(NULL),
  78   fQADiff(NULL),
  79   fMinMult(0),
  80   fMaxMult(10000000),
  81   fMinA(-1.0),
  82   fMaxA(-0.01),
  83   fMinB(0.01),
  84   fMaxB(1.0),
  85   fQA(kFALSE),
  86   fMCReactionPlaneAngle(0.),
  87   fCount(0),
  88   fNoOfLoops(1),
  89   fEllipticFlowValue(0.),
  90   fSigmaEllipticFlowValue(0.),
  91   fMultiplicityOfEvent(1000000000),
  92   fSigmaMultiplicityOfEvent(0),
  93   fMyTRandom3(NULL),
  94   fbAfterburnerOn(kFALSE)
  95 {
  96   // Constructor
  97   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()"<<endl;
  98 }
  99
 100 //________________________________________________________________________
 101 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed) :
 102   AliAnalysisTaskSE(name),
 103 //  fOutputFile(NULL),
 104   fEventMaker(NULL),
 105   fAnalysisType("ESD"),
 106   fCFManager1(NULL),
 107   fCFManager2(NULL),
 108   fQAInt(NULL),
 109   fQADiff(NULL),
 110   fMinMult(0),
 111   fMaxMult(10000000),
 112   fMinA(-1.0),
 113   fMaxA(-0.01),
 114   fMinB(0.01),
 115   fMaxB(1.0),
 116   fQA(on),
 117   fMCReactionPlaneAngle(0.),
 118   fCount(0),
 119   fNoOfLoops(1),
 120   fEllipticFlowValue(0.),
 121   fSigmaEllipticFlowValue(0.),
 122   fMultiplicityOfEvent(1000000000),
 123   fSigmaMultiplicityOfEvent(0),
 124   fMyTRandom3(NULL),
 125   fbAfterburnerOn(kFALSE)
 126 {
 127   // Constructor
 128   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)"<<endl;
 129   fMyTRandom3 = new TRandom3(iseed);
 130   gRandom->SetSeed(fMyTRandom3->Integer(65539));
 131
 132
 133   // Define output slots here
 134   // Define here the flow event output
 135   DefineOutput(1, AliFlowEventSimple::Class());
 136   if(on) {
 137     DefineOutput(2, TList::Class());
 138     DefineOutput(3, TList::Class()); }
 139   // and for testing open an output file
 140   //  fOutputFile = new TFile("FlowEvents.root","RECREATE");
 141
 142 }
 143
 144 //________________________________________________________________________
 145 AliAnalysisTaskFlowEvent::~AliAnalysisTaskFlowEvent()
 146 {
 147   //
 148   // Destructor
 149   //
 150   if (fMyTRandom3) delete fMyTRandom3;
 151   // objects in the output list are deleted
 152   // by the TSelector dtor (I hope)
 153
 154 }
 155
 156 //________________________________________________________________________
 157 void AliAnalysisTaskFlowEvent::UserCreateOutputObjects()
 158 {
 159   // Called at every worker node to initialize
 160   cout<<"AliAnalysisTaskFlowEvent::CreateOutputObjects()"<<endl;
 161
 162   if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMC0"  || fAnalysisType == "ESDMC1" || fAnalysisType == "MC")) {
 163     cout<<"WRONG ANALYSIS TYPE! only ESD, ESDMC0, ESDMC1, AOD and MC are allowed."<<endl;
 164     exit(1);
 165   }
 166
 167   // Flow Event maker
 168   fEventMaker = new AliFlowEventSimpleMaker();
 169 }
 170
 171 //________________________________________________________________________
 172 void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
 173 {
 174   // Main loop
 175   // Called for each event
 176   AliFlowEventSimple* fEvent = NULL;
 177   Double_t fRP = 0.; // the monte carlo reaction plane angle
 178   AliMCEvent* mcEvent = fMCEvent;
 179   // See if we can get Monte Carlo Information and if so get the reaction plane
 180
 181   if (mcEvent) {
 182
 183     //COCKTAIL with HIJING
 184     if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Cocktail Header")) { //returns 0 if matches
 185       AliGenCocktailEventHeader *headerC = dynamic_cast<AliGenCocktailEventHeader *> (mcEvent-> GenEventHeader());
 186       if (headerC) {
 187         TList *lhd = headerC->GetHeaders();
 188         if (lhd) {
 189           AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(0));
 190           if (hdh) {
 191             fRP = hdh->ReactionPlaneAngle();
 192             //cout<<"The reactionPlane from Hijing is: "<< fRP <<endl;
 193           }
 194         }
 195       }
 196       //else { cout<<"headerC is NULL"<<endl; }
 197     }
 198
 199     //GEVSIM
 200     else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"GeVSim header")) { //returns 0 if matches
 201       AliGenGeVSimEventHeader* headerG = (AliGenGeVSimEventHeader*)(mcEvent->GenEventHeader());
 202       if (headerG) {
 203         fRP = headerG->GetEventPlane();
 204         //cout<<"The reactionPlane from GeVSim is: "<< fRP <<endl;
 205       }
 206       //else { cout<<"headerG is NULL"<<endl; }
 207     }
 208
 209     //HIJING
 210     else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Hijing")) { //returns 0 if matches
 211       AliGenHijingEventHeader* headerH = (AliGenHijingEventHeader*)(mcEvent->GenEventHeader());
 212       if (headerH) {
 213         fRP = headerH->ReactionPlaneAngle();
 214         //cout<<"The reactionPlane from Hijing is: "<< fRP <<endl;
 215       }
 216       //else { cout<<"headerH is NULL"<<endl; }
 217     }
 218
 219     //EPOS
 220     else if (!strcmp(mcEvent->GenEventHeader()->GetName(),"EPOS")) {
 221       AliGenEposEventHeader* headerE = (AliGenEposEventHeader*)(mcEvent->GenEventHeader());
 222       if (headerE) {
 223         fRP = headerE->ReactionPlaneAngle();
 224         //cout<<"The reactionPlane from EPOS is: "<< fR <<endl;
 225       }
 226       //else { cout<<"headerE is NULL"<<endl; }
 227     }
 228
 229
 230   }
 231   //else {cout<<"No eventHandler!"<<endl; }
 232
 233
 234   fEventMaker->SetMCReactionPlaneAngle(fRP);
 235   //setting event cuts
 236   fEventMaker->SetMinMult(fMinMult);
 237   fEventMaker->SetMaxMult(fMaxMult);
 238   //setting ranges for eta subevents
 239   fEventMaker->SetSubeventEtaRange(fMinA,fMaxA,fMinB,fMaxB);
 240
 241   if (fEllipticFlowValue != 0.) {
 242     // set the value of the monte carlo event plane for the flow event
 243     cout << "settings for afterburner in TaskFlowEvent.cxx:" << endl;
 244     cout << "fCount = " << fCount << endl;
 245     cout << "fNoOfLoops = " << fNoOfLoops << endl;
 246     cout << "fEllipticFlowValue = " << fEllipticFlowValue << endl;
 247     cout << "fSigmaEllipticFlowValue = " << fSigmaEllipticFlowValue << endl;
 248     cout << "fMultiplicityOfEvent = " << fMultiplicityOfEvent << endl;
 249     cout << "fSigmaMultiplicityOfEvent = " << fSigmaMultiplicityOfEvent << endl;
 250
 251     Double_t xRPangle=0.;
 252     Double_t xNewFlowValue = 0.;
 253     Int_t nNewMultOfEvent = 100000000;
 254
 255     if (fbAfterburnerOn && fMyTRandom3) {
 256       xRPangle = TMath::TwoPi()*(fMyTRandom3->Rndm());
 257       xNewFlowValue = fMyTRandom3->Gaus(fEllipticFlowValue,fSigmaEllipticFlowValue);
 258       nNewMultOfEvent =  TMath::Nint(fMyTRandom3->Gaus(fMultiplicityOfEvent,fSigmaMultiplicityOfEvent));
 259     }
 260     else {
 261       cout << "no random generator pointer initialized " << endl;
 262     }
 263     cout << "xRPangle = " << xRPangle << endl;
 264     cout << "xNewFlowValue = " << xNewFlowValue << endl;
 265     cout << "nNewMultOfEvent = " << nNewMultOfEvent << endl;
 266     cout << "settings for after burner" << endl;
 267
 268     fEventMaker->SetMCReactionPlaneAngle(xRPangle);
 269     fEventMaker->SetNoOfLoops(fNoOfLoops);
 270     fEventMaker->SetEllipticFlowValue(xNewFlowValue);
 271     fEventMaker->SetMultiplicityOfEvent(nNewMultOfEvent);
 272     //end settings afterburner
 273   }
 274
 275
 276   // Fill the FlowEventSimple for MC input
 277   if (fAnalysisType == "MC") {
 278     if (!fCFManager1) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 279     if (!fCFManager2) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 280
 281     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 282     // This handler can return the current MC event
 283     if (!mcEvent) { Printf("ERROR: Could not retrieve MC event"); return;}
 284
 285     fCFManager1->SetMCEventInfo(mcEvent);
 286     fCFManager2->SetMCEventInfo(mcEvent);
 287
 288     // analysis
 289     Printf("Number of MC particles: %d", mcEvent->GetNumberOfTracks());
 290     fEvent = fEventMaker->FillTracks(mcEvent,fCFManager1,fCFManager2);
 291     // here we have the fEvent and want to make it available as an output stream
 292     // so no delete fEvent;
 293   }
 294   // Fill the FlowEventSimple for ESD input
 295   else if (fAnalysisType == "ESD") {
 296     if (!fCFManager1) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 297     if (!fCFManager2) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 298
 299     if (!fInputEvent) { Printf("ERROR: fInputEvent not available"); return;}
 300     //check the offline trigger (check if the event has the correct trigger)
 301     if (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected())
 302       {
 303         Printf("There are %d tracks in this event", fInputEvent->GetNumberOfTracks());
 304
 305         // analysis
 306         fEvent = fEventMaker->FillTracks(dynamic_cast <AliESDEvent*> (fInputEvent),fCFManager1,fCFManager2);
 307       }
 308   }
 309   // Fill the FlowEventSimple for ESD input combined with MC info
 310   else if (fAnalysisType == "ESDMC0" || fAnalysisType == "ESDMC1" ) {
 311     if (!fCFManager1) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 312     if (!fCFManager2) {cout << "ERROR: No pointer to correction framework cuts! " << endl; return; }
 313     if (!fInputEvent) { Printf("ERROR: fInputEvent not available"); return;}
 314     Printf("There are %d tracks in this event", fInputEvent->GetNumberOfTracks());
 315
 316     if (!mcEvent) {Printf("ERROR: Could not retrieve MC event"); return;}
 317
 318     fCFManager1->SetMCEventInfo(mcEvent);
 319     fCFManager2->SetMCEventInfo(mcEvent);
 320
 321
 322     if (fAnalysisType == "ESDMC0") {
 323       fEvent = fEventMaker->FillTracks(dynamic_cast <AliESDEvent*> (fInputEvent), mcEvent, fCFManager1, fCFManager2, 0); //0 = kine from ESD, 1 = kine from MC
 324     } else if (fAnalysisType == "ESDMC1") {
 325       fEvent = fEventMaker->FillTracks(dynamic_cast <AliESDEvent*> (fInputEvent), mcEvent, fCFManager1, fCFManager2, 1); //0 = kine from ESD, 1 = kine from MC
 326     }
 327   }
 328   // Fill the FlowEventSimple for AOD input
 329   else if (fAnalysisType == "AOD") {
 330     if (!fOutputAOD) {Printf("ERROR: fOutputAOD not available"); return;}
 331     Printf("There are %d tracks in this event", fOutputAOD->GetNumberOfTracks());
 332
 333     // analysis
 334     //For the moment don't use CF //AliFlowEventSimple* fEvent = fEventMaker->FillTracks(fOutputAOD,fCFManager1,fCFManager2);
 335     fEvent = fEventMaker->FillTracks(fOutputAOD);
 336   }
 337
 338   //fListHistos->Print();
 339   //  fOutputFile->WriteObject(fEvent,"myFlowEventSimple");
 340   PostData(1,fEvent);
 341   if (fQA) {
 342     PostData(2,fQAInt);
 343     PostData(3,fQADiff); }
 344 }
 345
 346 //________________________________________________________________________
 347 void AliAnalysisTaskFlowEvent::Terminate(Option_t *)
 348 {
 349   // Called once at the end of the query -- do not call in case of CAF
 350
 351 }
 352
 353