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