]> 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
Fixing T0 compilation
[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
dbf053f0 416 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
417
418}
419
420void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
421{
422
423 if(fUseGlobalSelection){
424 // no selection by the service task, we continue
425 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
426 PostData(1, fHistList);
427 return;
428 }
429
430 //
431 // Execute analysis for current event
432 //
433 AliESDEvent *fESD = 0;
434 if(fUseAODTrackInput){
435 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
436 if(!fAOD){
437 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
438 return;
439 }
440 // fethc the header
441 }
442 else{
443 // assume that the AOD is in the general output...
444 fAOD = AODEvent();
445 if(!fAOD){
446 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
447 return;
448 }
449 if(fDebug>0){
450 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
451 }
452 }
453
454 Bool_t selectEvent = false;
455 Bool_t physicsSelection = true; // fInputHandler->IsEventSelected();
456
457 if(fAOD){
458 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
459 if(vtxAOD->GetNContributors()>0){
460 Float_t zvtx = vtxAOD->GetZ();
461 Float_t yvtx = vtxAOD->GetY();
462 Float_t xvtx = vtxAOD->GetX();
463 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
464 if(TMath::Abs(zvtx)<8.&&r2<1.){
465 if(physicsSelection)selectEvent = true;
466 }
467 }
468 }
469 if(!selectEvent){
470 PostData(1, fHistList);
471 return;
472 }
473
474 fh1Trials->Fill("#sum{ntrials}",1);
475
476
477 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
478
479 // ==== General variables needed
480
481
482
483 // we simply fetch the tracks/mc particles as a list of AliVParticles
484
485 TList recParticles;
486 TList genParticles;
487
488 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
489 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
490 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
491 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
492
493 // find the jets....
494
495 vector<fastjet::PseudoJet> inputParticlesRec;
496 vector<fastjet::PseudoJet> inputParticlesRecRan;
497 for(int i = 0; i < recParticles.GetEntries(); i++){
498 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
499 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
500 // we talk total momentum here
501 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
502 jInp.set_user_index(i);
503 inputParticlesRec.push_back(jInp);
504
505 // the randomized input changes eta and phi, but keeps the p_T
506 Double_t pT = vp->Pt();
507 Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
508 Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
509
510
511 Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));
512 Double_t pZ = pT/TMath::Tan(theta);
513
514 Double_t pX = pT * TMath::Cos(phi);
515 Double_t pY = pT * TMath::Sin(phi);
516 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
517
518 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
519 jInpRan.set_user_index(i);
520 inputParticlesRecRan.push_back(jInpRan);
521
522 }
523
524 if(inputParticlesRec.size()==0){
525 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
526 PostData(1, fHistList);
527 return;
528 }
529
530 // run fast jet
531 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
532 vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
533 fastjet::ClusterSequence clustSeq(inputParticlesRec, jetDef);
534
535 inclusiveJets = clustSeq.inclusive_jets();
536 sortedJets = sorted_by_pt(inclusiveJets);
537
538 if(fDebug)Printf("%s:%d Found %d jets",(char*)__FILE__,__LINE__, inclusiveJets.size());
539
540 fh1NJetsRec->Fill(sortedJets.size());
541
542 // loop over all jets an fill information, first on is the leading jet
543
544 Int_t nRecOver = inclusiveJets.size();
545 Int_t nRec = inclusiveJets.size();
546 if(inclusiveJets.size()>0){
547 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
548 Float_t pt = leadingJet.Pt();
549
550 Int_t iCount = 0;
551 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
552 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
553 while(pt<ptCut&&iCount<nRec){
554 nRecOver--;
555 iCount++;
556 if(iCount<nRec){
557 pt = sortedJets[iCount].perp();
558 }
559 }
560 if(nRecOver<=0)break;
561 fh2NRecJetsPt->Fill(ptCut,nRecOver);
562 }
563 Float_t phi = leadingJet.Phi();
564 if(phi<0)phi+=TMath::Pi()*2.;
565 Float_t eta = leadingJet.Eta();
566 pt = leadingJet.Pt();
567
568 // correlation of leading jet with tracks
569 TIterator *recIter = recParticles.MakeIterator();
570 AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());
571
572 recIter->Reset();
573 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
574 Float_t tmpPt = tmpRecTrack->Pt();
575 // correlation
576 Float_t tmpPhi = tmpRecTrack->Phi();
577 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
578 Float_t dPhi = phi - tmpPhi;
579 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
580 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
581 // Float_t dEta = eta - tmpRecTrack->Eta();
582 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
583 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
584 }
585
586
587
588 for(int j = 0; j < nRec;j++){
589 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
590 Float_t tmpPt = tmpRec.Pt();
591 fh1PtJetsRecIn->Fill(tmpPt);
592 // Fill Spectra with constituents
593 vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
594 fh1NConstRec->Fill(constituents.size());
595 fh2NConstPt->Fill(tmpPt,constituents.size());
596 // loop over constiutents and fill spectrum
597 for(unsigned int ic = 0; ic < constituents.size();ic++){
598 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
599 fh1PtJetConstRec->Fill(part->Pt());
600 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
601 }
602
603 // correlation
604 Float_t tmpPhi = tmpRec.Phi();
605 Float_t tmpEta = tmpRec.Eta();
606 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
607
608 if(j==0){
609 fh1PtJetsLeadingRecIn->Fill(tmpPt);
610 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
611 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
612 fh1NConstLeadingRec->Fill(constituents.size());
613 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
614 continue;
615 }
616 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
617 fh2JetEtaPt->Fill(tmpEta,tmpPt);
618 Float_t dPhi = phi - tmpPhi;
619 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
620 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
621 Float_t dEta = eta - tmpRec.Eta();
622 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
623 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
624 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
94190d53 625 }
626 delete recIter;
dbf053f0 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*/