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