]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliAnalysisTaskJetCluster.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bug fix and small changed in the macros
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetCluster.cxx
CommitLineData
dbf053f0 1// **************************************
2// Task used for the correction of determiantion of reconstructed jet spectra
3// Compares input (gen) and output (rec) jets
4// *******************************************
5
6
7/**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9 * *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
12 * *
13 * Permission to use, copy, modify and distribute this software and its *
14 * documentation strictly for non-commercial purposes is hereby granted *
15 * without fee, provided that the above copyright notice appears in all *
16 * copies and that both the copyright notice and this permission notice *
17 * appear in the supporting documentation. The authors make no claims *
18 * about the suitability of this software for any purpose. It is *
19 * provided "as is" without express or implied warranty. *
20 **************************************************************************/
21
22
23#include <TROOT.h>
24#include <TRandom.h>
25#include <TSystem.h>
26#include <TInterpreter.h>
27#include <TChain.h>
28#include <TRandom.h>
29#include <TFile.h>
30#include <TKey.h>
31#include <TH1F.h>
32#include <TH2F.h>
33#include <TH3F.h>
34#include <TProfile.h>
35#include <TList.h>
36#include <TLorentzVector.h>
37#include <TClonesArray.h>
38#include "TDatabasePDG.h"
39
40#include "AliAnalysisTaskJetCluster.h"
41#include "AliAnalysisManager.h"
42#include "AliJetFinder.h"
43#include "AliJetHeader.h"
44#include "AliJetReader.h"
45#include "AliJetReaderHeader.h"
46#include "AliUA1JetHeaderV1.h"
47#include "AliESDEvent.h"
48#include "AliAODEvent.h"
49#include "AliAODHandler.h"
50#include "AliAODTrack.h"
51#include "AliAODJet.h"
52#include "AliAODMCParticle.h"
53#include "AliMCEventHandler.h"
54#include "AliMCEvent.h"
55#include "AliStack.h"
56#include "AliGenPythiaEventHeader.h"
57#include "AliJetKineReaderHeader.h"
58#include "AliGenCocktailEventHeader.h"
59#include "AliInputEventHandler.h"
60
61
62#include "fastjet/PseudoJet.hh"
63#include "fastjet/ClusterSequenceArea.hh"
64#include "fastjet/AreaDefinition.hh"
65#include "fastjet/JetDefinition.hh"
66// get info on how fastjet was configured
67#include "fastjet/config.h"
68
69
70ClassImp(AliAnalysisTaskJetCluster)
71
72AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
73 fAOD(0x0),
74 fUseAODTrackInput(kFALSE),
75 fUseAODMCInput(kFALSE),
76 fUseGlobalSelection(kFALSE),
77 fFilterMask(0),
78 fTrackTypeRec(kTrackUndef),
79 fTrackTypeGen(kTrackUndef),
80 fAvgTrials(1),
81 fExternalWeight(1),
82 fRecEtaWindow(0.5),
83 fRparam(1.0),
84 fAlgorithm(fastjet::kt_algorithm),
85 fStrategy(fastjet::Best),
86 fRecombScheme(fastjet::BIpt_scheme),
87 fAreaType(fastjet::active_area),
88 fh1Xsec(0x0),
89 fh1Trials(0x0),
90 fh1PtHard(0x0),
91 fh1PtHardNoW(0x0),
92 fh1PtHardTrials(0x0),
93 fh1NJetsRec(0x0),
94 fh1NConstRec(0x0),
95 fh1NConstLeadingRec(0x0),
96 fh1PtJetsRecIn(0x0),
97 fh1PtJetsLeadingRecIn(0x0),
98 fh1PtJetConstRec(0x0),
99 fh1PtJetConstLeadingRec(0x0),
100 fh1PtTracksRecIn(0x0),
101 fh1PtTracksLeadingRecIn(0x0),
102 fh1NJetsRecRan(0x0),
103 fh1NConstRecRan(0x0),
104 fh1PtJetsLeadingRecInRan(0x0),
105 fh1NConstLeadingRecRan(0x0),
106 fh1PtJetsRecInRan(0x0),
107 fh1PtTracksGenIn(0x0),
a0d08b6d 108 fh1Nch(0x0),
dbf053f0 109 fh2NRecJetsPt(0x0),
110 fh2NRecTracksPt(0x0),
111 fh2NConstPt(0x0),
112 fh2NConstLeadingPt(0x0),
113 fh2JetPhiEta(0x0),
114 fh2LeadingJetPhiEta(0x0),
115 fh2JetEtaPt(0x0),
116 fh2LeadingJetEtaPt(0x0),
117 fh2TrackEtaPt(0x0),
118 fh2LeadingTrackEtaPt(0x0),
119 fh2JetsLeadingPhiEta(0x0),
120 fh2JetsLeadingPhiPt(0x0),
121 fh2TracksLeadingPhiEta(0x0),
122 fh2TracksLeadingPhiPt(0x0),
123 fh2TracksLeadingJetPhiPt(0x0),
124 fh2JetsLeadingPhiPtW(0x0),
125 fh2TracksLeadingPhiPtW(0x0),
126 fh2TracksLeadingJetPhiPtW(0x0),
127 fh2NRecJetsPtRan(0x0),
128 fh2NConstPtRan(0x0),
129 fh2NConstLeadingPtRan(0x0),
a0d08b6d 130 fh2PtNch(0x0),
131 fh2PtNchRan(0x0),
132 fh2PtNchN(0x0),
133 fh2PtNchNRan(0x0),
dbf053f0 134 fh2TracksLeadingJetPhiPtRan(0x0),
135 fh2TracksLeadingJetPhiPtWRan(0x0),
136 fHistList(0x0)
137{
138
139}
140
141AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
142 AliAnalysisTaskSE(name),
143 fAOD(0x0),
144 fUseAODTrackInput(kFALSE),
145 fUseAODMCInput(kFALSE),
146 fUseGlobalSelection(kFALSE),
147 fFilterMask(0),
148 fTrackTypeRec(kTrackUndef),
149 fTrackTypeGen(kTrackUndef),
150 fAvgTrials(1),
151 fExternalWeight(1),
152 fRecEtaWindow(0.5),
153 fRparam(1.0),
154 fAlgorithm(fastjet::kt_algorithm),
155 fStrategy(fastjet::Best),
156 fRecombScheme(fastjet::BIpt_scheme),
157 fAreaType(fastjet::active_area),
158 fh1Xsec(0x0),
159 fh1Trials(0x0),
160 fh1PtHard(0x0),
161 fh1PtHardNoW(0x0),
162 fh1PtHardTrials(0x0),
163 fh1NJetsRec(0x0),
164 fh1NConstRec(0x0),
165 fh1NConstLeadingRec(0x0),
166 fh1PtJetsRecIn(0x0),
167 fh1PtJetsLeadingRecIn(0x0),
168 fh1PtJetConstRec(0x0),
169 fh1PtJetConstLeadingRec(0x0),
170 fh1PtTracksRecIn(0x0),
171 fh1PtTracksLeadingRecIn(0x0),
172 fh1NJetsRecRan(0x0),
173 fh1NConstRecRan(0x0),
174 fh1PtJetsLeadingRecInRan(0x0),
175 fh1NConstLeadingRecRan(0x0),
176 fh1PtJetsRecInRan(0x0),
177 fh1PtTracksGenIn(0x0),
a0d08b6d 178 fh1Nch(0x0),
dbf053f0 179 fh2NRecJetsPt(0x0),
180 fh2NRecTracksPt(0x0),
181 fh2NConstPt(0x0),
182 fh2NConstLeadingPt(0x0),
183 fh2JetPhiEta(0x0),
184 fh2LeadingJetPhiEta(0x0),
185 fh2JetEtaPt(0x0),
186 fh2LeadingJetEtaPt(0x0),
187 fh2TrackEtaPt(0x0),
188 fh2LeadingTrackEtaPt(0x0),
189 fh2JetsLeadingPhiEta(0x0),
190 fh2JetsLeadingPhiPt(0x0),
191 fh2TracksLeadingPhiEta(0x0),
192 fh2TracksLeadingPhiPt(0x0),
193 fh2TracksLeadingJetPhiPt(0x0),
194 fh2JetsLeadingPhiPtW(0x0),
195 fh2TracksLeadingPhiPtW(0x0),
196 fh2TracksLeadingJetPhiPtW(0x0),
197 fh2NRecJetsPtRan(0x0),
198 fh2NConstPtRan(0x0),
199 fh2NConstLeadingPtRan(0x0),
a0d08b6d 200 fh2PtNch(0x0),
201 fh2PtNchRan(0x0),
202 fh2PtNchN(0x0),
203 fh2PtNchNRan(0x0),
dbf053f0 204 fh2TracksLeadingJetPhiPtRan(0x0),
205 fh2TracksLeadingJetPhiPtWRan(0x0),
206 fHistList(0x0)
207{
208 DefineOutput(1, TList::Class());
209}
210
211
212
213Bool_t AliAnalysisTaskJetCluster::Notify()
214{
215 //
216 // Implemented Notify() to read the cross sections
217 // and number of trials from pyxsec.root
218 //
219 return kTRUE;
220}
221
222void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
223{
224
225 //
226 // Create the output container
227 //
228
229
230 // Connect the AOD
231
232
233 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
234
235 OpenFile(1);
236 if(!fHistList)fHistList = new TList();
237
238 Bool_t oldStatus = TH1::AddDirectoryStatus();
239 TH1::AddDirectory(kFALSE);
240
241 //
242 // Histogram
243
244 const Int_t nBinPt = 200;
245 Double_t binLimitsPt[nBinPt+1];
246 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
247 if(iPt == 0){
248 binLimitsPt[iPt] = 0.0;
249 }
250 else {// 1.0
251 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 0.5;
252 }
253 }
254
255 const Int_t nBinPhi = 90;
256 Double_t binLimitsPhi[nBinPhi+1];
257 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
258 if(iPhi==0){
259 binLimitsPhi[iPhi] = -1.*TMath::Pi();
260 }
261 else{
262 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
263 }
264 }
265
266
267
268 const Int_t nBinEta = 40;
269 Double_t binLimitsEta[nBinEta+1];
270 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
271 if(iEta==0){
272 binLimitsEta[iEta] = -2.0;
273 }
274 else{
275 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
276 }
277 }
278
a0d08b6d 279 const int nChMax = 100;
dbf053f0 280
281 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
282 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
283
284 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
285 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
286
287
288 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
289 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
290
291 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
292 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
293 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
294 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
295
296
297 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
298 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
299 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
300
301 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
302 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
303 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
304 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
305 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
306 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
307 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
308 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
309 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
a0d08b6d 310 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
dbf053f0 311
312 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
313 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
314 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
315 //
316
317
318 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
319 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
320 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
321 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
322
a0d08b6d 323 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
324 fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
325 fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
326 fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
327
328
329
dbf053f0 330 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
331 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
332 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
333 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
334
335 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
336 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
337 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
338 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
339
340 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
341 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
342 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
343 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
344
345
346
347 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
348 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
349 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
350 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
351 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
352 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
353 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
354 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
355 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
356 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
357 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
358 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
359
360 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
361 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
362 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
363 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
364
365 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
366 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
367 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
368 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
369
370
371
372 const Int_t saveLevel = 3; // large save level more histos
373 if(saveLevel>0){
374 fHistList->Add(fh1Xsec);
375 fHistList->Add(fh1Trials);
376
377 fHistList->Add(fh1NJetsRec);
378 fHistList->Add(fh1NConstRec);
379 fHistList->Add(fh1NConstLeadingRec);
380 fHistList->Add(fh1PtJetsRecIn);
381 fHistList->Add(fh1PtJetsLeadingRecIn);
382 fHistList->Add(fh1PtTracksRecIn);
383 fHistList->Add(fh1PtTracksLeadingRecIn);
384 fHistList->Add(fh1PtJetConstRec);
385 fHistList->Add(fh1PtJetConstLeadingRec);
386 fHistList->Add(fh1NJetsRecRan);
387 fHistList->Add(fh1NConstRecRan);
388 fHistList->Add(fh1PtJetsLeadingRecInRan);
389 fHistList->Add(fh1NConstLeadingRecRan);
390 fHistList->Add(fh1PtJetsRecInRan);
a0d08b6d 391 fHistList->Add(fh1Nch);
dbf053f0 392 fHistList->Add(fh2NRecJetsPt);
393 fHistList->Add(fh2NRecTracksPt);
394 fHistList->Add(fh2NConstPt);
395 fHistList->Add(fh2NConstLeadingPt);
a0d08b6d 396 fHistList->Add(fh2PtNch);
397 fHistList->Add(fh2PtNchRan);
398 fHistList->Add(fh2PtNchN);
399 fHistList->Add(fh2PtNchNRan);
dbf053f0 400 fHistList->Add(fh2JetPhiEta);
401 fHistList->Add(fh2LeadingJetPhiEta);
402 fHistList->Add(fh2JetEtaPt);
403 fHistList->Add(fh2LeadingJetEtaPt);
404 fHistList->Add(fh2TrackEtaPt);
405 fHistList->Add(fh2LeadingTrackEtaPt);
406 fHistList->Add(fh2JetsLeadingPhiEta );
407 fHistList->Add(fh2JetsLeadingPhiPt);
408 fHistList->Add(fh2TracksLeadingPhiEta);
409 fHistList->Add(fh2TracksLeadingPhiPt);
410 fHistList->Add(fh2TracksLeadingJetPhiPt);
411 fHistList->Add(fh2JetsLeadingPhiPtW);
412 fHistList->Add(fh2TracksLeadingPhiPtW);
413 fHistList->Add(fh2TracksLeadingJetPhiPtW);
414 fHistList->Add(fh2NRecJetsPtRan);
415 fHistList->Add(fh2NConstPtRan);
416 fHistList->Add(fh2NConstLeadingPtRan);
417 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
418 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
419 }
420
421 // =========== Switch on Sumw2 for all histos ===========
422 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
423 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
424 if (h1){
425 h1->Sumw2();
426 continue;
427 }
428 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
429 if(hn)hn->Sumw2();
430 }
431 TH1::AddDirectory(oldStatus);
432}
433
434void AliAnalysisTaskJetCluster::Init()
435{
436 //
437 // Initialization
438 //
439
dbf053f0 440 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
441
442}
443
444void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
445{
446
447 if(fUseGlobalSelection){
448 // no selection by the service task, we continue
449 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
450 PostData(1, fHistList);
451 return;
452 }
453
454 //
455 // Execute analysis for current event
456 //
457 AliESDEvent *fESD = 0;
458 if(fUseAODTrackInput){
459 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
460 if(!fAOD){
461 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
462 return;
463 }
464 // fethc the header
465 }
466 else{
467 // assume that the AOD is in the general output...
468 fAOD = AODEvent();
469 if(!fAOD){
470 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
471 return;
472 }
473 if(fDebug>0){
474 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
475 }
476 }
477
478 Bool_t selectEvent = false;
479 Bool_t physicsSelection = true; // fInputHandler->IsEventSelected();
480
481 if(fAOD){
482 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
483 if(vtxAOD->GetNContributors()>0){
484 Float_t zvtx = vtxAOD->GetZ();
485 Float_t yvtx = vtxAOD->GetY();
486 Float_t xvtx = vtxAOD->GetX();
487 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
488 if(TMath::Abs(zvtx)<8.&&r2<1.){
489 if(physicsSelection)selectEvent = true;
490 }
491 }
492 }
493 if(!selectEvent){
494 PostData(1, fHistList);
495 return;
496 }
497
498 fh1Trials->Fill("#sum{ntrials}",1);
499
500
501 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
502
503 // ==== General variables needed
504
505
506
507 // we simply fetch the tracks/mc particles as a list of AliVParticles
508
509 TList recParticles;
510 TList genParticles;
511
512 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
a0d08b6d 513 Float_t nCh = recParticles.GetEntries();
514 fh1Nch->Fill(nCh);
dbf053f0 515 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
516 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
517 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
518
519 // find the jets....
520
521 vector<fastjet::PseudoJet> inputParticlesRec;
522 vector<fastjet::PseudoJet> inputParticlesRecRan;
523 for(int i = 0; i < recParticles.GetEntries(); i++){
524 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
525 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
526 // we talk total momentum here
527 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
528 jInp.set_user_index(i);
529 inputParticlesRec.push_back(jInp);
530
531 // the randomized input changes eta and phi, but keeps the p_T
532 Double_t pT = vp->Pt();
533 Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
534 Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
535
536
537 Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));
538 Double_t pZ = pT/TMath::Tan(theta);
539
540 Double_t pX = pT * TMath::Cos(phi);
541 Double_t pY = pT * TMath::Sin(phi);
542 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
543
544 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
545 jInpRan.set_user_index(i);
546 inputParticlesRecRan.push_back(jInpRan);
547
548 }
549
550 if(inputParticlesRec.size()==0){
551 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
552 PostData(1, fHistList);
553 return;
554 }
555
556 // run fast jet
557 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
558 vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
559 fastjet::ClusterSequence clustSeq(inputParticlesRec, jetDef);
560
561 inclusiveJets = clustSeq.inclusive_jets();
562 sortedJets = sorted_by_pt(inclusiveJets);
563
564 if(fDebug)Printf("%s:%d Found %d jets",(char*)__FILE__,__LINE__, inclusiveJets.size());
565
566 fh1NJetsRec->Fill(sortedJets.size());
567
568 // loop over all jets an fill information, first on is the leading jet
569
570 Int_t nRecOver = inclusiveJets.size();
571 Int_t nRec = inclusiveJets.size();
572 if(inclusiveJets.size()>0){
573 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
574 Float_t pt = leadingJet.Pt();
575
576 Int_t iCount = 0;
577 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
578 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
579 while(pt<ptCut&&iCount<nRec){
580 nRecOver--;
581 iCount++;
582 if(iCount<nRec){
583 pt = sortedJets[iCount].perp();
584 }
585 }
586 if(nRecOver<=0)break;
587 fh2NRecJetsPt->Fill(ptCut,nRecOver);
588 }
589 Float_t phi = leadingJet.Phi();
590 if(phi<0)phi+=TMath::Pi()*2.;
591 Float_t eta = leadingJet.Eta();
592 pt = leadingJet.Pt();
593
594 // correlation of leading jet with tracks
595 TIterator *recIter = recParticles.MakeIterator();
596 AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());
597
598 recIter->Reset();
599 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
600 Float_t tmpPt = tmpRecTrack->Pt();
601 // correlation
602 Float_t tmpPhi = tmpRecTrack->Phi();
603 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
604 Float_t dPhi = phi - tmpPhi;
605 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
606 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
607 // Float_t dEta = eta - tmpRecTrack->Eta();
608 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
609 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
610 }
611
612
613
614 for(int j = 0; j < nRec;j++){
615 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
616 Float_t tmpPt = tmpRec.Pt();
617 fh1PtJetsRecIn->Fill(tmpPt);
618 // Fill Spectra with constituents
619 vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
620 fh1NConstRec->Fill(constituents.size());
a0d08b6d 621 fh2PtNch->Fill(nCh,tmpPt);
622 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
dbf053f0 623 fh2NConstPt->Fill(tmpPt,constituents.size());
624 // loop over constiutents and fill spectrum
625 for(unsigned int ic = 0; ic < constituents.size();ic++){
626 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
627 fh1PtJetConstRec->Fill(part->Pt());
628 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
629 }
630
631 // correlation
632 Float_t tmpPhi = tmpRec.Phi();
633 Float_t tmpEta = tmpRec.Eta();
634 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
635
636 if(j==0){
637 fh1PtJetsLeadingRecIn->Fill(tmpPt);
638 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
639 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
640 fh1NConstLeadingRec->Fill(constituents.size());
641 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
642 continue;
643 }
644 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
645 fh2JetEtaPt->Fill(tmpEta,tmpPt);
646 Float_t dPhi = phi - tmpPhi;
647 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
648 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
649 Float_t dEta = eta - tmpRec.Eta();
650 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
651 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
652 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
94190d53 653 }
654 delete recIter;
dbf053f0 655 }
656
657
658 // fill track information
659 Int_t nTrackOver = recParticles.GetSize();
660 // do the same for tracks and jets
661 if(nTrackOver>0){
662 TIterator *recIter = recParticles.MakeIterator();
663 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
664 Float_t pt = tmpRec->Pt();
665 // Printf("Leading track p_t %3.3E",pt);
666 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
667 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
668 while(pt<ptCut&&tmpRec){
669 nTrackOver--;
670 tmpRec = (AliVParticle*)(recIter->Next());
671 if(tmpRec){
672 pt = tmpRec->Pt();
673 }
674 }
675 if(nTrackOver<=0)break;
676 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
677 }
678
679 recIter->Reset();
680 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
681 Float_t phi = leading->Phi();
682 if(phi<0)phi+=TMath::Pi()*2.;
683 Float_t eta = leading->Eta();
684 pt = leading->Pt();
685 while((tmpRec = (AliVParticle*)(recIter->Next()))){
686 Float_t tmpPt = tmpRec->Pt();
687 Float_t tmpEta = tmpRec->Eta();
688 fh1PtTracksRecIn->Fill(tmpPt);
689 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
690 if(tmpRec==leading){
691 fh1PtTracksLeadingRecIn->Fill(tmpPt);
692 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
693 continue;
694 }
695 // correlation
696 Float_t tmpPhi = tmpRec->Phi();
697
698 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
699 Float_t dPhi = phi - tmpPhi;
700 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
701 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
702 Float_t dEta = eta - tmpRec->Eta();
703 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
704 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
705 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
706 }
707 delete recIter;
708 }
709
710 // find the random jets
711 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
712 fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
713
714 inclusiveJetsRan = clustSeqRan.inclusive_jets();
715 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
716
717 // fill the jet information from random track
718
719 fh1NJetsRecRan->Fill(sortedJetsRan.size());
720
721 // loop over all jets an fill information, first on is the leading jet
722
723 Int_t nRecOverRan = inclusiveJetsRan.size();
724 Int_t nRecRan = inclusiveJetsRan.size();
725 if(inclusiveJetsRan.size()>0){
726 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
727 Float_t pt = leadingJet.Pt();
728
729 Int_t iCount = 0;
730 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
731 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
732 while(pt<ptCut&&iCount<nRecRan){
733 nRecOverRan--;
734 iCount++;
735 if(iCount<nRecRan){
736 pt = sortedJetsRan[iCount].perp();
737 }
738 }
739 if(nRecOverRan<=0)break;
740 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
741 }
742 Float_t phi = leadingJet.Phi();
743 if(phi<0)phi+=TMath::Pi()*2.;
744 pt = leadingJet.Pt();
745
746 // correlation of leading jet with random tracks
747
748 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
749 {
750 Float_t tmpPt = inputParticlesRecRan[ip].perp();
751 // correlation
752 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
753 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
754 Float_t dPhi = phi - tmpPhi;
755 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
756 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
757 // Float_t dEta = eta - tmpRecTrack->Eta();
758 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
759 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
760 }
761
762
763
764 for(int j = 0; j < nRecRan;j++){
765 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
766 Float_t tmpPt = tmpRec.Pt();
767 fh1PtJetsRecInRan->Fill(tmpPt);
768 // Fill Spectra with constituents
769 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
770 fh1NConstRecRan->Fill(constituents.size());
771 fh2NConstPtRan->Fill(tmpPt,constituents.size());
a0d08b6d 772 fh2PtNchRan->Fill(nCh,tmpPt);
773 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
dbf053f0 774 // correlation
775 Float_t tmpPhi = tmpRec.Phi();
776 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
777
778 if(j==0){
779 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
780 fh1NConstLeadingRecRan->Fill(constituents.size());
781 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
782 continue;
783 }
784 }
785 }
786
787
788 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
789 PostData(1, fHistList);
790}
791
792void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
793{
794// Terminate analysis
795//
796 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
797}
798
799
800Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
801
802 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
803
804 Int_t iCount = 0;
805 if(type==kTrackAOD){
806 AliAODEvent *aod = 0;
807 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
808 else aod = AODEvent();
809 if(!aod){
810 return iCount;
811 }
812 for(int it = 0;it < aod->GetNumberOfTracks();++it){
813 AliAODTrack *tr = aod->GetTrack(it);
814 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
815 if(TMath::Abs(tr->Eta())>0.9)continue;
816 // if(tr->Pt()<0.3)continue;
817 list->Add(tr);
818 iCount++;
819 }
820 }
821 else if (type == kTrackKineAll||type == kTrackKineCharged){
822 AliMCEvent* mcEvent = MCEvent();
823 if(!mcEvent)return iCount;
824 // we want to have alivpartilces so use get track
825 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
826 if(!mcEvent->IsPhysicalPrimary(it))continue;
827 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
828 if(type == kTrackKineAll){
829 list->Add(part);
830 iCount++;
831 }
832 else if(type == kTrackKineCharged){
833 if(part->Particle()->GetPDG()->Charge()==0)continue;
834 list->Add(part);
835 iCount++;
836 }
837 }
838 }
839 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
840 AliAODEvent *aod = 0;
841 if(fUseAODMCInput)dynamic_cast<AliAODEvent*>(InputEvent());
842 else aod = AODEvent();
843 if(!aod)return iCount;
844 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
845 if(!tca)return iCount;
846 for(int it = 0;it < tca->GetEntriesFast();++it){
847 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
848 if(!part->IsPhysicalPrimary())continue;
849 if(type == kTrackAODMCAll){
850 list->Add(part);
851 iCount++;
852 }
853 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
854 if(part->Charge()==0)continue;
855 if(kTrackAODMCCharged){
856 list->Add(part);
857 }
858 else {
859 if(TMath::Abs(part->Eta())>0.9)continue;
860 list->Add(part);
861 }
862 iCount++;
863 }
864 }
865 }// AODMCparticle
866 list->Sort();
867 return iCount;
868
869}
870
871/*
872Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
873 for(int i = 0; i < particles.GetEntries(); i++){
874 AliVParticle *vp = (AliVParticle*)particles.At(i);
875 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
876 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
877 jInp.set_user_index(i);
878 inputParticles.push_back(jInp);
879 }
880
881 return 0;
882
883}
884*/