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 #include "AliFlowTrackCuts.h"
  66 #include "AliFlowEventCuts.h"
  67 #include "AliFlowCommonConstants.h"
  68 #include "AliAnalysisTaskFlowEvent.h"
  69
  70 #include "AliLog.h"
  71
  72 ClassImp(AliAnalysisTaskFlowEvent)
  73
  74 //________________________________________________________________________
  75 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent() :
  76   AliAnalysisTaskSE(),
  77   //  fOutputFile(NULL),
  78   fAnalysisType("ESD"),
  79   fRPType("Global"),
  80   fCFManager1(NULL),
  81   fCFManager2(NULL),
  82   fCutsEvent(NULL),
  83   fCutsRP(NULL),
  84   fCutsPOI(NULL),
  85   fQAInt(NULL),
  86   fQADiff(NULL),
  87   fMinMult(0),
  88   fMaxMult(10000000),
  89   fMinA(-1.0),
  90   fMaxA(-0.01),
  91   fMinB(0.01),
  92   fMaxB(1.0),
  93   fQA(kFALSE),
  94   fNbinsMult(10000),
  95   fNbinsPt(100),
  96   fNbinsPhi(100),
  97   fNbinsEta(200),
  98   fNbinsQ(500),
  99   fMultMin(0.),
 100   fMultMax(10000.),
 101   fPtMin(0.),
 102   fPtMax(10.),
 103   fPhiMin(0.),
 104   fPhiMax(TMath::TwoPi()),
 105   fEtaMin(-5.),
 106   fEtaMax(5.),
 107   fQMin(0.),
 108   fQMax(3.),
 109   fExcludedEtaMin(0.),
 110   fExcludedEtaMax(0.),
 111   fExcludedPhiMin(0.),
 112   fExcludedPhiMax(0.),
 113   fAfterburnerOn(kFALSE),
 114   fNonFlowNumberOfTrackClones(0),
 115   fV1(0.),
 116   fV2(0.),
 117   fV3(0.),
 118   fV4(0.),
 119   fMyTRandom3(NULL)
 120 {
 121   // Constructor
 122   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()"<<endl;
 123 }
 124
 125 //________________________________________________________________________
 126 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, TString RPtype, Bool_t on, UInt_t iseed) :
 127   AliAnalysisTaskSE(name),
 128   //  fOutputFile(NULL),
 129   fAnalysisType("ESD"),
 130   fRPType(RPtype),
 131   fCFManager1(NULL),
 132   fCFManager2(NULL),
 133   fCutsEvent(NULL),
 134   fCutsRP(NULL),
 135   fCutsPOI(NULL),
 136   fQAInt(NULL),
 137   fQADiff(NULL),
 138   fMinMult(0),
 139   fMaxMult(10000000),
 140   fMinA(-1.0),
 141   fMaxA(-0.01),
 142   fMinB(0.01),
 143   fMaxB(1.0),
 144   fQA(on),
 145   fNbinsMult(10000),
 146   fNbinsPt(100),
 147   fNbinsPhi(100),
 148   fNbinsEta(200),
 149   fNbinsQ(500),
 150   fMultMin(0.),
 151   fMultMax(10000.),
 152   fPtMin(0.),
 153   fPtMax(10.),
 154   fPhiMin(0.),
 155   fPhiMax(TMath::TwoPi()),
 156   fEtaMin(-5.),
 157   fEtaMax(5.),
 158   fQMin(0.),
 159   fQMax(3.),
 160   fExcludedEtaMin(0.),
 161   fExcludedEtaMax(0.),
 162   fExcludedPhiMin(0.),
 163   fExcludedPhiMax(0.),
 164   fAfterburnerOn(kFALSE),
 165   fNonFlowNumberOfTrackClones(0),
 166   fV1(0.),
 167   fV2(0.),
 168   fV3(0.),
 169   fV4(0.),
 170   fMyTRandom3(NULL)
 171 {
 172   // Constructor
 173   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)"<<endl;
 174   fMyTRandom3 = new TRandom3(iseed);
 175   gRandom->SetSeed(fMyTRandom3->Integer(65539));
 176
 177   //FMD input slot
 178   if (strcmp(RPtype,"FMD")==0) {
 179     DefineInput(1, TList::Class());
 180   }
 181
 182   // Define output slots here
 183   // Define here the flow event output
 184   DefineOutput(1, AliFlowEventSimple::Class());
 185   if(on)
 186   {
 187     DefineOutput(2, TList::Class());
 188     DefineOutput(3, TList::Class());
 189   }
 190   // and for testing open an output file
 191   //  fOutputFile = new TFile("FlowEvents.root","RECREATE");
 192
 193 }
 194
 195 //________________________________________________________________________
 196 AliAnalysisTaskFlowEvent::~AliAnalysisTaskFlowEvent()
 197 {
 198   //
 199   // Destructor
 200   //
 201   if (fMyTRandom3) delete fMyTRandom3;
 202   // objects in the output list are deleted
 203   // by the TSelector dtor (I hope)
 204
 205 }
 206
 207 //________________________________________________________________________
 208 void AliAnalysisTaskFlowEvent::UserCreateOutputObjects()
 209 {
 210   // Called at every worker node to initialize
 211   cout<<"AliAnalysisTaskFlowEvent::CreateOutputObjects()"<<endl;
 212
 213   if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD"  || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC" || fAnalysisType == "MK"))
 214   {
 215     AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD and MC are allowed.");
 216     exit(1);
 217   }
 218
 219   //set the common constants
 220   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 221   cc->SetNbinsMult(fNbinsMult);
 222   cc->SetNbinsPt(fNbinsPt);
 223   cc->SetNbinsPhi(fNbinsPhi);
 224   cc->SetNbinsEta(fNbinsEta);
 225   cc->SetNbinsQ(fNbinsQ);
 226   cc->SetMultMin(fMultMin);
 227   cc->SetMultMax(fMultMax);
 228   cc->SetPtMin(fPtMin);
 229   cc->SetPtMax(fPtMax);
 230   cc->SetPhiMin(fPhiMin);
 231   cc->SetPhiMax(fPhiMax);
 232   cc->SetEtaMin(fEtaMin);
 233   cc->SetEtaMax(fEtaMax);
 234   cc->SetQMin(fQMin);
 235   cc->SetQMax(fQMax);
 236
 237 }
 238
 239 //________________________________________________________________________
 240 void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
 241 {
 242   // Main loop
 243   // Called for each event
 244   AliFlowEvent* flowEvent = NULL;
 245   AliMCEvent*  mcEvent = MCEvent();                              // from TaskSE
 246   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
 247   AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
 248   AliMultiplicity* myTracklets = NULL;
 249   TH2F* histFMD = NULL;
 250
 251   if(GetNinputs()==2) {
 252     TList* FMDdata = dynamic_cast<TList*>(GetInputData(1));
 253     if(!FMDdata) {
 254       cout<<" No FMDdata "<<endl;
 255       exit(2);
 256     }
 257     histFMD = dynamic_cast<TH2F*>(FMDdata->FindObject("dNdetadphiHistogramTrVtx"));
 258     if (!histFMD) {
 259       cout<< "No histFMD"<<endl;
 260       exit(2);
 261     }
 262   }
 263
 264   //use the new and temporarily inclomplete way of doing things
 265   if (fAnalysisType == "MK")
 266   {
 267     if (!(fCutsRP&&fCutsPOI))
 268     {
 269       AliError("cuts not set");
 270       return;
 271     }
 272     if (fCutsEvent)
 273     {
 274       if (!fCutsEvent->IsSelected(InputEvent())) return;
 275     }
 276
 277     fCutsRP->SetMCevent( MCEvent() );
 278     fCutsPOI->SetMCevent( MCEvent() );
 279     flowEvent = new AliFlowEvent( InputEvent(), fCutsRP, fCutsPOI );
 280     if (myESD)
 281       flowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
 282   }
 283
 284   // Make the FlowEvent for MC input
 285   if (fAnalysisType == "MC")
 286   {
 287     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 288     // This handler can return the current MC event
 289     if (!(fCFManager1&&fCFManager2))
 290     {
 291       AliError("ERROR: No pointer to correction framework cuts! ");
 292       return;
 293     }
 294     if (!mcEvent)
 295     {
 296       AliError("ERROR: Could not retrieve MC event");
 297       return;
 298     }
 299
 300     fCFManager1->SetMCEventInfo(mcEvent);
 301     fCFManager2->SetMCEventInfo(mcEvent);
 302
 303     AliInfo(Form("Number of MC particles: %d", mcEvent->GetNumberOfTracks()));
 304
 305     //check multiplicity
 306     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtGenCuts,mcEvent))
 307     {
 308       AliWarning("Event does not pass multiplicity cuts"); return;
 309     }
 310     //make event
 311     flowEvent = new AliFlowEvent(mcEvent,fCFManager1,fCFManager2);
 312   }
 313
 314   // Make the FlowEvent for ESD input
 315   else if (fAnalysisType == "ESD")
 316   {
 317     if (!(fCFManager1&&fCFManager2))
 318     {
 319       AliError("ERROR: No pointer to correction framework cuts!");
 320       return;
 321     }
 322     if (!myESD)
 323     {
 324       AliError("ERROR: ESD not available");
 325       return;
 326     }
 327
 328     //check the offline trigger (check if the event has the correct trigger)
 329     AliInfo(Form("ESD has %d tracks", fInputEvent->GetNumberOfTracks()));
 330
 331     //check multiplicity
 332     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
 333     {
 334       AliWarning("Event does not pass multiplicity cuts"); return;
 335     }
 336
 337     //make event
 338     if (fRPType == "Global") {
 339       flowEvent = new AliFlowEvent(myESD,fCFManager1,fCFManager2);
 340     }
 341     if (fRPType == "TPCOnly") {
 342       flowEvent = new AliFlowEvent(myESD,fCFManager2,kFALSE);
 343     }
 344     if (fRPType == "TPCHybrid") {
 345       flowEvent = new AliFlowEvent(myESD,fCFManager2,kTRUE);
 346     }
 347     else if (fRPType == "Tracklet"){
 348       flowEvent = new AliFlowEvent(myESD,myTracklets,fCFManager2);
 349     }
 350     else if (fRPType == "FMD"){
 351       flowEvent = new AliFlowEvent(myESD,histFMD,fCFManager2);
 352     }
 353
 354     // if monte carlo event get reaction plane from monte carlo (depends on generator)
 355     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
 356     //set reference multiplicity, TODO: maybe move it to the constructor?
 357     flowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
 358   }
 359
 360   // Make the FlowEvent for ESD input combined with MC info
 361   else if (fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" )
 362   {
 363     if (!(fCFManager1&&fCFManager2))
 364     {
 365       AliError("ERROR: No pointer to correction framework cuts! ");
 366       return;
 367     }
 368     if (!myESD)
 369     {
 370       AliError("ERROR: ESD not available");
 371       return;
 372     }
 373     AliInfo(Form("There are %d tracks in this event", fInputEvent->GetNumberOfTracks()));
 374
 375     if (!mcEvent)
 376     {
 377       AliError("ERROR: Could not retrieve MC event");
 378       return;
 379     }
 380
 381     fCFManager1->SetMCEventInfo(mcEvent);
 382     fCFManager2->SetMCEventInfo(mcEvent);
 383
 384     //check multiplicity
 385     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
 386     {
 387       AliWarning("Event does not pass multiplicity cuts"); return;
 388     }
 389
 390     //make event
 391     if (fAnalysisType == "ESDMCkineESD")
 392     {
 393       flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kESDkine, fCFManager1, fCFManager2 );
 394     }
 395     else if (fAnalysisType == "ESDMCkineMC")
 396     {
 397       flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kMCkine, fCFManager1, fCFManager2 );
 398     }
 399     // if monte carlo event get reaction plane from monte carlo (depends on generator)
 400     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
 401     //set reference multiplicity, TODO: maybe move it to the constructor?
 402     flowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
 403   }
 404   // Make the FlowEventSimple for AOD input
 405   else if (fAnalysisType == "AOD")
 406   {
 407     if (!myAOD)
 408     {
 409       AliError("ERROR: AOD not available");
 410       return;
 411     }
 412     AliInfo(Form("AOD has %d tracks", myAOD->GetNumberOfTracks()));
 413     flowEvent = new AliFlowEvent(myAOD);
 414   }
 415
 416   //check final event cuts
 417   Int_t mult = flowEvent->NumberOfTracks();
 418   AliInfo(Form("FlowEvent has %i tracks",mult));
 419   if (mult<fMinMult || mult>fMaxMult)
 420   {
 421     AliWarning("FlowEvent cut on multiplicity"); return;
 422   }
 423
 424   //define dead zone
 425   flowEvent->DefineDeadZone(fExcludedEtaMin, fExcludedEtaMax, fExcludedPhiMin, fExcludedPhiMax );
 426
 427
 428   //////////////////////////////////////////////////////////////////////////////
 429   ///////////////////////////AFTERBURNER
 430   if (fAfterburnerOn)
 431   {
 432     //if reaction plane not set from elsewhere randomize it before adding flow
 433     if (!flowEvent->IsSetMCReactionPlaneAngle())
 434       flowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
 435
 436     flowEvent->AddFlow(fV1,fV2,fV3,fV4);     //add flow
 437     flowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
 438   }
 439   //////////////////////////////////////////////////////////////////////////////
 440
 441   //tag subEvents
 442   flowEvent->TagSubeventsInEta(fMinA,fMaxA,fMinB,fMaxB);
 443
 444   //fListHistos->Print();
 445   //fOutputFile->WriteObject(flowEvent,"myFlowEventSimple");
 446   PostData(1,flowEvent);
 447   if (fQA)
 448   {
 449     PostData(2,fQAInt);
 450     PostData(3,fQADiff);
 451   }
 452 }
 453
 454 //________________________________________________________________________
 455 void AliAnalysisTaskFlowEvent::Terminate(Option_t *)
 456 {
 457   // Called once at the end of the query -- do not call in case of CAF
 458 }
 459
 460