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 #include "TF1.h"
  37
  38 // ALICE Analysis Framework
  39 #include "AliAnalysisManager.h"
  40 #include "AliAnalysisTaskSE.h"
  41 #include "AliESDtrackCuts.h"
  42
  43 // ESD interface
  44 #include "AliESDEvent.h"
  45 #include "AliESDInputHandler.h"
  46
  47 // AOD interface
  48 #include "AliAODEvent.h"
  49 #include "AliAODInputHandler.h"
  50
  51 // Monte Carlo Event
  52 #include "AliMCEventHandler.h"
  53 #include "AliMCEvent.h"
  54
  55 // ALICE Correction Framework
  56 #include "AliCFManager.h"
  57
  58 // Interface to Event generators to get Reaction Plane Angle
  59 #include "AliGenCocktailEventHeader.h"
  60 #include "AliGenHijingEventHeader.h"
  61 #include "AliGenGeVSimEventHeader.h"
  62 #include "AliGenEposEventHeader.h"
  63
  64 // Interface to Load short life particles
  65 #include "TObjArray.h"
  66 #include "AliFlowCandidateTrack.h"
  67
  68 // Interface to make the Flow Event Simple used in the flow analysis methods
  69 #include "AliFlowEvent.h"
  70 #include "AliFlowTrackCuts.h"
  71 #include "AliFlowEventCuts.h"
  72 #include "AliFlowCommonConstants.h"
  73 #include "AliAnalysisTaskFlowEvent.h"
  74
  75 #include "AliLog.h"
  76
  77 ClassImp(AliAnalysisTaskFlowEvent)
  78
  79 //________________________________________________________________________
  80 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent() :
  81   AliAnalysisTaskSE(),
  82   //  fOutputFile(NULL),
  83   fAnalysisType("AUTOMATIC"),
  84   fRPType("Global"),
  85   fCFManager1(NULL),
  86   fCFManager2(NULL),
  87   fCutsEvent(NULL),
  88   fCutsRP(NULL),
  89   fCutsPOI(NULL),
  90   fQAInt(NULL),
  91   fQADiff(NULL),
  92   fMinMult(0),
  93   fMaxMult(10000000),
  94   fMinA(-1.0),
  95   fMaxA(-0.01),
  96   fMinB(0.01),
  97   fMaxB(1.0),
  98   fQA(kFALSE),
  99   fLoadCandidates(kFALSE),
 100   fNbinsMult(10000),
 101   fNbinsPt(100),
 102   fNbinsPhi(100),
 103   fNbinsEta(200),
 104   fNbinsQ(500),
 105   fMultMin(0.),
 106   fMultMax(10000.),
 107   fPtMin(0.),
 108   fPtMax(10.),
 109   fPhiMin(0.),
 110   fPhiMax(TMath::TwoPi()),
 111   fEtaMin(-5.),
 112   fEtaMax(5.),
 113   fQMin(0.),
 114   fQMax(3.),
 115   fExcludedEtaMin(0.),
 116   fExcludedEtaMax(0.),
 117   fExcludedPhiMin(0.),
 118   fExcludedPhiMax(0.),
 119   fAfterburnerOn(kFALSE),
 120   fNonFlowNumberOfTrackClones(0),
 121   fV1(0.),
 122   fV2(0.),
 123   fV3(0.),
 124   fV4(0.),
 125   fV2Function(0),
 126   fMyTRandom3(NULL)
 127 {
 128   // Constructor
 129   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()"<<endl;
 130 }
 131
 132 //________________________________________________________________________
 133 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates) :
 134   AliAnalysisTaskSE(name),
 135   //  fOutputFile(NULL),
 136   fAnalysisType("AUTOMATIC"),
 137   fRPType(RPtype),
 138   fCFManager1(NULL),
 139   fCFManager2(NULL),
 140   fCutsEvent(NULL),
 141   fCutsRP(NULL),
 142   fCutsPOI(NULL),
 143   fQAInt(NULL),
 144   fQADiff(NULL),
 145   fMinMult(0),
 146   fMaxMult(10000000),
 147   fMinA(-1.0),
 148   fMaxA(-0.01),
 149   fMinB(0.01),
 150   fMaxB(1.0),
 151   fQA(on),
 152   fLoadCandidates(bCandidates),
 153   fNbinsMult(10000),
 154   fNbinsPt(100),
 155   fNbinsPhi(100),
 156   fNbinsEta(200),
 157   fNbinsQ(500),
 158   fMultMin(0.),
 159   fMultMax(10000.),
 160   fPtMin(0.),
 161   fPtMax(10.),
 162   fPhiMin(0.),
 163   fPhiMax(TMath::TwoPi()),
 164   fEtaMin(-5.),
 165   fEtaMax(5.),
 166   fQMin(0.),
 167   fQMax(3.),
 168   fExcludedEtaMin(0.),
 169   fExcludedEtaMax(0.),
 170   fExcludedPhiMin(0.),
 171   fExcludedPhiMax(0.),
 172   fAfterburnerOn(kFALSE),
 173   fNonFlowNumberOfTrackClones(0),
 174   fV1(0.),
 175   fV2(0.),
 176   fV3(0.),
 177   fV4(0.),
 178   fV2Function(0),
 179   fMyTRandom3(NULL)
 180 {
 181   // Constructor
 182   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)"<<endl;
 183   fMyTRandom3 = new TRandom3(iseed);
 184   gRandom->SetSeed(fMyTRandom3->Integer(65539));
 185
 186   int availableINslot=1;
 187   //FMD input slot
 188   if (strcmp(RPtype,"FMD")==0) {
 189     DefineInput(availableINslot++, TList::Class());
 190   }
 191   //Candidates input slot
 192   if( fLoadCandidates )
 193     DefineInput(availableINslot, TObjArray::Class());
 194
 195   // Define output slots here
 196   // Define here the flow event output
 197   DefineOutput(1, AliFlowEventSimple::Class());
 198   if(on)
 199   {
 200     DefineOutput(2, TList::Class());
 201     DefineOutput(3, TList::Class());
 202   }
 203
 204   // and for testing open an output file
 205   //  fOutputFile = new TFile("FlowEvents.root","RECREATE");
 206
 207 }
 208
 209 //________________________________________________________________________
 210 AliAnalysisTaskFlowEvent::~AliAnalysisTaskFlowEvent()
 211 {
 212   //
 213   // Destructor
 214   //
 215   delete fMyTRandom3;
 216   // objects in the output list are deleted
 217   // by the TSelector dtor (I hope)
 218
 219 }
 220
 221 //________________________________________________________________________
 222 void AliAnalysisTaskFlowEvent::NotifyRun()
 223 {
 224   //at the beginning of each new run
 225         AliESDEvent* fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 226   if (!fESD) return;
 227
 228   Int_t run = fESD->GetRunNumber();
 229   AliInfo(Form("Stariting run #%i",run));
 230 }
 231
 232 //________________________________________________________________________
 233 void AliAnalysisTaskFlowEvent::UserCreateOutputObjects()
 234 {
 235   // Called at every worker node to initialize
 236   cout<<"AliAnalysisTaskFlowEvent::CreateOutputObjects()"<<endl;
 237
 238   if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD"  || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC" || fAnalysisType == "AUTOMATIC"))
 239   {
 240     AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD, MC and AUTOMATIC are allowed.");
 241     exit(1);
 242   }
 243
 244   //set the common constants
 245   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 246   cc->SetNbinsMult(fNbinsMult);
 247   cc->SetNbinsPt(fNbinsPt);
 248   cc->SetNbinsPhi(fNbinsPhi);
 249   cc->SetNbinsEta(fNbinsEta);
 250   cc->SetNbinsQ(fNbinsQ);
 251   cc->SetMultMin(fMultMin);
 252   cc->SetMultMax(fMultMax);
 253   cc->SetPtMin(fPtMin);
 254   cc->SetPtMax(fPtMax);
 255   cc->SetPhiMin(fPhiMin);
 256   cc->SetPhiMax(fPhiMax);
 257   cc->SetEtaMin(fEtaMin);
 258   cc->SetEtaMax(fEtaMax);
 259   cc->SetQMin(fQMin);
 260   cc->SetQMax(fQMax);
 261
 262 }
 263
 264 //________________________________________________________________________
 265 void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
 266 {
 267   // Main loop
 268   // Called for each event
 269   AliFlowEvent* flowEvent = NULL;
 270   AliMCEvent*  mcEvent = MCEvent();                              // from TaskSE
 271   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
 272   AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
 273   AliMultiplicity* myTracklets = NULL;
 274   AliESDPmdTrack* pmdtracks = NULL;//pmd
 275
 276   int availableINslot=1;
 277
 278   if (!(fCutsRP&&fCutsPOI&&fCutsEvent))
 279   {
 280     AliError("cuts not set");
 281     return;
 282   }
 283
 284   //use the new and temporarily incomplete way of doing things
 285   if (fAnalysisType == "AUTOMATIC")
 286   {
 287     //check event cuts
 288     if (!fCutsEvent->IsSelected(InputEvent())) return;
 289
 290     //first attach all possible information to the cuts
 291     fCutsRP->SetEvent( InputEvent(), MCEvent() );  //attach event
 292     fCutsPOI->SetEvent( InputEvent(), MCEvent() );
 293
 294     //then make the event
 295     flowEvent = new AliFlowEvent( fCutsRP, fCutsPOI );
 296     if (myESD)
 297       flowEvent->SetReferenceMultiplicity(fCutsEvent->GetReferenceMultiplicity(InputEvent()));
 298     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
 299   }
 300
 301   // Make the FlowEvent for MC input
 302   else if (fAnalysisType == "MC")
 303   {
 304     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 305     // This handler can return the current MC event
 306     if (!(fCFManager1&&fCFManager2))
 307     {
 308       AliError("ERROR: No pointer to correction framework cuts! ");
 309       return;
 310     }
 311     if (!mcEvent)
 312     {
 313       AliError("ERROR: Could not retrieve MC event");
 314       return;
 315     }
 316
 317     fCFManager1->SetMCEventInfo(mcEvent);
 318     fCFManager2->SetMCEventInfo(mcEvent);
 319
 320     AliInfo(Form("Number of MC particles: %d", mcEvent->GetNumberOfTracks()));
 321
 322     //check multiplicity
 323     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtGenCuts,mcEvent))
 324     {
 325       AliWarning("Event does not pass multiplicity cuts"); return;
 326     }
 327     //make event
 328     flowEvent = new AliFlowEvent(mcEvent,fCFManager1,fCFManager2);
 329   }
 330
 331   // Make the FlowEvent for ESD input
 332   else if (fAnalysisType == "ESD")
 333   {
 334     if (!(fCFManager1&&fCFManager2))
 335     {
 336       AliError("ERROR: No pointer to correction framework cuts!");
 337       return;
 338     }
 339     if (!myESD)
 340     {
 341       AliError("ERROR: ESD not available");
 342       return;
 343     }
 344
 345     //check the offline trigger (check if the event has the correct trigger)
 346     AliInfo(Form("ESD has %d tracks", fInputEvent->GetNumberOfTracks()));
 347
 348     //check multiplicity
 349     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
 350     {
 351       AliWarning("Event does not pass multiplicity cuts"); return;
 352     }
 353
 354     //make event
 355     if (fRPType == "Global") {
 356       flowEvent = new AliFlowEvent(myESD,fCFManager1,fCFManager2);
 357     }
 358     else if (fRPType == "TPCOnly") {
 359       flowEvent = new AliFlowEvent(myESD,fCFManager2,kFALSE);
 360     }
 361     else if (fRPType == "TPCHybrid") {
 362       flowEvent = new AliFlowEvent(myESD,fCFManager2,kTRUE);
 363     }
 364     else if (fRPType == "Tracklet"){
 365       flowEvent = new AliFlowEvent(myESD,myTracklets,fCFManager2);
 366     }
 367     else if (fRPType == "FMD"){
 368       TList* FMDdata = dynamic_cast<TList*>(GetInputData(availableINslot++));
 369       if(!FMDdata) {
 370         cout<<" No FMDdata "<<endl;
 371         return;
 372       }
 373       TH2F* histFMD = dynamic_cast<TH2F*>(FMDdata->FindObject("dNdetadphiHistogramTrVtx"));
 374       if (!histFMD) {
 375         cout<< "No histFMD"<<endl;
 376         return;
 377       }
 378       flowEvent = new AliFlowEvent(myESD,histFMD,fCFManager2);
 379     }
 380     else if (fRPType == "PMD"){
 381       flowEvent = new AliFlowEvent(myESD,pmdtracks,fCFManager2);
 382     }
 383     else return;
 384
 385     // if monte carlo event get reaction plane from monte carlo (depends on generator)
 386     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
 387     //set reference multiplicity, TODO: maybe move it to the constructor?
 388     flowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
 389   }
 390
 391   // Make the FlowEvent for ESD input combined with MC info
 392   else if (fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" )
 393   {
 394     if (!(fCFManager1&&fCFManager2))
 395     {
 396       AliError("ERROR: No pointer to correction framework cuts! ");
 397       return;
 398     }
 399     if (!myESD)
 400     {
 401       AliError("ERROR: ESD not available");
 402       return;
 403     }
 404     AliInfo(Form("There are %d tracks in this event", fInputEvent->GetNumberOfTracks()));
 405
 406     if (!mcEvent)
 407     {
 408       AliError("ERROR: Could not retrieve MC event");
 409       return;
 410     }
 411
 412     fCFManager1->SetMCEventInfo(mcEvent);
 413     fCFManager2->SetMCEventInfo(mcEvent);
 414
 415     //check multiplicity
 416     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
 417     {
 418       AliWarning("Event does not pass multiplicity cuts"); return;
 419     }
 420
 421     //make event
 422     if (fAnalysisType == "ESDMCkineESD")
 423     {
 424       flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kESDkine, fCFManager1, fCFManager2 );
 425     }
 426     else if (fAnalysisType == "ESDMCkineMC")
 427     {
 428       flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kMCkine, fCFManager1, fCFManager2 );
 429     }
 430     if (!flowEvent) return;
 431     // if monte carlo event get reaction plane from monte carlo (depends on generator)
 432     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
 433     //set reference multiplicity, TODO: maybe move it to the constructor?
 434     flowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
 435   }
 436   // Make the FlowEventSimple for AOD input
 437   else if (fAnalysisType == "AOD")
 438   {
 439     if (!myAOD)
 440     {
 441       AliError("ERROR: AOD not available");
 442       return;
 443     }
 444     AliInfo(Form("AOD has %d tracks", myAOD->GetNumberOfTracks()));
 445     flowEvent = new AliFlowEvent(myAOD);
 446   }
 447
 448   //inject candidates
 449   if(fLoadCandidates) {
 450     TObjArray* Candidates = dynamic_cast<TObjArray*>(GetInputData(availableINslot++));
 451     //if(Candidates->GetEntriesFast())
 452       //printf("I received %d candidates\n",Candidates->GetEntriesFast());
 453     if (Candidates)
 454     {
 455       for(int iCand=0; iCand!=Candidates->GetEntriesFast(); ++iCand ) {
 456         AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(Candidates->At(iCand));
 457         if (!cand) continue;
 458         cand->SetForPOISelection(kTRUE);
 459         cand->SetForRPSelection(kFALSE);
 460         //printf(" Ⱶ Checking at candidate %d with %d daughters\n",iCand,cand->GetNDaughters());
 461         for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
 462           //printf("    Ⱶ Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau) );
 463           for(int iRPs=0; iRPs!=flowEvent->NumberOfTracks(); ++iRPs ) {
 464             AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(flowEvent->GetTrack( iRPs ));
 465             if (!iRP) continue;
 466             if( !iRP->InRPSelection() )
 467               continue;
 468             if( cand->GetIDDaughter(iDau) == iRP->GetID() ) {
 469               //printf(" was in RP set");
 470               cand->SetDaughter( iDau, iRP );
 471               //temporarily untagging all daugters
 472               iRP->SetForRPSelection(kFALSE);
 473             }
 474           }
 475           //printf("\n");
 476         }
 477         flowEvent->AddTrack(cand);
 478       }
 479     }
 480   }
 481
 482   if (!flowEvent) return; //shuts up coverity
 483
 484   //check final event cuts
 485   Int_t mult = flowEvent->NumberOfTracks();
 486   AliInfo(Form("FlowEvent has %i tracks",mult));
 487   if (mult<fMinMult || mult>fMaxMult)
 488   {
 489     AliWarning("FlowEvent cut on multiplicity"); return;
 490   }
 491
 492   //define dead zone
 493   flowEvent->DefineDeadZone(fExcludedEtaMin, fExcludedEtaMax, fExcludedPhiMin, fExcludedPhiMax );
 494
 495
 496   //////////////////////////////////////////////////////////////////////////////
 497   ///////////////////////////AFTERBURNER
 498   if (fAfterburnerOn)
 499   {
 500     //if reaction plane not set from elsewhere randomize it before adding flow
 501     if (!flowEvent->IsSetMCReactionPlaneAngle())
 502       flowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
 503
 504     //flowEvent->AddFlow(fV1,fV2,fV3,fV4);     //add flow
 505     flowEvent->AddV2(fV2Function);
 506     flowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
 507   }
 508   //////////////////////////////////////////////////////////////////////////////
 509
 510   //tag subEvents
 511   flowEvent->TagSubeventsInEta(fMinA,fMaxA,fMinB,fMaxB);
 512
 513   //fListHistos->Print();
 514   //fOutputFile->WriteObject(flowEvent,"myFlowEventSimple");
 515   PostData(1,flowEvent);
 516   if (fQA)
 517   {
 518     PostData(2,fQAInt);
 519     PostData(3,fQADiff);
 520   }
 521 }
 522
 523 //________________________________________________________________________
 524 void AliAnalysisTaskFlowEvent::Terminate(Option_t *)
 525 {
 526   // Called once at the end of the query -- do not call in case of CAF
 527 }
 528