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