PWG/EMCAL/AliAnalysisTaskEmcal.cxx

   1 // $Id$
   2 //
   3 // Emcal base analysis task.
   4 //
   5 // Author: S.Aiola, M. Verweij
   6
   7 #include "AliAnalysisTaskEmcal.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TList.h>
  11 #include <TObject.h>
  12 #include <TH1F.h>
  13 #include <TProfile.h>
  14 #include <TSystem.h>
  15 #include <TFile.h>
  16 #include <TChain.h>
  17 #include <TKey.h>
  18
  19 #include "AliAODEvent.h"
  20 #include "AliAnalysisManager.h"
  21 #include "AliCentrality.h"
  22 #include "AliEMCALGeometry.h"
  23 #include "AliESDEvent.h"
  24 #include "AliEmcalParticle.h"
  25 #include "AliEventplane.h"
  26 #include "AliInputEventHandler.h"
  27 #include "AliLog.h"
  28 #include "AliMCParticle.h"
  29 #include "AliVCluster.h"
  30 #include "AliVEventHandler.h"
  31 #include "AliVParticle.h"
  32 #include "AliVCaloTrigger.h"
  33 #include "AliGenPythiaEventHeader.h"
  34 #include "AliAODMCHeader.h"
  35 #include "AliMCEvent.h"
  36 #include "AliAnalysisUtils.h"
  37 #include "AliEmcalTriggerPatchInfo.h"
  38
  39 #include "AliParticleContainer.h"
  40 #include "AliClusterContainer.h"
  41
  42 ClassImp(AliAnalysisTaskEmcal)
  43
  44 //________________________________________________________________________
  45 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() :
  46   AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
  47   fForceBeamType(kNA),
  48   fGeneralHistograms(kFALSE),
  49   fInitialized(kFALSE),
  50   fCreateHisto(kTRUE),
  51   fCaloCellsName(),
  52   fCaloTriggersName(),
  53   fCaloTriggerPatchInfoName(),
  54   fMinCent(-999),
  55   fMaxCent(-999),
  56   fMinVz(-999),
  57   fMaxVz(-999),
  58   fTrackPtCut(0),
  59   fMinNTrack(0),
  60   fUseAliAnaUtils(kFALSE),
  61   fAliAnalysisUtils(0x0),
  62   fOffTrigger(AliVEvent::kAny),
  63   fTrigClass(),
  64   fTriggerTypeSel(kND),
  65   fNbins(500),
  66   fMinBinPt(0),
  67   fMaxBinPt(250),
  68   fMinPtTrackInEmcal(0),
  69   fEventPlaneVsEmcal(-1),
  70   fMinEventPlane(-1e6),
  71   fMaxEventPlane(1e6),
  72   fCentEst("V0M"),
  73   fIsEmbedded(kFALSE),
  74   fIsPythia(kFALSE),
  75   fSelectPtHardBin(-999),
  76   fMinMCLabel(0),
  77   fMCLabelShift(0),
  78   fNcentBins(4),
  79   fGeom(0),
  80   fTracks(0),
  81   fCaloClusters(0),
  82   fCaloCells(0),
  83   fCaloTriggers(0),
  84   fTriggerPatchInfo(0),
  85   fCent(0),
  86   fCentBin(-1),
  87   fEPV0(-1.0),
  88   fEPV0A(-1.0),
  89   fEPV0C(-1.0),
  90   fNVertCont(0),
  91   fBeamType(kNA),
  92   fPythiaHeader(0),
  93   fPtHard(0),
  94   fPtHardBin(0),
  95   fNTrials(0),
  96   fParticleCollArray(),
  97   fClusterCollArray(),
  98   fMainTriggerPatch(0x0),
  99   fTriggerType(kND),
 100   fOutput(0),
 101   fHistTrialsAfterSel(0),
 102   fHistEventsAfterSel(0),
 103   fHistTrials(0),
 104   fHistXsection(0),
 105   fHistEvents(0),
 106   fHistPtHard(0),
 107   fHistCentrality(0),
 108   fHistZVertex(0),
 109   fHistEventPlane(0),
 110   fHistEventRejection(0)
 111 {
 112   // Default constructor.
 113
 114   fVertex[0] = 0;
 115   fVertex[1] = 0;
 116   fVertex[2] = 0;
 117
 118   fParticleCollArray.SetOwner(kTRUE);
 119   fClusterCollArray.SetOwner(kTRUE);
 120 }
 121
 122 //________________________________________________________________________
 123 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
 124   AliAnalysisTaskSE(name),
 125   fForceBeamType(kNA),
 126   fGeneralHistograms(kFALSE),
 127   fInitialized(kFALSE),
 128   fCreateHisto(histo),
 129   fCaloCellsName(),
 130   fCaloTriggersName(),
 131   fCaloTriggerPatchInfoName(),
 132   fMinCent(-999),
 133   fMaxCent(-999),
 134   fMinVz(-999),
 135   fMaxVz(-999),
 136   fTrackPtCut(0),
 137   fMinNTrack(0),
 138   fUseAliAnaUtils(kFALSE),
 139   fAliAnalysisUtils(0x0),
 140   fOffTrigger(AliVEvent::kAny),
 141   fTrigClass(),
 142   fTriggerTypeSel(kND),
 143   fNbins(500),
 144   fMinBinPt(0),
 145   fMaxBinPt(250),
 146   fMinPtTrackInEmcal(0),
 147   fEventPlaneVsEmcal(-1),
 148   fMinEventPlane(-1e6),
 149   fMaxEventPlane(1e6),
 150   fCentEst("V0M"),
 151   fIsEmbedded(kFALSE),
 152   fIsPythia(kFALSE),
 153   fSelectPtHardBin(-999),
 154   fMinMCLabel(0),
 155   fMCLabelShift(0),
 156   fNcentBins(4),
 157   fGeom(0),
 158   fTracks(0),
 159   fCaloClusters(0),
 160   fCaloCells(0),
 161   fCaloTriggers(0),
 162   fTriggerPatchInfo(0),
 163   fCent(0),
 164   fCentBin(-1),
 165   fEPV0(-1.0),
 166   fEPV0A(-1.0),
 167   fEPV0C(-1.0),
 168   fNVertCont(0),
 169   fBeamType(kNA),
 170   fPythiaHeader(0),
 171   fPtHard(0),
 172   fPtHardBin(0),
 173   fNTrials(0),
 174   fParticleCollArray(),
 175   fClusterCollArray(),
 176   fMainTriggerPatch(0x0),
 177   fTriggerType(kND),
 178   fOutput(0),
 179   fHistTrialsAfterSel(0),
 180   fHistEventsAfterSel(0),
 181   fHistTrials(0),
 182   fHistXsection(0),
 183   fHistEvents(0),
 184   fHistPtHard(0),
 185   fHistCentrality(0),
 186   fHistZVertex(0),
 187   fHistEventPlane(0),
 188   fHistEventRejection(0)
 189 {
 190   // Standard constructor.
 191
 192   fVertex[0] = 0;
 193   fVertex[1] = 0;
 194   fVertex[2] = 0;
 195
 196   fParticleCollArray.SetOwner(kTRUE);
 197   fClusterCollArray.SetOwner(kTRUE);
 198
 199   if (fCreateHisto) {
 200     DefineOutput(1, TList::Class());
 201   }
 202 }
 203
 204 //________________________________________________________________________
 205 AliAnalysisTaskEmcal::~AliAnalysisTaskEmcal()
 206 {
 207   // Destructor
 208 }
 209
 210 //________________________________________________________________________
 211 void AliAnalysisTaskEmcal::SetClusPtCut(Double_t cut, Int_t c)
 212 {
 213   AliClusterContainer *cont = GetClusterContainer(c);
 214   if (cont) cont->SetClusPtCut(cut);
 215   else AliError(Form("%s in SetClusPtCut(...): container %d not found",GetName(),c));
 216 }
 217
 218 //________________________________________________________________________
 219 void AliAnalysisTaskEmcal::SetClusTimeCut(Double_t min, Double_t max, Int_t c)
 220 {
 221   AliClusterContainer *cont = GetClusterContainer(c);
 222   if (cont) cont->SetClusTimeCut(min,max);
 223   else AliError(Form("%s in SetClusTimeCut(...): container %d not found",GetName(),c));
 224 }
 225
 226 //________________________________________________________________________
 227 void AliAnalysisTaskEmcal::SetTrackPtCut(Double_t cut, Int_t c)
 228 {
 229   AliParticleContainer *cont = GetParticleContainer(c);
 230   if (cont) cont->SetParticlePtCut(cut);
 231   else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
 232
 233   fTrackPtCut = cut;
 234 }
 235
 236 //________________________________________________________________________
 237 void AliAnalysisTaskEmcal::SetTrackEtaLimits(Double_t min, Double_t max, Int_t c)
 238 {
 239   AliParticleContainer *cont = GetParticleContainer(c);
 240   if (cont) cont->SetParticleEtaLimits(min,max);
 241   else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
 242 }
 243
 244 //________________________________________________________________________
 245 void AliAnalysisTaskEmcal::SetTrackPhiLimits(Double_t min, Double_t max, Int_t c)
 246 {
 247   AliParticleContainer *cont = GetParticleContainer(c);
 248   if (cont) cont->SetParticlePhiLimits(min,max);
 249   else AliError(Form("%s in SetTrackPhiLimits(...): container %d not found",GetName(),c));
 250 }
 251
 252 //________________________________________________________________________
 253 void AliAnalysisTaskEmcal::UserCreateOutputObjects()
 254 {
 255   // Create user output.
 256   if (!fCreateHisto)
 257     return;
 258
 259   OpenFile(1);
 260   fOutput = new TList();
 261   fOutput->SetOwner();
 262
 263   if (fForceBeamType == kpp)
 264     fNcentBins = 1;
 265
 266   if (!fGeneralHistograms)
 267     return;
 268
 269   if (fIsPythia) {
 270     fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
 271     fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
 272     fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
 273     fOutput->Add(fHistTrialsAfterSel);
 274
 275     fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
 276     fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
 277     fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
 278     fOutput->Add(fHistEventsAfterSel);
 279
 280     fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
 281     fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
 282     fHistTrials->GetYaxis()->SetTitle("trials");
 283     fOutput->Add(fHistTrials);
 284
 285     fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
 286     fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
 287     fHistXsection->GetYaxis()->SetTitle("xsection");
 288     fOutput->Add(fHistXsection);
 289
 290     fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
 291     fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
 292     fHistEvents->GetYaxis()->SetTitle("total events");
 293     fOutput->Add(fHistEvents);
 294
 295     const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 296     const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 297
 298     for (Int_t i = 1; i < 12; i++) {
 299       fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 300       fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 301
 302       fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 303       fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 304       fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 305     }
 306
 307     fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", fNbins*2, fMinBinPt, fMaxBinPt*4);
 308     fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
 309     fHistPtHard->GetYaxis()->SetTitle("counts");
 310     fOutput->Add(fHistPtHard);
 311   }
 312
 313   fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
 314   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
 315   fHistCentrality->GetYaxis()->SetTitle("counts");
 316   fOutput->Add(fHistCentrality);
 317
 318   fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
 319   fHistZVertex->GetXaxis()->SetTitle("z");
 320   fHistZVertex->GetYaxis()->SetTitle("counts");
 321   fOutput->Add(fHistZVertex);
 322
 323   fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
 324   fHistEventPlane->GetXaxis()->SetTitle("event plane");
 325   fHistEventPlane->GetYaxis()->SetTitle("counts");
 326   fOutput->Add(fHistEventPlane);
 327
 328   fHistEventRejection = new TH1I("fHistEventRejection","Reasons to reject event",20,0,20);
 329   fHistEventRejection->Fill("PhysSel",0);
 330   fHistEventRejection->Fill("trigger",0);
 331   fHistEventRejection->Fill("trigTypeSel",0);
 332   fHistEventRejection->Fill("Cent",0);
 333   fHistEventRejection->Fill("vertex contr.",0);
 334   fHistEventRejection->Fill("Vz",0);
 335   fHistEventRejection->Fill("trackInEmcal",0);
 336   fHistEventRejection->Fill("minNTrack",0);
 337   fHistEventRejection->Fill("VtxSel2013pA",0);
 338   fHistEventRejection->Fill("PileUp",0);
 339   fHistEventRejection->Fill("EvtPlane",0);
 340   fHistEventRejection->Fill("SelPtHardBin",0);
 341   fHistEventRejection->GetYaxis()->SetTitle("counts");
 342   fOutput->Add(fHistEventRejection);
 343
 344   PostData(1, fOutput);
 345 }
 346
 347 //________________________________________________________________________
 348 Bool_t AliAnalysisTaskEmcal::FillGeneralHistograms()
 349 {
 350   if (fIsPythia) {
 351     fHistEventsAfterSel->SetBinContent(fPtHardBin + 1, fHistEventsAfterSel->GetBinContent(fPtHardBin + 1) + 1);
 352     fHistTrialsAfterSel->SetBinContent(fPtHardBin + 1, fHistTrialsAfterSel->GetBinContent(fPtHardBin + 1) + fNTrials);
 353     fHistPtHard->Fill(fPtHard);
 354   }
 355
 356   fHistCentrality->Fill(fCent);
 357   fHistZVertex->Fill(fVertex[2]);
 358   fHistEventPlane->Fill(fEPV0);
 359
 360   return kTRUE;
 361 }
 362
 363 //________________________________________________________________________
 364 void AliAnalysisTaskEmcal::UserExec(Option_t *)
 365 {
 366   // Main loop, called for each event.
 367
 368   fMainTriggerPatch = NULL;
 369
 370   if (!fInitialized)
 371     ExecOnce();
 372
 373   if (!fInitialized)
 374     return;
 375
 376   if (!RetrieveEventObjects())
 377     return;
 378
 379   if (!IsEventSelected())
 380     return;
 381
 382   if (fGeneralHistograms && fCreateHisto) {
 383     if (!FillGeneralHistograms())
 384       return;
 385   }
 386
 387   if (!Run())
 388     return;
 389
 390   if (fCreateHisto) {
 391     if (!FillHistograms())
 392       return;
 393   }
 394
 395   if (fCreateHisto && fOutput) {
 396     // information for this iteration of the UserExec in the container
 397     PostData(1, fOutput);
 398   }
 399 }
 400
 401 //________________________________________________________________________
 402 Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Int_t c) const
 403 {
 404   // Return true if cluster is accepted.
 405
 406   if (!clus)
 407     return kFALSE;
 408
 409   AliClusterContainer *cont = GetClusterContainer(c);
 410   if (!cont) {
 411     AliError(Form("%s:Container %d not found",GetName(),c));
 412     return 0;
 413   }
 414
 415   return cont->AcceptCluster(clus);
 416
 417 }
 418
 419 //________________________________________________________________________
 420 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track, Int_t c) const
 421 {
 422   // Return true if track is accepted.
 423
 424   if (!track)
 425     return kFALSE;
 426
 427   AliParticleContainer *cont = GetParticleContainer(c);
 428   if (!cont) {
 429     AliError(Form("%s:Container %d not found",GetName(),c));
 430     return 0;
 431   }
 432
 433   return cont->AcceptParticle(track);
 434 }
 435
 436 //________________________________________________________________________
 437 Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
 438 {
 439   //
 440   // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 441   // Get the pt hard bin from the file path
 442   // This is to called in Notify and should provide the path to the AOD/ESD file
 443   // (Partially copied from AliAnalysisHelperJetTasks)
 444
 445   TString file(currFile);
 446   fXsec = 0;
 447   fTrials = 1;
 448
 449   if (file.Contains(".zip#")) {
 450     Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
 451     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 452     Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
 453     file.Replace(pos+1,pos2-pos1,"");
 454   }
 455   else {
 456     // not an archive take the basename....
 457     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 458   }
 459   AliDebug(1,Form("File name: %s",file.Data()));
 460
 461   // Get the pt hard bin
 462   TString strPthard(file);
 463
 464   strPthard.Remove(strPthard.Last('/'));
 465   strPthard.Remove(strPthard.Last('/'));
 466   if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
 467   strPthard.Remove(0,strPthard.Last('/')+1);
 468   if (strPthard.IsDec())
 469     pthard = strPthard.Atoi();
 470   else
 471     AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
 472
 473   // problem that we cannot really test the existance of a file in a archive so we have to live with open error message from root
 474   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
 475
 476   if (!fxsec) {
 477     // next trial fetch the histgram file
 478     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 479     if (!fxsec) {
 480         // not a severe condition but inciate that we have no information
 481       return kFALSE;
 482     }
 483     else{
 484       // find the tlist we want to be independtent of the name so use the Tkey
 485       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 486       if (!key) {
 487         fxsec->Close();
 488         return kFALSE;
 489       }
 490       TList *list = dynamic_cast<TList*>(key->ReadObj());
 491       if (!list) {
 492         fxsec->Close();
 493         return kFALSE;
 494       }
 495       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 496       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 497       fxsec->Close();
 498     }
 499   } // no tree pyxsec.root
 500   else {
 501     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 502     if (!xtree) {
 503       fxsec->Close();
 504       return kFALSE;
 505     }
 506     UInt_t   ntrials  = 0;
 507     Double_t  xsection  = 0;
 508     xtree->SetBranchAddress("xsection",&xsection);
 509     xtree->SetBranchAddress("ntrials",&ntrials);
 510     xtree->GetEntry(0);
 511     fTrials = ntrials;
 512     fXsec = xsection;
 513     fxsec->Close();
 514   }
 515   return kTRUE;
 516 }
 517
 518 //________________________________________________________________________
 519 Bool_t AliAnalysisTaskEmcal::UserNotify()
 520 {
 521   if (!fIsPythia || !fGeneralHistograms || !fCreateHisto)
 522     return kTRUE;
 523
 524   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 525   if (!tree) {
 526     AliError(Form("%s - UserNotify: No current tree!",GetName()));
 527     return kFALSE;
 528   }
 529
 530   Float_t xsection = 0;
 531   Float_t trials   = 0;
 532   Int_t   pthard   = 0;
 533
 534   TFile *curfile = tree->GetCurrentFile();
 535   if (!curfile) {
 536     AliError(Form("%s - UserNotify: No current file!",GetName()));
 537     return kFALSE;
 538   }
 539
 540   TChain *chain = dynamic_cast<TChain*>(tree);
 541   if (chain)
 542     tree = chain->GetTree();
 543
 544   Int_t nevents = tree->GetEntriesFast();
 545
 546   PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthard);
 547
 548   // TODO: Workaround
 549   if ((pthard < 0) || (pthard > 10))
 550     pthard = 0;
 551   fHistTrials->Fill(pthard, trials);
 552   fHistXsection->Fill(pthard, xsection);
 553   fHistEvents->Fill(pthard, nevents);
 554
 555   return kTRUE;
 556 }
 557
 558 //________________________________________________________________________
 559 void AliAnalysisTaskEmcal::ExecOnce()
 560 {
 561   // Init the analysis.
 562
 563   if (!InputEvent()) {
 564     AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
 565     return;
 566   }
 567
 568   fGeom = AliEMCALGeometry::GetInstance();
 569   if (!fGeom) {
 570     AliError(Form("%s: Can not create geometry", GetName()));
 571     return;
 572   }
 573
 574   if (fEventPlaneVsEmcal >= 0) {
 575     Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
 576     fMinEventPlane = ep - TMath::Pi() / 4;
 577     fMaxEventPlane = ep + TMath::Pi() / 4;
 578   }
 579
 580   //Load all requested track branches - each container knows name already
 581   for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
 582     AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
 583     cont->SetArray(InputEvent());
 584   }
 585
 586   if (fParticleCollArray.GetEntriesFast()>0) {
 587     fTracks = GetParticleArray(0);
 588     if (!fTracks) {
 589       AliError(Form("%s: Could not retrieve first track branch!", GetName()));
 590       return;
 591     }
 592   }
 593
 594   //Load all requested cluster branches - each container knows name already
 595   for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
 596     AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
 597     cont->SetArray(InputEvent());
 598   }
 599
 600   if (fClusterCollArray.GetEntriesFast()>0) {
 601     fCaloClusters = GetClusterArray(0);
 602     if (!fCaloClusters) {
 603       AliError(Form("%s: Could not retrieve first cluster branch!", GetName()));
 604       return;
 605     }
 606   }
 607
 608   if (!fCaloCellsName.IsNull() && !fCaloCells) {
 609     fCaloCells =  dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
 610     if (!fCaloCells) {
 611       AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
 612       return;
 613     }
 614   }
 615
 616   if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
 617     fCaloTriggers =  dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
 618     if (!fCaloTriggers) {
 619       AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
 620       return;
 621     }
 622   }
 623
 624   if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
 625     fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEmcalTriggerPatchInfo");
 626     if (!fTriggerPatchInfo) {
 627       AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
 628       return;
 629     }
 630
 631   }
 632
 633   fInitialized = kTRUE;
 634 }
 635
 636 //_____________________________________________________
 637 AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
 638 {
 639   // Get beam type : pp-AA-pA
 640   // ESDs have it directly, AODs get it from hardcoded run number ranges
 641
 642   if (fForceBeamType != kNA)
 643     return fForceBeamType;
 644
 645   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
 646   if (esd) {
 647     const AliESDRun *run = esd->GetESDRun();
 648     TString beamType = run->GetBeamType();
 649     if (beamType == "p-p")
 650       return kpp;
 651     else if (beamType == "A-A")
 652       return kAA;
 653     else if (beamType == "p-A")
 654       return kpA;
 655     else
 656       return kNA;
 657   } else {
 658     Int_t runNumber = InputEvent()->GetRunNumber();
 659     if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
 660         (runNumber >= 166529 && runNumber <= 170593))    // LHC11h
 661     {
 662       return kAA;
 663     }
 664     else if ((runNumber>=188365 && runNumber <= 188366) || // LHC12g
 665              (runNumber >= 195344 && runNumber <= 196608))  // LHC13b-f
 666     {
 667       return kpA;
 668     } else {
 669       return kpp;
 670     }
 671   }
 672 }
 673
 674 //_____________________________________________________
 675 AliAnalysisTaskEmcal::TriggerType AliAnalysisTaskEmcal::GetTriggerType()
 676 {
 677   // Get trigger type: kND, kJ1, kJ2
 678
 679   if (!fTriggerPatchInfo)
 680     return kND;
 681
 682   //number of patches in event
 683   Int_t nPatch = fTriggerPatchInfo->GetEntries();
 684
 685   //loop over patches to define trigger type of event
 686   Int_t nJ1 = 0;
 687   Int_t nJ2 = 0;
 688   AliEmcalTriggerPatchInfo *patch;
 689   for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
 690     patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
 691     if (patch->IsJetHigh()) nJ1++;
 692     if (patch->IsJetLow())  nJ2++;
 693   }
 694
 695   if (nJ1>0)
 696     return kJ1;
 697   else if (nJ2>0)
 698     return kJ2;
 699   else
 700     return kND;
 701
 702 }
 703
 704 //________________________________________________________________________
 705 Bool_t AliAnalysisTaskEmcal::IsEventSelected()
 706 {
 707   // Check if event is selected
 708
 709   if (fOffTrigger != AliVEvent::kAny) {
 710     UInt_t res = 0;
 711     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
 712     if (eev) {
 713       res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 714     } else {
 715       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
 716       if (aev) {
 717         res = aev->GetHeader()->GetOfflineTrigger();
 718       }
 719     }
 720     if ((res & fOffTrigger) == 0) {
 721       if (fGeneralHistograms)
 722         fHistEventRejection->Fill("PhysSel",1);
 723       return kFALSE;
 724     }
 725   }
 726
 727   if (!fTrigClass.IsNull()) {
 728     TString fired;
 729     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
 730     if (eev) {
 731       fired = eev->GetFiredTriggerClasses();
 732     } else {
 733       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
 734       if (aev) {
 735         fired = aev->GetFiredTriggerClasses();
 736       }
 737     }
 738     if (!fired.Contains("-B-")) {
 739       if (fGeneralHistograms)
 740         fHistEventRejection->Fill("trigger",1);
 741       return kFALSE;
 742     }
 743     TObjArray *arr = fTrigClass.Tokenize("|");
 744     if (!arr) {
 745       if (fGeneralHistograms)
 746         fHistEventRejection->Fill("trigger",1);
 747       return kFALSE;
 748     }
 749     Bool_t match = 0;
 750     for (Int_t i=0;i<arr->GetEntriesFast();++i) {
 751       TObject *obj = arr->At(i);
 752       if (!obj)
 753         continue;
 754       if (fired.Contains(obj->GetName())) {
 755         match = 1;
 756         break;
 757       }
 758     }
 759     delete arr;
 760     if (!match) {
 761       if (fGeneralHistograms)
 762         fHistEventRejection->Fill("trigger",1);
 763       return kFALSE;
 764     }
 765   }
 766
 767   if (fTriggerTypeSel != kND) {
 768     if (fTriggerType != fTriggerTypeSel) {
 769       if (fGeneralHistograms)
 770         fHistEventRejection->Fill("trigTypeSel",1);
 771       return kFALSE;
 772     }
 773   }
 774
 775
 776   if ((fMinCent != -999) && (fMaxCent != -999)) {
 777     if (fCent<fMinCent || fCent>fMaxCent) {
 778       if (fGeneralHistograms)
 779         fHistEventRejection->Fill("Cent",1);
 780       return kFALSE;
 781     }
 782   }
 783
 784   if ((fMinVz != -999) && (fMaxVz != -999)) {
 785     if (fNVertCont == 0 ) {
 786       if (fGeneralHistograms)
 787         fHistEventRejection->Fill("vertex contr.",1);
 788       return kFALSE;
 789     }
 790     Double_t vz = fVertex[2];
 791     if (vz<fMinVz || vz>fMaxVz) {
 792       if (fGeneralHistograms)
 793         fHistEventRejection->Fill("Vz",1);
 794       return kFALSE;
 795     }
 796   }
 797
 798   if (fMinPtTrackInEmcal > 0 && fGeom) {
 799     Bool_t trackInEmcalOk = kFALSE;
 800     Int_t ntracks = GetNParticles(0);
 801     for (Int_t i = 0; i < ntracks; i++) {
 802       AliVParticle *track = GetAcceptParticleFromArray(i,0);
 803       if (!track)
 804         continue;
 805
 806       Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 807       Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 808       Int_t runNumber = InputEvent()->GetRunNumber();
 809       if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
 810         phiMin = 1.4;
 811         phiMax = TMath::Pi();
 812       }
 813
 814       if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
 815         continue;
 816       if (track->Pt() > fMinPtTrackInEmcal) {
 817         trackInEmcalOk = kTRUE;
 818         break;
 819       }
 820     }
 821     if (!trackInEmcalOk) {
 822       if (fGeneralHistograms)
 823         fHistEventRejection->Fill("trackInEmcal",1);
 824       return kFALSE;
 825     }
 826   }
 827
 828   if (fMinNTrack > 0) {
 829     Int_t nTracksAcc = 0;
 830     Int_t ntracks = GetNParticles(0);
 831     for (Int_t i = 0; i < ntracks; i++) {
 832       AliVParticle *track = GetAcceptParticleFromArray(i,0);
 833       if (!track)
 834         continue;
 835       if (track->Pt() > fTrackPtCut) {
 836         nTracksAcc++;
 837         if (nTracksAcc>=fMinNTrack)
 838           break;
 839       }
 840     }
 841     if (nTracksAcc<fMinNTrack) {
 842       if (fGeneralHistograms)
 843         fHistEventRejection->Fill("minNTrack",1);
 844       return kFALSE;
 845     }
 846   }
 847
 848   if (fUseAliAnaUtils) {
 849     if (!fAliAnalysisUtils)
 850       fAliAnalysisUtils = new AliAnalysisUtils();
 851     fAliAnalysisUtils->SetMinVtxContr(2);
 852     fAliAnalysisUtils->SetMaxVtxZ(10.);
 853
 854     if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
 855       if (fGeneralHistograms)
 856         fHistEventRejection->Fill("VtxSel2013pA",1);
 857       return kFALSE;
 858     }
 859
 860     if (fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
 861       fHistEventRejection->Fill("PileUp",1);
 862       return kFALSE;
 863     }
 864   }
 865
 866   if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
 867       !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
 868       !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane))
 869     {
 870       if (fGeneralHistograms)
 871         fHistEventRejection->Fill("EvtPlane",1);
 872       return kFALSE;
 873     }
 874
 875   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)  {
 876       if (fGeneralHistograms)
 877         fHistEventRejection->Fill("SelPtHardBin",1);
 878       return kFALSE;
 879     }
 880
 881   return kTRUE;
 882 }
 883
 884 //________________________________________________________________________
 885 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
 886 {
 887   // Get array from event.
 888
 889   TClonesArray *arr = 0;
 890   TString sname(name);
 891   if (!sname.IsNull()) {
 892     arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
 893     if (!arr) {
 894       AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
 895       return 0;
 896     }
 897   } else {
 898     return 0;
 899   }
 900
 901   if (!clname)
 902     return arr;
 903
 904   TString objname(arr->GetClass()->GetName());
 905   TClass cls(objname);
 906   if (!cls.InheritsFrom(clname)) {
 907     AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
 908                     GetName(), cls.GetName(), name, clname));
 909     return 0;
 910   }
 911   return arr;
 912 }
 913
 914 //________________________________________________________________________
 915 Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
 916 {
 917   // Retrieve objects from event.
 918
 919   fVertex[0] = 0;
 920   fVertex[1] = 0;
 921   fVertex[2] = 0;
 922   fNVertCont = 0;
 923
 924   const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
 925   if (vert) {
 926     vert->GetXYZ(fVertex);
 927     fNVertCont = vert->GetNContributors();
 928   }
 929
 930   fBeamType = GetBeamType();
 931
 932   if (fBeamType == kAA || fBeamType == kpA ) {
 933     AliCentrality *aliCent = InputEvent()->GetCentrality();
 934     if (aliCent) {
 935       fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
 936       if (fNcentBins==4) {
 937         if      (fCent >=  0 && fCent <   10) fCentBin = 0;
 938         else if (fCent >= 10 && fCent <   30) fCentBin = 1;
 939         else if (fCent >= 30 && fCent <   50) fCentBin = 2;
 940         else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
 941         else {
 942           AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
 943           fCentBin = fNcentBins-1;
 944         }
 945       }
 946       else {
 947         Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
 948         fCentBin = TMath::FloorNint(fCent/centWidth);
 949         if (fCentBin>=fNcentBins) {
 950           AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
 951           fCentBin = fNcentBins-1;
 952         }
 953       }
 954     } else {
 955       AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
 956       fCentBin = 3;
 957     }
 958     AliEventplane *aliEP = InputEvent()->GetEventplane();
 959     if (aliEP) {
 960       fEPV0  = aliEP->GetEventplane("V0" ,InputEvent());
 961       fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
 962       fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
 963     } else {
 964       AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
 965     }
 966   } else {
 967     fCent = 99;
 968     fCentBin = 0;
 969   }
 970
 971   if (fIsPythia) {
 972
 973     if (MCEvent()) {
 974       fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 975       if (!fPythiaHeader) {
 976         // Check if AOD
 977         AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 978
 979         if (aodMCH) {
 980           for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
 981             fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
 982             if (fPythiaHeader) break;
 983           }
 984         }
 985       }
 986     }
 987
 988     if (fPythiaHeader) {
 989       fPtHard = fPythiaHeader->GetPtHard();
 990
 991       const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 992       const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 993       for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
 994         if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
 995           break;
 996       }
 997
 998       fNTrials = fPythiaHeader->Trials();
 999     }
1000   }
1001
1002   fTriggerType = GetTriggerType();
1003
1004   return kTRUE;
1005 }
1006
1007 //________________________________________________________________________
1008 AliParticleContainer* AliAnalysisTaskEmcal::AddParticleContainer(const char *n) {
1009
1010   // Add particle container
1011   // will be called in AddTask macro
1012
1013   TString tmp = TString(n);
1014   if (tmp.IsNull()) return 0;
1015
1016   AliParticleContainer *cont = 0x0;
1017   cont = new AliParticleContainer();
1018   cont->SetArrayName(n);
1019   TString contName = cont->GetArrayName();
1020
1021   fParticleCollArray.Add(cont);
1022
1023   return cont;
1024 }
1025
1026 //________________________________________________________________________
1027 AliClusterContainer* AliAnalysisTaskEmcal::AddClusterContainer(const char *n) {
1028
1029   // Add cluster container
1030   // will be called in AddTask macro
1031
1032   TString tmp = TString(n);
1033   if (tmp.IsNull()) return 0;
1034
1035   AliClusterContainer *cont = 0x0;
1036   cont = new AliClusterContainer();
1037   cont->SetArrayName(n);
1038
1039   fClusterCollArray.Add(cont);
1040
1041   return cont;
1042 }
1043
1044 //________________________________________________________________________
1045 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(Int_t i) const {
1046   // Get i^th particle container
1047
1048   if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1049   AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1050   return cont;
1051 }
1052
1053 //________________________________________________________________________
1054 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(Int_t i) const {
1055   // Get i^th cluster container
1056
1057   if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1058   AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1059   return cont;
1060 }
1061
1062 //________________________________________________________________________
1063 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(const char *name) const {
1064   // Get particle container with name
1065
1066   AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1067   return cont;
1068 }
1069
1070 //________________________________________________________________________
1071 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(const char *name) const {
1072   // Get cluster container with name
1073
1074   AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1075   return cont;
1076 }
1077
1078 //________________________________________________________________________
1079 TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const {
1080   // Get i^th TClonesArray with AliVParticle
1081
1082   AliParticleContainer *cont = GetParticleContainer(i);
1083   if (!cont) {
1084     AliError(Form("%s: Particle container %d not found",GetName(),i));
1085     return 0;
1086   }
1087   TString contName = cont->GetArrayName();
1088   return cont->GetArray();
1089 }
1090
1091 //________________________________________________________________________
1092 TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const {
1093   // Get i^th TClonesArray with AliVCluster
1094
1095   AliClusterContainer *cont = GetClusterContainer(i);
1096   if (!cont) {
1097     AliError(Form("%s:Cluster container %d not found",GetName(),i));
1098     return 0;
1099   }
1100   return cont->GetArray();
1101 }
1102
1103 //________________________________________________________________________
1104 AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const {
1105   // Get particle p if accepted from  container c
1106   // If particle not accepted return 0
1107
1108   AliParticleContainer *cont = GetParticleContainer(c);
1109   if (!cont) {
1110     AliError(Form("%s: Particle container %d not found",GetName(),c));
1111     return 0;
1112   }
1113   AliVParticle *vp = cont->GetAcceptParticle(p);
1114
1115   return vp;
1116 }
1117
1118 //________________________________________________________________________
1119 AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const {
1120   // Get particle p if accepted from  container c
1121   // If particle not accepted return 0
1122
1123   AliClusterContainer *cont = GetClusterContainer(c);
1124   if (!cont) {
1125     AliError(Form("%s: Cluster container %d not found",GetName(),c));
1126     return 0;
1127   }
1128   AliVCluster *vc = cont->GetAcceptCluster(cl);
1129
1130   return vc;
1131 }
1132
1133 //________________________________________________________________________
1134 Int_t AliAnalysisTaskEmcal::GetNParticles(Int_t i) const {
1135   // Get number of entries in particle array i
1136
1137   AliParticleContainer *cont = GetParticleContainer(i);
1138   if (!cont) {
1139     AliError(Form("%s: Particle container %d not found",GetName(),i));
1140     return 0;
1141   }
1142   return cont->GetNEntries();
1143 }
1144
1145 //________________________________________________________________________
1146 Int_t AliAnalysisTaskEmcal::GetNClusters(Int_t i) const {
1147   // Get number of entries in cluster array i
1148
1149   AliClusterContainer *cont = GetClusterContainer(i);
1150   if (!cont) {
1151     AliError(Form("%s: Cluster container %d not found",GetName(),i));
1152     return 0;
1153   }
1154   return cont->GetNEntries();
1155 }
1156
1157 //________________________________________________________________________
1158 AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch() {
1159   //get main trigger match; if not known yet, look for it and cache
1160
1161   if (fMainTriggerPatch)
1162     return fMainTriggerPatch;
1163
1164   if (!fTriggerPatchInfo) {
1165     AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1166     return 0;
1167   }
1168
1169   //number of patches in event
1170   Int_t nPatch = fTriggerPatchInfo->GetEntries();
1171
1172   //extract main trigger patch
1173   AliEmcalTriggerPatchInfo *patch;
1174   for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1175
1176     patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1177     if (patch->IsMainTrigger()) {
1178       fMainTriggerPatch = patch;
1179       break;
1180     }
1181   }
1182
1183   return fMainTriggerPatch;
1184 }