PWG/EMCAL/AliAnalysisTaskEmcal.cxx

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