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