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