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