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