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