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 <TProfile2D.h>
  11 #include <TH1F.h>
  12 #include <TH2F.h>
  13 #include <TCanvas.h>
  14 #include <TList.h>
  15 #include <TClonesArray.h>
  16 #include <TObject.h>
  17 #include <TMath.h>
  18 #include <TSystem.h>
  19 #include <TInterpreter.h>
  20 #include <TH1.h>
  21 #include "AliAnalysisTask.h"
  22 #include "AliCentrality.h"
  23 #include "AliStack.h"
  24 #include "AliESDEvent.h"
  25 #include "AliESDInputHandler.h"
  26 #include "AliAODEvent.h"
  27 #include "AliAODHandler.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliAnalysisTaskSE.h"
  30 #endif
  31
  32 #include <time.h>
  33 #include <TRandom3.h>
  34 #include "AliGenPythiaEventHeader.h"
  35 #include "AliAODMCHeader.h"
  36 #include "AliMCEvent.h"
  37 #include "AliLog.h"
  38 #include <AliEmcalJet.h>
  39 #include <AliPicoTrack.h>
  40 #include "AliVEventHandler.h"
  41 #include "AliVParticle.h"
  42 #include "AliAODMCParticle.h"
  43 #include "AliAnalysisUtils.h"
  44 #include "AliRhoParameter.h"
  45 #include "TVector3.h"
  46
  47 #include "AliAnalysisTaskChargedJetsPA.h"
  48 using std::min;
  49
  50 //TODO: Not accessing the particles when using MC
  51 //TODO: FillHistogram can be done better with virtual TH1(?)
  52 ClassImp(AliAnalysisTaskChargedJetsPA)
  53
  54 // ######################################################################################## DEFINE HISTOGRAMS
  55 void AliAnalysisTaskChargedJetsPA::Init()
  56 {
  57   #ifdef DEBUGMODE
  58     AliInfo("Creating histograms.");
  59   #endif
  60
  61   TH1D* tmpHisto = AddHistogram1D<TH1D>("hNumberEvents", "Number of events (0 = before cuts, 1 = after cuts)", "", 2, 0, 2, "stage","N^{Events}/cut");
  62   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
  63   tmpHisto->GetXaxis()->SetBinLabel(2, "After cuts");
  64
  65   tmpHisto = AddHistogram1D<TH1D>("hEventAcceptance", "Accepted events (0 = before cuts, 1 = after pile up, 2 = after vertex)", "", 3, 0, 3, "stage","N^{Events}/cut");
  66   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
  67   tmpHisto->GetXaxis()->SetBinLabel(2, "After pile up");
  68   tmpHisto->GetXaxis()->SetBinLabel(3, "After vertex");
  69
  70   tmpHisto = AddHistogram1D<TH1D>("hTrackAcceptance", "Accepted tracks (0 = before cuts, 1 = after eta, 2 = after pT)", "", 3, 0, 3, "stage","N^{Tracks}/cut");
  71   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
  72   tmpHisto->GetXaxis()->SetBinLabel(2, "After eta");
  73   tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}");
  74
  75   tmpHisto = AddHistogram1D<TH1D>("hJetAcceptance", "Accepted jets (0 = before cuts, 1 = after eta, 2 = after pT, 3 = after area)", "", 4, 0, 4, "stage","N^{Jets}/cut");
  76   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
  77   tmpHisto->GetXaxis()->SetBinLabel(2, "After eta");
  78   tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}");
  79   tmpHisto->GetXaxis()->SetBinLabel(4, "After area");
  80
  81   // NOTE: Jet histograms
  82   if (fAnalyzeJets)
  83   {
  84     // ######## Jet spectra
  85     AddHistogram1D<TH1D>("hRawJetPt", "Raw jets p_{T} distribution (before cuts)", "", 500, 0., 250., "p_{T} (GeV/c)", "dN^{Jets}/dp_{T}");
  86     AddHistogram2D<TH2D>("hJetPt", "Jets p_{T} distribution", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  87     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTImprovedCMS", "Jets p_{T} distribution, KT background (Improved CMS) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  88     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedPP", "Jets p_{T} distribution, pp background subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  89     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedExternal", "Jets p_{T} distribution, external bgrd. subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  90
  91
  92
  93     AddHistogram2D<TProfile2D>("hJetPtSubtractedRhoKTImprovedCMS", "Mean subtracted KT (CMS w/o signal) background from jets", "COLZ", 600, 0, 150, fNumberOfCentralityBins, 0, 100, "Jet p_{T}", "Centrality", "#rho mean");
  94
  95     if(fAnalyzeDeprecatedBackgrounds)
  96     {
  97       AddHistogram2D<TH2D>("hJetPtBgrdSubtractedTR", "Jets p_{T} distribution, TR background (Cone R=0.6 around jets excluded) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  98       AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTPbPb", "Jets p_{T} distribution, KT background (PbPb w/o ghosts) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
  99       AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTPbPbWithGhosts", "Jets p_{T} distribution, KT background (PbPb w/ ghosts) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
 100       AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTCMS", "Jets p_{T} distribution, KT background (CMS) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
 101       AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTMean", "Jets p_{T} distribution, KT background (Mean) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
 102       AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTTrackLike", "Jets p_{T} distribution, KT background (track-like) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
 103     }
 104
 105     // ######## Jet profiles
 106     if(fAnalyzeJetProfile)
 107     {
 108       AddHistogram2D<TH2D>("hJetProfile10GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 10 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
 109       AddHistogram2D<TH2D>("hJetProfile20GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 20 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
 110       AddHistogram2D<TH2D>("hJetProfile30GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 30 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
 111       AddHistogram2D<TH2D>("hJetProfile40GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 40 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
 112       AddHistogram2D<TH2D>("hJetProfile50GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 50 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
 113       AddHistogram2D<TH2D>("hJetProfile60GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 60 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
 114       AddHistogram2D<TH2D>("hJetProfile70GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 70 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
 115     }
 116     // ######## Jet stuff
 117     AddHistogram2D<TH2D>("hJetConstituentPt", "Jet constituents p_{T} distribution", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
 118     AddHistogram1D<TH1D>("hJetCountAll", "Number of Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
 119     AddHistogram1D<TH1D>("hJetCountAccepted", "Number of accepted Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
 120     AddHistogram2D<TH2D>("hJetCount", "Correlation jets/accepted jets", "", 200, 0., 200., 200, 0., 200., "N jets","N jets accepted", "d^{2}N^{Events}/dN^{Jets dN^{Jets, acc}}");
 121     AddHistogram1D<TH1D>("hLeadingJetPt", "Leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 122     AddHistogram1D<TH1D>("hSecondLeadingJetPt", "Second leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 123     AddHistogram1D<TH1D>("hCorrectedLeadingJetPt", "Corrected leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 124     AddHistogram1D<TH1D>("hCorrectedSecondLeadingJetPt", "Corrected second leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 125     AddHistogram1D<TH1D>("hJetDeltaPhi", "Jets combinatorial #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
 126     AddHistogram1D<TH1D>("hLeadingJetDeltaPhi", "1st and 2nd leading jet #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
 127
 128     // ########## Dijet stuff
 129     AddHistogram1D<TH1D>("hDijetLeadingJetPt", "Dijet leading jet p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 130     AddHistogram1D<TH1D>("hDijetConstituentsPt", "Dijet constituents p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 131     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}");
 132   }
 133
 134   // NOTE: Jet background histograms
 135   if (fAnalyzeBackground)
 136   {
 137     // ########## Default background estimates
 138     AddHistogram2D<TH2D>("hKTBackgroundImprovedCMS", "KT background density (Improved CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 139     AddHistogram2D<TH2D>("hKTBackgroundImprovedCMSExternal", "KT background density (Improved CMS approach from external task)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 140     AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMS", "Background fluctuations #delta p_{T} (KT, Improved CMS-like)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 141     AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMSPartialExclusion", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, partial jet exclusion)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 142     AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMSPartialExclusion_Signal", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, partial jet exclusion w/ 1/N_{sig} probability)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 143     AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMSFullExclusion", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, full leading jet exclusion)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 144     AddHistogram2D<TH2D>("hDeltaPtNoBackground", "Background fluctuations #delta p_{T} (No background)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 145     AddHistogram2D<TH2D>("hDeltaPtNoBackgroundNoEmptyCones", "Background fluctuations #delta p_{T} (No background, no empty cones)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 146
 147     AddHistogram1D<TProfile>("hKTMeanBackgroundImprovedCMS", "KT background mean (Improved CMS approach)", "", 100, 0, 100, "Centrality", "#rho mean");
 148
 149     AddHistogram2D<TH2D>("hDijetBackground", "Background density (dijets excluded)", "", 200, 0., 20., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 150     AddHistogram2D<TH2D>("hDijetBackgroundPerpendicular", "Background density (dijets excluded)", "", 200, 0., 20., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 151     AddHistogram2D<TH2D>("hDijetBackgroundPerpendicularDeltaPt", "Background fluctuations #delta p_{T} (dijets excluded)", "", 1201, -80.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 152
 153     if(fAnalyzeDeprecatedBackgrounds)
 154     {
 155       // ########## Different background estimates
 156       AddHistogram2D<TH2D>("hKTBackgroundPbPb", "KT background density (PbPb approach, no ghosts)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 157       AddHistogram2D<TH2D>("hKTBackgroundPbPbWithGhosts", "KT background density (PbPb approach w/ ghosts)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 158       AddHistogram2D<TH2D>("hKTBackgroundCMS", "KT background density (CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 159       AddHistogram2D<TH2D>("hKTBackgroundMean", "KT background density (Mean approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 160       AddHistogram2D<TH2D>("hKTBackgroundTrackLike", "KT background density (Track-like approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 161
 162       AddHistogram2D<TH2D>("hTRBackgroundNoExcl", "TR background density (No signal excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 163       AddHistogram2D<TH2D>("hTRBackgroundCone02", "TR background density (Cones R=0.2 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 164       AddHistogram2D<TH2D>("hTRBackgroundCone04", "TR background density (Cones R=0.4 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 165       AddHistogram2D<TH2D>("hTRBackgroundCone06", "TR background density (Cones R=0.6 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 166       AddHistogram2D<TH2D>("hTRBackgroundCone08", "TR background density (Cones R=0.8 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 167       AddHistogram2D<TH2D>("hTRBackgroundExact",  "TR background density (signal jets exactly excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 168
 169       // ########## Delta Pt
 170       AddHistogram2D<TH2D>("hDeltaPtKTPbPb", "Background fluctuations #delta p_{T} (KT, PbPb w/o ghosts)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 171       AddHistogram2D<TH2D>("hDeltaPtKTPbPbWithGhosts", "Background fluctuations #delta p_{T} (KT, PbPb w/ ghosts)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 172       AddHistogram2D<TH2D>("hDeltaPtKTCMS", "Background fluctuations #delta p_{T} (KT, CMS-like)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 173       AddHistogram2D<TH2D>("hDeltaPtKTMean", "Background fluctuations #delta p_{T} (KT, Mean)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 174       AddHistogram2D<TH2D>("hDeltaPtKTTrackLike", "Background fluctuations #delta p_{T} (KT, track-like)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 175       AddHistogram2D<TH2D>("hDeltaPtTR", "Background fluctuations #delta p_{T} (TR, cone R=0.6)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
 176
 177       // ########## Profiles for background means vs. centrality
 178       AddHistogram1D<TProfile>("hKTMeanBackgroundPbPb", "KT background mean (PbPb approach w/o ghosts)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
 179       AddHistogram1D<TProfile>("hKTMeanBackgroundPbPbWithGhosts", "KT background mean (PbPb approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
 180       AddHistogram1D<TProfile>("hKTMeanBackgroundCMS", "KT background mean (CMS approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
 181       AddHistogram1D<TProfile>("hKTMeanBackgroundMean", "KT background mean (Mean approach)", "",  fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
 182       AddHistogram1D<TProfile>("hKTMeanBackgroundTPC", "KT background mean (Track-like approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
 183       AddHistogram1D<TProfile>("hTRMeanBackground", "TR background mean", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
 184     }
 185   }
 186   // NOTE: Jet constituent correlations
 187   if(fAnalyzeMassCorrelation)
 188   {
 189     AddHistogram1D<TH1D>("hJetMassFromConstituents", "Jet mass by mass of charged constituents", "", 200, 0., 200., "N jets","dN^{Jets}/dN^{mass}");
 190     AddHistogram1D<TH1D>("hJetMass", "Jet mass from fastjet", "", 200, 0., 200., "N jets","dN^{Jets}/dN^{mass}");
 191
 192     AddHistogram2D<TH2D>("hJetPtVsMass", "Correlation jet pt/summed constituent jet", "", 400, 0., 100., 400, 0., 100., "p_{T}","Mass", "d^{2}N}/dN^{p_{T}dM");
 193     AddHistogram2D<TH2D>("hJetPtVsJetMass", "Correlation jet pt/summed jet", "", 400, 0., 100., 400, 0., 100., "p_{T}","Mass", "d^{2}N}/dN^{p_{T}dM");
 194
 195     AddHistogram2D<TH2D>("hJetPtVsProtonCount", "Correlation jet pt/amount of proton in jet", "", 400, 0., 100., 20, 0., 1., "p_{T}","Proton percentage", "d^{2}N}/dN^{p_{T}dPercentage");
 196     AddHistogram2D<TH2D>("hJetPtVsPionCount", "Correlation jet pt/amount of pion in jet", "", 400, 0., 100., 20, 0., 1., "p_{T}","Pion percentage", "d^{2}N}/dN^{p_{T}dPercentage");
 197     AddHistogram2D<TH2D>("hJetPtVsKaonCount", "Correlation jet pt/amount of kaon in jet", "", 400, 0., 100., 20, 0., 1., "p_{T}","Kaon percentage", "d^{2}N}/dN^{p_{T}dPercentage");
 198     AddHistogram2D<TH2D>("hJetPtVsElectronCount", "Correlation jet pt/amount of proton in jet", "", 400, 0., 100., 20, 0., 1., "p_{T}","Electron percentage", "d^{2}N}/dN^{p_{T}dPercentage");
 199     AddHistogram2D<TH2D>("hJetPtVsOthersCount", "Correlation jet pt/amount of other particles in jet", "", 400, 0., 100., 20, 0., 1., "p_{T}","Others percentage", "d^{2}N}/dN^{p_{T}dPercentage");
 200     AddHistogram1D<TH1D>("hJetConstituentProtonCount", "Proton count in jets", "", 100, 0., 100., "N protons","dN^{Jets}/dN^{protons}");
 201     AddHistogram1D<TH1D>("hJetConstituentPionCount", "Pion count in jets", "", 100, 0., 100., "N pions","dN^{Jets}/dN^{pions}");
 202     AddHistogram1D<TH1D>("hJetConstituentKaonCount", "Kaon count in jets", "", 100, 0., 100., "N kaons","dN^{Jets}/dN^{kaons}");
 203     AddHistogram1D<TH1D>("hJetConstituentElectronCount", "Electron count in jets", "", 100, 0., 100., "N electrons","dN^{Jets}/dN^{electrons}");
 204     AddHistogram1D<TH1D>("hJetConstituentOthersCount", "Others count in jets", "", 100, 0., 100., "N others","dN^{Jets}/dN^{others}");
 205
 206     AddHistogram2D<TH2D>("hJetPtVsMass_6_14", "Correlation jet pt/summed constituent jet", "", 400, 0., 100., 400, 0., 100., "p_{T}","Mass", "d^{2}N}/dN^{p_{T}dM");
 207     AddHistogram2D<TH2D>("hJetPtVsJetMass_6_14", "Correlation jet pt/summed jet", "", 400, 0., 100., 400, 0., 100., "p_{T}","Mass", "d^{2}N}/dN^{p_{T}dM");
 208     AddHistogram2D<TH2D>("hJetPtVsProtonCount_6_14", "Correlation jet pt/amount of proton in jet", "", 400, 0., 100., 8, 6., 14., "p_{T}","Proton count", "d^{2}N}/dN^{p_{T}dN");
 209     AddHistogram2D<TH2D>("hJetPtVsPionCount_6_14", "Correlation jet pt/amount of pion in jet", "", 400, 0., 100., 8, 6., 14., "p_{T}","Pion count", "d^{2}N}/dN^{p_{T}dN");
 210     AddHistogram2D<TH2D>("hJetPtVsKaonCount_6_14", "Correlation jet pt/amount of kaon in jet", "", 400, 0., 100., 8, 6., 14., "p_{T}","Kaon count", "d^{2}N}/dN^{p_{T}dN");
 211     AddHistogram2D<TH2D>("hJetPtVsElectronCount_6_14", "Correlation jet pt/amount of proton in jet", "", 400, 0., 100., 8, 6., 14., "p_{T}","Electron count", "d^{2}N}/dN^{p_{T}dN");
 212     AddHistogram2D<TH2D>("hJetPtVsOthersCount_6_14", "Correlation jet pt/amount of other particles in jet", "", 400, 0., 100., 8, 6., 14., "p_{T}","Others count", "d^{2}N}/dN^{p_{T}dN");
 213     AddHistogram1D<TH1D>("hJetConstituentProtonCount_6_14", "Proton count in jets", "", 100, 0., 100., "N protons","dN^{Jets}/dN^{protons}");
 214     AddHistogram1D<TH1D>("hJetConstituentPionCount_6_14", "Pion count in jets", "", 100, 0., 100., "N pions","dN^{Jets}/dN^{pions}");
 215     AddHistogram1D<TH1D>("hJetConstituentKaonCount_6_14", "Kaon count in jets", "", 100, 0., 100., "N kaons","dN^{Jets}/dN^{kaons}");
 216     AddHistogram1D<TH1D>("hJetConstituentElectronCount_6_14", "Electron count in jets", "", 100, 0., 100., "N electrons","dN^{Jets}/dN^{electrons}");
 217     AddHistogram1D<TH1D>("hJetConstituentOthersCount_6_14", "Others count in jets", "", 100, 0., 100., "N others","dN^{Jets}/dN^{others}");
 218
 219     AddHistogram2D<TH2D>("hJetPtVsMass_2_X", "Correlation jet pt/summed constituent jet", "", 400, 0., 100., 400, 0., 100., "p_{T}","Mass", "d^{2}N}/dN^{p_{T}dM");
 220     AddHistogram2D<TH2D>("hJetPtVsJetMass_2_X", "Correlation jet pt/summed jet", "", 400, 0., 100., 400, 0., 100., "p_{T}","Mass", "d^{2}N}/dN^{p_{T}dM");
 221     AddHistogram2D<TH2D>("hJetPtVsProtonCount_2_X", "Correlation jet pt/amount of proton in jet", "", 400, 0., 100., 68, 2., 70., "p_{T}","Proton count", "d^{2}N}/dN^{p_{T}dN");
 222     AddHistogram2D<TH2D>("hJetPtVsPionCount_2_X", "Correlation jet pt/amount of pion in jet", "", 400, 0., 100., 68, 2., 70., "p_{T}","Pion count", "d^{2}N}/dN^{p_{T}dN");
 223     AddHistogram2D<TH2D>("hJetPtVsKaonCount_2_X", "Correlation jet pt/amount of kaon in jet", "", 400, 0., 100., 68, 2., 70., "p_{T}","Kaon count", "d^{2}N}/dN^{p_{T}dN");
 224     AddHistogram2D<TH2D>("hJetPtVsElectronCount_2_X", "Correlation jet pt/amount of proton in jet", "", 400, 0., 100., 68, 2., 70., "p_{T}","Electron count", "d^{2}N}/dN^{p_{T}dN");
 225     AddHistogram2D<TH2D>("hJetPtVsOthersCount_2_X", "Correlation jet pt/amount of other particles in jet", "", 400, 0., 100., 68, 2., 70., "p_{T}","Others count", "d^{2}N}/dN^{p_{T}dN");
 226     AddHistogram1D<TH1D>("hJetConstituentProtonCount_2_X", "Proton count in jets", "", 100, 0., 100., "N protons","dN^{Jets}/dN^{protons}");
 227     AddHistogram1D<TH1D>("hJetConstituentPionCount_2_X", "Pion count in jets", "", 100, 0., 100., "N pions","dN^{Jets}/dN^{pions}");
 228     AddHistogram1D<TH1D>("hJetConstituentKaonCount_2_X", "Kaon count in jets", "", 100, 0., 100., "N kaons","dN^{Jets}/dN^{kaons}");
 229     AddHistogram1D<TH1D>("hJetConstituentElectronCount_2_X", "Electron count in jets", "", 100, 0., 100., "N electrons","dN^{Jets}/dN^{electrons}");
 230     AddHistogram1D<TH1D>("hJetConstituentOthersCount_2_X", "Others count in jets", "", 100, 0., 100., "N others","dN^{Jets}/dN^{others}");
 231
 232     AddHistogram2D<TH2D>("hJetPtVsMass_2_7", "Correlation jet pt/summed constituent jet", "", 400, 0., 100., 400, 0., 100., "p_{T}","Mass", "d^{2}N}/dN^{p_{T}dM");
 233     AddHistogram2D<TH2D>("hJetPtVsJetMass_2_7", "Correlation jet pt/summed jet", "", 400, 0., 100., 400, 0., 100., "p_{T}","Mass", "d^{2}N}/dN^{p_{T}dM");
 234     AddHistogram2D<TH2D>("hJetPtVsProtonCount_2_7", "Correlation jet pt/amount of proton in jet", "", 400, 0., 100., 6, 2., 8., "p_{T}","Proton count", "d^{2}N}/dN^{p_{T}dN");
 235     AddHistogram2D<TH2D>("hJetPtVsPionCount_2_7", "Correlation jet pt/amount of pion in jet", "", 400, 0., 100., 6, 2., 8., "p_{T}","Pion count", "d^{2}N}/dN^{p_{T}dN");
 236     AddHistogram2D<TH2D>("hJetPtVsKaonCount_2_7", "Correlation jet pt/amount of kaon in jet", "", 400, 0., 100., 6, 2., 8., "p_{T}","Kaon count", "d^{2}N}/dN^{p_{T}dN");
 237     AddHistogram2D<TH2D>("hJetPtVsElectronCount_2_7", "Correlation jet pt/amount of proton in jet", "", 400, 0., 100., 6, 2., 8., "p_{T}","Electron count", "d^{2}N}/dN^{p_{T}dN");
 238     AddHistogram2D<TH2D>("hJetPtVsOthersCount_2_7", "Correlation jet pt/amount of other particles in jet", "", 400, 0., 100., 6, 2., 8., "p_{T}","Others count", "d^{2}N}/dN^{p_{T}dN");
 239     AddHistogram1D<TH1D>("hJetConstituentProtonCount_2_7", "Proton count in jets", "", 100, 0., 100., "N protons","dN^{Jets}/dN^{protons}");
 240     AddHistogram1D<TH1D>("hJetConstituentPionCount_2_7", "Pion count in jets", "", 100, 0., 100., "N pions","dN^{Jets}/dN^{pions}");
 241     AddHistogram1D<TH1D>("hJetConstituentKaonCount_2_7", "Kaon count in jets", "", 100, 0., 100., "N kaons","dN^{Jets}/dN^{kaons}");
 242     AddHistogram1D<TH1D>("hJetConstituentElectronCount_2_7", "Electron count in jets", "", 100, 0., 100., "N electrons","dN^{Jets}/dN^{electrons}");
 243     AddHistogram1D<TH1D>("hJetConstituentOthersCount_2_7", "Others count in jets", "", 100, 0., 100., "N others","dN^{Jets}/dN^{others}");
 244
 245   }
 246
 247   // NOTE: Track & Cluster & QA histograms
 248   if (fAnalyzeQA)
 249   {
 250     AddHistogram1D<TH1D>("hVertexX", "X distribution of the vertex", "", 2000, -1., 1., "#Delta x(cm)","dN^{Events}/dx");
 251     AddHistogram1D<TH1D>("hVertexY", "Y distribution of the vertex", "", 2000, -1., 1., "#Delta y(cm)","dN^{Events}/dy");
 252     AddHistogram2D<TH2D>("hVertexXY", "XY distribution of the vertex", "COLZ", 500, -1., 1., 500, -1., 1.,"#Delta x(cm)", "#Delta y(cm)","dN^{Events}/dxdy");
 253     AddHistogram1D<TH1D>("hVertexZBeforeVertexCut", "Z distribution of the vertex (before std. vertex cut)", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz");
 254     AddHistogram1D<TH1D>("hVertexZAfterVertexCut", "Z distribution of the vertex (after std. vertex cut)", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz");
 255     AddHistogram1D<TH1D>("hVertexR", "R distribution of the vertex", "", 100, 0., 1., "#Delta r(cm)","dN^{Events}/dr");
 256     AddHistogram1D<TH1D>("hCentralityV0M", "Centrality distribution V0M", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
 257     AddHistogram1D<TH1D>("hCentralityV0A", "Centrality distribution V0A", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
 258     AddHistogram1D<TH1D>("hCentralityV0C", "Centrality distribution V0C", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
 259     AddHistogram1D<TH1D>("hCentralityZNA", "Centrality distribution ZNA", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
 260     AddHistogram1D<TH1D>("hCentrality", Form("Centrality distribution %s", fCentralityType.Data()), "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
 261
 262
 263     AddHistogram2D<TH2D>("hTrackCountAcc", "Number of tracks in acceptance vs. centrality", "LEGO2", 750, 0., 750., fNumberOfCentralityBins, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
 264     AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
 265     AddHistogram2D<TH2D>("hTrackPtNegEta", "Tracks p_{T} distribution (negative #eta)", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
 266     AddHistogram2D<TH2D>("hTrackPtPosEta", "Tracks p_{T} distribution (positive #eta)", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
 267     AddHistogram1D<TH1D>("hTrackCharge", "Charge", "", 11, -5, 5, "Charge (e)","dN^{Tracks}/dq");
 268     AddHistogram1D<TH1D>("hTrackPhi", "Track #phi distribution", "", 360, 0, TMath::TwoPi(), "#phi","dN^{Tracks}/d#phi");
 269     AddHistogram2D<TH2D>("hTrackPhiEta", "Track angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -2.5, 2.5, "#phi","#eta","dN^{Tracks}/(d#phi d#eta)");
 270
 271     AddHistogram2D<TH2D>("hTrackPhiPtCut", "Track #phi distribution for different pT cuts", "LEGO2", 360, 0, TMath::TwoPi(), 20, 0, 20, "#phi", "p_{T} lower cut", "dN^{Tracks}/d#phi dp_{T}");
 272     AddHistogram2D<TH2D>("hTrackPhiTrackType", "Track #phi distribution for different track types", "LEGO2", 360, 0, TMath::TwoPi(), 3, 0, 3, "#phi", "Label", "dN^{Tracks}/d#phi");
 273     AddHistogram2D<TH2D>("hTrackEta", "Track #eta distribution", "COLZ", 180, -fTrackEtaWindow, +fTrackEtaWindow, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Tracks}/d#eta");
 274     if (fAnalyzeJets)
 275     {
 276       // ######## Jet QA
 277       AddHistogram1D<TH1D>("hRawJetArea", "Jets area distribution w/o area cut", "", 200, 0., 2., "Area","dN^{Jets}/dA");
 278       AddHistogram1D<TH1D>("hJetArea", "Jets area distribution", "", 200, 0., 2., "Area","dN^{Jets}/dA");
 279       AddHistogram2D<TH2D>("hRawJetPhiEta", "Raw Jets angular distribution w/o #eta cut", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dN^{Jets}/(d#phi d#eta)");
 280       AddHistogram2D<TH2D>("hJetEta", "Jets #eta distribution", "COLZ", 180, -fTrackEtaWindow, +fTrackEtaWindow, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta");
 281       AddHistogram2D<TH2D>("hJetEta2GeVTracks", "Jets #eta distribution, track p_{T} > 2 GeV", "COLZ", 180, -fTrackEtaWindow, +fTrackEtaWindow, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta");
 282       AddHistogram2D<TH2D>("hJetEta4GeVTracks", "Jets #eta distribution, track p_{T} > 4 GeV", "COLZ", 180, -fTrackEtaWindow, +fTrackEtaWindow, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta");
 283       AddHistogram2D<TH2D>("hJetPhiEta", "Jets angular distribution", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dN^{Jets}/(d#phi d#eta)");
 284       AddHistogram2D<TH2D>("hJetPtVsConstituentCount", "Jets number of constituents vs. jet p_{T}", "COLZ", 400, 0., 200., 100, 0., 100., "p_{T}","N^{Tracks}","dN^{Jets}/(dp_{T} dN^{tracks})");
 285     }
 286   }
 287
 288   // NOTE: Pythia histograms
 289   if (fAnalyzePythia)
 290   {
 291     AddHistogram1D<TH1D>("hPythiaPtHard", "Pythia p_{T} hard distribution", "", 2000, 0, 400, "p_{T} hard","dN^{Events}/dp_{T,hard}");
 292   }
 293
 294   // register Histograms
 295   for (Int_t i = 0; i < fHistCount; i++)
 296   {
 297     fOutputList->Add(fHistList->At(i));
 298   }
 299
 300   PostData(1,fOutputList); // important for merging
 301
 302 }
 303
 304 //________________________________________________________________________
 305 AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName) : AliAnalysisTaskSE(name), fOutputList(0), fAnalyzeJets(1), fAnalyzeJetProfile(1), fAnalyzeQA(1), fAnalyzeBackground(1), fAnalyzeDeprecatedBackgrounds(1), fAnalyzePythia(0), fAnalyzeMassCorrelation(0), fHasTracks(0), fHasJets(0), fHasBackgroundJets(0), fIsKinematics(0), fUseDefaultVertexCut(1), fUsePileUpCut(1), fSetCentralityToOne(0), fPartialAnalysisNParts(1), fPartialAnalysisIndex(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fNumPtHardBins(11), fUsePtHardBin(-1), fRhoTaskName(), fNcoll(6.88348), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fTRBackgroundConeRadius(0.6), fNumberRandCones(8), fNumberExcludedJets(-1), fDijetMaxAngleDeviation(10.0), fPhysicalJetRadius(0.6), fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetPt(0.15), fMinJetArea(0.5), fMinBackgroundJetPt(0.0), fMinDijetLeadingPt(10.0), fNumberOfCentralityBins(20), 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), fIsDEBUG(0), fEventCounter(0)
 306 {
 307   #ifdef DEBUGMODE
 308     AliInfo("Calling constructor.");
 309   #endif
 310
 311   // Every instance of this task gets his own number
 312   static Int_t instance = 0;
 313   fTaskInstanceCounter = instance;
 314   instance++;
 315
 316   fTrackArrayName = new TString(trackArrayName);
 317   if (fTrackArrayName->Contains("MCParticles") || fTrackArrayName->Contains("mcparticles"))
 318     fIsKinematics = kTRUE;
 319
 320   fJetArrayName = new TString(jetArrayName);
 321   if (strcmp(fJetArrayName->Data(),"") == 0)
 322   {
 323     fAnalyzeJets = kFALSE;
 324     fAnalyzeJetProfile = kFALSE;
 325   }
 326   else
 327     fAnalyzeJets = kTRUE;
 328
 329   fBackgroundJetArrayName = new TString(backgroundJetArrayName);
 330   if (strcmp(fBackgroundJetArrayName->Data(),"") == 0)
 331     fAnalyzeBackground = kFALSE;
 332   else
 333     fAnalyzeBackground = kTRUE;
 334
 335   DefineOutput(1, TList::Class());
 336
 337   fHistList = new TList();
 338
 339   for(Int_t i=0;i<1024;i++)
 340     fSignalJets[i] = NULL;
 341
 342
 343   #ifdef DEBUGMODE
 344     AliInfo("Constructor done.");
 345   #endif
 346
 347 }
 348
 349 //________________________________________________________________________
 350 inline Double_t AliAnalysisTaskChargedJetsPA::GetConePt(Double_t eta, Double_t phi, Double_t radius)
 351 {
 352   Double_t tmpConePt = 0.0;
 353
 354   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
 355   {
 356     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
 357     if (IsTrackInAcceptance(tmpTrack))
 358       if(IsTrackInCone(tmpTrack, eta, phi, radius))
 359         tmpConePt = tmpConePt + tmpTrack->Pt();
 360   }
 361   return tmpConePt;
 362 }
 363
 364 //________________________________________________________________________
 365 inline Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedConePt(Double_t eta, Double_t phi, Double_t radius, Double_t background)
 366 {
 367   Double_t tmpConePt = 0.0;
 368
 369   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
 370   {
 371     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
 372     if (IsTrackInAcceptance(tmpTrack))
 373       if(IsTrackInCone(tmpTrack, eta, phi, radius))
 374         tmpConePt = tmpConePt + tmpTrack->Pt();
 375   }
 376   tmpConePt -= background*radius*radius*TMath::Pi(); // subtract background
 377   return tmpConePt;
 378 }
 379
 380
 381 //________________________________________________________________________
 382 inline Double_t AliAnalysisTaskChargedJetsPA::GetPtHard()
 383 {
 384   #ifdef DEBUGMODE
 385     AliInfo("Starting GetPtHard.");
 386   #endif
 387   AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 388   if (MCEvent())
 389     if (!pythiaHeader)
 390     {
 391       // Check if AOD
 392       AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 393
 394       if (aodMCH)
 395       {
 396         for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++)
 397         {
 398           pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
 399           if (pythiaHeader) break;
 400         }
 401       }
 402     }
 403
 404   #ifdef DEBUGMODE
 405     AliInfo("Ending GetPtHard.");
 406   #endif
 407   if (pythiaHeader)
 408     return pythiaHeader->GetPtHard();
 409
 410   AliWarning(Form("In task %s: GetPtHard() failed!", GetName()));
 411   return -1.0;
 412 }
 413
 414
 415 //________________________________________________________________________
 416 inline Double_t AliAnalysisTaskChargedJetsPA::GetPythiaTrials()
 417 {
 418   #ifdef DEBUGMODE
 419     AliInfo("Starting GetPythiaTrials.");
 420   #endif
 421   AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 422   if (MCEvent())
 423     if (!pythiaHeader)
 424     {
 425       // Check if AOD
 426       AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 427
 428       if (aodMCH)
 429       {
 430         for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++)
 431         {
 432           pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
 433           if (pythiaHeader) break;
 434         }
 435       }
 436     }
 437
 438   #ifdef DEBUGMODE
 439     AliInfo("Ending GetPythiaTrials.");
 440   #endif
 441   if (pythiaHeader)
 442     return pythiaHeader->Trials();
 443
 444   AliWarning(Form("In task %s: GetPythiaTrials() failed!", GetName()));
 445   return -1.0;
 446 }
 447
 448
 449
 450 //________________________________________________________________________
 451 inline Int_t AliAnalysisTaskChargedJetsPA::GetPtHardBin()
 452 {
 453   #ifdef DEBUGMODE
 454     AliInfo("Starting GetPtHardBin.");
 455   #endif
 456   // ########## PT HARD BIN EDGES
 457   const Int_t kPtHardLowerEdges[] =  { 0, 5,11,21,36,57, 84,117,152,191,234};
 458   const Int_t kPtHardHigherEdges[] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 459
 460   Int_t tmpPtHardBin = 0;
 461   Double_t tmpPtHard = GetPtHard();
 462
 463   for (tmpPtHardBin = 0; tmpPtHardBin <= fNumPtHardBins; tmpPtHardBin++)
 464     if (tmpPtHard >= kPtHardLowerEdges[tmpPtHardBin] && tmpPtHard < kPtHardHigherEdges[tmpPtHardBin])
 465       break;
 466
 467   #ifdef DEBUGMODE
 468     AliInfo("Ending GetPtHardBin.");
 469   #endif
 470   return tmpPtHardBin;
 471 }
 472
 473 //________________________________________________________________________
 474 Double_t AliAnalysisTaskChargedJetsPA::GetExternalRho()
 475 {
 476   // Get rho from event.
 477   AliRhoParameter *rho = 0;
 478   if (!fRhoTaskName.IsNull()) {
 479     rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoTaskName.Data()));
 480     if (!rho) {
 481       AliWarning(Form("%s: Could not retrieve rho with name %s!", GetName(), fRhoTaskName.Data()));
 482       return 0;
 483     }
 484   }
 485   else
 486     return 0;
 487
 488   return (rho->GetVal());
 489 }
 490
 491
 492 //________________________________________________________________________
 493 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius)
 494 {
 495   // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
 496   Double_t trackPhi = 0.0;
 497   if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
 498     trackPhi = track->Phi() - TMath::TwoPi();
 499   else if (track->Phi() < (phi+radius - TMath::TwoPi()))
 500     trackPhi = track->Phi() + TMath::TwoPi();
 501   else
 502     trackPhi = track->Phi();
 503
 504   if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
 505     return kTRUE;
 506
 507   return kFALSE;
 508 }
 509
 510 //________________________________________________________________________
 511 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInJet(AliEmcalJet* jet, Int_t trackIndex)
 512 {
 513   for (Int_t i = 0; i < jet->GetNumberOfTracks(); ++i)
 514   {
 515     Int_t jetTrack = jet->TrackAt(i);
 516     if (jetTrack == trackIndex)
 517       return kTRUE;
 518   }
 519   return kFALSE;
 520 }
 521
 522 //________________________________________________________________________
 523 inline Bool_t AliAnalysisTaskChargedJetsPA::IsJetOverlapping(AliEmcalJet* jet1, AliEmcalJet* jet2)
 524 {
 525   for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
 526   {
 527     Int_t jet1Track = jet1->TrackAt(i);
 528     for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
 529     {
 530       Int_t jet2Track = jet2->TrackAt(j);
 531       if (jet1Track == jet2Track)
 532         return kTRUE;
 533     }
 534   }
 535   return kFALSE;
 536 }
 537
 538 //________________________________________________________________________
 539 inline Bool_t AliAnalysisTaskChargedJetsPA::IsEventInAcceptance(AliVEvent* event)
 540 {
 541   if (!event)
 542     return kFALSE;
 543
 544   FillHistogram("hEventAcceptance", 0.5); // number of events before manual cuts
 545   if(fUsePileUpCut)
 546     if(fHelperClass->IsPileUpEvent(event))
 547       return kFALSE;
 548
 549   FillHistogram("hEventAcceptance", 1.5); // number of events after pileup cuts
 550
 551   if(fAnalyzeQA)
 552     FillHistogram("hVertexZBeforeVertexCut",event->GetPrimaryVertex()->GetZ());
 553
 554   if(fUseDefaultVertexCut)
 555   {
 556     if(!fHelperClass->IsVertexSelected2013pA(event))
 557       return kFALSE;
 558   }
 559   else // Failsafe vertex cut
 560   {
 561     if(TMath::Abs(event->GetPrimaryVertex()->GetZ()) > 10.0)
 562       return kFALSE;
 563   }
 564
 565   if(fAnalyzeQA)
 566     FillHistogram("hVertexZAfterVertexCut",event->GetPrimaryVertex()->GetZ());
 567
 568   FillHistogram("hEventAcceptance", 2.5); // number of events after vertex cut
 569
 570   return kTRUE;
 571 }
 572
 573 //________________________________________________________________________
 574 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInAcceptance(AliVParticle* track)
 575 {
 576   FillHistogram("hTrackAcceptance", 0.5);
 577   if (track != 0)
 578   {
 579     if(fIsKinematics)
 580     {
 581       // TODO: Only working for AOD MC
 582       if((!track->Charge()) || (!(static_cast<AliAODMCParticle*>(track))->IsPhysicalPrimary()) )
 583         return kFALSE;
 584     }
 585     if (TMath::Abs(track->Eta()) <= fTrackEtaWindow)
 586     {
 587       FillHistogram("hTrackAcceptance", 1.5);
 588       if (track->Pt() >= fMinTrackPt)
 589       {
 590         FillHistogram("hTrackAcceptance", 2.5);
 591         return kTRUE;
 592       }
 593     }
 594   }
 595   return kFALSE;
 596 }
 597
 598 //________________________________________________________________________
 599 inline Bool_t AliAnalysisTaskChargedJetsPA::IsBackgroundJetInAcceptance(AliEmcalJet *jet)
 600 {
 601   if (jet != 0)
 602     if (TMath::Abs(jet->Eta()) <= fBackgroundJetEtaWindow)
 603       if (jet->Pt() >= fMinBackgroundJetPt)
 604         return kTRUE;
 605
 606   return kFALSE;
 607 }
 608
 609 //________________________________________________________________________
 610 inline Bool_t AliAnalysisTaskChargedJetsPA::IsSignalJetInAcceptance(AliEmcalJet *jet)
 611 {
 612   FillHistogram("hJetAcceptance", 0.5);
 613   if (jet != 0)
 614     if (TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow)
 615     {
 616       FillHistogram("hJetAcceptance", 1.5);
 617       if (jet->Pt() >= fMinJetPt)
 618       {
 619         FillHistogram("hJetAcceptance", 2.5);
 620         if (jet->Area() >= fMinJetArea)
 621         {
 622           FillHistogram("hJetAcceptance", 3.5);
 623           return kTRUE;
 624         }
 625       }
 626     }
 627   return kFALSE;
 628 }
 629
 630 //________________________________________________________________________
 631 inline Bool_t AliAnalysisTaskChargedJetsPA::IsDijet(AliEmcalJet *jet1, AliEmcalJet *jet2)
 632 {
 633   // Output from GetDeltaPhi is < pi in any case
 634   if ((jet1 != 0) && (jet2 != 0))
 635     if((TMath::Pi() - GetDeltaPhi(jet1->Phi(),jet2->Phi())) < fDijetMaxAngleDeviation)
 636       if ((jet1->Pt() > fMinDijetLeadingPt) && (jet2->Pt() > fMinDijetLeadingPt)) //TODO: Introduce recoil jet?
 637         return kTRUE;
 638
 639   return kFALSE;
 640 }
 641
 642 //________________________________________________________________________
 643 void AliAnalysisTaskChargedJetsPA::ExecOnce()
 644 {
 645   #ifdef DEBUGMODE
 646     AliInfo("Starting ExecOnce.");
 647   #endif
 648   fInitialized = kTRUE;
 649
 650   // Check for track array
 651   if (strcmp(fTrackArrayName->Data(), "") != 0)
 652   {
 653     fTrackArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrayName->Data()));
 654     fHasTracks = kTRUE;
 655     if (!fTrackArray)
 656     {
 657       AliWarning(Form("%s: Could not retrieve tracks %s! This is OK, if tracks are not demanded.", GetName(), fTrackArrayName->Data()));
 658       fHasTracks = kFALSE;
 659     }
 660     else
 661     {
 662       TClass *cl = fTrackArray->GetClass();
 663       if (!cl->GetBaseClass("AliVParticle"))
 664       {
 665         AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName->Data()));
 666         fTrackArray = 0;
 667         fHasTracks = kFALSE;
 668       }
 669     }
 670   }
 671
 672   // Check for jet array
 673   if (strcmp(fJetArrayName->Data(), "") != 0)
 674   {
 675     fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrayName->Data()));
 676     fHasJets = kTRUE;
 677
 678     if (!fJetArray)
 679     {
 680       AliWarning(Form("%s: Could not retrieve jets %s! This is OK, if jets are not demanded.", GetName(), fJetArrayName->Data()));
 681       fHasJets = kFALSE;
 682     }
 683     else
 684     {
 685       if (!fJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
 686       {
 687         AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName->Data()));
 688         fJetArray = 0;
 689         fHasJets = kFALSE;
 690       }
 691     }
 692   }
 693
 694   // Check for background object
 695   if (strcmp(fBackgroundJetArrayName->Data(), "") != 0)
 696   {
 697     fBackgroundJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fBackgroundJetArrayName->Data()));
 698     fHasBackgroundJets = kTRUE;
 699     if (!fBackgroundJetArray)
 700     {
 701       AliInfo(Form("%s: Could not retrieve background jets %s! This is OK, if background is not demanded.", GetName(), fBackgroundJetArrayName->Data()));
 702       fHasBackgroundJets = kFALSE;
 703     }
 704   }
 705
 706   // Look, if initialization is OK
 707   if (!fHasTracks && fAnalyzeBackground)
 708   {
 709     AliError(Form("%s: Tracks NOT successfully casted although demanded! Deactivating background analysis",GetName()));
 710     fAnalyzeBackground = kFALSE;
 711   }
 712   if ((!fHasJets && fAnalyzeJets) || (!fHasJets && fAnalyzeBackground))
 713   {
 714     AliError(Form("%s: Jets NOT successfully casted although demanded!  Deactivating jet- and background analysis",GetName()));
 715     fAnalyzeJets = kFALSE;
 716     fAnalyzeBackground = kFALSE;
 717   }
 718   if (!fHasBackgroundJets && fAnalyzeBackground)
 719   {
 720     AliError(Form("%s: Background NOT successfully casted although demanded!  Deactivating background analysis",GetName()));
 721     fAnalyzeBackground = kFALSE;
 722   }
 723
 724   // Initialize helper class (for vertex selection & pile up correction)
 725   fHelperClass = new AliAnalysisUtils();
 726   fHelperClass->SetCutOnZVertexSPD(kFALSE);
 727   // Histogram init
 728   Init();
 729
 730   #ifdef DEBUGMODE
 731     AliInfo("ExecOnce done.");
 732   #endif
 733
 734 }
 735
 736 //________________________________________________________________________
 737 void AliAnalysisTaskChargedJetsPA::GetSignalJets()
 738 {
 739   // Reset vars
 740   fFirstLeadingJet = NULL;
 741   fSecondLeadingJet = NULL;
 742   fNumberSignalJets = 0;
 743
 744   TList tmpJets;
 745   for (Int_t i = 0; i < fJetArray->GetEntries(); i++)
 746   {
 747     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(i));
 748     if (!jet)
 749     {
 750       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 751       continue;
 752     }
 753     if (!IsSignalJetInAcceptance(jet))
 754       continue;
 755
 756     for (Int_t j = 0; j <= tmpJets.GetEntries(); j++)
 757     {
 758       if (j>tmpJets.GetEntries()-1) // When passed last item add the jet at the end
 759       {
 760         tmpJets.Add(jet);
 761         break;
 762       }
 763
 764       AliEmcalJet* listJet = static_cast<AliEmcalJet*>(tmpJets.At(j));
 765
 766       if(jet->Pt() < listJet->Pt()) // Insert jet before that one in list if pt smaller
 767       {
 768         tmpJets.AddAt(jet, j);
 769         break;
 770       }
 771     }
 772   }
 773
 774   for (Int_t i = 0; i < tmpJets.GetEntries(); i++)
 775   {
 776     AliEmcalJet* jet = static_cast<AliEmcalJet*>(tmpJets.At(i));
 777     fSignalJets[fNumberSignalJets] = jet;
 778     fNumberSignalJets++;
 779   }
 780
 781   if (fNumberSignalJets > 0)
 782     fFirstLeadingJet  = static_cast<AliEmcalJet*>(tmpJets.At(0));
 783   if (fNumberSignalJets > 1)
 784     fSecondLeadingJet = static_cast<AliEmcalJet*>(tmpJets.At(1));
 785
 786 }
 787
 788 //________________________________________________________________________
 789 Int_t AliAnalysisTaskChargedJetsPA::GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray, Bool_t isSignalJets)
 790 {
 791 // Writes first two leading jets into already registered array jetIDArray
 792
 793   if (!jetArray)
 794   {
 795     AliError("Could not get the jet array to get leading jets from it!");
 796     return 0;
 797   }
 798
 799   Float_t maxJetPts[] = {0, 0};
 800   jetIDArray[0] = -1;
 801   jetIDArray[1] = -1;
 802
 803   Int_t jetCount = jetArray->GetEntries();
 804   Int_t jetCountAccepted = 0;
 805
 806   for (Int_t i = 0; i < jetCount; i++)
 807   {
 808     AliEmcalJet* jet = static_cast<AliEmcalJet*>(jetArray->At(i));
 809     if (!jet)
 810     {
 811       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 812       continue;
 813     }
 814
 815     if(isSignalJets)
 816     {
 817       if (!IsSignalJetInAcceptance(jet)) continue;
 818     }
 819     else
 820     {
 821       if (!IsBackgroundJetInAcceptance(jet)) continue;
 822     }
 823
 824     if (jet->Pt() > maxJetPts[0])
 825     {
 826       maxJetPts[1] = maxJetPts[0];
 827       jetIDArray[1] = jetIDArray[0];
 828       maxJetPts[0] = jet->Pt();
 829       jetIDArray[0] = i;
 830     }
 831     else if (jet->Pt() > maxJetPts[1])
 832     {
 833       maxJetPts[1] = jet->Pt();
 834       jetIDArray[1] = i;
 835     }
 836     jetCountAccepted++;
 837   }
 838   return jetCountAccepted;
 839 }
 840
 841
 842 //________________________________________________________________________
 843 Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Double_t background)
 844 {
 845   #ifdef DEBUGMODE
 846     AliInfo("Getting corrected jet spectra.");
 847   #endif
 848
 849   if(!jet)
 850   {
 851     AliError("Jet pointer passed to GetCorrectedJetPt() not valid!");
 852     return -1.0;
 853   }
 854
 855   Double_t correctedPt = -1.0;
 856
 857   // if the passed background is not valid, do not subtract it
 858   if(background < 0)
 859     background = 0;
 860
 861   // Subtract background
 862   correctedPt = jet->Pt() - background * jet->Area();
 863
 864   #ifdef DEBUGMODE
 865     AliInfo("Got corrected jet spectra.");
 866   #endif
 867
 868   return correctedPt;
 869 }
 870
 871
 872
 873 //________________________________________________________________________
 874 Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t leadingJetExclusionProbability)
 875 {
 876   #ifdef DEBUGMODE
 877     AliInfo("Getting Delta Pt.");
 878   #endif
 879
 880   // Define an invalid delta pt
 881   Double_t deltaPt = -10000.0;
 882
 883   // Define eta range
 884   Double_t etaMin, etaMax;
 885   etaMin = -(fTrackEtaWindow-fRandConeRadius);
 886   etaMax = +(fTrackEtaWindow-fRandConeRadius);
 887
 888   // Define random cone
 889   Bool_t coneValid = kTRUE;
 890   Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
 891   Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
 892
 893   // if there is a jet, check for overlap if demanded
 894   if(leadingJetExclusionProbability)
 895   {
 896     AliEmcalJet* tmpLeading = dynamic_cast<AliEmcalJet*>(fJetArray->At(0));
 897     // Get leading jet (regardless of pT)
 898     for (Int_t i = 1; i<fJetArray->GetEntries(); i++)
 899     {
 900       AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
 901       // if jet is in acceptance and higher, take as new leading
 902       if (tmpJet)
 903         if ((TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow) && (tmpJet->Area() >= fMinJetArea))
 904           if((!tmpLeading) || (tmpJet->Pt() > tmpLeading->Pt()))
 905             tmpLeading = tmpJet;
 906     }
 907     if(tmpLeading)
 908     {
 909       Double_t excludedJetPhi = tmpLeading->Phi();
 910       Double_t excludedJetEta = tmpLeading->Eta();
 911       Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, excludedJetPhi);
 912
 913       // Check, if cone has overlap with jet
 914       if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(tmpRandConeEta-excludedJetEta)*TMath::Abs(tmpRandConeEta-excludedJetEta) <= fRandConeRadius*fRandConeRadius)
 915       {
 916         // Define probability to exclude the RC
 917         Double_t probability = leadingJetExclusionProbability;
 918
 919         // Only exclude cone with a given probability
 920         if (fRandom->Rndm()<=probability)
 921           coneValid = kFALSE;
 922       }
 923     }
 924   }
 925
 926
 927   // Get the cones' pt and calculate delta pt
 928   if (coneValid)
 929     deltaPt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius) - (rho*fRandConeRadius*fRandConeRadius*TMath::Pi());
 930
 931   return deltaPt;
 932   #ifdef DEBUGMODE
 933     AliInfo("Got Delta Pt.");
 934   #endif
 935 }
 936
 937 //________________________________________________________________________
 938 void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensityAll(Int_t numberExcludeLeadingJets, Double_t& rhoPbPb, Double_t& rhoPbPbWithGhosts, Double_t& rhoCMS, Double_t& rhoImprovedCMS, Double_t& rhoMean, Double_t& rhoTrackLike)
 939 {
 940   #ifdef DEBUGMODE
 941     AliInfo("Getting ALL KT background density.");
 942   #endif
 943
 944   static Double_t tmpRhoPbPb[1024];
 945   static Double_t tmpRhoPbPbWithGhosts[1024];
 946   static Double_t tmpRhoMean[1024];
 947   static Double_t tmpRhoCMS[1024];
 948   static Double_t tmpRhoImprovedCMS[1024];
 949   Double_t tmpCoveredArea = 0.0;
 950   Double_t tmpSummedArea = 0.0;
 951   Double_t tmpPtTrackLike = 0.0;
 952   Double_t tmpAreaTrackLike = 0.0;
 953
 954   // Setting invalid values
 955   rhoPbPb = 0.0;
 956   rhoPbPbWithGhosts = 0.0;
 957   rhoCMS = 0.0;
 958   rhoImprovedCMS = 0.0;
 959   rhoMean = 0.0;
 960   rhoTrackLike = 0.0;
 961
 962   Int_t rhoPbPbJetCount = 0;
 963   Int_t rhoPbPbWithGhostsJetCount = 0;
 964   Int_t rhoCMSJetCount = 0;
 965   Int_t rhoImprovedCMSJetCount = 0;
 966   Int_t rhoMeanJetCount = 0;
 967
 968
 969   // Find 2 leading KT jets for the original PbPb approach
 970   Int_t leadingKTJets[]   = {-1, -1};
 971   GetLeadingJets(fBackgroundJetArray, &leadingKTJets[0], kFALSE);
 972
 973   // Exclude UP TO numberExcludeLeadingJets
 974   if(numberExcludeLeadingJets==-1)
 975     numberExcludeLeadingJets = fNumberSignalJets;
 976   if (fNumberSignalJets < numberExcludeLeadingJets)
 977     numberExcludeLeadingJets = fNumberSignalJets;
 978
 979   for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++)
 980   {
 981     AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
 982
 983     if (!backgroundJet)
 984     {
 985       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 986       continue;
 987     }
 988
 989     tmpSummedArea += backgroundJet->Area();
 990     if(backgroundJet->Pt() > 0.150)
 991       tmpCoveredArea += backgroundJet->Area();
 992
 993     if (!IsBackgroundJetInAcceptance(backgroundJet))
 994       continue;
 995
 996     // Search for overlap with signal jets
 997     Bool_t isOverlapping = kFALSE;
 998     for(Int_t j=0;j<numberExcludeLeadingJets;j++)
 999     {
1000       AliEmcalJet* signalJet = fSignalJets[j];
1001
1002       if(signalJet->Pt() >= 5.0)
1003         if(IsJetOverlapping(signalJet, backgroundJet))
1004         {
1005           isOverlapping = kTRUE;
1006           break;
1007         }
1008     }
1009
1010     Double_t tmpRho = 0.0;
1011     if(backgroundJet->Area())
1012       tmpRho = backgroundJet->Pt() / backgroundJet->Area();
1013
1014     // PbPb approach (take ghosts into account)
1015     if ((i != leadingKTJets[0]) && (i != leadingKTJets[1]))
1016     {
1017       tmpRhoPbPbWithGhosts[rhoPbPbWithGhostsJetCount] = tmpRho;
1018       rhoPbPbWithGhostsJetCount++;
1019     }
1020
1021     if(backgroundJet->Pt() > 0.150)
1022     {
1023       // CMS approach: don't take ghosts into acount
1024       tmpRhoCMS[rhoCMSJetCount] = tmpRho;
1025       rhoCMSJetCount++;
1026
1027       // Improved CMS approach: like CMS but excluding signal
1028       if(!isOverlapping)
1029       {
1030         tmpRhoImprovedCMS[rhoImprovedCMSJetCount] = tmpRho;
1031         rhoImprovedCMSJetCount++;
1032       }
1033
1034       // PbPb w/o ghosts approach (just neglect ghosts)
1035       if ((i != leadingKTJets[0]) && (i != leadingKTJets[1]))
1036       {
1037         tmpRhoPbPb[rhoPbPbJetCount] = tmpRho;
1038         rhoPbPbJetCount++;
1039       }
1040     }
1041
1042     // (no overlap with signal jets)
1043     if(!isOverlapping)
1044     {
1045       // Mean approach
1046       tmpRhoMean[rhoMeanJetCount] = tmpRho;
1047       rhoMeanJetCount++;
1048
1049       // Track like approach approach
1050       tmpPtTrackLike += backgroundJet->Pt();
1051       tmpAreaTrackLike += backgroundJet->Area();
1052     }
1053
1054   }
1055
1056   if (tmpAreaTrackLike > 0)
1057     rhoTrackLike = tmpPtTrackLike/tmpAreaTrackLike;
1058   if (rhoPbPbJetCount > 0)
1059     rhoPbPb = TMath::Median(rhoPbPbJetCount, tmpRhoPbPb);
1060   if (rhoPbPbWithGhostsJetCount > 0)
1061     rhoPbPbWithGhosts = TMath::Median(rhoPbPbWithGhostsJetCount, tmpRhoPbPbWithGhosts);
1062   if (rhoCMSJetCount > 0)
1063     rhoCMS = TMath::Median(rhoCMSJetCount, tmpRhoCMS) * tmpCoveredArea/tmpSummedArea;
1064   if (rhoImprovedCMSJetCount > 0)
1065   {
1066     rhoImprovedCMS = TMath::Median(rhoImprovedCMSJetCount, tmpRhoImprovedCMS) * tmpCoveredArea/tmpSummedArea;
1067   }
1068   if (rhoMeanJetCount > 0)
1069     rhoMean = TMath::Mean(rhoMeanJetCount, tmpRhoMean);
1070
1071   #ifdef DEBUGMODE
1072     AliInfo("Got ALL KT background density.");
1073   #endif
1074 }
1075
1076 //________________________________________________________________________
1077 void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoImprovedCMS)
1078 {
1079   #ifdef DEBUGMODE
1080     AliInfo("Getting KT background density.");
1081   #endif
1082
1083   static Double_t tmpRhoImprovedCMS[1024];
1084   Double_t tmpCoveredArea = 0.0;
1085   Double_t tmpSummedArea = 0.0;
1086
1087   // Setting invalid values
1088   rhoImprovedCMS = 0.0;
1089
1090   Int_t rhoImprovedCMSJetCount = 0;
1091
1092   // Exclude UP TO numberExcludeLeadingJets
1093   if(numberExcludeLeadingJets==-1)
1094     numberExcludeLeadingJets = fNumberSignalJets;
1095   if (fNumberSignalJets < numberExcludeLeadingJets)
1096     numberExcludeLeadingJets = fNumberSignalJets;
1097
1098   for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++)
1099   {
1100     AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
1101
1102     if (!backgroundJet)
1103     {
1104       AliError(Form("%s: Could not receive jet %d", GetName(), i));
1105       continue;
1106     }
1107
1108     // Search for overlap with signal jets
1109     Bool_t isOverlapping = kFALSE;
1110     for(Int_t j=0;j<numberExcludeLeadingJets;j++)
1111     {
1112       AliEmcalJet* signalJet = fSignalJets[j];
1113       if(signalJet->Pt() >= 5.0)
1114         if(IsJetOverlapping(signalJet, backgroundJet))
1115         {
1116           isOverlapping = kTRUE;
1117           break;
1118         }
1119     }
1120
1121     tmpSummedArea += backgroundJet->Area();
1122     if(backgroundJet->Pt() > 0.150)
1123       tmpCoveredArea += backgroundJet->Area();
1124
1125     if (!IsBackgroundJetInAcceptance(backgroundJet))
1126       continue;
1127
1128     Double_t tmpRho = backgroundJet->Pt() / backgroundJet->Area();
1129
1130     if(backgroundJet->Pt() > 0.150)
1131       if(!isOverlapping)
1132       {
1133         tmpRhoImprovedCMS[rhoImprovedCMSJetCount] = tmpRho;
1134         rhoImprovedCMSJetCount++;
1135       }
1136   }
1137
1138   if (rhoImprovedCMSJetCount > 0)
1139   {
1140     rhoImprovedCMS = TMath::Median(rhoImprovedCMSJetCount, tmpRhoImprovedCMS) * tmpCoveredArea/tmpSummedArea;
1141   }
1142   #ifdef DEBUGMODE
1143     AliInfo("Got KT background density.");
1144   #endif
1145 }
1146
1147
1148 //________________________________________________________________________
1149 void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoNoExclusion, Double_t& rhoConeExclusion02, Double_t& rhoConeExclusion04, Double_t& rhoConeExclusion06, Double_t& rhoConeExclusion08, Double_t& rhoExactExclusion)
1150 {
1151   #ifdef DEBUGMODE
1152     AliInfo("Getting TR background density.");
1153   #endif
1154
1155   Double_t summedTracksPtCone04 = 0.0;
1156   Double_t summedTracksPtCone02 = 0.0;
1157   Double_t summedTracksPtCone06 = 0.0;
1158   Double_t summedTracksPtCone08 = 0.0;
1159   Double_t summedTracksPtWithinJets = 0.0;
1160   Double_t summedTracksPt = 0.0;
1161
1162   // Setting invalid values
1163   rhoNoExclusion = 0.0;
1164   rhoConeExclusion02 = 0.0;
1165   rhoConeExclusion04 = 0.0;
1166   rhoConeExclusion06 = 0.0;
1167   rhoConeExclusion08 = 0.0;
1168   rhoExactExclusion  = 0.0;
1169
1170   // Exclude UP TO numberExcludeLeadingJets
1171   if(numberExcludeLeadingJets==-1)
1172     numberExcludeLeadingJets = fNumberSignalJets;
1173   if (fNumberSignalJets < numberExcludeLeadingJets)
1174     numberExcludeLeadingJets = fNumberSignalJets;
1175
1176   Int_t fSignalJetCount5GeV = 0;
1177   for(Int_t j=0;j<numberExcludeLeadingJets;j++)
1178   {
1179     AliEmcalJet* signalJet = fSignalJets[j];
1180     if(signalJet->Pt() < 5.0)
1181       continue;
1182     fSignalJetCount5GeV++;
1183   }
1184
1185   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
1186   {
1187     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
1188     Bool_t trackWithinJet = kFALSE; Bool_t trackWithin02Cone = kFALSE; Bool_t trackWithin04Cone = kFALSE; Bool_t trackWithin06Cone = kFALSE; Bool_t trackWithin08Cone = kFALSE;
1189
1190     if (IsTrackInAcceptance(tmpTrack))
1191     {
1192       // Check if tracks overlaps with jet
1193       for(Int_t j=0;j<numberExcludeLeadingJets;j++)
1194       {
1195         AliEmcalJet* signalJet = fSignalJets[j];
1196
1197         if(signalJet->Pt() < 5.0)
1198           continue;
1199
1200         // Exact jet exclusion
1201         if (IsTrackInJet(signalJet, i))
1202           trackWithinJet = kTRUE;
1203
1204         // Cone exclusions
1205         if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.2))
1206         {
1207           trackWithin02Cone = kTRUE;
1208           trackWithin04Cone = kTRUE;
1209           trackWithin06Cone = kTRUE;
1210           trackWithin08Cone = kTRUE;
1211           break;
1212         }
1213         else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.4))
1214         {
1215           trackWithin04Cone = kTRUE;
1216           trackWithin06Cone = kTRUE;
1217           trackWithin08Cone = kTRUE;
1218         }
1219         else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.6))
1220         {
1221           trackWithin06Cone = kTRUE;
1222           trackWithin08Cone = kTRUE;
1223         }
1224         else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.8))
1225         {
1226           trackWithin08Cone = kTRUE;
1227         }
1228       }
1229
1230       if(!trackWithin08Cone)
1231       {
1232         summedTracksPtCone08 += tmpTrack->Pt();
1233       }
1234       if(!trackWithin06Cone)
1235       {
1236         summedTracksPtCone06 += tmpTrack->Pt();
1237       }
1238       if(!trackWithin04Cone)
1239       {
1240         summedTracksPtCone04 += tmpTrack->Pt();
1241       }
1242       if(!trackWithin02Cone)
1243       {
1244         summedTracksPtCone02 += tmpTrack->Pt();
1245       }
1246       if(!trackWithinJet)
1247       {
1248         summedTracksPtWithinJets += tmpTrack->Pt();
1249       }
1250       summedTracksPt += tmpTrack->Pt();
1251
1252     }
1253   }
1254
1255   // Calculate the correct area where the tracks were taking from
1256
1257   Double_t tmpFullTPCArea = (2.0*fTrackEtaWindow) * TMath::TwoPi();
1258   Double_t tmpAreaCone02     = tmpFullTPCArea;
1259   Double_t tmpAreaCone04     = tmpFullTPCArea;
1260   Double_t tmpAreaCone06     = tmpFullTPCArea;
1261   Double_t tmpAreaCone08     = tmpFullTPCArea;
1262   Double_t tmpAreaWithinJets = tmpFullTPCArea;
1263   std::vector<Double_t> tmpEtas(fSignalJetCount5GeV);
1264   std::vector<Double_t> tmpPhis(fSignalJetCount5GeV);
1265
1266   Int_t iSignal = 0;
1267   for(Int_t i=0;i<numberExcludeLeadingJets;i++)
1268   {
1269     AliEmcalJet* tmpJet = fSignalJets[i];
1270
1271     if(tmpJet->Pt() < 5.0)
1272       continue;
1273
1274     tmpEtas[iSignal] = tmpJet->Eta();
1275     tmpPhis[iSignal] = tmpJet->Phi();
1276     tmpAreaWithinJets -= tmpJet->Area();
1277
1278     iSignal++;
1279   }
1280
1281   tmpAreaCone02 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fSignalJetCount5GeV, tmpEtas, tmpPhis, 0.2, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi());
1282   tmpAreaCone04 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fSignalJetCount5GeV, tmpEtas, tmpPhis, 0.4, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi());
1283   tmpAreaCone06 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fSignalJetCount5GeV, tmpEtas, tmpPhis, 0.6, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi());
1284   tmpAreaCone08 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fSignalJetCount5GeV, tmpEtas, tmpPhis, 0.8, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi());
1285
1286   rhoConeExclusion02 = summedTracksPtCone02/tmpAreaCone02;
1287   rhoConeExclusion04 = summedTracksPtCone04/tmpAreaCone04;
1288   rhoConeExclusion06 = summedTracksPtCone06/tmpAreaCone06;
1289   rhoConeExclusion08 = summedTracksPtCone08/tmpAreaCone08;
1290   rhoExactExclusion  = summedTracksPtWithinJets/tmpAreaWithinJets;
1291   rhoNoExclusion     = summedTracksPt/tmpFullTPCArea;
1292
1293
1294   #ifdef DEBUGMODE
1295     AliInfo("Got TR background density.");
1296   #endif
1297 }
1298
1299 //________________________________________________________________________
1300 void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular)
1301 {
1302   #ifdef DEBUGMODE
1303     AliInfo("Getting TR background density.");
1304   #endif
1305
1306   // Setting invalid values
1307   Double_t summedTracksPt = 0.0;
1308   rhoMean = 0.0;
1309   area = -1.0;
1310
1311   Double_t tmpRadius = 0.0;
1312   if (doSearchPerpendicular)
1313     tmpRadius = 0.4*TMath::Pi(); // exclude 90 degrees around jets
1314   else
1315     tmpRadius = 0.8;
1316
1317   numberExcludeLeadingJets = 2; // dijet is excluded here in any case
1318
1319
1320
1321   if (!fTrackArray || !fJetArray)
1322   {
1323     AliError("Could not get the track/jet array to calculate track rho!");
1324     return;
1325   }
1326
1327   Int_t   trackCount = fTrackArray->GetEntries();
1328   Int_t   trackCountAccepted = 0;
1329   for (Int_t i = 0; i < trackCount; i++)
1330   {
1331     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
1332     if (IsTrackInAcceptance(tmpTrack))
1333     {
1334       if (IsTrackInCone(tmpTrack, excludeJet1->Eta(), excludeJet1->Phi(), tmpRadius))
1335         continue;
1336
1337       if (numberExcludeLeadingJets > 1)
1338         if (IsTrackInCone(tmpTrack, excludeJet2->Eta(), excludeJet2->Phi(), tmpRadius))
1339           continue;
1340
1341         // Add track pt to array
1342         summedTracksPt = summedTracksPt + tmpTrack->Pt();
1343         trackCountAccepted++;
1344     }
1345   }
1346
1347   if (trackCountAccepted > 0)
1348   {
1349     Double_t tmpArea = 2.0*fTrackEtaWindow*TMath::TwoPi()  - 2*(tmpRadius*tmpRadius * TMath::Pi()); //TPC area - excluding jet area
1350     rhoMean = summedTracksPt/tmpArea;
1351     area = tmpArea;
1352   }
1353
1354   #ifdef DEBUGMODE
1355     AliInfo("Got TR background density.");
1356   #endif
1357 }
1358
1359 //________________________________________________________________________
1360 void AliAnalysisTaskChargedJetsPA::GetPPBackgroundDensity(Double_t& background, AliEmcalJet* jet)
1361 {
1362   // This is the background that was used for the pp 7 TeV ALICE paper
1363   Double_t jetMom[3] = { jet->Px(), jet->Py(), jet->Pz() };
1364   TVector3 jet3mom1(jetMom);
1365   TVector3 jet3mom2(jetMom);
1366   background = 0;
1367
1368   jet3mom1.RotateZ(TMath::Pi());
1369   jet3mom2.RotateZ(-TMath::Pi());
1370
1371   for (int i = 0; i < fTrackArray->GetEntries(); i++)
1372   {
1373     AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i));
1374
1375     Double_t trackMom[3] = { track->Px(), track->Py(), track->Pz() };
1376     TVector3 track3mom(trackMom);
1377
1378     Double_t dR1 = jet3mom1.DeltaR(track3mom);
1379     Double_t dR2 = jet3mom2.DeltaR(track3mom);
1380
1381     if (dR1 <= fSignalJetRadius || dR2 <= fSignalJetRadius)
1382       background += track3mom.Pt();
1383   }
1384
1385   background /= (2 * TMath::Pi() * fSignalJetRadius * fSignalJetRadius);
1386 }
1387
1388 //________________________________________________________________________
1389 void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
1390 {
1391   #ifdef DEBUGMODE
1392     AliInfo("Starting Calculate().");
1393   #endif
1394   ////////////////////// NOTE: initialization & casting
1395
1396   fEventCounter++;
1397
1398   // Check, if analysis should be done in pt hard bins
1399   if(fUsePtHardBin != -1)
1400     if(GetPtHardBin() != fUsePtHardBin)
1401       return;
1402
1403   // This is to take only every Nth event
1404   if((fEventCounter+fPartialAnalysisIndex) % fPartialAnalysisNParts != 0)
1405     return;
1406
1407   FillHistogram("hNumberEvents",0.5);
1408
1409   if(!IsEventInAcceptance(event))
1410     return;
1411
1412   FillHistogram("hNumberEvents",1.5);
1413
1414   #ifdef DEBUGMODE
1415     AliInfo("Calculate()::Init done.");
1416   #endif
1417
1418   ////////////////////// NOTE: Get Centrality, (Leading)Signal jets and Background
1419
1420   // Get centrality
1421   AliCentrality* tmpCentrality = NULL;
1422   tmpCentrality = event->GetCentrality();
1423   Double_t centralityPercentile = -1.0;
1424   Double_t centralityPercentileV0A = 0.0;
1425   Double_t centralityPercentileV0C = 0.0;
1426   Double_t centralityPercentileV0M = 0.0;
1427   Double_t centralityPercentileZNA = 0.0;
1428   if (tmpCentrality != NULL)
1429   {
1430     centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
1431     centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A");
1432     centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C");
1433     centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M");
1434     centralityPercentileZNA = tmpCentrality->GetCentralityPercentile("ZNA");
1435   }
1436
1437   if((centralityPercentile < 0.0) || (centralityPercentile > 100.0))
1438     AliWarning(Form("Centrality value not valid (c=%E)",centralityPercentile));
1439
1440   if(fSetCentralityToOne)
1441     centralityPercentile = 1.0;
1442
1443   // Get jets
1444   if (fAnalyzeBackground || fAnalyzeJets)
1445     GetSignalJets();
1446
1447   // Get background estimates
1448   Double_t              backgroundKTImprovedCMS = -1.0;
1449   Double_t              backgroundKTImprovedCMSExternal = -1.0;
1450   Double_t              backgroundDijet = -1.0;
1451   Double_t              backgroundDijetPerpendicular = -1.0;
1452   Double_t              deltaPtDijetPerpendicular = -10000.0;
1453   Double_t              backgroundKTPbPb = -1.0;
1454   Double_t              backgroundKTPbPbWithGhosts = -1.0;
1455   Double_t              backgroundKTCMS = -1.0;
1456   Double_t              backgroundKTMean = -1.0;
1457   Double_t              backgroundKTTrackLike = -1.0;
1458   Double_t              backgroundTRNoExcl = -1.0;
1459   Double_t              backgroundTRCone02 = -1.0;
1460   Double_t              backgroundTRCone04 = -1.0;
1461   Double_t              backgroundTRCone06 = -1.0;
1462   Double_t              backgroundTRCone08 = -1.0;
1463   Double_t              backgroundTRExact  = -1.0;
1464   Double_t              backgroundPP       = -1.0;
1465   // Calculate background for different jet exclusions
1466
1467   if (fAnalyzeBackground)
1468   {
1469
1470     if(fAnalyzeDeprecatedBackgrounds)
1471       GetKTBackgroundDensityAll    (fNumberExcludedJets, backgroundKTPbPb, backgroundKTPbPbWithGhosts, backgroundKTCMS, backgroundKTImprovedCMS, backgroundKTMean, backgroundKTTrackLike);
1472     else
1473       GetKTBackgroundDensity       (fNumberExcludedJets, backgroundKTImprovedCMS);
1474
1475     if(fAnalyzeDeprecatedBackgrounds)
1476       GetTRBackgroundDensity    (fNumberExcludedJets, backgroundTRNoExcl, backgroundTRCone02, backgroundTRCone04, backgroundTRCone06, backgroundTRCone08, backgroundTRExact);
1477
1478     backgroundKTImprovedCMSExternal = GetExternalRho();
1479   }
1480
1481   #ifdef DEBUGMODE
1482     AliInfo("Calculate()::Centrality&SignalJets&Background-Calculation done.");
1483   #endif
1484
1485   if (fAnalyzeQA)
1486   {
1487     FillHistogram("hVertexX",event->GetPrimaryVertex()->GetX());
1488     FillHistogram("hVertexY",event->GetPrimaryVertex()->GetY());
1489     FillHistogram("hVertexXY",event->GetPrimaryVertex()->GetX(), event->GetPrimaryVertex()->GetY());
1490     FillHistogram("hVertexR",TMath::Sqrt(event->GetPrimaryVertex()->GetX()*event->GetPrimaryVertex()->GetX() + event->GetPrimaryVertex()->GetY()*event->GetPrimaryVertex()->GetY()));
1491     FillHistogram("hCentralityV0M",centralityPercentileV0M);
1492     FillHistogram("hCentralityV0A",centralityPercentileV0A);
1493     FillHistogram("hCentralityV0C",centralityPercentileV0C);
1494     FillHistogram("hCentralityZNA",centralityPercentileZNA);
1495     FillHistogram("hCentrality",centralityPercentile);
1496
1497     Int_t trackCountAcc = 0;
1498     Int_t nTracks = fTrackArray->GetEntries();
1499     for (Int_t i = 0; i < nTracks; i++)
1500     {
1501       AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i));
1502
1503       if (track != 0)
1504         if (track->Pt() >= fMinTrackPt)
1505           FillHistogram("hTrackPhiEta", track->Phi(),track->Eta(), 1);
1506
1507       if (IsTrackInAcceptance(track))
1508       {
1509         FillHistogram("hTrackPt", track->Pt(), centralityPercentile);
1510         if(track->Eta() >= 0)
1511           FillHistogram("hTrackPtPosEta", track->Pt(), centralityPercentile);
1512         else
1513           FillHistogram("hTrackPtNegEta", track->Pt(), centralityPercentile);
1514
1515         FillHistogram("hTrackEta", track->Eta(), centralityPercentile);
1516         FillHistogram("hTrackPhi", track->Phi());
1517
1518         if(static_cast<AliPicoTrack*>(track))
1519           FillHistogram("hTrackPhiTrackType", track->Phi(), (static_cast<AliPicoTrack*>(track))->GetTrackType());
1520
1521         for(Int_t j=0;j<20;j++)
1522           if(track->Pt() > j)
1523             FillHistogram("hTrackPhiPtCut", track->Phi(), track->Pt());
1524
1525         FillHistogram("hTrackCharge", track->Charge());
1526         trackCountAcc++;
1527       }
1528     }
1529     FillHistogram("hTrackCountAcc", trackCountAcc, centralityPercentile);
1530
1531
1532 /*
1533     // This is code that is run locally to get some special events
1534     TFile* fileOutput = new TFile("SpecialEvents.root", "UPDATE");
1535     if(fSecondLeadingJet&&(fSecondLeadingJet->Pt()>10.))
1536     {
1537       cout << "Event found\n";
1538       TH2D* tmpEvent = new TH2D(Form("Event%lu", fEventCounter), "", 40, -0.9, 0.9, 40, 0., TMath::TwoPi());
1539       tmpEvent->GetXaxis()->SetTitle("#eta");
1540       tmpEvent->GetYaxis()->SetTitle("#phi");
1541       tmpEvent->SetOption("LEGO2");
1542       tmpEvent->Sumw2();
1543
1544       for (Int_t i = 0; i < nTracks; i++)
1545       {
1546         AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i));
1547
1548         if (IsTrackInAcceptance(track))
1549         {
1550           tmpEvent->Fill(track->Eta(), track->Phi(), track->Pt());
1551         }
1552       }
1553       tmpEvent->Write(0, kOverwrite);
1554     }
1555     fileOutput->Close();
1556 */
1557   }
1558   #ifdef DEBUGMODE
1559     AliInfo("Calculate()::QA done.");
1560   #endif
1561
1562   ////////////////////// NOTE: Jet analysis and calculations
1563
1564   if (fAnalyzeJets)
1565   {
1566     for (Int_t i = 0; i<fJetArray->GetEntries(); i++)
1567     {
1568       AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
1569       if (!tmpJet)
1570         continue;
1571
1572       FillHistogram("hRawJetPt", tmpJet->Pt());
1573       if (tmpJet->Pt() >= 5.0)
1574       {
1575       // ### RAW JET ANALYSIS
1576         if (tmpJet->Area() >= fMinJetArea)
1577           FillHistogram("hRawJetPhiEta", tmpJet->Phi(), tmpJet->Eta());
1578         if (TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow)
1579           FillHistogram("hRawJetArea", tmpJet->Area());
1580       }
1581
1582       if(IsSignalJetInAcceptance(tmpJet))
1583       {
1584       // ### SIGNAL JET ANALYSIS
1585         // Jet spectra
1586         FillHistogram("hJetPt", tmpJet->Pt(), centralityPercentile);
1587         FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile);
1588         FillHistogram("hJetPtBgrdSubtractedExternal", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal), centralityPercentile);
1589         FillHistogram("hJetPtSubtractedRhoKTImprovedCMS", tmpJet->Pt(), centralityPercentile, backgroundKTImprovedCMS);
1590
1591         // pp background
1592         GetPPBackgroundDensity(backgroundPP, tmpJet);
1593         FillHistogram("hJetPtBgrdSubtractedPP", GetCorrectedJetPt(tmpJet, backgroundPP), centralityPercentile);
1594
1595         if(fAnalyzeDeprecatedBackgrounds)
1596         {
1597           FillHistogram("hJetPtBgrdSubtractedTR", GetCorrectedJetPt(tmpJet, backgroundTRCone06), centralityPercentile);
1598           FillHistogram("hJetPtBgrdSubtractedKTPbPb", GetCorrectedJetPt(tmpJet, backgroundKTPbPb), centralityPercentile);
1599           FillHistogram("hJetPtBgrdSubtractedKTPbPbWithGhosts", GetCorrectedJetPt(tmpJet, backgroundKTPbPbWithGhosts), centralityPercentile);
1600           FillHistogram("hJetPtBgrdSubtractedKTCMS", GetCorrectedJetPt(tmpJet, backgroundKTCMS), centralityPercentile);
1601           FillHistogram("hJetPtBgrdSubtractedKTMean", GetCorrectedJetPt(tmpJet, backgroundKTMean), centralityPercentile);
1602           FillHistogram("hJetPtBgrdSubtractedKTTrackLike", GetCorrectedJetPt(tmpJet, backgroundKTTrackLike), centralityPercentile);
1603         }
1604
1605         // Jet profile analysis
1606         if(TMath::Abs(tmpJet->Eta()) <= 0.3)
1607         {
1608           if(tmpJet->Pt()>=70.0)
1609           {
1610             FillHistogram("hJetProfile70GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1611             FillHistogram("hJetProfile70GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1612             FillHistogram("hJetProfile70GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1613             FillHistogram("hJetProfile70GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1614             FillHistogram("hJetProfile70GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1615             FillHistogram("hJetProfile70GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1616             FillHistogram("hJetProfile70GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1617             FillHistogram("hJetProfile70GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1618             FillHistogram("hJetProfile70GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1619             FillHistogram("hJetProfile70GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1620             FillHistogram("hJetProfile70GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1621             FillHistogram("hJetProfile70GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1622           }
1623           else if(GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal)>=60.0)
1624           {
1625             FillHistogram("hJetProfile60GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1626             FillHistogram("hJetProfile60GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1627             FillHistogram("hJetProfile60GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1628             FillHistogram("hJetProfile60GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1629             FillHistogram("hJetProfile60GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1630             FillHistogram("hJetProfile60GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1631             FillHistogram("hJetProfile60GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1632             FillHistogram("hJetProfile60GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1633             FillHistogram("hJetProfile60GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1634             FillHistogram("hJetProfile60GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1635             FillHistogram("hJetProfile60GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1636             FillHistogram("hJetProfile60GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1637           }
1638           else if(GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal)>=50.0)
1639           {
1640             FillHistogram("hJetProfile50GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1641             FillHistogram("hJetProfile50GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1642             FillHistogram("hJetProfile50GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1643             FillHistogram("hJetProfile50GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1644             FillHistogram("hJetProfile50GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1645             FillHistogram("hJetProfile50GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1646             FillHistogram("hJetProfile50GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1647             FillHistogram("hJetProfile50GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1648             FillHistogram("hJetProfile50GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1649             FillHistogram("hJetProfile50GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1650             FillHistogram("hJetProfile50GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1651             FillHistogram("hJetProfile50GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1652           }
1653           else if(GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal)>=40.0)
1654           {
1655             FillHistogram("hJetProfile40GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1656             FillHistogram("hJetProfile40GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1657             FillHistogram("hJetProfile40GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1658             FillHistogram("hJetProfile40GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1659             FillHistogram("hJetProfile40GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1660             FillHistogram("hJetProfile40GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1661             FillHistogram("hJetProfile40GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1662             FillHistogram("hJetProfile40GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1663             FillHistogram("hJetProfile40GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1664             FillHistogram("hJetProfile40GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1665             FillHistogram("hJetProfile40GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1666             FillHistogram("hJetProfile40GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1667           }
1668           else if(GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal)>=30.0)
1669           {
1670             FillHistogram("hJetProfile30GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1671             FillHistogram("hJetProfile30GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1672             FillHistogram("hJetProfile30GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1673             FillHistogram("hJetProfile30GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1674             FillHistogram("hJetProfile30GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1675             FillHistogram("hJetProfile30GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1676             FillHistogram("hJetProfile30GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1677             FillHistogram("hJetProfile30GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1678             FillHistogram("hJetProfile30GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1679             FillHistogram("hJetProfile30GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1680             FillHistogram("hJetProfile30GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1681             FillHistogram("hJetProfile30GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1682           }
1683           else if(GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal)>=20.0)
1684           {
1685             FillHistogram("hJetProfile20GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1686             FillHistogram("hJetProfile20GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1687             FillHistogram("hJetProfile20GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1688             FillHistogram("hJetProfile20GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1689             FillHistogram("hJetProfile20GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1690             FillHistogram("hJetProfile20GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1691             FillHistogram("hJetProfile20GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1692             FillHistogram("hJetProfile20GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1693             FillHistogram("hJetProfile20GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1694             FillHistogram("hJetProfile20GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1695             FillHistogram("hJetProfile20GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1696             FillHistogram("hJetProfile20GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1697           }
1698           else if(GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal)>=10.0)
1699           {
1700             FillHistogram("hJetProfile10GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1701             FillHistogram("hJetProfile10GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1702             FillHistogram("hJetProfile10GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1703             FillHistogram("hJetProfile10GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1704             FillHistogram("hJetProfile10GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1705             FillHistogram("hJetProfile10GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1706             FillHistogram("hJetProfile10GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1707             FillHistogram("hJetProfile10GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1708             FillHistogram("hJetProfile10GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1709             FillHistogram("hJetProfile10GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1710             FillHistogram("hJetProfile10GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1711             FillHistogram("hJetProfile10GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, backgroundKTImprovedCMSExternal))/GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMSExternal));
1712           }
1713         }
1714         FillHistogram("hJetPtVsConstituentCount", tmpJet->Pt(),tmpJet->GetNumberOfTracks());
1715
1716         if((fAnalyzeQA) && (tmpJet->Pt() >= 5.0))
1717         {
1718           Double_t lowestTrackPt = 1e99;
1719           Double_t highestTrackPt = 0.0;
1720           for(Int_t j=0; j<tmpJet->GetNumberOfTracks(); j++)
1721           {
1722             FillHistogram("hJetConstituentPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
1723             // Find the lowest pT of a track in the jet
1724             if (tmpJet->TrackAt(j, fTrackArray)->Pt() < lowestTrackPt)
1725               lowestTrackPt = tmpJet->TrackAt(j, fTrackArray)->Pt();
1726             if (tmpJet->TrackAt(j, fTrackArray)->Pt() > highestTrackPt)
1727               highestTrackPt = tmpJet->TrackAt(j, fTrackArray)->Pt();
1728           }
1729           FillHistogram("hJetArea", tmpJet->Area());
1730           // Signal jet vs. signal jet - "Combinatorial"
1731           for (Int_t j = i+1; j<fNumberSignalJets; j++)
1732             FillHistogram("hJetDeltaPhi", GetDeltaPhi(tmpJet->Phi(), fSignalJets[j]->Phi()));
1733
1734           FillHistogram("hJetPhiEta", tmpJet->Phi(),tmpJet->Eta());
1735           FillHistogram("hJetEta", tmpJet->Eta(), centralityPercentile);
1736
1737           if(lowestTrackPt>=2.0)
1738             FillHistogram("hJetEta2GeVTracks", tmpJet->Eta(), centralityPercentile);
1739           if(lowestTrackPt>=4.0)
1740             FillHistogram("hJetEta4GeVTracks", tmpJet->Eta(), centralityPercentile);
1741         }
1742
1743
1744         // Jet constituent mass analyses (for PYTHIA)
1745         if (fAnalyzeMassCorrelation)
1746         {
1747           Double_t jetMass = 0.0;
1748           Int_t constituentCount = tmpJet->GetNumberOfTracks();
1749           Int_t electronCount = 0;
1750           Int_t pionCount = 0;
1751           Int_t protonCount = 0;
1752           Int_t kaonCount = 0;
1753           Int_t othersCount = 0;
1754
1755           for(Int_t j=0; j<tmpJet->GetNumberOfTracks(); j++)
1756           {
1757             AliVParticle* tmpParticle = tmpJet->TrackAt(j, fTrackArray);
1758
1759             jetMass += tmpParticle->M();
1760             if(TMath::Abs(tmpParticle->PdgCode()) == 11) // electron, positron
1761               electronCount++;
1762             else if(TMath::Abs(tmpParticle->PdgCode()) == 211) // p-, p+
1763               pionCount++;
1764             else if(TMath::Abs(tmpParticle->PdgCode()) == 2212) // p, pbar
1765               protonCount++;
1766             else if(TMath::Abs(tmpParticle->PdgCode()) == 321) // kaon+,kaon-
1767               kaonCount++;
1768             else
1769               othersCount++;
1770           }
1771
1772           FillHistogram("hJetMassFromConstituents", jetMass);
1773           FillHistogram("hJetMass", tmpJet->M());
1774
1775           FillHistogram("hJetPtVsMass", tmpJet->Pt(), jetMass);
1776           FillHistogram("hJetPtVsJetMass", tmpJet->Pt(), tmpJet->M());
1777
1778           FillHistogram("hJetPtVsProtonCount", tmpJet->Pt(), protonCount/(static_cast<Double_t>(constituentCount)));
1779           FillHistogram("hJetPtVsPionCount", tmpJet->Pt(), pionCount/(static_cast<Double_t>(constituentCount)));
1780           FillHistogram("hJetPtVsKaonCount", tmpJet->Pt(), kaonCount/(static_cast<Double_t>(constituentCount)));
1781           FillHistogram("hJetPtVsElectronCount", tmpJet->Pt(), electronCount/(static_cast<Double_t>(constituentCount)));
1782           FillHistogram("hJetPtVsOthersCount", tmpJet->Pt(), othersCount/(static_cast<Double_t>(constituentCount)));
1783           FillHistogram("hJetConstituentProtonCount", protonCount);
1784           FillHistogram("hJetConstituentPionCount", pionCount);
1785           FillHistogram("hJetConstituentKaonCount", kaonCount);
1786           FillHistogram("hJetConstituentElectronCount", electronCount);
1787           FillHistogram("hJetConstituentOthersCount", othersCount);
1788
1789           // Results for the "normal" jet (avoiding single particle jets)
1790           if((constituentCount>=6) && (constituentCount<=14))
1791           {
1792             FillHistogram("hJetPtVsMass_6_14", tmpJet->Pt(), jetMass);
1793             FillHistogram("hJetPtVsJetMass_6_14", tmpJet->Pt(), tmpJet->M());
1794
1795             FillHistogram("hJetPtVsProtonCount_6_14", tmpJet->Pt(), protonCount);
1796             FillHistogram("hJetPtVsPionCount_6_14", tmpJet->Pt(), pionCount);
1797             FillHistogram("hJetPtVsKaonCount_6_14", tmpJet->Pt(), kaonCount);
1798             FillHistogram("hJetPtVsElectronCount_6_14", tmpJet->Pt(), electronCount);
1799             FillHistogram("hJetPtVsOthersCount_6_14", tmpJet->Pt(), othersCount);
1800             FillHistogram("hJetConstituentProtonCount_6_14", protonCount);
1801             FillHistogram("hJetConstituentPionCount_6_14", pionCount);
1802             FillHistogram("hJetConstituentKaonCount_6_14", kaonCount);
1803             FillHistogram("hJetConstituentElectronCount_6_14", electronCount);
1804             FillHistogram("hJetConstituentOthersCount_6_14", othersCount);
1805           }
1806           if((constituentCount>=2))
1807           {
1808             FillHistogram("hJetPtVsMass_2_X", tmpJet->Pt(), jetMass);
1809             FillHistogram("hJetPtVsJetMass_2_X", tmpJet->Pt(), tmpJet->M());
1810
1811             FillHistogram("hJetPtVsProtonCount_2_X", tmpJet->Pt(), protonCount);
1812             FillHistogram("hJetPtVsPionCount_2_X", tmpJet->Pt(), pionCount);
1813             FillHistogram("hJetPtVsKaonCount_2_X", tmpJet->Pt(), kaonCount);
1814             FillHistogram("hJetPtVsElectronCount_2_X", tmpJet->Pt(), electronCount);
1815             FillHistogram("hJetPtVsOthersCount_2_X", tmpJet->Pt(), othersCount);
1816             FillHistogram("hJetConstituentProtonCount_2_X", protonCount);
1817             FillHistogram("hJetConstituentPionCount_2_X", pionCount);
1818             FillHistogram("hJetConstituentKaonCount_2_X", kaonCount);
1819             FillHistogram("hJetConstituentElectronCount_2_X", electronCount);
1820             FillHistogram("hJetConstituentOthersCount_2_X", othersCount);
1821           }
1822           if((constituentCount>=2) && (constituentCount<=7))
1823           {
1824             FillHistogram("hJetPtVsMass_2_7", tmpJet->Pt(), jetMass);
1825             FillHistogram("hJetPtVsJetMass_2_7", tmpJet->Pt(), tmpJet->M());
1826
1827             FillHistogram("hJetPtVsProtonCount_2_7", tmpJet->Pt(), protonCount);
1828             FillHistogram("hJetPtVsPionCount_2_7", tmpJet->Pt(), pionCount);
1829             FillHistogram("hJetPtVsKaonCount_2_7", tmpJet->Pt(), kaonCount);
1830             FillHistogram("hJetPtVsElectronCount_2_7", tmpJet->Pt(), electronCount);
1831             FillHistogram("hJetPtVsOthersCount_2_7", tmpJet->Pt(), othersCount);
1832             FillHistogram("hJetConstituentProtonCount_2_7", protonCount);
1833             FillHistogram("hJetConstituentPionCount_2_7", pionCount);
1834             FillHistogram("hJetConstituentKaonCount_2_7", kaonCount);
1835             FillHistogram("hJetConstituentElectronCount_2_7", electronCount);
1836             FillHistogram("hJetConstituentOthersCount_2_7", othersCount);
1837           }
1838         }
1839       }
1840     }
1841
1842     // ### DIJETS
1843     if(fNumberSignalJets >= 2)
1844     {
1845       FillHistogram("hLeadingJetDeltaPhi", GetDeltaPhi(fFirstLeadingJet->Phi(), fSecondLeadingJet->Phi()));
1846
1847       if (IsDijet(fFirstLeadingJet, fSecondLeadingJet))
1848       {
1849         FillHistogram("hDijetConstituentsPt", fFirstLeadingJet->Pt());
1850         FillHistogram("hDijetConstituentsPt", fSecondLeadingJet->Pt());
1851
1852         FillHistogram("hDijetLeadingJetPt", fFirstLeadingJet->Pt());
1853         FillHistogram("hDijetPtCorrelation", fFirstLeadingJet->Pt(), fSecondLeadingJet->Pt());
1854         Double_t dummyArea = 0;
1855         GetTRBackgroundDensity (2, backgroundDijet, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kFALSE);
1856         GetTRBackgroundDensity (2, backgroundDijetPerpendicular, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kTRUE);
1857         deltaPtDijetPerpendicular = GetDeltaPt(backgroundDijetPerpendicular);
1858       }
1859     }
1860
1861     // ### SOME JET PLOTS
1862     FillHistogram("hJetCountAll", fJetArray->GetEntries());
1863     FillHistogram("hJetCountAccepted", fNumberSignalJets);
1864     FillHistogram("hJetCount", fJetArray->GetEntries(), fNumberSignalJets);
1865     if (fFirstLeadingJet)
1866     {
1867       FillHistogram("hLeadingJetPt", fFirstLeadingJet->Pt());
1868       FillHistogram("hCorrectedLeadingJetPt", GetCorrectedJetPt(fFirstLeadingJet,backgroundKTImprovedCMS));
1869     }
1870     if (fSecondLeadingJet)
1871     {
1872       FillHistogram("hSecondLeadingJetPt", fSecondLeadingJet->Pt());
1873       FillHistogram("hCorrectedSecondLeadingJetPt", GetCorrectedJetPt(fSecondLeadingJet,backgroundKTImprovedCMS));
1874     }
1875   } //endif AnalyzeJets
1876
1877   #ifdef DEBUGMODE
1878     AliInfo("Calculate()::Jets done.");
1879   #endif
1880   ////////////////////// NOTE: Background analysis
1881
1882   if (fAnalyzeBackground)
1883   {
1884     // Calculate background in centrality classes
1885     FillHistogram("hKTBackgroundImprovedCMS", backgroundKTImprovedCMS, centralityPercentile);
1886
1887     FillHistogram("hKTBackgroundImprovedCMSExternal", backgroundKTImprovedCMSExternal, centralityPercentile);
1888
1889     FillHistogram("hKTMeanBackgroundImprovedCMS", centralityPercentile, backgroundKTImprovedCMS);
1890
1891     // In case of dijets -> look at the background
1892     if (backgroundDijet >= 0)
1893       FillHistogram("hDijetBackground", backgroundDijet, centralityPercentile);
1894     if (backgroundDijetPerpendicular >= 0)
1895       FillHistogram("hDijetBackgroundPerpendicular", backgroundDijetPerpendicular, centralityPercentile);
1896     if (deltaPtDijetPerpendicular > -10000.)
1897       FillHistogram("hDijetBackgroundPerpendicularDeltaPt", deltaPtDijetPerpendicular, centralityPercentile);
1898
1899
1900     if(fAnalyzeDeprecatedBackgrounds)
1901     {
1902       FillHistogram("hKTBackgroundPbPb", backgroundKTPbPb, centralityPercentile);
1903       FillHistogram("hKTBackgroundPbPbWithGhosts", backgroundKTPbPbWithGhosts, centralityPercentile);
1904       FillHistogram("hKTBackgroundCMS", backgroundKTCMS, centralityPercentile);
1905       FillHistogram("hKTBackgroundMean", backgroundKTMean, centralityPercentile);
1906       FillHistogram("hKTBackgroundTrackLike", backgroundKTTrackLike, centralityPercentile);
1907
1908       FillHistogram("hTRBackgroundNoExcl", backgroundTRNoExcl, centralityPercentile);
1909       FillHistogram("hTRBackgroundCone02", backgroundTRCone02, centralityPercentile);
1910       FillHistogram("hTRBackgroundCone04", backgroundTRCone04, centralityPercentile);
1911       FillHistogram("hTRBackgroundCone06", backgroundTRCone06, centralityPercentile);
1912       FillHistogram("hTRBackgroundCone08", backgroundTRCone08, centralityPercentile);
1913       FillHistogram("hTRBackgroundExact", backgroundTRExact, centralityPercentile);
1914
1915       // Calculate background profiles in terms of centrality
1916       FillHistogram("hKTMeanBackgroundPbPb", centralityPercentile,  backgroundKTPbPb);
1917       FillHistogram("hKTMeanBackgroundPbPbWithGhosts", centralityPercentile,  backgroundKTPbPbWithGhosts);
1918       FillHistogram("hKTMeanBackgroundCMS", centralityPercentile, backgroundKTCMS);
1919       FillHistogram("hKTMeanBackgroundMean", centralityPercentile, backgroundKTMean);
1920       FillHistogram("hKTMeanBackgroundTPC", centralityPercentile, backgroundKTTrackLike);
1921       FillHistogram("hTRMeanBackground", centralityPercentile,  backgroundTRCone06);
1922     }
1923
1924
1925     // Calculate the delta pt
1926     Double_t tmpDeltaPtNoBackground = GetDeltaPt(0.0);
1927     Double_t tmpDeltaPtKTImprovedCMS = GetDeltaPt(backgroundKTImprovedCMS);
1928
1929     Double_t tmpDeltaPtKTImprovedCMSPartialExclusion = 0.0;
1930     if(fNcoll)
1931       tmpDeltaPtKTImprovedCMSPartialExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0/fNcoll);
1932     else
1933       tmpDeltaPtKTImprovedCMSPartialExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0);
1934
1935     Double_t tmpDeltaPtKTImprovedCMSPartialExclusion_Signal = 0.0;
1936     if(fNumberSignalJets)
1937       tmpDeltaPtKTImprovedCMSPartialExclusion_Signal = GetDeltaPt(backgroundKTImprovedCMS, 1.0/fNumberSignalJets);
1938     else
1939       tmpDeltaPtKTImprovedCMSPartialExclusion_Signal = GetDeltaPt(backgroundKTImprovedCMS, 1.0);
1940
1941     Double_t tmpDeltaPtKTImprovedCMSFullExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0);
1942
1943     Double_t tmpDeltaPtKTPbPb = 0;
1944     Double_t tmpDeltaPtKTPbPbWithGhosts = 0;
1945     Double_t tmpDeltaPtKTCMS = 0;
1946     Double_t tmpDeltaPtKTMean = 0;
1947     Double_t tmpDeltaPtKTTrackLike = 0;
1948     Double_t tmpDeltaPtTR = 0;
1949
1950     if(fAnalyzeDeprecatedBackgrounds)
1951     {
1952       tmpDeltaPtKTPbPb = GetDeltaPt(backgroundKTPbPb);
1953       tmpDeltaPtKTPbPbWithGhosts = GetDeltaPt(backgroundKTPbPbWithGhosts);
1954       tmpDeltaPtKTCMS = GetDeltaPt(backgroundKTCMS);
1955       tmpDeltaPtKTMean = GetDeltaPt(backgroundKTMean);
1956       tmpDeltaPtKTTrackLike = GetDeltaPt(backgroundKTTrackLike);
1957       tmpDeltaPtTR = GetDeltaPt(backgroundTRCone06);
1958     }
1959
1960     // If valid, fill the delta pt histograms
1961
1962     if(tmpDeltaPtKTImprovedCMS > -10000.0)
1963       FillHistogram("hDeltaPtKTImprovedCMS", tmpDeltaPtKTImprovedCMS, centralityPercentile);
1964     if(tmpDeltaPtKTImprovedCMSPartialExclusion > -10000.0)
1965       FillHistogram("hDeltaPtKTImprovedCMSPartialExclusion", tmpDeltaPtKTImprovedCMSPartialExclusion, centralityPercentile);
1966     if(tmpDeltaPtKTImprovedCMSPartialExclusion_Signal > -10000.0)
1967       FillHistogram("hDeltaPtKTImprovedCMSPartialExclusion_Signal", tmpDeltaPtKTImprovedCMSPartialExclusion_Signal, centralityPercentile);
1968     if(tmpDeltaPtKTImprovedCMSFullExclusion > -10000.0)
1969       FillHistogram("hDeltaPtKTImprovedCMSFullExclusion", tmpDeltaPtKTImprovedCMSFullExclusion, centralityPercentile);
1970
1971     if(tmpDeltaPtNoBackground > 0.000001)
1972       FillHistogram("hDeltaPtNoBackgroundNoEmptyCones", tmpDeltaPtNoBackground, centralityPercentile);
1973     else if(tmpDeltaPtNoBackground > -10000.0)
1974       FillHistogram("hDeltaPtNoBackground", tmpDeltaPtNoBackground, centralityPercentile);
1975
1976
1977     if(fAnalyzeDeprecatedBackgrounds)
1978     {
1979       if(tmpDeltaPtKTPbPb > -10000.0)
1980         FillHistogram("hDeltaPtKTPbPb", tmpDeltaPtKTPbPb, centralityPercentile);
1981       if(tmpDeltaPtKTPbPbWithGhosts > -10000.0)
1982         FillHistogram("hDeltaPtKTPbPbWithGhosts", tmpDeltaPtKTPbPbWithGhosts, centralityPercentile);
1983       if(tmpDeltaPtKTCMS > -10000.0)
1984         FillHistogram("hDeltaPtKTCMS", tmpDeltaPtKTCMS, centralityPercentile);
1985       if(tmpDeltaPtKTMean > -10000.0)
1986         FillHistogram("hDeltaPtKTMean", tmpDeltaPtKTMean, centralityPercentile);
1987       if(tmpDeltaPtKTTrackLike > -10000.0)
1988         FillHistogram("hDeltaPtKTTrackLike", tmpDeltaPtKTTrackLike, centralityPercentile);
1989
1990       if(tmpDeltaPtTR > -10000.0)
1991         FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile);
1992     }
1993   }
1994
1995   #ifdef DEBUGMODE
1996     AliInfo("Calculate()::Background done.");
1997   #endif
1998
1999   ////////////////////// NOTE: Pythia histograms
2000   if(fAnalyzePythia)
2001   {
2002     FillHistogram("hPythiaPtHard", GetPtHard());
2003     FillHistogram("hPythiaNTrials", GetPtHardBin()+0.1, GetPythiaTrials());
2004     FillHistogram("hPythiaXSection", GetPtHardBin()+0.1, fCrossSection);
2005
2006         #ifdef DEBUGMODE
2007       AliInfo("Calculate()::Pythia done.");
2008     #endif
2009   }
2010   #ifdef DEBUGMODE
2011     AliInfo("Calculate() done.");
2012   #endif
2013 }
2014
2015 //________________________________________________________________________
2016 Bool_t AliAnalysisTaskChargedJetsPA::UserNotify()
2017 {
2018   // Implemented Notify() to read the cross sections
2019   // and number of trials from pyxsec.root
2020   //
2021   #ifdef DEBUGMODE
2022     AliInfo("UserNotify started.");
2023   #endif
2024
2025   if(fAnalyzePythia)
2026   {
2027     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
2028     TFile *currFile = tree->GetCurrentFile();
2029
2030     TString file(currFile->GetName());
2031
2032     if(file.Contains("root_archive.zip#")){
2033       Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
2034       Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
2035       file.Replace(pos+1,20,"");
2036     }
2037     else {
2038       // not an archive take the basename....
2039       file.ReplaceAll(gSystem->BaseName(file.Data()),"");
2040     }
2041
2042     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
2043     if(!fxsec){
2044       // next trial fetch the histgram file
2045       fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
2046       if(!fxsec){
2047           // not a severe condition but inciate that we have no information
2048         return kFALSE;
2049       }
2050       else{
2051         // find the tlist we want to be independtent of the name so use the Tkey
2052         TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
2053         if(!key){
2054           fxsec->Close();
2055           return kFALSE;
2056         }
2057         TList *list = dynamic_cast<TList*>(key->ReadObj());
2058         if(!list){
2059           fxsec->Close();
2060           return kFALSE;
2061         }
2062         fCrossSection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
2063         fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
2064         fxsec->Close();
2065       }
2066     } // no tree pyxsec.root
2067     else {
2068       TTree *xtree = (TTree*)fxsec->Get("Xsection");
2069       if(!xtree){
2070         fxsec->Close();
2071         return kFALSE;
2072       }
2073       UInt_t   ntrials  = 0;
2074       Double_t  xsection  = 0;
2075       xtree->SetBranchAddress("xsection",&xsection);
2076       xtree->SetBranchAddress("ntrials",&ntrials);
2077       xtree->GetEntry(0);
2078       fTrials = ntrials;
2079       fCrossSection = xsection;
2080       fxsec->Close();
2081     }
2082   }
2083   #ifdef DEBUGMODE
2084     AliInfo("UserNotify ended.");
2085   #endif
2086   return kTRUE;
2087 }
2088
2089
2090 //________________________________________________________________________
2091 inline Double_t AliAnalysisTaskChargedJetsPA::EtaToTheta(Double_t arg)
2092   {return 2.*atan(exp(-arg));}
2093 //________________________________________________________________________
2094 inline Double_t AliAnalysisTaskChargedJetsPA::ThetaToEta(Double_t arg)
2095 {
2096   if ((arg > TMath::Pi()) || (arg < 0.0))
2097   {
2098     AliError(Form("ThetaToEta got wrong input! (%f)", arg));
2099     return 0.0;
2100   }
2101   return -log(tan(arg/2.));
2102 }
2103 //________________________________________________________________________
2104 inline Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPhi(Double_t phi1, Double_t phi2)
2105   {return min(TMath::Abs(phi1-phi2),TMath::TwoPi() - TMath::Abs(phi1-phi2));}
2106
2107 //________________________________________________________________________
2108 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)
2109 {
2110   const Int_t kTests = 1000;
2111   Int_t hits = 0;
2112   TRandom3 randomGen(0);
2113
2114   // Loop over kTests-many tests
2115   for (Int_t i=0; i<kTests; i++)
2116   {
2117     //Choose random position in rectangle for the tester
2118     Double_t tmpTestX = randomGen.Uniform(rPosXmin, rPosXmax);
2119     Double_t tmpTestY = randomGen.Uniform(rPosYmin, rPosYmax);
2120
2121     //Check, if tester is in circle. If yes, increment circle counter.
2122     Double_t tmpDistance = TMath::Sqrt( (tmpTestX - cPosX)*(tmpTestX - cPosX) + (tmpTestY - cPosY)*(tmpTestY - cPosY) );
2123     if(tmpDistance < cRadius)
2124       hits++;
2125   }
2126
2127   // return ratio
2128   return (static_cast<Double_t>(hits)/static_cast<Double_t>(kTests));
2129 }
2130
2131 //________________________________________________________________________
2132 Double_t AliAnalysisTaskChargedJetsPA::MCGetOverlapMultipleCirclesRectancle(Int_t numCircles, std::vector<Double_t> cPosX, std::vector<Double_t> cPosY, Double_t cRadius, Double_t rPosXmin, Double_t rPosXmax, Double_t rPosYmin, Double_t rPosYmax)
2133 {
2134
2135   const Int_t kTests = 1000;
2136   Int_t hits = 0;
2137   TRandom3 randomGen(0);
2138
2139   // Loop over kTests-many tests
2140   for (Int_t i=0; i<kTests; i++)
2141   {
2142     //Choose random position in rectangle for the tester
2143     Double_t tmpTestX = randomGen.Uniform(rPosXmin, rPosXmax);
2144     Double_t tmpTestY = randomGen.Uniform(rPosYmin, rPosYmax);
2145
2146     //Check, if tester is in one of the circles. If yes, increment circle counter.
2147     for(Int_t j=0; j<numCircles; j++)
2148     {
2149       Double_t tmpDistance = TMath::Sqrt( (tmpTestX - cPosX[j])*(tmpTestX - cPosX[j]) + (tmpTestY - cPosY[j])*(tmpTestY - cPosY[j]) );
2150       if(tmpDistance < cRadius)
2151       {
2152         hits++;
2153         break;
2154       }
2155     }
2156   }
2157
2158   // return ratio
2159   return (static_cast<Double_t>(hits)/static_cast<Double_t>(kTests));
2160
2161 }
2162
2163 //________________________________________________________________________
2164 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x)
2165 {
2166   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(key)));
2167   if(!tmpHist)
2168   {
2169     AliError(Form("Cannot find histogram <%s> ",key)) ;
2170     return;
2171   }
2172
2173   tmpHist->Fill(x);
2174 }
2175
2176 //________________________________________________________________________
2177 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y)
2178 {
2179   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(key)));
2180   if(!tmpHist)
2181   {
2182     AliError(Form("Cannot find histogram <%s> ",key));
2183     return;
2184   }
2185
2186   if (tmpHist->IsA()->GetBaseClass("TH1"))
2187     static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
2188   else if (tmpHist->IsA()->GetBaseClass("TH2"))
2189     static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
2190 }
2191
2192 //________________________________________________________________________
2193 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y, Double_t add)
2194 {
2195   TH2* tmpHist = static_cast<TH2*>(fOutputList->FindObject(GetHistoName(key)));
2196   if(!tmpHist)
2197   {
2198     AliError(Form("Cannot find histogram <%s> ",key));
2199     return;
2200   }
2201
2202   tmpHist->Fill(x,y,add);
2203 }
2204 //________________________________________________________________________
2205 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)
2206 {
2207   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax);
2208
2209   tmpHist->GetXaxis()->SetTitle(xTitle);
2210   tmpHist->GetYaxis()->SetTitle(yTitle);
2211   tmpHist->SetOption(options);
2212   tmpHist->SetMarkerStyle(kFullCircle);
2213   tmpHist->Sumw2();
2214
2215   fHistList->Add(tmpHist);
2216   fHistCount++;
2217
2218   return tmpHist;
2219 }
2220
2221 //________________________________________________________________________
2222 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)
2223 {
2224   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax, yBins, yMin, yMax);
2225   tmpHist->GetXaxis()->SetTitle(xTitle);
2226   tmpHist->GetYaxis()->SetTitle(yTitle);
2227   tmpHist->GetZaxis()->SetTitle(zTitle);
2228   tmpHist->SetOption(options);
2229   tmpHist->SetMarkerStyle(kFullCircle);
2230   tmpHist->Sumw2();
2231
2232   fHistList->Add(tmpHist);
2233   fHistCount++;
2234
2235   return tmpHist;
2236 }
2237
2238 //________________________________________________________________________
2239 void AliAnalysisTaskChargedJetsPA::Terminate(Option_t *)
2240 {
2241   PostData(1, fOutputList);
2242
2243   // Mandatory
2244   fOutputList = dynamic_cast<TList*> (GetOutputData(1)); // '1' refers to the output slot
2245   if (!fOutputList) {
2246     printf("ERROR: Output list not available\n");
2247     return;
2248   }
2249 }
2250
2251 //________________________________________________________________________
2252 AliAnalysisTaskChargedJetsPA::~AliAnalysisTaskChargedJetsPA()
2253 {
2254   // Destructor. Clean-up the output list, but not the histograms that are put inside
2255   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
2256   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
2257     delete fOutputList;
2258   }
2259 }
2260
2261 //________________________________________________________________________
2262 void AliAnalysisTaskChargedJetsPA::UserCreateOutputObjects()
2263 {
2264   // called once to create user defined output objects like histograms, plots etc.
2265   // and to put it on the output list.
2266   // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
2267
2268   fRandom = new TRandom3(0);
2269
2270
2271   fOutputList = new TList();
2272   fOutputList->SetOwner(); // otherwise it produces leaks in merging
2273
2274   // NOTE: Pythia histograms
2275   AddHistogram1D<TProfile>("hPythiaXSection", "Pythia cross section distribution", "", fNumPtHardBins+1, 0, fNumPtHardBins+1, "p_{T} hard bin","dN^{Events}/dp_{T,hard}");
2276   AddHistogram1D<TH1D>("hPythiaNTrials", "Pythia trials (no correction for manual cuts)", "", fNumPtHardBins+1, 0, fNumPtHardBins+1, "p_{T} hard bin", "Trials");
2277
2278
2279   PostData(1, fOutputList);
2280 }
2281
2282 //________________________________________________________________________
2283 void AliAnalysisTaskChargedJetsPA::UserExec(Option_t *)
2284 {
2285   #ifdef DEBUGMODE
2286     AliInfo("UserExec() started.");
2287   #endif
2288
2289   if (!InputEvent())
2290   {
2291     AliError("??? Event pointer == 0 ???");
2292     return;
2293   }
2294
2295   if (!fInitialized)
2296     ExecOnce(); // Get tracks, jets, background from arrays if not already given + Init Histos
2297
2298   Calculate(InputEvent());
2299
2300   PostData(1, fOutputList);
2301 }