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