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