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     fXsection = fPythiaHeader->GetXsection();
 357     if(fXsection>0.) {
 358       fHistXsection->Fill(fPtHardBin+1, fXsection);
 359       fHistTrials->Fill(fPtHardBin+1, fPythiaHeader->Trials());
 360     }
 361   }
 362
 363   fHistCentrality->Fill(fCent);
 364   fHistZVertex->Fill(fVertex[2]);
 365   fHistEventPlane->Fill(fEPV0);
 366
 367   return kTRUE;
 368 }
 369
 370 //________________________________________________________________________
 371 void AliAnalysisTaskEmcal::UserExec(Option_t *)
 372 {
 373   // Main loop, called for each event.
 374
 375   fMainTriggerPatch = NULL;
 376
 377   if (!fInitialized)
 378     ExecOnce();
 379
 380   if (!fInitialized)
 381     return;
 382
 383   if (!RetrieveEventObjects())
 384     return;
 385
 386   if (!IsEventSelected())
 387     return;
 388
 389   if (fGeneralHistograms && fCreateHisto) {
 390     if (!FillGeneralHistograms())
 391       return;
 392   }
 393
 394   if (!Run())
 395     return;
 396
 397   if (fCreateHisto) {
 398     if (!FillHistograms())
 399       return;
 400   }
 401
 402   if (fCreateHisto && fOutput) {
 403     // information for this iteration of the UserExec in the container
 404     PostData(1, fOutput);
 405   }
 406 }
 407
 408 //________________________________________________________________________
 409 Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Int_t c) const
 410 {
 411   // Return true if cluster is accepted.
 412
 413   if (!clus)
 414     return kFALSE;
 415
 416   AliClusterContainer *cont = GetClusterContainer(c);
 417   if (!cont) {
 418     AliError(Form("%s:Container %d not found",GetName(),c));
 419     return 0;
 420   }
 421
 422   return cont->AcceptCluster(clus);
 423 }
 424
 425 //________________________________________________________________________
 426 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track, Int_t c) const
 427 {
 428   // Return true if track is accepted.
 429
 430   if (!track)
 431     return kFALSE;
 432
 433   AliParticleContainer *cont = GetParticleContainer(c);
 434   if (!cont) {
 435     AliError(Form("%s:Container %d not found",GetName(),c));
 436     return 0;
 437   }
 438
 439   return cont->AcceptParticle(track);
 440 }
 441
 442 //________________________________________________________________________
 443 Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
 444 {
 445   //
 446   // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 447   // Get the pt hard bin from the file path
 448   // This is to called in Notify and should provide the path to the AOD/ESD file
 449   // (Partially copied from AliAnalysisHelperJetTasks)
 450
 451   TString file(currFile);
 452   fXsec = 0;
 453   fTrials = 1;
 454
 455   if (file.Contains(".zip#")) {
 456     Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
 457     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 458     Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
 459     file.Replace(pos+1,pos2-pos1,"");
 460   } else {
 461     // not an archive take the basename....
 462     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 463   }
 464   AliDebug(1,Form("File name: %s",file.Data()));
 465
 466   // Get the pt hard bin
 467   TString strPthard(file);
 468
 469   strPthard.Remove(strPthard.Last('/'));
 470   strPthard.Remove(strPthard.Last('/'));
 471   if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
 472   strPthard.Remove(0,strPthard.Last('/')+1);
 473   if (strPthard.IsDec())
 474     pthard = strPthard.Atoi();
 475   else
 476     AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
 477
 478   // 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
 479   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
 480
 481   if (!fxsec) {
 482     // next trial fetch the histgram file
 483     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 484     if (!fxsec) {
 485         // not a severe condition but inciate that we have no information
 486       return kFALSE;
 487     } else {
 488       // find the tlist we want to be independtent of the name so use the Tkey
 489       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 490       if (!key) {
 491         fxsec->Close();
 492         return kFALSE;
 493       }
 494       TList *list = dynamic_cast<TList*>(key->ReadObj());
 495       if (!list) {
 496         fxsec->Close();
 497         return kFALSE;
 498       }
 499       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 500       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 501       fxsec->Close();
 502     }
 503   } else { // no tree pyxsec.root
 504     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 505     if (!xtree) {
 506       fxsec->Close();
 507       return kFALSE;
 508     }
 509     UInt_t   ntrials  = 0;
 510     Double_t  xsection  = 0;
 511     xtree->SetBranchAddress("xsection",&xsection);
 512     xtree->SetBranchAddress("ntrials",&ntrials);
 513     xtree->GetEntry(0);
 514     fTrials = ntrials;
 515     fXsec = xsection;
 516     fxsec->Close();
 517   }
 518   return kTRUE;
 519 }
 520
 521 //________________________________________________________________________
 522 Bool_t AliAnalysisTaskEmcal::UserNotify()
 523 {
 524   // Called when file changes.
 525
 526   if (!fIsPythia || !fGeneralHistograms || !fCreateHisto)
 527     return kTRUE;
 528
 529   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 530   if (!tree) {
 531     AliError(Form("%s - UserNotify: No current tree!",GetName()));
 532     return kFALSE;
 533   }
 534
 535   Float_t trials   = 0;
 536   Int_t   pthard   = 0;
 537
 538   TFile *curfile = tree->GetCurrentFile();
 539   if (!curfile) {
 540     AliError(Form("%s - UserNotify: No current file!",GetName()));
 541     return kFALSE;
 542   }
 543
 544   TChain *chain = dynamic_cast<TChain*>(tree);
 545   if (chain)
 546     tree = chain->GetTree();
 547
 548   Int_t nevents = tree->GetEntriesFast();
 549
 550   PythiaInfoFromFile(curfile->GetName(), fXsection, trials, pthard);
 551
 552   // TODO: Workaround
 553   if ((pthard < 0) || (pthard > 10))
 554     pthard = 0;
 555   fHistTrials->Fill(pthard, trials);
 556   fHistXsection->Fill(pthard, fXsection);
 557   fHistEvents->Fill(pthard, nevents);
 558
 559   return kTRUE;
 560 }
 561
 562 //________________________________________________________________________
 563 void AliAnalysisTaskEmcal::ExecOnce()
 564 {
 565   // Init the analysis.
 566
 567   if (!InputEvent()) {
 568     AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
 569     return;
 570   }
 571
 572   fGeom = AliEMCALGeometry::GetInstance();
 573   if (!fGeom) {
 574     AliError(Form("%s: Can not create geometry", GetName()));
 575     return;
 576   }
 577
 578   if (fEventPlaneVsEmcal >= 0) {
 579     Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
 580     fMinEventPlane = ep - TMath::Pi() / 4;
 581     fMaxEventPlane = ep + TMath::Pi() / 4;
 582   }
 583
 584   //Load all requested track branches - each container knows name already
 585   for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
 586     AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
 587     cont->SetArray(InputEvent());
 588   }
 589
 590   if (fParticleCollArray.GetEntriesFast()>0) {
 591     fTracks = GetParticleArray(0);
 592     if (!fTracks) {
 593       AliError(Form("%s: Could not retrieve first track branch!", GetName()));
 594       return;
 595     }
 596   }
 597
 598   //Load all requested cluster branches - each container knows name already
 599   for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
 600     AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
 601     cont->SetArray(InputEvent());
 602   }
 603
 604   if (fClusterCollArray.GetEntriesFast()>0) {
 605     fCaloClusters = GetClusterArray(0);
 606     if (!fCaloClusters) {
 607       AliError(Form("%s: Could not retrieve first cluster branch!", GetName()));
 608       return;
 609     }
 610   }
 611
 612   if (!fCaloCellsName.IsNull() && !fCaloCells) {
 613     fCaloCells =  dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
 614     if (!fCaloCells) {
 615       AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
 616       return;
 617     }
 618   }
 619
 620   if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
 621     fCaloTriggers =  dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
 622     if (!fCaloTriggers) {
 623       AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
 624       return;
 625     }
 626   }
 627
 628   if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
 629     fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEmcalTriggerPatchInfo");
 630     if (!fTriggerPatchInfo) {
 631       AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
 632       return;
 633     }
 634
 635   }
 636
 637   fInitialized = kTRUE;
 638 }
 639
 640 //_____________________________________________________
 641 AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
 642 {
 643   // Get beam type : pp-AA-pA
 644   // ESDs have it directly, AODs get it from hardcoded run number ranges
 645
 646   if (fForceBeamType != kNA)
 647     return fForceBeamType;
 648
 649   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
 650   if (esd) {
 651     const AliESDRun *run = esd->GetESDRun();
 652     TString beamType = run->GetBeamType();
 653     if (beamType == "p-p")
 654       return kpp;
 655     else if (beamType == "A-A")
 656       return kAA;
 657     else if (beamType == "p-A")
 658       return kpA;
 659     else
 660       return kNA;
 661   } else {
 662     Int_t runNumber = InputEvent()->GetRunNumber();
 663     if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
 664         (runNumber >= 166529 && runNumber <= 170593)) {  // LHC11h
 665       return kAA;
 666     } else if ((runNumber>=188365 && runNumber <= 188366) ||   // LHC12g
 667                (runNumber >= 195344 && runNumber <= 196608)) { // LHC13b-f
 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 Bool_t AliAnalysisTaskEmcal::IsEventSelected()
 706 {
 707   // Check if event is selected
 708
 709   if (fOffTrigger != AliVEvent::kAny) {
 710     UInt_t res = 0;
 711     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
 712     if (eev) {
 713       res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 714     } else {
 715       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
 716       if (aev) {
 717         res = aev->GetHeader()->GetOfflineTrigger();
 718       }
 719     }
 720     if ((res & fOffTrigger) == 0) {
 721       if (fGeneralHistograms)
 722         fHistEventRejection->Fill("PhysSel",1);
 723       return kFALSE;
 724     }
 725   }
 726
 727   if (!fTrigClass.IsNull()) {
 728     TString fired;
 729     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
 730     if (eev) {
 731       fired = eev->GetFiredTriggerClasses();
 732     } else {
 733       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
 734       if (aev) {
 735         fired = aev->GetFiredTriggerClasses();
 736       }
 737     }
 738     if (!fired.Contains("-B-")) {
 739       if (fGeneralHistograms)
 740         fHistEventRejection->Fill("trigger",1);
 741       return kFALSE;
 742     }
 743     TObjArray *arr = fTrigClass.Tokenize("|");
 744     if (!arr) {
 745       if (fGeneralHistograms)
 746         fHistEventRejection->Fill("trigger",1);
 747       return kFALSE;
 748     }
 749     Bool_t match = 0;
 750     for (Int_t i=0;i<arr->GetEntriesFast();++i) {
 751       TObject *obj = arr->At(i);
 752       if (!obj)
 753         continue;
 754       if (fired.Contains(obj->GetName())) {
 755         match = 1;
 756         break;
 757       }
 758     }
 759     delete arr;
 760     if (!match) {
 761       if (fGeneralHistograms)
 762         fHistEventRejection->Fill("trigger",1);
 763       return kFALSE;
 764     }
 765   }
 766
 767   if (fTriggerTypeSel != kND) {
 768     if (fTriggerType != fTriggerTypeSel) {
 769       if (fGeneralHistograms)
 770         fHistEventRejection->Fill("trigTypeSel",1);
 771       return kFALSE;
 772     }
 773   }
 774
 775   if ((fMinCent != -999) && (fMaxCent != -999)) {
 776     if (fCent<fMinCent || fCent>fMaxCent) {
 777       if (fGeneralHistograms)
 778         fHistEventRejection->Fill("Cent",1);
 779       return kFALSE;
 780     }
 781   }
 782
 783   if (fUseAliAnaUtils) {
 784     if (!fAliAnalysisUtils)
 785       fAliAnalysisUtils = new AliAnalysisUtils();
 786     fAliAnalysisUtils->SetMinVtxContr(2);
 787     fAliAnalysisUtils->SetMaxVtxZ(999);
 788     if(fMinVz<-10.) fMinVz = -10.;
 789     if(fMinVz>10.)  fMaxVz = 10.;
 790
 791     if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
 792       if (fGeneralHistograms)
 793         fHistEventRejection->Fill("VtxSel2013pA",1);
 794       return kFALSE;
 795     }
 796
 797     if (fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
 798       fHistEventRejection->Fill("PileUp",1);
 799       return kFALSE;
 800     }
 801   }
 802
 803   if ((fMinVz != -999) && (fMaxVz != -999)) {
 804     if (fNVertCont == 0 ) {
 805       if (fGeneralHistograms)
 806         fHistEventRejection->Fill("vertex contr.",1);
 807       return kFALSE;
 808     }
 809     Double_t vz = fVertex[2];
 810     if (vz<fMinVz || vz>fMaxVz) {
 811       if (fGeneralHistograms)
 812         fHistEventRejection->Fill("Vz",1);
 813       return kFALSE;
 814     }
 815   }
 816
 817   if (fMinPtTrackInEmcal > 0 && fGeom) {
 818     Bool_t trackInEmcalOk = kFALSE;
 819     Int_t ntracks = GetNParticles(0);
 820     for (Int_t i = 0; i < ntracks; i++) {
 821       AliVParticle *track = GetAcceptParticleFromArray(i,0);
 822       if (!track)
 823         continue;
 824
 825       Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 826       Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 827       Int_t runNumber = InputEvent()->GetRunNumber();
 828       if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
 829         phiMin = 1.4;
 830         phiMax = TMath::Pi();
 831       }
 832
 833       if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
 834         continue;
 835       if (track->Pt() > fMinPtTrackInEmcal) {
 836         trackInEmcalOk = kTRUE;
 837         break;
 838       }
 839     }
 840     if (!trackInEmcalOk) {
 841       if (fGeneralHistograms)
 842         fHistEventRejection->Fill("trackInEmcal",1);
 843       return kFALSE;
 844     }
 845   }
 846
 847   if (fMinNTrack > 0) {
 848     Int_t nTracksAcc = 0;
 849     Int_t ntracks = GetNParticles(0);
 850     for (Int_t i = 0; i < ntracks; i++) {
 851       AliVParticle *track = GetAcceptParticleFromArray(i,0);
 852       if (!track)
 853         continue;
 854       if (track->Pt() > fTrackPtCut) {
 855         nTracksAcc++;
 856         if (nTracksAcc>=fMinNTrack)
 857           break;
 858       }
 859     }
 860     if (nTracksAcc<fMinNTrack) {
 861       if (fGeneralHistograms)
 862         fHistEventRejection->Fill("minNTrack",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       } else {
 947         Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
 948         fCentBin = TMath::FloorNint(fCent/centWidth);
 949         if (fCentBin>=fNcentBins) {
 950           AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
 951           fCentBin = fNcentBins-1;
 952         }
 953       }
 954     } else {
 955       AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
 956       fCentBin = 3;
 957     }
 958     AliEventplane *aliEP = InputEvent()->GetEventplane();
 959     if (aliEP) {
 960       fEPV0  = aliEP->GetEventplane("V0" ,InputEvent());
 961       fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
 962       fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
 963     } else {
 964       AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
 965     }
 966   } else {
 967     fCent = 99;
 968     fCentBin = 0;
 969   }
 970
 971   if (fIsPythia) {
 972
 973     if (MCEvent()) {
 974       fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 975       if (!fPythiaHeader) {
 976         // Check if AOD
 977         AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 978
 979         if (aodMCH) {
 980           for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
 981             fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
 982             if (fPythiaHeader) break;
 983           }
 984         }
 985       }
 986     }
 987
 988     if (fPythiaHeader) {
 989       fPtHard = fPythiaHeader->GetPtHard();
 990
 991       const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 992       const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 993       for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
 994         if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
 995           break;
 996       }
 997
 998       fNTrials = fPythiaHeader->Trials();
 999     }
1000   }
1001
1002   fTriggerType = GetTriggerType();
1003
1004   return kTRUE;
1005 }
1006
1007 //________________________________________________________________________
1008 AliParticleContainer* AliAnalysisTaskEmcal::AddParticleContainer(const char *n)
1009 {
1010   // Add particle container
1011   // will be called in AddTask macro
1012
1013   TString tmp = TString(n);
1014   if (tmp.IsNull()) return 0;
1015
1016   AliParticleContainer *cont = 0x0;
1017   cont = new AliParticleContainer();
1018   cont->SetArrayName(n);
1019   TString contName = cont->GetArrayName();
1020
1021   fParticleCollArray.Add(cont);
1022
1023   return cont;
1024 }
1025
1026 //________________________________________________________________________
1027 AliClusterContainer* AliAnalysisTaskEmcal::AddClusterContainer(const char *n)
1028 {
1029   // Add cluster container
1030   // will be called in AddTask macro
1031
1032   TString tmp = TString(n);
1033   if (tmp.IsNull()) return 0;
1034
1035   AliClusterContainer *cont = 0x0;
1036   cont = new AliClusterContainer();
1037   cont->SetArrayName(n);
1038
1039   fClusterCollArray.Add(cont);
1040
1041   return cont;
1042 }
1043
1044 //________________________________________________________________________
1045 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(Int_t i) const
1046 {
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 {
1057   // Get i^th cluster container
1058
1059   if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1060   AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1061   return cont;
1062 }
1063
1064 //________________________________________________________________________
1065 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(const char *name) const
1066 {
1067   // Get particle container with name
1068
1069   AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1070   return cont;
1071 }
1072
1073 //________________________________________________________________________
1074 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(const char *name) const
1075 {
1076   // Get cluster container with name
1077
1078   AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1079   return cont;
1080 }
1081
1082 //________________________________________________________________________
1083 TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const
1084 {
1085   // Get i^th TClonesArray with AliVParticle
1086
1087   AliParticleContainer *cont = GetParticleContainer(i);
1088   if (!cont) {
1089     AliError(Form("%s: Particle container %d not found",GetName(),i));
1090     return 0;
1091   }
1092   TString contName = cont->GetArrayName();
1093   return cont->GetArray();
1094 }
1095
1096 //________________________________________________________________________
1097 TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const
1098 {
1099   // Get i^th TClonesArray with AliVCluster
1100
1101   AliClusterContainer *cont = GetClusterContainer(i);
1102   if (!cont) {
1103     AliError(Form("%s:Cluster container %d not found",GetName(),i));
1104     return 0;
1105   }
1106   return cont->GetArray();
1107 }
1108
1109 //________________________________________________________________________
1110 AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const
1111 {
1112   // Get particle p if accepted from  container c
1113   // If particle not accepted return 0
1114
1115   AliParticleContainer *cont = GetParticleContainer(c);
1116   if (!cont) {
1117     AliError(Form("%s: Particle container %d not found",GetName(),c));
1118     return 0;
1119   }
1120   AliVParticle *vp = cont->GetAcceptParticle(p);
1121
1122   return vp;
1123 }
1124
1125 //________________________________________________________________________
1126 AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const
1127 {
1128   // Get particle p if accepted from  container c
1129   // If particle not accepted return 0
1130
1131   AliClusterContainer *cont = GetClusterContainer(c);
1132   if (!cont) {
1133     AliError(Form("%s: Cluster container %d not found",GetName(),c));
1134     return 0;
1135   }
1136   AliVCluster *vc = cont->GetAcceptCluster(cl);
1137
1138   return vc;
1139 }
1140
1141 //________________________________________________________________________
1142 Int_t AliAnalysisTaskEmcal::GetNParticles(Int_t i) const
1143 {
1144   // Get number of entries in particle array i
1145
1146   AliParticleContainer *cont = GetParticleContainer(i);
1147   if (!cont) {
1148     AliError(Form("%s: Particle container %d not found",GetName(),i));
1149     return 0;
1150   }
1151   return cont->GetNEntries();
1152 }
1153
1154 //________________________________________________________________________
1155 Int_t AliAnalysisTaskEmcal::GetNClusters(Int_t i) const
1156 {
1157   // Get number of entries in cluster array i
1158
1159   AliClusterContainer *cont = GetClusterContainer(i);
1160   if (!cont) {
1161     AliError(Form("%s: Cluster container %d not found",GetName(),i));
1162     return 0;
1163   }
1164   return cont->GetNEntries();
1165 }
1166
1167 //________________________________________________________________________
1168 AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch()
1169 {
1170   //get main trigger match; if not known yet, look for it and cache
1171
1172   if (fMainTriggerPatch)
1173     return fMainTriggerPatch;
1174
1175   if (!fTriggerPatchInfo) {
1176     AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1177     return 0;
1178   }
1179
1180   //number of patches in event
1181   Int_t nPatch = fTriggerPatchInfo->GetEntries();
1182
1183   //extract main trigger patch
1184   AliEmcalTriggerPatchInfo *patch;
1185   for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1186
1187     patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1188     if (patch->IsMainTrigger()) {
1189       fMainTriggerPatch = patch;
1190       break;
1191     }
1192   }
1193
1194   return fMainTriggerPatch;
1195 }
1196
1197 //________________________________________________________________________
1198 void AliAnalysisTaskEmcal::AddObjectToEvent(TObject *obj)
1199 {
1200   // Add object to event
1201
1202   if (!(InputEvent()->FindListObject(obj->GetName()))) {
1203     InputEvent()->AddObject(obj);
1204   } else {
1205     AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1206   }
1207 }