PWG/EMCAL/AliAnalysisTaskEmcal.cxx

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