]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - JETAN/AliAnalysisTaskJetCluster.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fixes for event selection
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetCluster.cxx
... / ...
CommitLineData
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),
108 fh1Nch(0x0),
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),
130 fh2PtNch(0x0),
131 fh2PtNchRan(0x0),
132 fh2PtNchN(0x0),
133 fh2PtNchNRan(0x0),
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),
178 fh1Nch(0x0),
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),
200 fh2PtNch(0x0),
201 fh2PtNchRan(0x0),
202 fh2PtNchN(0x0),
203 fh2PtNchNRan(0x0),
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
279 const int nChMax = 100;
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);
310 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
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
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
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);
391 fHistList->Add(fh1Nch);
392 fHistList->Add(fh2NRecJetsPt);
393 fHistList->Add(fh2NRecTracksPt);
394 fHistList->Add(fh2NConstPt);
395 fHistList->Add(fh2NConstLeadingPt);
396 fHistList->Add(fh2PtNch);
397 fHistList->Add(fh2PtNchRan);
398 fHistList->Add(fh2PtNchN);
399 fHistList->Add(fh2PtNchNRan);
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
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 = (fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
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);
513 Float_t nCh = recParticles.GetEntries();
514 fh1Nch->Fill(nCh);
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 fh1NJetsRec->Fill(sortedJets.size());
565
566 // loop over all jets an fill information, first on is the leading jet
567
568 Int_t nRecOver = inclusiveJets.size();
569 Int_t nRec = inclusiveJets.size();
570 if(inclusiveJets.size()>0){
571 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
572 Float_t pt = leadingJet.Pt();
573
574 Int_t iCount = 0;
575 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
576 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
577 while(pt<ptCut&&iCount<nRec){
578 nRecOver--;
579 iCount++;
580 if(iCount<nRec){
581 pt = sortedJets[iCount].perp();
582 }
583 }
584 if(nRecOver<=0)break;
585 fh2NRecJetsPt->Fill(ptCut,nRecOver);
586 }
587 Float_t phi = leadingJet.Phi();
588 if(phi<0)phi+=TMath::Pi()*2.;
589 Float_t eta = leadingJet.Eta();
590 pt = leadingJet.Pt();
591
592 // correlation of leading jet with tracks
593 TIterator *recIter = recParticles.MakeIterator();
594 AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());
595
596 recIter->Reset();
597 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
598 Float_t tmpPt = tmpRecTrack->Pt();
599 // correlation
600 Float_t tmpPhi = tmpRecTrack->Phi();
601 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
602 Float_t dPhi = phi - tmpPhi;
603 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
604 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
605 // Float_t dEta = eta - tmpRecTrack->Eta();
606 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
607 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
608 }
609
610
611
612 for(int j = 0; j < nRec;j++){
613 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
614 Float_t tmpPt = tmpRec.Pt();
615 fh1PtJetsRecIn->Fill(tmpPt);
616 // Fill Spectra with constituents
617 vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
618 fh1NConstRec->Fill(constituents.size());
619 fh2PtNch->Fill(nCh,tmpPt);
620 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
621 fh2NConstPt->Fill(tmpPt,constituents.size());
622 // loop over constiutents and fill spectrum
623 for(unsigned int ic = 0; ic < constituents.size();ic++){
624 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
625 fh1PtJetConstRec->Fill(part->Pt());
626 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
627 }
628
629 // correlation
630 Float_t tmpPhi = tmpRec.Phi();
631 Float_t tmpEta = tmpRec.Eta();
632 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
633
634 if(j==0){
635 fh1PtJetsLeadingRecIn->Fill(tmpPt);
636 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
637 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
638 fh1NConstLeadingRec->Fill(constituents.size());
639 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
640 continue;
641 }
642 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
643 fh2JetEtaPt->Fill(tmpEta,tmpPt);
644 Float_t dPhi = phi - tmpPhi;
645 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
646 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
647 Float_t dEta = eta - tmpRec.Eta();
648 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
649 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
650 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
651 }
652 delete recIter;
653 }
654
655
656 // fill track information
657 Int_t nTrackOver = recParticles.GetSize();
658 // do the same for tracks and jets
659 if(nTrackOver>0){
660 TIterator *recIter = recParticles.MakeIterator();
661 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
662 Float_t pt = tmpRec->Pt();
663 // Printf("Leading track p_t %3.3E",pt);
664 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
665 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
666 while(pt<ptCut&&tmpRec){
667 nTrackOver--;
668 tmpRec = (AliVParticle*)(recIter->Next());
669 if(tmpRec){
670 pt = tmpRec->Pt();
671 }
672 }
673 if(nTrackOver<=0)break;
674 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
675 }
676
677 recIter->Reset();
678 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
679 Float_t phi = leading->Phi();
680 if(phi<0)phi+=TMath::Pi()*2.;
681 Float_t eta = leading->Eta();
682 pt = leading->Pt();
683 while((tmpRec = (AliVParticle*)(recIter->Next()))){
684 Float_t tmpPt = tmpRec->Pt();
685 Float_t tmpEta = tmpRec->Eta();
686 fh1PtTracksRecIn->Fill(tmpPt);
687 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
688 if(tmpRec==leading){
689 fh1PtTracksLeadingRecIn->Fill(tmpPt);
690 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
691 continue;
692 }
693 // correlation
694 Float_t tmpPhi = tmpRec->Phi();
695
696 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
697 Float_t dPhi = phi - tmpPhi;
698 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
699 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
700 Float_t dEta = eta - tmpRec->Eta();
701 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
702 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
703 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
704 }
705 delete recIter;
706 }
707
708 // find the random jets
709 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
710 fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
711
712 inclusiveJetsRan = clustSeqRan.inclusive_jets();
713 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
714
715 // fill the jet information from random track
716
717 fh1NJetsRecRan->Fill(sortedJetsRan.size());
718
719 // loop over all jets an fill information, first on is the leading jet
720
721 Int_t nRecOverRan = inclusiveJetsRan.size();
722 Int_t nRecRan = inclusiveJetsRan.size();
723 if(inclusiveJetsRan.size()>0){
724 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
725 Float_t pt = leadingJet.Pt();
726
727 Int_t iCount = 0;
728 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
729 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
730 while(pt<ptCut&&iCount<nRecRan){
731 nRecOverRan--;
732 iCount++;
733 if(iCount<nRecRan){
734 pt = sortedJetsRan[iCount].perp();
735 }
736 }
737 if(nRecOverRan<=0)break;
738 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
739 }
740 Float_t phi = leadingJet.Phi();
741 if(phi<0)phi+=TMath::Pi()*2.;
742 pt = leadingJet.Pt();
743
744 // correlation of leading jet with random tracks
745
746 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
747 {
748 Float_t tmpPt = inputParticlesRecRan[ip].perp();
749 // correlation
750 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
751 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
752 Float_t dPhi = phi - tmpPhi;
753 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
754 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
755 // Float_t dEta = eta - tmpRecTrack->Eta();
756 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
757 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
758 }
759
760
761
762 for(int j = 0; j < nRecRan;j++){
763 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
764 Float_t tmpPt = tmpRec.Pt();
765 fh1PtJetsRecInRan->Fill(tmpPt);
766 // Fill Spectra with constituents
767 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
768 fh1NConstRecRan->Fill(constituents.size());
769 fh2NConstPtRan->Fill(tmpPt,constituents.size());
770 fh2PtNchRan->Fill(nCh,tmpPt);
771 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
772 // correlation
773 Float_t tmpPhi = tmpRec.Phi();
774 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
775
776 if(j==0){
777 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
778 fh1NConstLeadingRecRan->Fill(constituents.size());
779 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
780 continue;
781 }
782 }
783 }
784
785
786 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
787 PostData(1, fHistList);
788}
789
790void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
791{
792// Terminate analysis
793//
794 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
795}
796
797
798Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
799
800 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
801
802 Int_t iCount = 0;
803 if(type==kTrackAOD){
804 AliAODEvent *aod = 0;
805 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
806 else aod = AODEvent();
807 if(!aod){
808 return iCount;
809 }
810 for(int it = 0;it < aod->GetNumberOfTracks();++it){
811 AliAODTrack *tr = aod->GetTrack(it);
812 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
813 if(TMath::Abs(tr->Eta())>0.9)continue;
814 // if(tr->Pt()<0.3)continue;
815 list->Add(tr);
816 iCount++;
817 }
818 }
819 else if (type == kTrackKineAll||type == kTrackKineCharged){
820 AliMCEvent* mcEvent = MCEvent();
821 if(!mcEvent)return iCount;
822 // we want to have alivpartilces so use get track
823 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
824 if(!mcEvent->IsPhysicalPrimary(it))continue;
825 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
826 if(type == kTrackKineAll){
827 list->Add(part);
828 iCount++;
829 }
830 else if(type == kTrackKineCharged){
831 if(part->Particle()->GetPDG()->Charge()==0)continue;
832 list->Add(part);
833 iCount++;
834 }
835 }
836 }
837 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
838 AliAODEvent *aod = 0;
839 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
840 else aod = AODEvent();
841 if(!aod)return iCount;
842 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
843 if(!tca)return iCount;
844 for(int it = 0;it < tca->GetEntriesFast();++it){
845 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
846 if(!part->IsPhysicalPrimary())continue;
847 if(type == kTrackAODMCAll){
848 list->Add(part);
849 iCount++;
850 }
851 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
852 if(part->Charge()==0)continue;
853 if(kTrackAODMCCharged){
854 list->Add(part);
855 }
856 else {
857 if(TMath::Abs(part->Eta())>0.9)continue;
858 list->Add(part);
859 }
860 iCount++;
861 }
862 }
863 }// AODMCparticle
864 list->Sort();
865 return iCount;
866
867}
868
869/*
870Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
871 for(int i = 0; i < particles.GetEntries(); i++){
872 AliVParticle *vp = (AliVParticle*)particles.At(i);
873 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
874 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
875 jInp.set_user_index(i);
876 inputParticles.push_back(jInp);
877 }
878
879 return 0;
880
881}
882*/