Fixed minor mem leak in the d'tor and added protection
[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 <TRefArray.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 "AliESDEvent.h"
46#include "AliAODEvent.h"
47#include "AliAODHandler.h"
48#include "AliAODTrack.h"
49#include "AliAODJet.h"
50#include "AliAODMCParticle.h"
51#include "AliMCEventHandler.h"
52#include "AliMCEvent.h"
53#include "AliStack.h"
54#include "AliGenPythiaEventHeader.h"
55#include "AliJetKineReaderHeader.h"
56#include "AliGenCocktailEventHeader.h"
57#include "AliInputEventHandler.h"
58#include "AliAODJetEventBackground.h"
59
60#include "fastjet/PseudoJet.hh"
61#include "fastjet/ClusterSequenceArea.hh"
62#include "fastjet/AreaDefinition.hh"
63#include "fastjet/JetDefinition.hh"
64// get info on how fastjet was configured
65#include "fastjet/config.h"
66
67
68ClassImp(AliAnalysisTaskJetCluster)
69
70AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
71 delete fRef;
72}
73
74AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
75 fAOD(0x0),
76 fAODExtension(0x0),
77 fRef(new TRefArray),
78 fUseAODTrackInput(kFALSE),
79 fUseAODMCInput(kFALSE),
80 fUseGlobalSelection(kFALSE),
81 fUseBackgroundCalc(kFALSE),
82 fFilterMask(0),
83 fTrackTypeRec(kTrackUndef),
84 fTrackTypeGen(kTrackUndef),
85 fNSkipLeadingRan(0),
86 fAvgTrials(1),
87 fExternalWeight(1),
88 fRecEtaWindow(0.5),
89 fTrackPtCut(0.),
90 fJetOutputMinPt(1),
91 fNonStdBranch(""),
92 fNonStdFile(""),
93 fRparam(1.0),
94 fAlgorithm(fastjet::kt_algorithm),
95 fStrategy(fastjet::Best),
96 fRecombScheme(fastjet::BIpt_scheme),
97 fAreaType(fastjet::active_area),
98 fGhostArea(0.01),
99 fActiveAreaRepeats(1),
100 fGhostEtamax(1.5),
101 fh1Xsec(0x0),
102 fh1Trials(0x0),
103 fh1PtHard(0x0),
104 fh1PtHardNoW(0x0),
105 fh1PtHardTrials(0x0),
106 fh1NJetsRec(0x0),
107 fh1NConstRec(0x0),
108 fh1NConstLeadingRec(0x0),
109 fh1PtJetsRecIn(0x0),
110 fh1PtJetsLeadingRecIn(0x0),
111 fh1PtJetConstRec(0x0),
112 fh1PtJetConstLeadingRec(0x0),
113 fh1PtTracksRecIn(0x0),
114 fh1PtTracksLeadingRecIn(0x0),
115 fh1NJetsRecRan(0x0),
116 fh1NConstRecRan(0x0),
117 fh1PtJetsLeadingRecInRan(0x0),
118 fh1NConstLeadingRecRan(0x0),
119 fh1PtJetsRecInRan(0x0),
120 fh1PtTracksGenIn(0x0),
121 fh1Nch(0x0),
122 fh2NRecJetsPt(0x0),
123 fh2NRecTracksPt(0x0),
124 fh2NConstPt(0x0),
125 fh2NConstLeadingPt(0x0),
126 fh2JetPhiEta(0x0),
127 fh2LeadingJetPhiEta(0x0),
128 fh2JetEtaPt(0x0),
129 fh2LeadingJetEtaPt(0x0),
130 fh2TrackEtaPt(0x0),
131 fh2LeadingTrackEtaPt(0x0),
132 fh2JetsLeadingPhiEta(0x0),
133 fh2JetsLeadingPhiPt(0x0),
134 fh2TracksLeadingPhiEta(0x0),
135 fh2TracksLeadingPhiPt(0x0),
136 fh2TracksLeadingJetPhiPt(0x0),
137 fh2JetsLeadingPhiPtW(0x0),
138 fh2TracksLeadingPhiPtW(0x0),
139 fh2TracksLeadingJetPhiPtW(0x0),
140 fh2NRecJetsPtRan(0x0),
141 fh2NConstPtRan(0x0),
142 fh2NConstLeadingPtRan(0x0),
143 fh2PtNch(0x0),
144 fh2PtNchRan(0x0),
145 fh2PtNchN(0x0),
146 fh2PtNchNRan(0x0),
147 fh2TracksLeadingJetPhiPtRan(0x0),
148 fh2TracksLeadingJetPhiPtWRan(0x0),
149 fHistList(0x0)
150{
151 for(int i = 0;i<3;i++){
152 fh1BiARandomCones[i] = 0;
153 fh1BiARandomConesRan[i] = 0;
154 }
155}
156
157AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
158 AliAnalysisTaskSE(name),
159 fAOD(0x0),
160 fAODExtension(0x0),
161 fRef(new TRefArray),
162 fUseAODTrackInput(kFALSE),
163 fUseAODMCInput(kFALSE),
164 fUseGlobalSelection(kFALSE),
165 fUseBackgroundCalc(kFALSE),
166 fFilterMask(0),
167 fTrackTypeRec(kTrackUndef),
168 fTrackTypeGen(kTrackUndef),
169 fNSkipLeadingRan(0),
170 fAvgTrials(1),
171 fExternalWeight(1),
172 fRecEtaWindow(0.5),
173 fTrackPtCut(0.),
174 fJetOutputMinPt(1),
175 fNonStdBranch(""),
176 fNonStdFile(""),
177 fRparam(1.0),
178 fAlgorithm(fastjet::kt_algorithm),
179 fStrategy(fastjet::Best),
180 fRecombScheme(fastjet::BIpt_scheme),
181 fAreaType(fastjet::active_area),
182 fGhostArea(0.01),
183 fActiveAreaRepeats(1),
184 fGhostEtamax(1.5),
185 fh1Xsec(0x0),
186 fh1Trials(0x0),
187 fh1PtHard(0x0),
188 fh1PtHardNoW(0x0),
189 fh1PtHardTrials(0x0),
190 fh1NJetsRec(0x0),
191 fh1NConstRec(0x0),
192 fh1NConstLeadingRec(0x0),
193 fh1PtJetsRecIn(0x0),
194 fh1PtJetsLeadingRecIn(0x0),
195 fh1PtJetConstRec(0x0),
196 fh1PtJetConstLeadingRec(0x0),
197 fh1PtTracksRecIn(0x0),
198 fh1PtTracksLeadingRecIn(0x0),
199 fh1NJetsRecRan(0x0),
200 fh1NConstRecRan(0x0),
201 fh1PtJetsLeadingRecInRan(0x0),
202 fh1NConstLeadingRecRan(0x0),
203 fh1PtJetsRecInRan(0x0),
204 fh1PtTracksGenIn(0x0),
205 fh1Nch(0x0),
206 fh2NRecJetsPt(0x0),
207 fh2NRecTracksPt(0x0),
208 fh2NConstPt(0x0),
209 fh2NConstLeadingPt(0x0),
210 fh2JetPhiEta(0x0),
211 fh2LeadingJetPhiEta(0x0),
212 fh2JetEtaPt(0x0),
213 fh2LeadingJetEtaPt(0x0),
214 fh2TrackEtaPt(0x0),
215 fh2LeadingTrackEtaPt(0x0),
216 fh2JetsLeadingPhiEta(0x0),
217 fh2JetsLeadingPhiPt(0x0),
218 fh2TracksLeadingPhiEta(0x0),
219 fh2TracksLeadingPhiPt(0x0),
220 fh2TracksLeadingJetPhiPt(0x0),
221 fh2JetsLeadingPhiPtW(0x0),
222 fh2TracksLeadingPhiPtW(0x0),
223 fh2TracksLeadingJetPhiPtW(0x0),
224 fh2NRecJetsPtRan(0x0),
225 fh2NConstPtRan(0x0),
226 fh2NConstLeadingPtRan(0x0),
227 fh2PtNch(0x0),
228 fh2PtNchRan(0x0),
229 fh2PtNchN(0x0),
230 fh2PtNchNRan(0x0),
231 fh2TracksLeadingJetPhiPtRan(0x0),
232 fh2TracksLeadingJetPhiPtWRan(0x0),
233 fHistList(0x0)
234{
235 for(int i = 0;i<3;i++){
236 fh1BiARandomCones[i] = 0;
237 fh1BiARandomConesRan[i] = 0;
238 }
239 DefineOutput(1, TList::Class());
240}
241
242
243
244Bool_t AliAnalysisTaskJetCluster::Notify()
245{
246 //
247 // Implemented Notify() to read the cross sections
248 // and number of trials from pyxsec.root
249 //
250 return kTRUE;
251}
252
253void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
254{
255
256 //
257 // Create the output container
258 //
259
260
261
262
263 // Connect the AOD
264
265
266 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
267
268
269
270 if(fNonStdBranch.Length()!=0)
271 {
272 // only create the output branch if we have a name
273 // Create a new branch for jets...
274 // -> cleared in the UserExec....
275 // here we can also have the case that the brnaches are written to a separate file
276
277 TClonesArray *tca = new TClonesArray("AliAODJet", 0);
278 tca->SetName(fNonStdBranch.Data());
279 AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
280
281
282 TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
283 tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
284 AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
285
286
287 if(fUseBackgroundCalc){
288 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
289 AliAODJetEventBackground* evBkg = new AliAODJetEventBackground();
290 evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
291 AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data());
292 }
293 }
294
295 if(fNonStdFile.Length()!=0){
296 //
297 // case that we have an AOD extension we need to fetch the jets from the extended output
298 // we identifay the extension aod event by looking for the branchname
299 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
300 TObjArray* extArray = aodH->GetExtensions();
301 if (extArray) {
302 TIter next(extArray);
303 while ((fAODExtension=(AliAODExtension*)next())){
304 TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data());
305 if(fDebug>10){
306 Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__);
307 fAODExtension->GetAOD()->Dump();
308 }
309 if(obj){
310 if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data());
311 break;
312 }
313 fAODExtension = 0;
314 }
315 }
316 }
317 }
318
319
320 OpenFile(1);
321 if(!fHistList)fHistList = new TList();
322 fHistList->SetOwner();
323
324 Bool_t oldStatus = TH1::AddDirectoryStatus();
325 TH1::AddDirectory(kFALSE);
326
327 //
328 // Histogram
329
330 const Int_t nBinPt = 200;
331 Double_t binLimitsPt[nBinPt+1];
332 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
333 if(iPt == 0){
334 binLimitsPt[iPt] = 0.0;
335 }
336 else {// 1.0
337 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 0.5;
338 }
339 }
340
341 const Int_t nBinPhi = 90;
342 Double_t binLimitsPhi[nBinPhi+1];
343 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
344 if(iPhi==0){
345 binLimitsPhi[iPhi] = -1.*TMath::Pi();
346 }
347 else{
348 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
349 }
350 }
351
352
353
354 const Int_t nBinEta = 40;
355 Double_t binLimitsEta[nBinEta+1];
356 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
357 if(iEta==0){
358 binLimitsEta[iEta] = -2.0;
359 }
360 else{
361 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
362 }
363 }
364
365 const int nChMax = 100;
366
367 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
368 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
369
370 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
371 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
372
373
374 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
375 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
376
377 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
378 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
379 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
380 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
381
382
383 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
384 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
385 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
386
387 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
388 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
389 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
390 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
391 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
392 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
393 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
394 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
395 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
396 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
397
398
399 for(int i = 0;i<3;i++){
400 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
401 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
402 }
403
404 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
405 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
406 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
407 //
408
409
410 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
411 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
412 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
413 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
414
415 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
416 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);
417 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);
418 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);
419
420
421
422 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
423 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
424 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
425 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
426
427 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
428 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
429 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
430 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
431
432 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
433 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
434 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
435 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
436
437
438
439 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
440 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
441 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
442 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
443 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
444 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
445 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
446 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
447 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
448 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
449 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
450 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
451
452 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
453 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
454 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
455 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
456
457 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
458 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
459 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
460 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
461
462
463
464 const Int_t saveLevel = 3; // large save level more histos
465 if(saveLevel>0){
466 fHistList->Add(fh1Xsec);
467 fHistList->Add(fh1Trials);
468
469 fHistList->Add(fh1NJetsRec);
470 fHistList->Add(fh1NConstRec);
471 fHistList->Add(fh1NConstLeadingRec);
472 fHistList->Add(fh1PtJetsRecIn);
473 fHistList->Add(fh1PtJetsLeadingRecIn);
474 fHistList->Add(fh1PtTracksRecIn);
475 fHistList->Add(fh1PtTracksLeadingRecIn);
476 fHistList->Add(fh1PtJetConstRec);
477 fHistList->Add(fh1PtJetConstLeadingRec);
478 fHistList->Add(fh1NJetsRecRan);
479 fHistList->Add(fh1NConstRecRan);
480 fHistList->Add(fh1PtJetsLeadingRecInRan);
481 fHistList->Add(fh1NConstLeadingRecRan);
482 fHistList->Add(fh1PtJetsRecInRan);
483 fHistList->Add(fh1Nch);
484 for(int i = 0;i<3;i++){
485 fHistList->Add(fh1BiARandomCones[i]);
486 fHistList->Add(fh1BiARandomConesRan[i]);
487 }
488 fHistList->Add(fh2NRecJetsPt);
489 fHistList->Add(fh2NRecTracksPt);
490 fHistList->Add(fh2NConstPt);
491 fHistList->Add(fh2NConstLeadingPt);
492 fHistList->Add(fh2PtNch);
493 fHistList->Add(fh2PtNchRan);
494 fHistList->Add(fh2PtNchN);
495 fHistList->Add(fh2PtNchNRan);
496 fHistList->Add(fh2JetPhiEta);
497 fHistList->Add(fh2LeadingJetPhiEta);
498 fHistList->Add(fh2JetEtaPt);
499 fHistList->Add(fh2LeadingJetEtaPt);
500 fHistList->Add(fh2TrackEtaPt);
501 fHistList->Add(fh2LeadingTrackEtaPt);
502 fHistList->Add(fh2JetsLeadingPhiEta );
503 fHistList->Add(fh2JetsLeadingPhiPt);
504 fHistList->Add(fh2TracksLeadingPhiEta);
505 fHistList->Add(fh2TracksLeadingPhiPt);
506 fHistList->Add(fh2TracksLeadingJetPhiPt);
507 fHistList->Add(fh2JetsLeadingPhiPtW);
508 fHistList->Add(fh2TracksLeadingPhiPtW);
509 fHistList->Add(fh2TracksLeadingJetPhiPtW);
510 fHistList->Add(fh2NRecJetsPtRan);
511 fHistList->Add(fh2NConstPtRan);
512 fHistList->Add(fh2NConstLeadingPtRan);
513 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
514 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
515 }
516
517 // =========== Switch on Sumw2 for all histos ===========
518 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
519 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
520 if (h1){
521 h1->Sumw2();
522 continue;
523 }
524 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
525 if(hn)hn->Sumw2();
526 }
527 TH1::AddDirectory(oldStatus);
528}
529
530void AliAnalysisTaskJetCluster::Init()
531{
532 //
533 // Initialization
534 //
535
536 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
537
538}
539
540void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
541{
542
543 if(fUseGlobalSelection){
544 // no selection by the service task, we continue
545 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
546 PostData(1, fHistList);
547 return;
548 }
549
550 // handle and reset the output jet branch
551 // only need this once
552 TClonesArray* jarray = 0;
553 TClonesArray* jarrayran = 0;
554 AliAODJetEventBackground* evBkg = 0;
555 if(fNonStdBranch.Length()!=0) {
556 if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
557 if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
558 if(jarray)jarray->Delete(); // this is our responsibility, clear before filling again
559 if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
560 if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
561 if(jarrayran)jarrayran->Delete(); // this is our responsibility, clear before filling again
562
563 if(fUseBackgroundCalc){
564 if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
565 if(!evBkg) evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
566 if(evBkg)evBkg->Reset();
567 }
568 }
569
570
571
572 //
573 // Execute analysis for current event
574 //
575 AliESDEvent *fESD = 0;
576 if(fUseAODTrackInput){
577 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
578 if(!fAOD){
579 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
580 return;
581 }
582 // fethc the header
583 }
584 else{
585 // assume that the AOD is in the general output...
586 fAOD = AODEvent();
587 if(!fAOD){
588 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
589 return;
590 }
591 if(fDebug>0){
592 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
593 }
594 }
595
596 Bool_t selectEvent = false;
597 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
598
599 if(fAOD){
600 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
601 TString vtxTitle(vtxAOD->GetTitle());
602
603 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
604 Float_t zvtx = vtxAOD->GetZ();
605 Float_t yvtx = vtxAOD->GetY();
606 Float_t xvtx = vtxAOD->GetX();
607 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
608 if(TMath::Abs(zvtx)<8.&&r2<1.){
609 if(physicsSelection)selectEvent = true;
610 }
611 }
612 }
613 if(!selectEvent){
614 PostData(1, fHistList);
615 return;
616 }
617
618 fh1Trials->Fill("#sum{ntrials}",1);
619
620
621 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
622
623 // ==== General variables needed
624
625
626
627 // we simply fetch the tracks/mc particles as a list of AliVParticles
628
629 TList recParticles;
630 TList genParticles;
631
632 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
633 Float_t nCh = recParticles.GetEntries();
634 fh1Nch->Fill(nCh);
635 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
636 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
637 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
638
639 // find the jets....
640
641 vector<fastjet::PseudoJet> inputParticlesRec;
642 vector<fastjet::PseudoJet> inputParticlesRecRan;
643
644
645 // create a random jet within the acceptance
646 const float rRandomCone2 = 0.4*0.4;
647 float etaRandomCone = gRandom->Rndm()-0.5; // +- 0.5
648 float phiRandomCone = gRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
649 float ptRandomCone = 0;
650 float ptRandomConeRan = 0;
651
652 for(int i = 0; i < recParticles.GetEntries(); i++){
653 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
654 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
655 // we talk total momentum here
656 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
657 jInp.set_user_index(i);
658 inputParticlesRec.push_back(jInp);
659
660 if(evBkg){
661 float deta = vp->Eta()-etaRandomCone;
662 float dphi = vp->Phi()-phiRandomCone;
663 if(dphi>TMath::Pi())dphi -= 2.*TMath::Pi();
664 if(dphi<-TMath::Pi())dphi += 2.*TMath::Pi();
665 float dr2 = dphi*dphi+deta*deta;
666 if(dr2<rRandomCone2){
667 // within random jet
668 ptRandomCone += vp->Pt();
669 }
670 }
671
672 // the randomized input changes eta and phi, but keeps the p_T
673 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
674 Double_t pT = vp->Pt();
675 Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
676 Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
677
678 Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));
679 Double_t pZ = pT/TMath::Tan(theta);
680
681 Double_t pX = pT * TMath::Cos(phi);
682 Double_t pY = pT * TMath::Sin(phi);
683 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
684 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
685 jInpRan.set_user_index(i);
686 inputParticlesRecRan.push_back(jInpRan);
687
688 if(evBkg){
689 float deta = eta-etaRandomCone;
690 float dphi = phi-phiRandomCone;
691 if(dphi>TMath::Pi())dphi -= 2.*TMath::Pi();
692 if(dphi<-TMath::Pi())dphi += 2.*TMath::Pi();
693 float dr2 = dphi*dphi+deta*deta;
694 if(dr2<rRandomCone2){
695 // within random jet
696 ptRandomConeRan += pT;
697 }
698 }
699 }
700
701 // fill the tref array, only needed when we write out jets
702 if(jarray){
703 if(i == 0){
704 fRef->Delete(); // make sure to delete before placement new...
705 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
706 }
707 fRef->Add(vp);
708 }
709 }
710
711 if(inputParticlesRec.size()==0){
712 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
713 PostData(1, fHistList);
714 return;
715 }
716
717 // run fast jet
718 // employ setters for these...
719
720
721 // now create the object that holds info about ghosts
722 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
723 // reduce CPU time...
724 fGhostArea = 0.5;
725 fActiveAreaRepeats = 0;
726 }
727 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
728 fastjet::AreaType areaType = fastjet::active_area;
729 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
730 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
731 vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
732 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
733
734 //range where to compute background
735 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
736 phiMin = 0;
737 phiMax = 2*TMath::Pi();
738 rapMax = fGhostEtamax - fRparam;
739 rapMin = - fGhostEtamax + fRparam;
740 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
741
742
743
744 inclusiveJets = clustSeq.inclusive_jets();
745 sortedJets = sorted_by_pt(inclusiveJets);
746
747 fh1NJetsRec->Fill(sortedJets.size());
748
749 // loop over all jets an fill information, first on is the leading jet
750
751 Int_t nRecOver = inclusiveJets.size();
752 Int_t nRec = inclusiveJets.size();
753 if(inclusiveJets.size()>0){
754 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
755 Float_t pt = leadingJet.Pt();
756 Int_t nAodOutJets = 0;
757 Int_t nAodOutTracks = 0;
758 AliAODJet *aodOutJet = 0;
759
760 Int_t iCount = 0;
761 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
762 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
763 while(pt<ptCut&&iCount<nRec){
764 nRecOver--;
765 iCount++;
766 if(iCount<nRec){
767 pt = sortedJets[iCount].perp();
768 }
769 }
770 if(nRecOver<=0)break;
771 fh2NRecJetsPt->Fill(ptCut,nRecOver);
772 }
773 Float_t phi = leadingJet.Phi();
774 if(phi<0)phi+=TMath::Pi()*2.;
775 Float_t eta = leadingJet.Eta();
776 pt = leadingJet.Pt();
777
778 // correlation of leading jet with tracks
779 TIterator *recIter = recParticles.MakeIterator();
780 recIter->Reset();
781 AliVParticle *tmpRecTrack = 0;
782 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
783 Float_t tmpPt = tmpRecTrack->Pt();
784 // correlation
785 Float_t tmpPhi = tmpRecTrack->Phi();
786 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
787 Float_t dPhi = phi - tmpPhi;
788 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
789 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
790 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
791 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
792 }
793
794
795
796 for(int j = 0; j < nRec;j++){
797 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
798 aodOutJet = 0;
799 nAodOutTracks = 0;
800 Float_t tmpPt = tmpRec.Pt();
801 fh1PtJetsRecIn->Fill(tmpPt);
802 // Fill Spectra with constituents
803 vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
804 fh1NConstRec->Fill(constituents.size());
805 fh2PtNch->Fill(nCh,tmpPt);
806 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
807 fh2NConstPt->Fill(tmpPt,constituents.size());
808 // loop over constiutents and fill spectrum
809
810 // Add the jet information and the track references to the output AOD
811
812 if(tmpPt>fJetOutputMinPt&&jarray){
813 aodOutJet = new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
814 Double_t area=clustSeq.area(sortedJets[j]);
815
816 aodOutJet->SetEffArea(area,0);
817 }
818
819
820 for(unsigned int ic = 0; ic < constituents.size();ic++){
821 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
822 fh1PtJetConstRec->Fill(part->Pt());
823 if(aodOutJet){
824 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
825 }
826 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
827 }
828
829
830
831
832
833
834
835
836 // correlation
837 Float_t tmpPhi = tmpRec.Phi();
838 Float_t tmpEta = tmpRec.Eta();
839 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
840
841 if(j==0){
842 fh1PtJetsLeadingRecIn->Fill(tmpPt);
843 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
844 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
845 fh1NConstLeadingRec->Fill(constituents.size());
846 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
847 continue;
848 }
849 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
850 fh2JetEtaPt->Fill(tmpEta,tmpPt);
851 Float_t dPhi = phi - tmpPhi;
852 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
853 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
854 Float_t dEta = eta - tmpRec.Eta();
855 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
856 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
857 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
858 }
859 delete recIter;
860
861 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
862
863 if(evBkg){
864 vector<fastjet::PseudoJet> jets2=sortedJets;
865 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
866 Double_t bkg1=0;
867 Double_t sigma1=0.;
868 Double_t meanarea1=0.;
869 Double_t bkg2=0;
870 Double_t sigma2=0.;
871 Double_t meanarea2=0.;
872
873 float areaRandomCone = rRandomCone2 *TMath::Pi();
874 clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
875 evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
876 fh1BiARandomCones[0]->Fill(ptRandomCone-(bkg1*areaRandomCone));
877 fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
878
879 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
880 evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
881 fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
882 fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
883 }
884 }
885
886
887
888
889
890 // fill track information
891 Int_t nTrackOver = recParticles.GetSize();
892 // do the same for tracks and jets
893
894 if(nTrackOver>0){
895 TIterator *recIter = recParticles.MakeIterator();
896 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
897 Float_t pt = tmpRec->Pt();
898
899 // Printf("Leading track p_t %3.3E",pt);
900 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
901 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
902 while(pt<ptCut&&tmpRec){
903 nTrackOver--;
904 tmpRec = (AliVParticle*)(recIter->Next());
905 if(tmpRec){
906 pt = tmpRec->Pt();
907 }
908 }
909 if(nTrackOver<=0)break;
910 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
911 }
912
913 recIter->Reset();
914 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
915 Float_t phi = leading->Phi();
916 if(phi<0)phi+=TMath::Pi()*2.;
917 Float_t eta = leading->Eta();
918 pt = leading->Pt();
919 while((tmpRec = (AliVParticle*)(recIter->Next()))){
920 Float_t tmpPt = tmpRec->Pt();
921 Float_t tmpEta = tmpRec->Eta();
922 fh1PtTracksRecIn->Fill(tmpPt);
923 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
924 if(tmpRec==leading){
925 fh1PtTracksLeadingRecIn->Fill(tmpPt);
926 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
927 continue;
928 }
929 // correlation
930 Float_t tmpPhi = tmpRec->Phi();
931
932 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
933 Float_t dPhi = phi - tmpPhi;
934 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
935 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
936 Float_t dEta = eta - tmpRec->Eta();
937 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
938 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
939 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
940 }
941 delete recIter;
942 }
943
944 // find the random jets
945 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
946 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
947
948 inclusiveJetsRan = clustSeqRan.inclusive_jets();
949 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
950
951 // fill the jet information from random track
952
953 fh1NJetsRecRan->Fill(sortedJetsRan.size());
954
955 // loop over all jets an fill information, first on is the leading jet
956
957 Int_t nRecOverRan = inclusiveJetsRan.size();
958 Int_t nRecRan = inclusiveJetsRan.size();
959 if(inclusiveJetsRan.size()>0){
960 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
961 Float_t pt = leadingJet.Pt();
962
963 Int_t iCount = 0;
964 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
965 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
966 while(pt<ptCut&&iCount<nRecRan){
967 nRecOverRan--;
968 iCount++;
969 if(iCount<nRecRan){
970 pt = sortedJetsRan[iCount].perp();
971 }
972 }
973 if(nRecOverRan<=0)break;
974 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
975 }
976 Float_t phi = leadingJet.Phi();
977 if(phi<0)phi+=TMath::Pi()*2.;
978 pt = leadingJet.Pt();
979
980 // correlation of leading jet with random tracks
981
982 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
983 {
984 Float_t tmpPt = inputParticlesRecRan[ip].perp();
985 // correlation
986 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
987 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
988 Float_t dPhi = phi - tmpPhi;
989 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
990 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
991 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
992 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
993 }
994
995 Int_t nAodOutJetsRan = 0;
996 AliAODJet *aodOutJetRan = 0;
997 for(int j = 0; j < nRecRan;j++){
998 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
999 Float_t tmpPt = tmpRec.Pt();
1000 fh1PtJetsRecInRan->Fill(tmpPt);
1001 // Fill Spectra with constituents
1002 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1003 fh1NConstRecRan->Fill(constituents.size());
1004 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1005 fh2PtNchRan->Fill(nCh,tmpPt);
1006 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1007
1008
1009 if(tmpPt>fJetOutputMinPt&&jarrayran){
1010 aodOutJetRan = new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1011 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1012
1013 aodOutJetRan->SetEffArea(arearan,0); }
1014
1015
1016
1017
1018 for(unsigned int ic = 0; ic < constituents.size();ic++){
1019 // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1020 // fh1PtJetConstRec->Fill(part->Pt());
1021 if(aodOutJetRan){
1022 aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
1023 }}
1024
1025
1026
1027
1028 // correlation
1029 Float_t tmpPhi = tmpRec.Phi();
1030 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1031
1032 if(j==0){
1033 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1034 fh1NConstLeadingRecRan->Fill(constituents.size());
1035 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1036 continue;
1037 }
1038 }
1039
1040
1041 if(evBkg){
1042 Double_t bkg3=0.;
1043 Double_t sigma3=0.;
1044 Double_t meanarea3=0.;
1045 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1046 evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
1047 float areaRandomCone = rRandomCone2 *TMath::Pi();
1048 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1049 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1050 }
1051
1052
1053
1054 }
1055
1056
1057 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1058 PostData(1, fHistList);
1059}
1060
1061void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1062{
1063// Terminate analysis
1064//
1065 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1066}
1067
1068
1069Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1070
1071 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1072
1073 Int_t iCount = 0;
1074 if(type==kTrackAOD){
1075 AliAODEvent *aod = 0;
1076 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1077 else aod = AODEvent();
1078 if(!aod){
1079 return iCount;
1080 }
1081 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1082 AliAODTrack *tr = aod->GetTrack(it);
1083 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1084 if(TMath::Abs(tr->Eta())>0.9)continue;
1085 // if(tr->Pt()<0.3)continue;
1086 list->Add(tr);
1087 iCount++;
1088 }
1089 }
1090 else if (type == kTrackKineAll||type == kTrackKineCharged){
1091 AliMCEvent* mcEvent = MCEvent();
1092 if(!mcEvent)return iCount;
1093 // we want to have alivpartilces so use get track
1094 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1095 if(!mcEvent->IsPhysicalPrimary(it))continue;
1096 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1097 if(type == kTrackKineAll){
1098 list->Add(part);
1099 iCount++;
1100 }
1101 else if(type == kTrackKineCharged){
1102 if(part->Particle()->GetPDG()->Charge()==0)continue;
1103 list->Add(part);
1104 iCount++;
1105 }
1106 }
1107 }
1108 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1109 AliAODEvent *aod = 0;
1110 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1111 else aod = AODEvent();
1112 if(!aod)return iCount;
1113 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1114 if(!tca)return iCount;
1115 for(int it = 0;it < tca->GetEntriesFast();++it){
1116 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1117 if(!part->IsPhysicalPrimary())continue;
1118 if(type == kTrackAODMCAll){
1119 list->Add(part);
1120 iCount++;
1121 }
1122 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1123 if(part->Charge()==0)continue;
1124 if(kTrackAODMCCharged){
1125 list->Add(part);
1126 }
1127 else {
1128 if(TMath::Abs(part->Eta())>0.9)continue;
1129 list->Add(part);
1130 }
1131 iCount++;
1132 }
1133 }
1134 }// AODMCparticle
1135 list->Sort();
1136 return iCount;
1137
1138}
1139
1140/*
1141Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1142 for(int i = 0; i < particles.GetEntries(); i++){
1143 AliVParticle *vp = (AliVParticle*)particles.At(i);
1144 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1145 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1146 jInp.set_user_index(i);
1147 inputParticles.push_back(jInp);
1148 }
1149
1150 return 0;
1151
1152}
1153*/