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