PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx

   1 #ifndef ALIANALYSISTASKSE_H
   2 #include <Riostream.h>
   3 #include <TROOT.h>
   4 #include <TFile.h>
   5 #include <TCint.h>
   6 #include <TChain.h>
   7 #include <TTree.h>
   8 #include <TKey.h>
   9 #include <TProfile.h>
  10 #include <TH1F.h>
  11 #include <TH2F.h>
  12 #include <TCanvas.h>
  13 #include <TList.h>
  14 #include <TClonesArray.h>
  15 #include <TObject.h>
  16 #include <TMath.h>
  17 #include <TSystem.h>
  18 #include <TInterpreter.h>
  19 #include <TH1.h>
  20 #include "AliAnalysisTask.h"
  21 #include "AliCentrality.h"
  22 #include "AliStack.h"
  23 #include "AliESDEvent.h"
  24 #include "AliESDInputHandler.h"
  25 #include "AliAODEvent.h"
  26 #include "AliAODHandler.h"
  27 #include "AliAnalysisManager.h"
  28 #include "AliAnalysisTaskSE.h"
  29 #endif
  30
  31 #include <TRandom3.h>
  32 #include "AliGenPythiaEventHeader.h"
  33 #include "AliMCEvent.h"
  34 #include "AliLog.h"
  35 #include <AliEmcalJet.h>
  36 #include "AliVEventHandler.h"
  37 #include "AliVParticle.h"
  38 #include "AliAnalysisUtils.h"
  39
  40
  41 #include "AliAnalysisTaskChargedJetsPA.h"
  42
  43 //TODO: Not accessing the particles when using MC
  44 //TODO: FillHistogram can be done better with virtual TH1(?)
  45 ClassImp(AliAnalysisTaskChargedJetsPA)
  46
  47 // ######################################################################################## DEFINE HISTOGRAMS
  48 void AliAnalysisTaskChargedJetsPA::Init()
  49 {
  50   #ifdef DEBUGMODE
  51     AliInfo("Creating histograms.");
  52   #endif
  53
  54   AddHistogram1D<TH1D>("hNumberEvents", "Number of events (0 = before, 1 = after vertex cuts)", "", 2, 0, 2, "#Delta z(cm)","N^{Events}/cut");
  55
  56   // NOTE: Jet histograms
  57   if (fAnalyzeJets)
  58   {
  59     // ######## Jet spectra
  60     AddHistogram2D<TH2D>("hJetPt", "Jets p_{T} distribution", "", 1000, 0., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  61     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedRC", "Jets p_{T} distribution, RC background subtracted", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  62     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKT", "Jets p_{T} distribution, KT background subtracted", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  63     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedTR", "Jets p_{T} distribution, TR background subtracted", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  64     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedRCNoEtaCorr", "Jets p_{T} distribution, RC background subtracted (no #eta correction)", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  65     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTNoEtaCorr", "Jets p_{T} distribution, KT background subtracted (no #eta correction)", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  66     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedTRNoEtaCorr", "Jets p_{T} distribution, TR background subtracted (no #eta correction)", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  67
  68     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedRCEtaBins", "Jets p_{T} distribution, RC background (in #eta bin) subtracted (no #eta correction)", "", 500, -50., 200., 5, -0.5, 0.5, "p_{T} (GeV/c)","#eta","dN^{Jets}/dp_{T}d#eta");
  69     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTEtaBins", "Jets p_{T} distribution, KT background (in #eta bin) subtracted (no #eta correction)", "", 500, -50., 200., 5, -0.5, 0.5, "p_{T} (GeV/c)","#eta","dN^{Jets}/dp_{T}d#eta");
  70     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedTREtaBins", "Jets p_{T} distribution, TR background (in #eta bin) subtracted (no #eta correction)", "", 500, -50., 200., 5, -0.5, 0.5, "p_{T} (GeV/c)","#eta","dN^{Jets}/dp_{T}d#eta");
  71
  72     // ######## Jet stuff
  73     AddHistogram1D<TH1D>("hJetCountAll", "Number of Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
  74     AddHistogram1D<TH1D>("hJetCountAccepted", "Number of accepted Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
  75     AddHistogram1D<TH1D>("hLeadingJetPt", "Leading jet p_{T}", "", 500,  0, 100, "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
  76     AddHistogram1D<TH1D>("hSecondLeadingJetPt", "Second Leading jet p_{T}", "", 500,  0, 100, "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
  77     AddHistogram1D<TH1D>("hJetDeltaPhi", "Jets combinatorial #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
  78     AddHistogram1D<TH1D>("hLeadingJetDeltaPhi", "1st and 2nd leading jet #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
  79   }
  80
  81   // NOTE: Jet background histograms
  82   if (fAnalyzeBackground)
  83   {
  84     // ########## Different background estimates
  85     AddHistogram2D<TH2D>("hRCBackground", "RC background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
  86     AddHistogram2D<TH2D>("hKTBackground", "KT background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
  87     AddHistogram2D<TH2D>("hTRBackground", "TR background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
  88
  89     // ########## Delta Pt
  90     AddHistogram2D<TH2D>("hDeltaPtKT", "Background fluctuations #delta p_{T} (KT)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  91     AddHistogram2D<TH2D>("hDeltaPtRC", "Background fluctuations #delta p_{T} (RC)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  92     AddHistogram2D<TH2D>("hDeltaPtTR", "Background fluctuations #delta p_{T} (TR)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  93     AddHistogram2D<TH2D>("hDeltaPtKTNoEtaCorr", "Background fluctuations #delta p_{T} (KT, no #eta correction)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  94     AddHistogram2D<TH2D>("hDeltaPtRCNoEtaCorr", "Background fluctuations #delta p_{T} (RC, no #eta correction)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  95     AddHistogram2D<TH2D>("hDeltaPtTRNoEtaCorr", "Background fluctuations #delta p_{T} (TR, no #eta correction)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  96     AddHistogram2D<TH2D>("hDeltaPtKTNoEtaCorrNoExcl", "Background fluctuations #delta p_{T} (KT, no #eta correction, no leading jet correction)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  97     AddHistogram2D<TH2D>("hDeltaPtRCNoEtaCorrNoExcl", "Background fluctuations #delta p_{T} (RC, no #eta correction, no leading jet correction)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  98     AddHistogram2D<TH2D>("hDeltaPtTRNoEtaCorrNoExcl", "Background fluctuations #delta p_{T} (TR, no #eta correction, no leading jet correction)", "", 600, -40., 80., 5, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
  99
 100     // ########## Min bias background in eta bins
 101     AddHistogram2D<TH2D>("hRCBackgroundEtaBins", "RC background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta");
 102     AddHistogram2D<TH2D>("hKTBackgroundEtaBins", "KT background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta");
 103     AddHistogram2D<TH2D>("hTRBackgroundEtaBins", "TR background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta");
 104
 105     AddHistogram2D<TH2D>("hRCBackgroundEtaBinsCorrected", "RC background density (2 leading jets excluded, #eta-corrected)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta");
 106     AddHistogram2D<TH2D>("hKTBackgroundEtaBinsCorrected", "KT background density (2 leading jets excluded, #eta-corrected)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta");
 107     AddHistogram2D<TH2D>("hTRBackgroundEtaBinsCorrected", "TR background density (2 leading jets excluded, #eta-corrected)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta");
 108
 109     // ########## Dijet stuff
 110     AddHistogram1D<TH1D>("hDijetLeadingJetPt", "Dijet leading jet p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 111     AddHistogram1D<TH1D>("hDijetConstituentsPt", "Dijet constituents p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 112     AddHistogram2D<TH2D>("hDijetPtCorrelation", "Dijet constituents p_{T} correlation", "COLZ", 500, 5., 100., 500, 5., 100., "1st leading jet p_{T} (GeV/c)","2nd leading jet p_{T} (GeV/c)","dN^{Dijets}/d^{2}p_{T}");
 113     AddHistogram2D<TH2D>("hDijetBackground", "Background density (dijets excluded)", "", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 114     AddHistogram2D<TH2D>("hDijetBackgroundPerpendicular", "Background density (dijets excluded)", "", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 115   }
 116
 117   // NOTE: Pythia histograms
 118   if (fAnalyzePythia)
 119   {
 120     AddHistogram1D<TH1D>("hPythiaPtHard", "Pythia p_{T} hard distribution", "", 2000, 0, 400, "p_{T} hard","dN^{Events}/dp_{T,hard}");
 121     AddHistogram1D<TProfile>("hPythiaXSection", "Pythia cross section distribution", "", fNumPtHardBins+2, -1, fNumPtHardBins+1, "p_{T} hard bin","dN^{Events}/dp_{T,hard}");
 122     AddHistogram1D<TH1D>("hPythiaNTrials", "Pythia trials (no correction for manual cuts)", "", fNumPtHardBins+2, -1, fNumPtHardBins+1, "p_{T} hard bin", "Trials");
 123   }
 124
 125   // register Histograms
 126   for (Int_t i = 0; i < fHistCount; i++)
 127   {
 128     fOutputList->Add(fHistList->At(i));
 129   }
 130
 131   PostData(1,fOutputList); // important for merging
 132
 133 }
 134
 135 //________________________________________________________________________
 136 AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName) : AliAnalysisTaskSE(name), fOutputList(0), fAnalyzeJets(1), fAnalyzeBackground(1), fAnalyzePythia(0), fHasTracks(0), fHasJets(0), fHasBackgroundJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fNumPtHardBins(11), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fTRBackgroundConeRadius(0.4), fNumberRandCones(8), fNumberExcludedJets(2), fDijetMaxAngleDeviation(10.0), fJetKTEtaCorrection(0), fJetRCEtaCorrection(0), fJetTREtaCorrection(0), fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinJetPt(1.0), fMinJetArea(0.4), fMinBackgroundJetPt(0.15), fMinDijetLeadingPt(10.0), fCentralityType("V0A"), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0),  fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0)
 137 {
 138   #ifdef DEBUGMODE
 139     AliInfo("Calling constructor.");
 140   #endif
 141
 142   // Every instance of this task gets his own number
 143   static Int_t instance = 0;
 144   fTaskInstanceCounter = instance;
 145   instance++;
 146
 147   fTrackArrayName = new TString(trackArrayName);
 148   if (fTrackArrayName->Contains("MCParticles")) //TODO: Not working for now
 149     fIsMC = kTRUE;
 150
 151   fJetArrayName = new TString(jetArrayName);
 152   if (strcmp(fJetArrayName->Data(),"") == 0)
 153     fAnalyzeJets = kFALSE;
 154   else
 155     fAnalyzeJets = kTRUE;
 156
 157   fBackgroundJetArrayName = new TString(backgroundJetArrayName);
 158   if (strcmp(fBackgroundJetArrayName->Data(),"") == 0)
 159     fAnalyzeBackground = kFALSE;
 160   else
 161     fAnalyzeBackground = kTRUE;
 162
 163   DefineOutput(1, TList::Class());
 164
 165   fHistList = new TList();
 166
 167   #ifdef DEBUGMODE
 168     AliInfo("Constructor done.");
 169   #endif
 170
 171 }
 172
 173 //________________________________________________________________________
 174 inline Double_t AliAnalysisTaskChargedJetsPA::GetConePt(Double_t eta, Double_t phi, Double_t radius)
 175 {
 176   Double_t tmpConePt = 0.0;
 177
 178   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
 179   {
 180     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
 181     if (IsTrackInAcceptance(tmpTrack))
 182       if(IsTrackInCone(tmpTrack, eta, phi, radius))
 183         tmpConePt = tmpConePt + tmpTrack->Pt();
 184   }
 185   return tmpConePt;
 186 }
 187
 188
 189 //________________________________________________________________________
 190 inline Double_t AliAnalysisTaskChargedJetsPA::GetPtHard()
 191 {
 192   Double_t tmpPtHard = -1.0;
 193
 194   if (!MCEvent())
 195     AliError("MCEvent not accessible although demanded!");
 196   else
 197   {
 198     AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 199     if (!pythiaHeader)
 200       AliError("Pythia Header not accessible!");
 201     else
 202       tmpPtHard = pythiaHeader->GetPtHard();
 203   }
 204   return tmpPtHard;
 205 }
 206
 207
 208 //________________________________________________________________________
 209 inline Int_t AliAnalysisTaskChargedJetsPA::GetPtHardBin()
 210 {
 211   // ########## PT HARD BIN EDGES
 212   const Int_t kPtHardLowerEdges[] =  { 0, 5,11,21,36,57, 84,117,152,191,234};
 213   const Int_t kPtHardHigherEdges[] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 214
 215   Int_t tmpPtHardBin = 0;
 216   Double_t tmpPtHard = GetPtHard();
 217
 218   for (tmpPtHardBin = 0; tmpPtHardBin <= fNumPtHardBins; tmpPtHardBin++)
 219     if (tmpPtHard >= kPtHardLowerEdges[tmpPtHardBin] && tmpPtHard < kPtHardHigherEdges[tmpPtHardBin])
 220       break;
 221
 222   return tmpPtHardBin;
 223 }
 224
 225
 226 //________________________________________________________________________
 227 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius)
 228 {
 229   // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
 230   Double_t trackPhi = 0.0;
 231   if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
 232     trackPhi = track->Phi() - TMath::TwoPi();
 233   else if (track->Phi() < (phi+radius - TMath::TwoPi()))
 234     trackPhi = track->Phi() + TMath::TwoPi();
 235   else
 236     trackPhi = track->Phi();
 237
 238   if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
 239     return kTRUE;
 240
 241   return kFALSE;
 242 }
 243
 244 //________________________________________________________________________
 245 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInAcceptance(AliVParticle* track)
 246 {
 247   if (track != 0)
 248     if (TMath::Abs(track->Eta()) <= fTrackEtaWindow)
 249       if (track->Pt() >= fMinTrackPt)
 250         return kTRUE;
 251
 252   return kFALSE;
 253 }
 254
 255 //________________________________________________________________________
 256 inline Bool_t AliAnalysisTaskChargedJetsPA::IsBackgroundJetInAcceptance(AliEmcalJet *jet)
 257 {
 258   if (jet != 0)
 259     if (TMath::Abs(jet->Eta()) <= fBackgroundJetEtaWindow)
 260       if (jet->Pt() >= fMinBackgroundJetPt)
 261         return kTRUE;
 262
 263   return kFALSE;
 264 }
 265
 266 //________________________________________________________________________
 267 inline Bool_t AliAnalysisTaskChargedJetsPA::IsSignalJetInAcceptance(AliEmcalJet *jet)
 268 {
 269   if (jet != 0)
 270     if (TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow)
 271       if (jet->Pt() >= fMinJetPt)
 272         if (jet->Area() >= fMinJetArea)
 273           return kTRUE;
 274
 275   return kFALSE;
 276 }
 277
 278 //________________________________________________________________________
 279 inline Bool_t AliAnalysisTaskChargedJetsPA::IsDijet(AliEmcalJet *jet1, AliEmcalJet *jet2)
 280 {
 281   // Output from GetDeltaPhi is < pi in any case
 282   if ((jet1 != 0) && (jet2 != 0))
 283     if((TMath::Pi() - GetDeltaPhi(jet1->Phi(),jet2->Phi())) < fDijetMaxAngleDeviation)
 284       if ((jet1->Pt() > fMinDijetLeadingPt) && (jet2->Pt() > fMinDijetLeadingPt)) //TODO: Introduce recoil jet?
 285         return kTRUE;
 286
 287   return kFALSE;
 288 }
 289
 290 //________________________________________________________________________
 291 void AliAnalysisTaskChargedJetsPA::ExecOnce()
 292 {
 293   #ifdef DEBUGMODE
 294     AliInfo("Starting ExecOnce.");
 295   #endif
 296   fInitialized = kTRUE;
 297
 298   // Check for track array
 299   if (strcmp(fTrackArrayName->Data(), "") != 0)
 300   {
 301     fTrackArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrayName->Data()));
 302     fHasTracks = kTRUE;
 303     if (!fTrackArray)
 304     {
 305       AliWarning(Form("%s: Could not retrieve tracks %s! This is OK, if tracks are not demanded.", GetName(), fTrackArrayName->Data()));
 306       fHasTracks = kFALSE;
 307     }
 308     else
 309     {
 310       TClass *cl = fTrackArray->GetClass();
 311       if (!cl->GetBaseClass("AliVParticle"))
 312       {
 313         AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName->Data()));
 314         fTrackArray = 0;
 315         fHasTracks = kFALSE;
 316       }
 317     }
 318   }
 319
 320   // Check for jet array
 321   if (strcmp(fJetArrayName->Data(), "") != 0)
 322   {
 323     fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrayName->Data()));
 324     fHasJets = kTRUE;
 325
 326     if (!fJetArray)
 327     {
 328       AliWarning(Form("%s: Could not retrieve jets %s! This is OK, if jets are not demanded.", GetName(), fJetArrayName->Data()));
 329       fHasJets = kFALSE;
 330     }
 331     else
 332     {
 333       if (!fJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
 334       {
 335         AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName->Data()));
 336         fJetArray = 0;
 337         fHasJets = kFALSE;
 338       }
 339     }
 340   }
 341
 342   // Check for background object
 343   if (strcmp(fBackgroundJetArrayName->Data(), "") != 0)
 344   {
 345     fBackgroundJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fBackgroundJetArrayName->Data()));
 346     fHasBackgroundJets = kTRUE;
 347     if (!fBackgroundJetArray)
 348     {
 349       AliInfo(Form("%s: Could not retrieve background jets %s! This is OK, if background is not demanded.", GetName(), fBackgroundJetArrayName->Data()));
 350       fHasBackgroundJets = kFALSE;
 351     }
 352   }
 353
 354   // Look, if initialization is OK
 355   if (!fHasTracks && fAnalyzeBackground)
 356   {
 357     AliError(Form("%s: Tracks NOT successfully casted although demanded! Deactivating background analysis",GetName()));
 358     fAnalyzeBackground = kFALSE;
 359   }
 360   if ((!fHasJets && fAnalyzeJets) || (!fHasJets && fAnalyzeBackground))
 361   {
 362     AliError(Form("%s: Jets NOT successfully casted although demanded!  Deactivating jet- and background analysis",GetName()));
 363     fAnalyzeJets = kFALSE;
 364     fAnalyzeBackground = kFALSE;
 365   }
 366   if (!fHasBackgroundJets && fAnalyzeBackground)
 367   {
 368     AliError(Form("%s: Background NOT successfully casted although demanded!  Deactivating background analysis",GetName()));
 369     fAnalyzeBackground = kFALSE;
 370   }
 371
 372   // Initialize helper class (for vertex selection)
 373   fHelperClass = new AliAnalysisUtils();
 374
 375   // Histogram init
 376   Init();
 377
 378   #ifdef DEBUGMODE
 379     AliInfo("ExecOnce done.");
 380   #endif
 381
 382 }
 383
 384 //________________________________________________________________________
 385 void AliAnalysisTaskChargedJetsPA::GetSignalJets()
 386 {
 387   // Reset vars
 388   fFirstLeadingJet = NULL;
 389   fSecondLeadingJet = NULL;
 390   fNumberSignalJets = 0;
 391
 392   Float_t maxJetPts[] = {0, 0};
 393   Int_t jetIDArray[]   = {-1, -1};
 394   Int_t jetCount = fJetArray->GetEntries();
 395
 396   // Go through all jets and save signal jets and the two leading ones
 397   for (Int_t i = 0; i < jetCount; i++)
 398   {
 399     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(i));
 400     if (!jet)
 401     {
 402       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 403       continue;
 404     }
 405
 406     if (!IsSignalJetInAcceptance(jet)) continue;
 407
 408     if (jet->Pt() > maxJetPts[0])
 409     {
 410       maxJetPts[1] = maxJetPts[0];
 411       jetIDArray[1] = jetIDArray[0];
 412       maxJetPts[0] = jet->Pt();
 413       jetIDArray[0] = i;
 414     }
 415     else if (jet->Pt() > maxJetPts[1])
 416     {
 417       maxJetPts[1] = jet->Pt();
 418       jetIDArray[1] = i;
 419     }
 420     fSignalJets[fNumberSignalJets] = jet;
 421     fNumberSignalJets++;
 422   }
 423
 424   if (fNumberSignalJets > 0)
 425     fFirstLeadingJet  = static_cast<AliEmcalJet*>(fJetArray->At(jetIDArray[0]));
 426   if (fNumberSignalJets > 1)
 427     fSecondLeadingJet = static_cast<AliEmcalJet*>(fJetArray->At(jetIDArray[1]));
 428
 429 }
 430
 431 //________________________________________________________________________
 432 Int_t AliAnalysisTaskChargedJetsPA::GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray, Bool_t isSignalJets)
 433 {
 434 // Writes first two leading jets into already registered array jetIDArray
 435
 436   if (!jetArray)
 437   {
 438     AliError("Could not get the jet array to get leading jets from it!");
 439     return 0;
 440   }
 441
 442   Float_t maxJetPts[] = {0, 0};
 443   jetIDArray[0] = -1;
 444   jetIDArray[1] = -1;
 445
 446   Int_t jetCount = jetArray->GetEntries();
 447   Int_t jetCountAccepted = 0;
 448
 449   for (Int_t i = 0; i < jetCount; i++)
 450   {
 451     AliEmcalJet* jet = static_cast<AliEmcalJet*>(jetArray->At(i));
 452     if (!jet)
 453     {
 454       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 455       continue;
 456     }
 457
 458     if(isSignalJets)
 459     {
 460       if (!IsSignalJetInAcceptance(jet)) continue;
 461     }
 462     else
 463     {
 464       if (!IsBackgroundJetInAcceptance(jet)) continue;
 465     }
 466
 467     if (jet->Pt() > maxJetPts[0])
 468     {
 469       maxJetPts[1] = maxJetPts[0];
 470       jetIDArray[1] = jetIDArray[0];
 471       maxJetPts[0] = jet->Pt();
 472       jetIDArray[0] = i;
 473     }
 474     else if (jet->Pt() > maxJetPts[1])
 475     {
 476       maxJetPts[1] = jet->Pt();
 477       jetIDArray[1] = i;
 478     }
 479     jetCountAccepted++;
 480   }
 481   return jetCountAccepted;
 482 }
 483
 484 //________________________________________________________________________
 485 Double_t AliAnalysisTaskChargedJetsPA::GetBackgroundEtaCorrFactor(EtaCorrectionMode mode, Double_t eta)
 486 {
 487   if ((eta>=-0.5) && (eta<-0.3))
 488     return GetBackgroundEtaBinCorrFactor(mode, 1);
 489   else if ((eta>=-0.3) && (eta<-0.1))
 490     return GetBackgroundEtaBinCorrFactor(mode, 2);
 491   else if ((eta>=-0.1) && (eta<+0.1))
 492     return GetBackgroundEtaBinCorrFactor(mode, 3);
 493   else if ((eta>=+0.1) && (eta<+0.3))
 494     return GetBackgroundEtaBinCorrFactor(mode, 4);
 495   else if ((eta>=+0.3) && (eta<=+0.5))
 496     return GetBackgroundEtaBinCorrFactor(mode, 5);
 497   else
 498     AliError(Form("Wrong eta value! Eta=%1.4f", eta));
 499
 500   return 1.0;
 501 }
 502
 503 //________________________________________________________________________
 504 Double_t AliAnalysisTaskChargedJetsPA::GetBackgroundEtaBinCorrFactor(EtaCorrectionMode mode, Int_t eta)
 505 {
 506   Double_t corrFactor = 1.0;
 507
 508   if((eta < 1) || (eta>5))
 509   {
 510     AliError("Wrong eta bin!");
 511     return corrFactor;
 512   }
 513
 514   if((mode == kKTEtaCorrection) && fJetKTEtaCorrection)
 515     corrFactor = fJetKTEtaCorrection->GetBinContent(eta);
 516   else if((mode == kRCEtaCorrection) && fJetRCEtaCorrection)
 517     corrFactor = fJetRCEtaCorrection->GetBinContent(eta);
 518   else if((mode == kTREtaCorrection) && fJetTREtaCorrection)
 519     corrFactor = fJetTREtaCorrection->GetBinContent(eta);
 520   else
 521     corrFactor = 1.0;
 522
 523   return corrFactor;
 524 }
 525
 526 //________________________________________________________________________
 527 Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Double_t background, EtaCorrectionMode mode)
 528 {
 529   #ifdef DEBUGMODE
 530     AliInfo("Getting corrected jet spectra.");
 531   #endif
 532
 533   if(!jet)
 534   {
 535     AliError("Jet pointer passed to GetCorrectedJet() not valid!");
 536     return -1.0;
 537   }
 538
 539   Double_t correctedPt = -1.0;
 540
 541   // if the passed background is not valid, do not subtract it
 542   if(background < 0)
 543     background = 0;
 544
 545   // Get Eta corrected background
 546   // The correction factor is to take the eta dependence of the background into account and to correct for it.
 547   Double_t tmpCorrectedBackground = background * GetBackgroundEtaCorrFactor(mode, jet->Eta());
 548
 549   // Subtract background
 550   correctedPt = jet->Pt() - tmpCorrectedBackground * jet->Area();
 551
 552   #ifdef DEBUGMODE
 553     AliInfo("Got corrected jet spectra.");
 554   #endif
 555
 556   return correctedPt;
 557 }
 558
 559
 560
 561 //________________________________________________________________________
 562 void AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t& deltaPt, Double_t rho, EtaCorrectionMode mode, Bool_t leadingJetExclusion)
 563 {
 564   #ifdef DEBUGMODE
 565     AliInfo("Getting Delta Pt.");
 566   #endif
 567
 568   // Define an invalid delta pt
 569   deltaPt = -10000.0;
 570
 571   // Define eta range
 572   Double_t etaMin, etaMax;
 573   etaMin = -(fTrackEtaWindow-fRandConeRadius);
 574   etaMax = +(fTrackEtaWindow-fRandConeRadius);
 575
 576   // Define random cone
 577   Bool_t coneValid = kTRUE;
 578   Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
 579   Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
 580
 581   // Apply eta correction on background if demanded
 582   rho *= GetBackgroundEtaCorrFactor(mode, tmpRandConeEta);
 583
 584   AliEmcalJet* tmpJet = fFirstLeadingJet;
 585   // if there is a jet, check for overlap if demanded
 586   if(tmpJet && leadingJetExclusion)
 587   {
 588     Double_t excludedJetPhi = tmpJet->Phi();
 589     Double_t excludedJetEta = tmpJet->Eta();
 590     Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, excludedJetPhi);
 591
 592     // Check, if cone has overlap with jet
 593     if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(tmpRandConeEta-excludedJetEta)*TMath::Abs(tmpRandConeEta-excludedJetEta) <= fRandConeRadius*fRandConeRadius)
 594     {
 595       // Define probability to exclude the RC
 596       Double_t probability = 1 - (fNumberSignalJets-1)/fNumberSignalJets;
 597
 598       // Only exclude cone with a given probability
 599       if (fRandom->Rndm()<=probability)
 600         coneValid = kFALSE;
 601     }
 602   }
 603
 604   // Get the cones' pt and calculate delta pt
 605   if (coneValid)
 606     deltaPt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius) - (rho*fRandConeRadius*fRandConeRadius*TMath::Pi());
 607
 608   #ifdef DEBUGMODE
 609     AliInfo("Got Delta Pt.");
 610   #endif
 611 }
 612
 613 //________________________________________________________________________
 614 void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMedian, Double_t& areaMean, Double_t etaMin, Double_t etaMax)
 615 {
 616   #ifdef DEBUGMODE
 617     AliInfo("Getting KT background density.");
 618   #endif
 619
 620   // static declaration. Advantage: more speed. Disadvantage: Problematic for events with more than 1024 jets :)
 621   static Double_t tmpRhos[1024];
 622   static Double_t tmpAreas[1024];
 623   Int_t maxJetIds[]   = {-1, -1}; // Indices for excludes jets (up to two)
 624
 625   // Setting invalid values
 626   rhoMedian = -1.0;
 627   areaMean= -1.0;
 628
 629   // Exclude UP TO numberExcludeLeadingJets
 630   Int_t numberBgrdJets = GetLeadingJets(fBackgroundJetArray, &maxJetIds[0], kFALSE);
 631   if (numberBgrdJets < numberExcludeLeadingJets)
 632     numberExcludeLeadingJets = numberBgrdJets;
 633   if ((etaMin == 0) && (etaMax == 0))
 634   {
 635     etaMin = -fBackgroundJetEtaWindow;
 636     etaMax = +fBackgroundJetEtaWindow;
 637   }
 638
 639   Int_t jetCountAccepted = 0;
 640   Int_t jetCount = fBackgroundJetArray->GetEntries();
 641
 642   for (Int_t i = 0; i < jetCount; i++)
 643   {
 644     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
 645     if (!jet)
 646     {
 647       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 648       continue;
 649     }
 650
 651     // exclude leading jets
 652     if (numberExcludeLeadingJets > 0)
 653       if (i == maxJetIds[0])
 654         continue;
 655     if (numberExcludeLeadingJets > 1)
 656       if (i == maxJetIds[1])
 657         continue;
 658
 659
 660
 661     if (!IsBackgroundJetInAcceptance(jet))
 662       continue;
 663     if (!((jet->Eta() >= etaMin) && (jet->Eta() < etaMax)))
 664       continue;
 665
 666
 667     tmpRhos[jetCountAccepted] = jet->Pt() / jet->Area();
 668     tmpAreas[jetCountAccepted] = jet->Area();
 669     jetCountAccepted++;
 670   }
 671
 672   if (jetCountAccepted > 0)
 673   {
 674     rhoMedian = TMath::Median(jetCountAccepted, tmpRhos);
 675     areaMean   = TMath::Mean(jetCountAccepted, tmpAreas);
 676   }
 677   #ifdef DEBUGMODE
 678     AliInfo("Got KT background density.");
 679   #endif
 680 }
 681
 682
 683 //________________________________________________________________________
 684 Int_t AliAnalysisTaskChargedJetsPA::GetRCBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& rhoMedian, Double_t etaMin, Double_t etaMax, Int_t numberRandCones)
 685 {
 686   #ifdef DEBUGMODE
 687     AliInfo("Getting RC background density.");
 688   #endif
 689
 690   if(numberRandCones == 0)
 691     numberRandCones = fNumberRandCones;
 692
 693   std::vector<AliEmcalJet> tmpCones(numberRandCones);
 694
 695   // Setting invalid values
 696   rhoMean = -1.0;
 697   rhoMedian = -1.0;
 698
 699   // Exclude UP TO numberExcludeLeadingJets
 700   if (fNumberSignalJets < 2)
 701     numberExcludeLeadingJets = fNumberSignalJets;
 702
 703   // Search given amount of RCs
 704   Int_t numAcceptedRCs = 0;
 705   for(Int_t i=0;i<numberRandCones;i++)
 706   {
 707     Double_t tmpRandConeEta = 0.0;
 708     Double_t tmpRandConePhi = 0.0;
 709     Double_t excludedJetEta = 0.0;
 710     Double_t excludedJetPhi = 0.0;
 711
 712     // Search random cone in acceptance with no overlap with already excluded jets (leading jets and random cones)
 713     Bool_t coneValid = kTRUE;
 714
 715     // Set the random cone position
 716     if ((etaMin == 0) && (etaMax == 0))
 717       tmpRandConeEta = (fTrackEtaWindow-fRandConeRadius)*(2.0*fRandom->Rndm()-1.0); // full RC is in acceptance
 718     else
 719       tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
 720
 721     tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
 722
 723     // Go through all excluded leading jets and check if there's an overlap
 724
 725     for(Int_t j=0;j<numberExcludeLeadingJets;j++)
 726     {
 727       AliEmcalJet* tmpJet = NULL;
 728
 729       if (j==0)
 730         tmpJet = fFirstLeadingJet;
 731       else if (j==1)
 732         tmpJet = fSecondLeadingJet;
 733       else
 734         AliFatal("Trying to exclude more than 2 jets in RC background -- not implemented.");
 735
 736       excludedJetPhi = tmpJet->Phi();
 737       excludedJetEta = tmpJet->Eta();
 738       Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, excludedJetPhi);
 739
 740       if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(tmpRandConeEta-excludedJetEta)*TMath::Abs(tmpRandConeEta-excludedJetEta) <= fRandConeRadius*fRandConeRadius)
 741       {
 742         coneValid = kFALSE;
 743         break;
 744       }
 745     }
 746
 747     // RC is accepted, so save it
 748     if(coneValid)
 749     {
 750       AliEmcalJet tmpJet(GetConePt(tmpRandConeEta, tmpRandConePhi, fRandConeRadius), tmpRandConeEta, tmpRandConePhi, 0.0);
 751       tmpCones[numAcceptedRCs] = tmpJet;
 752       numAcceptedRCs++;
 753     }
 754   }
 755
 756   // Calculate Rho and the mean from the RCs (no excluded jets are considered!)
 757   if(numAcceptedRCs > 0)
 758   {
 759     std::vector<Double_t> tmpRho(numAcceptedRCs);
 760     for (Int_t i=0; i<numAcceptedRCs;i++)
 761       tmpRho[i] = tmpCones[i].Pt()/(fRandConeRadius*fRandConeRadius*TMath::Pi());
 762
 763     rhoMean = TMath::Mean(tmpRho.begin(), tmpRho.end());
 764     rhoMedian = 0.0; // NOT IMPLEMENTED because TMath::Median is not working with iterators
 765   }
 766
 767   #ifdef DEBUGMODE
 768     AliInfo("Got RC background density.");
 769   #endif
 770   return numAcceptedRCs;
 771 }
 772
 773 //________________________________________________________________________
 774 void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, Double_t etaMin, Double_t etaMax)
 775 {
 776   #ifdef DEBUGMODE
 777     AliInfo("Getting TR background density.");
 778   #endif
 779
 780   Double_t summedTracksPt = 0.0;
 781
 782   if ((etaMin == 0) && (etaMax == 0))
 783   {
 784     etaMin = -fTrackEtaWindow;
 785     etaMax = +fTrackEtaWindow;
 786   }
 787
 788   // Setting invalid values
 789   rhoMean = -1.0;
 790   area = -1.0;
 791   // Exclude UP TO numberExcludeLeadingJets
 792   if (fNumberSignalJets < 2)
 793     numberExcludeLeadingJets = fNumberSignalJets;
 794
 795
 796   Int_t   trackCount = fTrackArray->GetEntries();
 797   Int_t   trackCountAccepted = 0;
 798   for (Int_t i = 0; i < trackCount; i++)
 799   {
 800     Bool_t  trackValid = kTRUE;
 801     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
 802     if (IsTrackInAcceptance(tmpTrack))
 803       if ((tmpTrack->Eta() >= etaMin) && (tmpTrack->Eta() < etaMax))
 804       {
 805         for (Int_t j = 0; j < numberExcludeLeadingJets; j++)
 806         {
 807           AliEmcalJet* tmpJet = NULL;
 808           if (j==0)
 809             tmpJet = fFirstLeadingJet;
 810           else if (j==1)
 811             tmpJet = fSecondLeadingJet;
 812           else
 813             AliFatal("Trying to exclude more than 2 jets in track background -- not implemented.");
 814
 815           if (IsTrackInCone(tmpTrack, tmpJet->Eta(), tmpJet->Phi(), fTRBackgroundConeRadius))
 816           {
 817             trackValid = kFALSE;
 818             break;
 819           }
 820         }
 821         if (trackValid)
 822         {
 823           // Add track pt to array
 824           summedTracksPt = summedTracksPt + tmpTrack->Pt();
 825           trackCountAccepted++;
 826         }
 827       }
 828   }
 829
 830   if (trackCountAccepted > 0)
 831   {
 832     Double_t tmpArea = 0.0;
 833
 834     tmpArea = (2.0*fTrackEtaWindow) * TMath::TwoPi() * (etaMax-etaMin)/(2.0*fTrackEtaWindow); // area of the used eta strip
 835
 836     // Now: exclude the part of the leading jet that is in the strip
 837     if (numberExcludeLeadingJets == 2)
 838       tmpArea = tmpArea*(1.0-MCGetOverlapCircleRectancle(fFirstLeadingJet->Eta(), fFirstLeadingJet->Phi(), fTRBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi()) -MCGetOverlapCircleRectancle(fSecondLeadingJet->Eta(), fSecondLeadingJet->Phi(), fTRBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi()));
 839     else if (numberExcludeLeadingJets == 1)
 840       tmpArea = tmpArea*(1.0-MCGetOverlapCircleRectancle(fFirstLeadingJet->Eta(), fFirstLeadingJet->Phi(), fTRBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi()));
 841
 842     rhoMean = summedTracksPt/tmpArea;
 843     area  = tmpArea;
 844   }
 845
 846   #ifdef DEBUGMODE
 847     AliInfo("Got TR background density.");
 848   #endif
 849 }
 850
 851 //________________________________________________________________________
 852 void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular)
 853 {
 854   #ifdef DEBUGMODE
 855     AliInfo("Getting TR background density.");
 856   #endif
 857
 858   // Setting invalid values
 859   Double_t summedTracksPt = 0.0;
 860   rhoMean = -1.0;
 861   area = -1.0;
 862
 863   Double_t tmpRadius = 0.0;
 864   if (doSearchPerpendicular)
 865     tmpRadius = 0.5*TMath::Pi(); // exclude 90 degrees around jets
 866   else
 867     tmpRadius = fSignalJetRadius;
 868
 869   numberExcludeLeadingJets = 2; // dijet is excluded here in any case
 870
 871
 872
 873   if (!fTrackArray || !fJetArray)
 874   {
 875     AliError("Could not get the track/jet array to calculate track rho!");
 876     return;
 877   }
 878
 879   Int_t   trackCount = fTrackArray->GetEntries();
 880   Int_t   trackCountAccepted = 0;
 881   for (Int_t i = 0; i < trackCount; i++)
 882   {
 883     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
 884     if (IsTrackInAcceptance(tmpTrack))
 885     {
 886       if (IsTrackInCone(tmpTrack, excludeJet1->Eta(), excludeJet1->Phi(), tmpRadius))
 887         continue;
 888
 889       if (numberExcludeLeadingJets > 1)
 890         if (IsTrackInCone(tmpTrack, excludeJet2->Eta(), excludeJet2->Phi(), tmpRadius))
 891           continue;
 892
 893         // Add track pt to array
 894         summedTracksPt = summedTracksPt + tmpTrack->Pt();
 895         trackCountAccepted++;
 896     }
 897   }
 898
 899   if (trackCountAccepted > 0)
 900   {
 901     Double_t tmpArea = 2.0*fTrackEtaWindow*TMath::TwoPi()  - 2*(tmpRadius*tmpRadius * TMath::Pi()); //TPC area - excluding jet area
 902     rhoMean = summedTracksPt/tmpArea;
 903     area = tmpArea;
 904   }
 905
 906   #ifdef DEBUGMODE
 907     AliInfo("Got TR background density.");
 908   #endif
 909 }
 910
 911 //________________________________________________________________________
 912 void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
 913 {
 914   #ifdef DEBUGMODE
 915     AliInfo("Starting Calculate().");
 916   #endif
 917   ////////////////////// NOTE: initialization & casting
 918
 919   // Additional cuts
 920   FillHistogram("hNumberEvents", 0.5); // number of events before manual cuts
 921
 922   if(!fHelperClass->IsVertexSelected2013pA(event))
 923     return;
 924
 925   FillHistogram("hNumberEvents", 1.5); // number of events after manual cuts
 926
 927   #ifdef DEBUGMODE
 928     AliInfo("Calculate()::Init done.");
 929   #endif
 930
 931   ////////////////////// NOTE: Get Centrality, (Leading)Signal jets and Background
 932
 933   // Get centrality
 934   AliCentrality* tmpCentrality = NULL;
 935   tmpCentrality = event->GetCentrality();
 936   Double_t centralityPercentile = 0.0;
 937   if (tmpCentrality != NULL)
 938     centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
 939
 940   // Get jets
 941   if (fAnalyzeBackground || fAnalyzeJets)
 942     GetSignalJets();
 943
 944   // Get background estimates
 945   Double_t              backgroundKTMedian = -1.0;
 946   Double_t              backgroundRCMean = -1.0;
 947   Double_t              backgroundRCMedian = -1.0;
 948   Double_t              backgroundTRMean = -1.0;
 949   Double_t              backgroundKTAreaMean = -1.0;
 950   Double_t              backgroundTRAreaMean = -1.0;
 951   Double_t              dijetBackground = -1.0;
 952   Double_t              dijetBackgroundPerpendicular = -1.0;
 953
 954   if (fAnalyzeBackground)
 955   {
 956     GetRCBackgroundDensity    (fNumberExcludedJets, backgroundRCMean, backgroundRCMedian);
 957     GetTRBackgroundDensity    (fNumberExcludedJets, backgroundTRMean, backgroundTRAreaMean);
 958     GetKTBackgroundDensity    (fNumberExcludedJets, backgroundKTMedian, backgroundKTAreaMean);
 959   }
 960
 961   #ifdef DEBUGMODE
 962     AliInfo("Calculate()::Centrality&SignalJets&Background-Calculation done.");
 963   #endif
 964
 965   ////////////////////// NOTE: Jet analysis and calculations
 966
 967   if (fAnalyzeJets)
 968   {
 969     // ### SIGNAL JET ANALYSIS
 970     for (Int_t i = 0; i<fNumberSignalJets; i++)
 971     {
 972       AliEmcalJet* tmpJet = fSignalJets[i];
 973
 974       // Jet spectra
 975       FillHistogram("hJetPt", tmpJet->Pt(), centralityPercentile);
 976       FillHistogram("hJetPtBgrdSubtractedRC", GetCorrectedJetPt(tmpJet, backgroundRCMean, kRCEtaCorrection), centralityPercentile);
 977       FillHistogram("hJetPtBgrdSubtractedKT", GetCorrectedJetPt(tmpJet, backgroundKTMedian, kKTEtaCorrection), centralityPercentile);
 978       FillHistogram("hJetPtBgrdSubtractedTR", GetCorrectedJetPt(tmpJet, backgroundTRMean, kTREtaCorrection), centralityPercentile);
 979
 980       FillHistogram("hJetPtBgrdSubtractedRCNoEtaCorr", GetCorrectedJetPt(tmpJet, backgroundRCMean, kNoEtaCorrection), centralityPercentile);
 981       FillHistogram("hJetPtBgrdSubtractedKTNoEtaCorr", GetCorrectedJetPt(tmpJet, backgroundKTMedian, kNoEtaCorrection), centralityPercentile);
 982       FillHistogram("hJetPtBgrdSubtractedTRNoEtaCorr", GetCorrectedJetPt(tmpJet, backgroundTRMean, kNoEtaCorrection), centralityPercentile);
 983
 984       // Signal jet vs. signal jet - "Combinatorial"
 985       for (Int_t j = i+1; j<fNumberSignalJets; j++)
 986         FillHistogram("hJetDeltaPhi", GetDeltaPhi(tmpJet->Phi(), fSignalJets[j]->Phi()));
 987     }
 988
 989     // ### DIJETS
 990     if(fNumberSignalJets >= 2)
 991     {
 992       FillHistogram("hLeadingJetDeltaPhi", GetDeltaPhi(fFirstLeadingJet->Phi(), fSecondLeadingJet->Phi()));
 993
 994       if (IsDijet(fFirstLeadingJet, fSecondLeadingJet))
 995       {
 996         FillHistogram("hDijetConstituentsPt", fFirstLeadingJet->Pt());
 997         FillHistogram("hDijetConstituentsPt", fSecondLeadingJet->Pt());
 998
 999         FillHistogram("hDijetLeadingJetPt", fFirstLeadingJet->Pt());
1000         FillHistogram("hDijetPtCorrelation", fFirstLeadingJet->Pt(), fSecondLeadingJet->Pt());
1001         Double_t dummyArea = 0;
1002         GetTRBackgroundDensity (2, dijetBackground, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kFALSE);
1003         GetTRBackgroundDensity (2, dijetBackgroundPerpendicular, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kTRUE);
1004       }
1005     }
1006
1007     // ### SOME JET PLOTS
1008     FillHistogram("hJetCountAll", fJetArray->GetEntries());
1009     FillHistogram("hJetCountAccepted", fNumberSignalJets);
1010     if (fFirstLeadingJet)
1011       FillHistogram("hLeadingJetPt", fFirstLeadingJet->Pt());
1012     if (fSecondLeadingJet)
1013       FillHistogram("hSecondLeadingJetPt", fSecondLeadingJet->Pt());
1014
1015   } //endif AnalyzeJets
1016
1017   #ifdef DEBUGMODE
1018     AliInfo("Calculate()::Jets done.");
1019   #endif
1020   ////////////////////// NOTE: Background analysis
1021
1022   if (fAnalyzeBackground)
1023   {
1024     // Calculate (non-eta corrected) background in centrality classes
1025     FillHistogram("hRCBackground", backgroundRCMean, centralityPercentile);
1026     FillHistogram("hTRBackground", backgroundTRMean, centralityPercentile);
1027     FillHistogram("hKTBackground", backgroundKTMedian, centralityPercentile);
1028
1029     // Calculate backgrounds in eta bins
1030     for(Int_t i=0;i<5;i++)
1031     {
1032       Double_t dummy = 0.0;
1033       Double_t tmpKTRho = 0.0;
1034       Double_t tmpRCRho = 0.0;
1035       Double_t tmpTRRho = 0.0;
1036       Double_t tmpKTRhoCorr = 0.0;
1037       Double_t tmpRCRhoCorr = 0.0;
1038       Double_t tmpTRRhoCorr = 0.0;
1039
1040       Double_t etaMin = -(fTrackEtaWindow-fSignalJetRadius) + 2*(fTrackEtaWindow-fSignalJetRadius)/5 *  i;
1041       Double_t etaMax = -(fTrackEtaWindow-fSignalJetRadius) + 2*(fTrackEtaWindow-fSignalJetRadius)/5 * (i+1);
1042
1043       // Calculate backgrounds
1044       GetKTBackgroundDensity    (fNumberExcludedJets, tmpKTRho, dummy, etaMin, etaMax);
1045       GetTRBackgroundDensity    (fNumberExcludedJets, tmpTRRho, dummy, etaMin, etaMax);
1046       GetRCBackgroundDensity    (fNumberExcludedJets, tmpRCRho, dummy, etaMin, etaMax);
1047
1048       // Add eta-correction
1049       tmpKTRhoCorr = tmpKTRho * GetBackgroundEtaCorrFactor(kKTEtaCorrection, (etaMin+etaMax)/2.0);
1050       tmpTRRhoCorr = tmpTRRho * GetBackgroundEtaCorrFactor(kTREtaCorrection, (etaMin+etaMax)/2.0);
1051       tmpRCRhoCorr = tmpRCRho * GetBackgroundEtaCorrFactor(kRCEtaCorrection, (etaMin+etaMax)/2.0);
1052
1053       FillHistogram("hRCBackgroundEtaBins", tmpRCRho, (etaMin+etaMax)/2.0);
1054       FillHistogram("hTRBackgroundEtaBins", tmpTRRho, (etaMin+etaMax)/2.0);
1055       FillHistogram("hKTBackgroundEtaBins", tmpKTRho, (etaMin+etaMax)/2.0);
1056       FillHistogram("hRCBackgroundEtaBinsCorrected", tmpRCRhoCorr, (etaMin+etaMax)/2.0);
1057       FillHistogram("hTRBackgroundEtaBinsCorrected", tmpTRRhoCorr, (etaMin+etaMax)/2.0);
1058       FillHistogram("hKTBackgroundEtaBinsCorrected", tmpKTRhoCorr, (etaMin+etaMax)/2.0);
1059
1060       for (Int_t j = 0; j<fNumberSignalJets; j++)
1061       {
1062         AliEmcalJet* tmpJet = fSignalJets[j];
1063         FillHistogram("hJetPtBgrdSubtractedRCEtaBins", GetCorrectedJetPt(tmpJet, tmpRCRho, kNoEtaCorrection), (etaMin+etaMax)/2.0);
1064         FillHistogram("hJetPtBgrdSubtractedKTEtaBins", GetCorrectedJetPt(tmpJet, tmpKTRho, kNoEtaCorrection), (etaMin+etaMax)/2.0);
1065         FillHistogram("hJetPtBgrdSubtractedTREtaBins", GetCorrectedJetPt(tmpJet, tmpTRRho, kNoEtaCorrection), (etaMin+etaMax)/2.0);
1066       }
1067     }
1068
1069     // In case of dijets -> look at the background
1070     if (dijetBackground >= 0)
1071       FillHistogram("hDijetBackground", dijetBackground, centralityPercentile);
1072     if (dijetBackgroundPerpendicular >= 0)
1073       FillHistogram("hDijetBackgroundPerpendicular", dijetBackgroundPerpendicular, centralityPercentile);
1074
1075
1076     // Calculate the delta pt
1077     Double_t tmpDeltaPtKT = 0.0;
1078     Double_t tmpDeltaPtRC = 0.0;
1079     Double_t tmpDeltaPtTR = 0.0;
1080     Double_t tmpDeltaPtKTNoEta = 0.0;
1081     Double_t tmpDeltaPtRCNoEta = 0.0;
1082     Double_t tmpDeltaPtTRNoEta = 0.0;
1083     Double_t tmpDeltaPtKTNoEtaNoExcl = 0.0;
1084     Double_t tmpDeltaPtRCNoEtaNoExcl = 0.0;
1085     Double_t tmpDeltaPtTRNoEtaNoExcl = 0.0;
1086     GetDeltaPt(tmpDeltaPtKT, backgroundKTMedian, kKTEtaCorrection);
1087     GetDeltaPt(tmpDeltaPtRC, backgroundRCMean, kRCEtaCorrection);
1088     GetDeltaPt(tmpDeltaPtTR, backgroundTRMean, kTREtaCorrection);
1089     GetDeltaPt(tmpDeltaPtKTNoEta, backgroundKTMedian, kNoEtaCorrection);
1090     GetDeltaPt(tmpDeltaPtRCNoEta, backgroundRCMean, kNoEtaCorrection);
1091     GetDeltaPt(tmpDeltaPtTRNoEta, backgroundTRMean, kNoEtaCorrection);
1092     GetDeltaPt(tmpDeltaPtKTNoEtaNoExcl, backgroundKTMedian, kNoEtaCorrection, kFALSE);
1093     GetDeltaPt(tmpDeltaPtRCNoEtaNoExcl, backgroundRCMean, kNoEtaCorrection, kFALSE);
1094     GetDeltaPt(tmpDeltaPtTRNoEtaNoExcl, backgroundTRMean, kNoEtaCorrection, kFALSE);
1095
1096     // If valid, fill the delta pt histograms
1097     if(tmpDeltaPtKT > -10000.0)
1098       FillHistogram("hDeltaPtKT", tmpDeltaPtKT, centralityPercentile);
1099     if(tmpDeltaPtRC > -10000.0)
1100       FillHistogram("hDeltaPtRC", tmpDeltaPtRC, centralityPercentile);
1101     if(tmpDeltaPtTR > -10000.0)
1102       FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile);
1103     if(tmpDeltaPtKTNoEta > -10000.0)
1104       FillHistogram("hDeltaPtKTNoEtaCorr", tmpDeltaPtKTNoEta, centralityPercentile);
1105     if(tmpDeltaPtRCNoEta > -10000.0)
1106       FillHistogram("hDeltaPtRCNoEtaCorr", tmpDeltaPtRCNoEta, centralityPercentile);
1107     if(tmpDeltaPtTRNoEta > -10000.0)
1108       FillHistogram("hDeltaPtTRNoEtaCorr", tmpDeltaPtTRNoEta, centralityPercentile);
1109     if(tmpDeltaPtKTNoEtaNoExcl > -10000.0)
1110       FillHistogram("hDeltaPtKTNoEtaCorrNoExcl", tmpDeltaPtKTNoEtaNoExcl, centralityPercentile);
1111     if(tmpDeltaPtRCNoEtaNoExcl > -10000.0)
1112       FillHistogram("hDeltaPtRCNoEtaCorrNoExcl", tmpDeltaPtRCNoEtaNoExcl, centralityPercentile);
1113     if(tmpDeltaPtTRNoEtaNoExcl > -10000.0)
1114       FillHistogram("hDeltaPtTRNoEtaCorrNoExcl", tmpDeltaPtTRNoEtaNoExcl, centralityPercentile);
1115
1116   }
1117
1118   #ifdef DEBUGMODE
1119     AliInfo("Calculate()::Background done.");
1120   #endif
1121
1122   ////////////////////// NOTE: Pythia histograms
1123   if(fAnalyzePythia)
1124   {
1125     FillHistogram("hPythiaPtHard", GetPtHard());
1126     FillHistogram("hPythiaNTrials", GetPtHardBin()-0.1, fTrials);
1127     FillHistogram("hPythiaXSection", GetPtHardBin()-0.1, fCrossSection);
1128
1129     #ifdef DEBUGMODE
1130       AliInfo("Calculate()::Pythia done.");
1131     #endif
1132   }
1133   #ifdef DEBUGMODE
1134     AliInfo("Calculate() done.");
1135   #endif
1136 }
1137
1138 //________________________________________________________________________
1139 Bool_t AliAnalysisTaskChargedJetsPA::Notify()
1140 {
1141   // Implemented Notify() to read the cross sections
1142   // and number of trials from pyxsec.root
1143   //
1144   #ifdef DEBUGMODE
1145     AliInfo("Notify started.");
1146   #endif
1147
1148   if(fAnalyzePythia)
1149   {
1150     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1151     TFile *currFile = tree->GetCurrentFile();
1152
1153     TString file(currFile->GetName());
1154
1155     if(file.Contains("root_archive.zip#")){
1156       Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1157       Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1158       file.Replace(pos+1,20,"");
1159     }
1160     else {
1161       // not an archive take the basename....
1162       file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1163     }
1164
1165     TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
1166     if(!fxsec){
1167       // next trial fetch the histgram file
1168       fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1169       if(!fxsec){
1170           // not a severe condition but inciate that we have no information
1171         return kFALSE;
1172       }
1173       else{
1174         // find the tlist we want to be independtent of the name so use the Tkey
1175         TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1176         if(!key){
1177           fxsec->Close();
1178           return kFALSE;
1179         }
1180         TList *list = dynamic_cast<TList*>(key->ReadObj());
1181         if(!list){
1182           fxsec->Close();
1183           return kFALSE;
1184         }
1185         fCrossSection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1186         fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1187         fxsec->Close();
1188       }
1189     } // no tree pyxsec.root
1190     else {
1191       TTree *xtree = (TTree*)fxsec->Get("Xsection");
1192       if(!xtree){
1193         fxsec->Close();
1194         return kFALSE;
1195       }
1196       UInt_t   ntrials  = 0;
1197       Double_t  xsection  = 0;
1198       xtree->SetBranchAddress("xsection",&xsection);
1199       xtree->SetBranchAddress("ntrials",&ntrials);
1200       xtree->GetEntry(0);
1201       fTrials = ntrials;
1202       fCrossSection = xsection;
1203       fxsec->Close();
1204     }
1205     #ifdef DEBUGMODE
1206       AliInfo("Notify ended.");
1207     #endif
1208   }
1209   return kTRUE;
1210 }
1211
1212 //________________________________________________________________________
1213 void AliAnalysisTaskChargedJetsPA::SetKTEtaCorrectionFactors(TH1D* histo)
1214 {
1215   // COPY given histogram
1216   fJetKTEtaCorrection = new TH1D(*histo);
1217
1218   if (!fJetKTEtaCorrection)
1219     AliError(Form("Setting the correction factors with %s (%s) failed! You won't get eta-corrected spectra!", histo->GetName(), histo->IsA()->GetName()));
1220
1221   // Look, if given histogram is compatible with given code
1222   if (fJetKTEtaCorrection->GetXaxis()->GetNbins() != 5)
1223     AliError(Form("Setting the correction factors failed, because the given histogram is not compatible! You need nbinX=5 (currently:%d)",fJetKTEtaCorrection->GetXaxis()->GetNbins()));
1224 }
1225
1226 //________________________________________________________________________
1227 void AliAnalysisTaskChargedJetsPA::SetRCEtaCorrectionFactors(TH1D* histo)
1228 {
1229   // COPY given histogram
1230   fJetRCEtaCorrection = new TH1D(*histo);
1231
1232   if (!fJetRCEtaCorrection)
1233     AliError(Form("Setting the correction factors with %s (%s) failed! You won't get eta-corrected spectra!", histo->GetName(), histo->IsA()->GetName()));
1234
1235   // Look, if given histogram is compatible with given code
1236   if (fJetRCEtaCorrection->GetXaxis()->GetNbins() != 5)
1237     AliError(Form("Setting the correction factors failed, because the given histogram is not compatible! You need nbinX=5 (currently:%d)",fJetRCEtaCorrection->GetXaxis()->GetNbins()));
1238 }
1239
1240 //________________________________________________________________________
1241 void AliAnalysisTaskChargedJetsPA::SetTREtaCorrectionFactors(TH1D* histo)
1242 {
1243   // COPY given histogram
1244   fJetTREtaCorrection = new TH1D(*histo);
1245
1246   if (!fJetTREtaCorrection)
1247     AliError(Form("Setting the correction factors with %s (%s) failed! You won't get eta-corrected spectra!", histo->GetName(), histo->IsA()->GetName()));
1248
1249   // Look, if given histogram is compatible with given code
1250   if (fJetTREtaCorrection->GetXaxis()->GetNbins() != 5)
1251     AliError(Form("Setting the correction factors failed, because the given histogram is not compatible! You need nbinX=5 (currently:%d)",fJetTREtaCorrection->GetXaxis()->GetNbins()));
1252 }
1253
1254 //________________________________________________________________________
1255 inline Double_t AliAnalysisTaskChargedJetsPA::EtaToTheta(Double_t arg)
1256   {return 2.*atan(exp(-arg));}
1257 //________________________________________________________________________
1258 inline Double_t AliAnalysisTaskChargedJetsPA::ThetaToEta(Double_t arg)
1259 {
1260   if ((arg > TMath::Pi()) || (arg < 0.0))
1261   {
1262     AliError(Form("ThetaToEta got wrong input! (%f)", arg));
1263     return 0.0;
1264   }
1265   return -log(tan(arg/2.));
1266 }
1267 //________________________________________________________________________
1268 inline Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPhi(Double_t phi1, Double_t phi2)
1269   {return min(TMath::Abs(phi1-phi2),TMath::TwoPi() - TMath::Abs(phi1-phi2));}
1270
1271 //________________________________________________________________________
1272 Double_t AliAnalysisTaskChargedJetsPA::MCGetOverlapCircleRectancle(Double_t cPosX, Double_t cPosY, Double_t cRadius, Double_t rPosXmin, Double_t rPosXmax, Double_t rPosYmin, Double_t rPosYmax)
1273 {
1274   const Int_t kTests = 1000;
1275   Int_t hits = 0;
1276   TRandom3 randomGen(0);
1277
1278   // Loop over kTests-many tests
1279   for (Int_t i=0; i<kTests; i++)
1280   {
1281     //Choose random position in rectangle for the tester
1282     Double_t tmpTestX = randomGen.Uniform(rPosXmin, rPosXmax);
1283     Double_t tmpTestY = randomGen.Uniform(rPosYmin, rPosYmax);
1284
1285     //Check, if tester is in circle. If yes, increment circle counter.
1286     Double_t tmpDistance = TMath::Sqrt( (tmpTestX - cPosX)*(tmpTestX - cPosX) + (tmpTestY - cPosY)*(tmpTestY - cPosY) );
1287     if(tmpDistance < cRadius)
1288       hits++;
1289   }
1290
1291   // return ratio
1292   return (static_cast<Double_t>(hits)/static_cast<Double_t>(kTests));
1293 }
1294
1295 //________________________________________________________________________
1296 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x)
1297 {
1298   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(key)));
1299   if(!tmpHist)
1300   {
1301     AliWarning(Form("Cannot find histogram <%s> ",key)) ;
1302     return;
1303   }
1304
1305   tmpHist->Fill(x);
1306 }
1307
1308 //________________________________________________________________________
1309 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y)
1310 {
1311   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(key)));
1312   if(!tmpHist)
1313   {
1314     AliWarning(Form("Cannot find histogram <%s> ",key));
1315     return;
1316   }
1317
1318   if (tmpHist->IsA()->GetBaseClass("TH1"))
1319     static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1320   else if (tmpHist->IsA()->GetBaseClass("TH2"))
1321     static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1322 }
1323
1324 //________________________________________________________________________
1325 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y, Double_t add)
1326 {
1327   TH2* tmpHist = static_cast<TH2*>(fOutputList->FindObject(GetHistoName(key)));
1328   if(!tmpHist)
1329   {
1330     AliWarning(Form("Cannot find histogram <%s> ",key));
1331     return;
1332   }
1333
1334   tmpHist->Fill(x,y,add);
1335 }
1336 //________________________________________________________________________
1337 template <class T> T* AliAnalysisTaskChargedJetsPA::AddHistogram1D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, const char* xTitle, const char* yTitle)
1338 {
1339   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax);
1340
1341   tmpHist->GetXaxis()->SetTitle(xTitle);
1342   tmpHist->GetYaxis()->SetTitle(yTitle);
1343   tmpHist->SetOption(options);
1344   tmpHist->SetMarkerStyle(kFullCircle);
1345   tmpHist->Sumw2();
1346
1347   fHistList->Add(tmpHist);
1348   fHistCount++;
1349
1350   return tmpHist;
1351 }
1352
1353 //________________________________________________________________________
1354 template <class T> T* AliAnalysisTaskChargedJetsPA::AddHistogram2D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, const char* xTitle, const char* yTitle, const char* zTitle)
1355 {
1356   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax, yBins, yMin, yMax);
1357   tmpHist->GetXaxis()->SetTitle(xTitle);
1358   tmpHist->GetYaxis()->SetTitle(yTitle);
1359   tmpHist->GetZaxis()->SetTitle(zTitle);
1360   tmpHist->SetOption(options);
1361   tmpHist->SetMarkerStyle(kFullCircle);
1362   tmpHist->Sumw2();
1363
1364   fHistList->Add(tmpHist);
1365   fHistCount++;
1366
1367   return tmpHist;
1368 }
1369
1370 //________________________________________________________________________
1371 void AliAnalysisTaskChargedJetsPA::Terminate(Option_t *)
1372 {
1373   PostData(1, fOutputList);
1374
1375   // Mandatory
1376   fOutputList = dynamic_cast<TList*> (GetOutputData(1)); // '1' refers to the output slot
1377   if (!fOutputList) {
1378     printf("ERROR: Output list not available\n");
1379     return;
1380   }
1381 }
1382
1383 //________________________________________________________________________
1384 AliAnalysisTaskChargedJetsPA::~AliAnalysisTaskChargedJetsPA()
1385 {
1386   // Destructor. Clean-up the output list, but not the histograms that are put inside
1387   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
1388   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
1389     delete fOutputList;
1390   }
1391 }
1392
1393 //________________________________________________________________________
1394 void AliAnalysisTaskChargedJetsPA::UserCreateOutputObjects()
1395 {
1396   // called once to create user defined output objects like histograms, plots etc.
1397   // and to put it on the output list.
1398   // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
1399
1400   fRandom = new TRandom3(0);
1401
1402   fOutputList = new TList();
1403   fOutputList->SetOwner(); // otherwise it produces leaks in merging
1404
1405   PostData(1, fOutputList);
1406 }
1407
1408 //________________________________________________________________________
1409 void AliAnalysisTaskChargedJetsPA::UserExec(Option_t *)
1410 {
1411   #ifdef DEBUGMODE
1412     AliInfo("UserExec() started.");
1413   #endif
1414
1415   if (!InputEvent())
1416   {
1417     AliError("??? Event pointer == 0 ???");
1418     return;
1419   }
1420
1421   if (!fInitialized)
1422     ExecOnce(); // Get tracks, jets, background from arrays if not already given + Init Histos
1423
1424   Calculate(InputEvent());
1425
1426   PostData(1, fOutputList);
1427 }