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 "TH2F.h"
  34 #include "TRandom3.h"
  35 #include "TTimeStamp.h"
  36
  37 // ALICE Analysis Framework
  38 #include "AliAnalysisManager.h"
  39 #include "AliAnalysisTaskSE.h"
  40 #include "AliESDtrackCuts.h"
  41
  42 // ESD interface
  43 #include "AliESDEvent.h"
  44 #include "AliESDInputHandler.h"
  45
  46 // AOD interface
  47 #include "AliAODEvent.h"
  48 #include "AliAODInputHandler.h"
  49
  50 // Monte Carlo Event
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCEvent.h"
  53
  54 // ALICE Correction Framework
  55 #include "AliCFManager.h"
  56
  57 // Interface to Event generators to get Reaction Plane Angle
  58 #include "AliGenCocktailEventHeader.h"
  59 #include "AliGenHijingEventHeader.h"
  60 #include "AliGenGeVSimEventHeader.h"
  61 #include "AliGenEposEventHeader.h"
  62
  63 // Interface to make the Flow Event Simple used in the flow analysis methods
  64 #include "AliFlowEvent.h"
  65
  66 #include "AliAnalysisTaskFlowEvent.h"
  67
  68 #include "AliLog.h"
  69
  70 ClassImp(AliAnalysisTaskFlowEvent)
  71
  72 //________________________________________________________________________
  73 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent() :
  74   AliAnalysisTaskSE(),
  75   //  fOutputFile(NULL),
  76   fAnalysisType("ESD"),
  77   fRPType("Global"),
  78   fCFManager1(NULL),
  79   fCFManager2(NULL),
  80   fQAInt(NULL),
  81   fQADiff(NULL),
  82   fMinMult(0),
  83   fMaxMult(10000000),
  84   fMinA(-1.0),
  85   fMaxA(-0.01),
  86   fMinB(0.01),
  87   fMaxB(1.0),
  88   fQA(kFALSE),
  89   fMCReactionPlaneAngle(0.),
  90   fCount(0),
  91   fNoOfLoops(1),
  92   fEllipticFlowValue(0.),
  93   fSigmaEllipticFlowValue(0.),
  94   fMultiplicityOfEvent(1000000000),
  95   fSigmaMultiplicityOfEvent(0),
  96   fMyTRandom3(NULL),
  97   fbAfterburnerOn(kFALSE)
  98 {
  99   // Constructor
 100   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()"<<endl;
 101 }
 102
 103 //________________________________________________________________________
 104 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, TString RPtype, Bool_t on, UInt_t iseed) :
 105   AliAnalysisTaskSE(name),
 106   //  fOutputFile(NULL),
 107   fAnalysisType("ESD"),
 108   fRPType(RPtype),
 109   fCFManager1(NULL),
 110   fCFManager2(NULL),
 111   fQAInt(NULL),
 112   fQADiff(NULL),
 113   fMinMult(0),
 114   fMaxMult(10000000),
 115   fMinA(-1.0),
 116   fMaxA(-0.01),
 117   fMinB(0.01),
 118   fMaxB(1.0),
 119   fQA(on),
 120   fMCReactionPlaneAngle(0.),
 121   fCount(0),
 122   fNoOfLoops(1),
 123   fEllipticFlowValue(0.),
 124   fSigmaEllipticFlowValue(0.),
 125   fMultiplicityOfEvent(1000000000),
 126   fSigmaMultiplicityOfEvent(0),
 127   fMyTRandom3(NULL),
 128   fbAfterburnerOn(kFALSE)
 129 {
 130   // Constructor
 131   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)"<<endl;
 132   fMyTRandom3 = new TRandom3(iseed);
 133   gRandom->SetSeed(fMyTRandom3->Integer(65539));
 134
 135   //FMD input slot
 136   if (strcmp(RPtype,"FMD")==0) {
 137     DefineInput(1, TList::Class());
 138   }
 139
 140   // Define output slots here
 141   // Define here the flow event output
 142   DefineOutput(1, AliFlowEventSimple::Class());
 143   if(on)
 144   {
 145     DefineOutput(2, TList::Class());
 146     DefineOutput(3, TList::Class());
 147   }
 148   // and for testing open an output file
 149   //  fOutputFile = new TFile("FlowEvents.root","RECREATE");
 150
 151 }
 152
 153 //________________________________________________________________________
 154 AliAnalysisTaskFlowEvent::~AliAnalysisTaskFlowEvent()
 155 {
 156   //
 157   // Destructor
 158   //
 159   if (fMyTRandom3) delete fMyTRandom3;
 160   // objects in the output list are deleted
 161   // by the TSelector dtor (I hope)
 162
 163 }
 164
 165 //________________________________________________________________________
 166 void AliAnalysisTaskFlowEvent::UserCreateOutputObjects()
 167 {
 168   // Called at every worker node to initialize
 169   cout<<"AliAnalysisTaskFlowEvent::CreateOutputObjects()"<<endl;
 170
 171   if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD"  || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC"))
 172   {
 173     AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD and MC are allowed.");
 174     exit(1);
 175   }
 176 }
 177
 178 //________________________________________________________________________
 179 void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
 180 {
 181   // Main loop
 182   // Called for each event
 183   AliFlowEvent* flowEvent = NULL;
 184   AliMCEvent*  mcEvent = MCEvent();                              // from TaskSE
 185   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
 186   AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
 187   AliMultiplicity* myTracklets = NULL;
 188   TH2F* histFMD = NULL;
 189
 190   if(GetNinputs()==2) {
 191     TList* FMDdata = dynamic_cast<TList*>(GetInputData(1));
 192     if(!FMDdata) {
 193       cout<<" No FMDdata "<<endl;
 194       exit(2);
 195     }
 196     histFMD = dynamic_cast<TH2F*>(FMDdata->FindObject("dNdetadphiHistogramTrVtx"));
 197     if (!histFMD) {
 198       cout<< "No histFMD"<<endl;
 199       exit(2);
 200     }
 201   }
 202
 203
 204   // Make the FlowEvent for MC input
 205   if (fAnalysisType == "MC")
 206   {
 207     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 208     // This handler can return the current MC event
 209     if (!(fCFManager1&&fCFManager2))
 210     {
 211       AliError("ERROR: No pointer to correction framework cuts! ");
 212       return;
 213     }
 214     if (!mcEvent)
 215     {
 216       AliError("ERROR: Could not retrieve MC event");
 217       return;
 218     }
 219
 220     fCFManager1->SetMCEventInfo(mcEvent);
 221     fCFManager2->SetMCEventInfo(mcEvent);
 222
 223     AliInfo(Form("Number of MC particles: %d", mcEvent->GetNumberOfTracks()));
 224
 225     //check multiplicity
 226     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtGenCuts,mcEvent))
 227     {
 228       AliWarning("Event does not pass multiplicity cuts"); return;
 229     }
 230     //make event
 231     flowEvent = new AliFlowEvent(mcEvent,fCFManager1,fCFManager2);
 232   }
 233
 234   // Make the FlowEvent for ESD input
 235   else if (fAnalysisType == "ESD")
 236   {
 237     if (!(fCFManager1&&fCFManager2))
 238     {
 239       AliError("ERROR: No pointer to correction framework cuts!");
 240       return;
 241     }
 242     if (!myESD)
 243     {
 244       AliError("ERROR: ESD not available");
 245       return;
 246     }
 247
 248     //check the offline trigger (check if the event has the correct trigger)
 249     AliInfo(Form("ESD has %d tracks", fInputEvent->GetNumberOfTracks()));
 250
 251     //check multiplicity
 252     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
 253     {
 254       AliWarning("Event does not pass multiplicity cuts"); return;
 255     }
 256
 257     //make event
 258     if (fRPType == "Global") {
 259       flowEvent = new AliFlowEvent(myESD,fCFManager1,fCFManager2);
 260     }
 261     else if (fRPType == "Tracklet"){
 262       flowEvent = new AliFlowEvent(myESD,myTracklets,fCFManager2);
 263     }
 264     else if (fRPType == "FMD"){
 265       flowEvent = new AliFlowEvent(myESD,histFMD,fCFManager2);
 266     }
 267
 268     // if monte carlo event get reaction plane from monte carlo (depends on generator)
 269     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
 270     //set reference multiplicity, TODO: maybe move it to the constructor?
 271     flowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
 272   }
 273
 274   // Make the FlowEvent for ESD input combined with MC info
 275   else if (fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" )
 276   {
 277     if (!(fCFManager1&&fCFManager2))
 278     {
 279       AliError("ERROR: No pointer to correction framework cuts! ");
 280       return;
 281     }
 282     if (!myESD)
 283     {
 284       AliError("ERROR: ESD not available");
 285       return;
 286     }
 287     AliInfo(Form("There are %d tracks in this event", fInputEvent->GetNumberOfTracks()));
 288
 289     if (!mcEvent)
 290     {
 291       AliError("ERROR: Could not retrieve MC event");
 292       return;
 293     }
 294
 295     fCFManager1->SetMCEventInfo(mcEvent);
 296     fCFManager2->SetMCEventInfo(mcEvent);
 297
 298     //check multiplicity
 299     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
 300     {
 301       AliWarning("Event does not pass multiplicity cuts"); return;
 302     }
 303
 304     //make event
 305     if (fAnalysisType == "ESDMCkineESD")
 306     {
 307       flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kESDkine, fCFManager1, fCFManager2 );
 308     }
 309     else if (fAnalysisType == "ESDMCkineMC")
 310     {
 311       flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kMCkine, fCFManager1, fCFManager2 );
 312     }
 313     // if monte carlo event get reaction plane from monte carlo (depends on generator)
 314     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
 315     //set reference multiplicity, TODO: maybe move it to the constructor?
 316     flowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
 317   }
 318   // Make the FlowEventSimple for AOD input
 319   else if (fAnalysisType == "AOD")
 320   {
 321     if (!myAOD)
 322     {
 323       AliError("ERROR: AOD not available");
 324       return;
 325     }
 326     AliInfo(Form("AOD has %d tracks", myAOD->GetNumberOfTracks()));
 327     flowEvent = new AliFlowEvent(myAOD);
 328   }
 329
 330   //check final event cuts
 331   Int_t mult = flowEvent->NumberOfTracks();
 332   AliInfo(Form("FlowEvent has %i tracks",mult));
 333   if (mult<fMinMult || mult>fMaxMult)
 334   {
 335     AliWarning("FlowEvent cut on multiplicity"); return;
 336   }
 337
 338   //tag subEvents
 339   flowEvent->TagSubeventsInEta(fMinA,fMaxA,fMinB,fMaxB);
 340
 341   //fListHistos->Print();
 342   //fOutputFile->WriteObject(flowEvent,"myFlowEventSimple");
 343   PostData(1,flowEvent);
 344   if (fQA)
 345   {
 346     PostData(2,fQAInt);
 347     PostData(3,fQADiff);
 348   }
 349 }
 350
 351 //________________________________________________________________________
 352 void AliAnalysisTaskFlowEvent::Terminate(Option_t *)
 353 {
 354   // Called once at the end of the query -- do not call in case of CAF
 355 }
 356
 357