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