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 ((fMinVz != -999) && (fMaxVz != -999)) {
 778     if (fNVertCont == 0 ) {
 779       if (fGeneralHistograms)
 780         fHistEventRejection->Fill("vertex contr.",1);
 781       return kFALSE;
 782     }
 783     Double_t vz = fVertex[2];
 784     if (vz<fMinVz || vz>fMaxVz) {
 785       if (fGeneralHistograms)
 786         fHistEventRejection->Fill("Vz",1);
 787       return kFALSE;
 788     }
 789   }
 790
 791   if (fMinPtTrackInEmcal > 0 && fGeom) {
 792     Bool_t trackInEmcalOk = kFALSE;
 793     Int_t ntracks = GetNParticles(0);
 794     for (Int_t i = 0; i < ntracks; i++) {
 795       AliVParticle *track = GetAcceptParticleFromArray(i,0);
 796       if (!track)
 797         continue;
 798
 799       Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 800       Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 801       Int_t runNumber = InputEvent()->GetRunNumber();
 802       if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
 803         phiMin = 1.4;
 804         phiMax = TMath::Pi();
 805       }
 806
 807       if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
 808         continue;
 809       if (track->Pt() > fMinPtTrackInEmcal) {
 810         trackInEmcalOk = kTRUE;
 811         break;
 812       }
 813     }
 814     if (!trackInEmcalOk) {
 815       if (fGeneralHistograms)
 816         fHistEventRejection->Fill("trackInEmcal",1);
 817       return kFALSE;
 818     }
 819   }
 820
 821   if (fMinNTrack > 0) {
 822     Int_t nTracksAcc = 0;
 823     Int_t ntracks = GetNParticles(0);
 824     for (Int_t i = 0; i < ntracks; i++) {
 825       AliVParticle *track = GetAcceptParticleFromArray(i,0);
 826       if (!track)
 827         continue;
 828       if (track->Pt() > fTrackPtCut) {
 829         nTracksAcc++;
 830         if (nTracksAcc>=fMinNTrack)
 831           break;
 832       }
 833     }
 834     if (nTracksAcc<fMinNTrack) {
 835       if (fGeneralHistograms)
 836         fHistEventRejection->Fill("minNTrack",1);
 837       return kFALSE;
 838     }
 839   }
 840
 841   if (fUseAliAnaUtils) {
 842     if (!fAliAnalysisUtils)
 843       fAliAnalysisUtils = new AliAnalysisUtils();
 844     fAliAnalysisUtils->SetMinVtxContr(2);
 845     fAliAnalysisUtils->SetMaxVtxZ(10.);
 846
 847     if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
 848       if (fGeneralHistograms)
 849         fHistEventRejection->Fill("VtxSel2013pA",1);
 850       return kFALSE;
 851     }
 852
 853     if (fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
 854       fHistEventRejection->Fill("PileUp",1);
 855       return kFALSE;
 856     }
 857   }
 858
 859   if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
 860       !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
 861       !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane))
 862     {
 863       if (fGeneralHistograms)
 864         fHistEventRejection->Fill("EvtPlane",1);
 865       return kFALSE;
 866     }
 867
 868   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)  {
 869       if (fGeneralHistograms)
 870         fHistEventRejection->Fill("SelPtHardBin",1);
 871       return kFALSE;
 872     }
 873
 874   return kTRUE;
 875 }
 876
 877 //________________________________________________________________________
 878 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
 879 {
 880   // Get array from event.
 881
 882   TClonesArray *arr = 0;
 883   TString sname(name);
 884   if (!sname.IsNull()) {
 885     arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
 886     if (!arr) {
 887       AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
 888       return 0;
 889     }
 890   } else {
 891     return 0;
 892   }
 893
 894   if (!clname)
 895     return arr;
 896
 897   TString objname(arr->GetClass()->GetName());
 898   TClass cls(objname);
 899   if (!cls.InheritsFrom(clname)) {
 900     AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
 901                     GetName(), cls.GetName(), name, clname));
 902     return 0;
 903   }
 904   return arr;
 905 }
 906
 907 //________________________________________________________________________
 908 Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
 909 {
 910   // Retrieve objects from event.
 911
 912   fVertex[0] = 0;
 913   fVertex[1] = 0;
 914   fVertex[2] = 0;
 915   fNVertCont = 0;
 916
 917   const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
 918   if (vert) {
 919     vert->GetXYZ(fVertex);
 920     fNVertCont = vert->GetNContributors();
 921   }
 922
 923   fBeamType = GetBeamType();
 924
 925   if (fBeamType == kAA || fBeamType == kpA ) {
 926     AliCentrality *aliCent = InputEvent()->GetCentrality();
 927     if (aliCent) {
 928       fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
 929       if (fNcentBins==4) {
 930         if      (fCent >=  0 && fCent <   10) fCentBin = 0;
 931         else if (fCent >= 10 && fCent <   30) fCentBin = 1;
 932         else if (fCent >= 30 && fCent <   50) fCentBin = 2;
 933         else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
 934         else {
 935           AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
 936           fCentBin = fNcentBins-1;
 937         }
 938       } else {
 939         Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
 940         fCentBin = TMath::FloorNint(fCent/centWidth);
 941         if (fCentBin>=fNcentBins) {
 942           AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
 943           fCentBin = fNcentBins-1;
 944         }
 945       }
 946     } else {
 947       AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
 948       fCentBin = 3;
 949     }
 950     AliEventplane *aliEP = InputEvent()->GetEventplane();
 951     if (aliEP) {
 952       fEPV0  = aliEP->GetEventplane("V0" ,InputEvent());
 953       fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
 954       fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
 955     } else {
 956       AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
 957     }
 958   } else {
 959     fCent = 99;
 960     fCentBin = 0;
 961   }
 962
 963   if (fIsPythia) {
 964
 965     if (MCEvent()) {
 966       fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 967       if (!fPythiaHeader) {
 968         // Check if AOD
 969         AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 970
 971         if (aodMCH) {
 972           for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
 973             fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
 974             if (fPythiaHeader) break;
 975           }
 976         }
 977       }
 978     }
 979
 980     if (fPythiaHeader) {
 981       fPtHard = fPythiaHeader->GetPtHard();
 982
 983       const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 984       const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 985       for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
 986         if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
 987           break;
 988       }
 989
 990       fNTrials = fPythiaHeader->Trials();
 991     }
 992   }
 993
 994   fTriggerType = GetTriggerType();
 995
 996   return kTRUE;
 997 }
 998
 999 //________________________________________________________________________
1000 AliParticleContainer* AliAnalysisTaskEmcal::AddParticleContainer(const char *n)
1001 {
1002   // Add particle container
1003   // will be called in AddTask macro
1004
1005   TString tmp = TString(n);
1006   if (tmp.IsNull()) return 0;
1007
1008   AliParticleContainer *cont = 0x0;
1009   cont = new AliParticleContainer();
1010   cont->SetArrayName(n);
1011   TString contName = cont->GetArrayName();
1012
1013   fParticleCollArray.Add(cont);
1014
1015   return cont;
1016 }
1017
1018 //________________________________________________________________________
1019 AliClusterContainer* AliAnalysisTaskEmcal::AddClusterContainer(const char *n)
1020 {
1021   // Add cluster container
1022   // will be called in AddTask macro
1023
1024   TString tmp = TString(n);
1025   if (tmp.IsNull()) return 0;
1026
1027   AliClusterContainer *cont = 0x0;
1028   cont = new AliClusterContainer();
1029   cont->SetArrayName(n);
1030
1031   fClusterCollArray.Add(cont);
1032
1033   return cont;
1034 }
1035
1036 //________________________________________________________________________
1037 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(Int_t i) const
1038 {
1039   // Get i^th particle container
1040
1041   if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1042   AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1043   return cont;
1044 }
1045
1046 //________________________________________________________________________
1047 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(Int_t i) const
1048 {
1049   // Get i^th cluster container
1050
1051   if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1052   AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1053   return cont;
1054 }
1055
1056 //________________________________________________________________________
1057 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(const char *name) const
1058 {
1059   // Get particle container with name
1060
1061   AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1062   return cont;
1063 }
1064
1065 //________________________________________________________________________
1066 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(const char *name) const
1067 {
1068   // Get cluster container with name
1069
1070   AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1071   return cont;
1072 }
1073
1074 //________________________________________________________________________
1075 TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const
1076 {
1077   // Get i^th TClonesArray with AliVParticle
1078
1079   AliParticleContainer *cont = GetParticleContainer(i);
1080   if (!cont) {
1081     AliError(Form("%s: Particle container %d not found",GetName(),i));
1082     return 0;
1083   }
1084   TString contName = cont->GetArrayName();
1085   return cont->GetArray();
1086 }
1087
1088 //________________________________________________________________________
1089 TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const
1090 {
1091   // Get i^th TClonesArray with AliVCluster
1092
1093   AliClusterContainer *cont = GetClusterContainer(i);
1094   if (!cont) {
1095     AliError(Form("%s:Cluster container %d not found",GetName(),i));
1096     return 0;
1097   }
1098   return cont->GetArray();
1099 }
1100
1101 //________________________________________________________________________
1102 AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const
1103 {
1104   // Get particle p if accepted from  container c
1105   // If particle not accepted return 0
1106
1107   AliParticleContainer *cont = GetParticleContainer(c);
1108   if (!cont) {
1109     AliError(Form("%s: Particle container %d not found",GetName(),c));
1110     return 0;
1111   }
1112   AliVParticle *vp = cont->GetAcceptParticle(p);
1113
1114   return vp;
1115 }
1116
1117 //________________________________________________________________________
1118 AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const
1119 {
1120   // Get particle p if accepted from  container c
1121   // If particle not accepted return 0
1122
1123   AliClusterContainer *cont = GetClusterContainer(c);
1124   if (!cont) {
1125     AliError(Form("%s: Cluster container %d not found",GetName(),c));
1126     return 0;
1127   }
1128   AliVCluster *vc = cont->GetAcceptCluster(cl);
1129
1130   return vc;
1131 }
1132
1133 //________________________________________________________________________
1134 Int_t AliAnalysisTaskEmcal::GetNParticles(Int_t i) const
1135 {
1136   // Get number of entries in particle array i
1137
1138   AliParticleContainer *cont = GetParticleContainer(i);
1139   if (!cont) {
1140     AliError(Form("%s: Particle container %d not found",GetName(),i));
1141     return 0;
1142   }
1143   return cont->GetNEntries();
1144 }
1145
1146 //________________________________________________________________________
1147 Int_t AliAnalysisTaskEmcal::GetNClusters(Int_t i) const
1148 {
1149   // Get number of entries in cluster array i
1150
1151   AliClusterContainer *cont = GetClusterContainer(i);
1152   if (!cont) {
1153     AliError(Form("%s: Cluster container %d not found",GetName(),i));
1154     return 0;
1155   }
1156   return cont->GetNEntries();
1157 }
1158
1159 //________________________________________________________________________
1160 AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch()
1161 {
1162   //get main trigger match; if not known yet, look for it and cache
1163
1164   if (fMainTriggerPatch)
1165     return fMainTriggerPatch;
1166
1167   if (!fTriggerPatchInfo) {
1168     AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1169     return 0;
1170   }
1171
1172   //number of patches in event
1173   Int_t nPatch = fTriggerPatchInfo->GetEntries();
1174
1175   //extract main trigger patch
1176   AliEmcalTriggerPatchInfo *patch;
1177   for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1178
1179     patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1180     if (patch->IsMainTrigger()) {
1181       fMainTriggerPatch = patch;
1182       break;
1183     }
1184   }
1185
1186   return fMainTriggerPatch;
1187 }
1188
1189 //________________________________________________________________________
1190 void AliAnalysisTaskEmcal::AddObjectToEvent(TObject *obj)
1191 {
1192   // Add object to event
1193
1194   if (!(InputEvent()->FindListObject(obj->GetName()))) {
1195     InputEvent()->AddObject(obj);
1196   } else {
1197     AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1198   }
1199 }