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