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