PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx

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