PWG/EMCAL/AliAnalysisTaskEmcal.cxx

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