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