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