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