]> 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
Added PostData in UserCreateOutputObjects
[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>
d7ca0e14 24#include <TRandom3.h>
dbf053f0 25#include <TSystem.h>
26#include <TInterpreter.h>
27#include <TChain.h>
54424110 28#include <TRefArray.h>
dbf053f0 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"
dbf053f0 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"
d6e66a82 58#include "AliAODJetEventBackground.h"
dbf053f0 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
54424110 70AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
71 delete fRef;
d7ca0e14 72 delete fRandom;
8f85b773 73
74 if(fTCAJetsOut)fTCAJetsOut->Delete();
75 delete fTCAJetsOut;
76 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
77 delete fTCAJetsOutRan;
78 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
79 delete fTCARandomConesOut;
80 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
81 delete fTCARandomConesOutRan;
82 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
83 delete fAODJetBackgroundOut;
54424110 84}
85
dbf053f0 86AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
87 fAOD(0x0),
54424110 88 fAODExtension(0x0),
89 fRef(new TRefArray),
dbf053f0 90 fUseAODTrackInput(kFALSE),
91 fUseAODMCInput(kFALSE),
92 fUseGlobalSelection(kFALSE),
9fd9f55b 93 fUseBackgroundCalc(kFALSE),
dbf053f0 94 fFilterMask(0),
95 fTrackTypeRec(kTrackUndef),
bd80a748 96 fTrackTypeGen(kTrackUndef),
97 fNSkipLeadingRan(0),
c4856527 98 fNRandomCones(0),
bd80a748 99 fAvgTrials(1),
8a320ab9 100 fExternalWeight(1),
101 fTrackEtaWindow(0.9),
dbf053f0 102 fRecEtaWindow(0.5),
54424110 103 fTrackPtCut(0.),
104 fJetOutputMinPt(1),
0341d0d8 105 fJetTriggerPtCut(0),
aae7993b 106 fCentCutUp(0),
107 fCentCutLo(0),
54424110 108 fNonStdBranch(""),
0341d0d8 109 fBackgroundBranch(""),
54424110 110 fNonStdFile(""),
dbf053f0 111 fRparam(1.0),
112 fAlgorithm(fastjet::kt_algorithm),
113 fStrategy(fastjet::Best),
114 fRecombScheme(fastjet::BIpt_scheme),
115 fAreaType(fastjet::active_area),
d6e66a82 116 fGhostArea(0.01),
117 fActiveAreaRepeats(1),
118 fGhostEtamax(1.5),
8f85b773 119 fTCAJetsOut(0x0),
120 fTCAJetsOutRan(0x0),
121 fTCARandomConesOut(0x0),
122 fTCARandomConesOutRan(0x0),
123 fAODJetBackgroundOut(0x0),
d7ca0e14 124 fRandom(0),
dbf053f0 125 fh1Xsec(0x0),
126 fh1Trials(0x0),
127 fh1PtHard(0x0),
128 fh1PtHardNoW(0x0),
129 fh1PtHardTrials(0x0),
130 fh1NJetsRec(0x0),
131 fh1NConstRec(0x0),
132 fh1NConstLeadingRec(0x0),
133 fh1PtJetsRecIn(0x0),
134 fh1PtJetsLeadingRecIn(0x0),
135 fh1PtJetConstRec(0x0),
136 fh1PtJetConstLeadingRec(0x0),
137 fh1PtTracksRecIn(0x0),
138 fh1PtTracksLeadingRecIn(0x0),
139 fh1NJetsRecRan(0x0),
140 fh1NConstRecRan(0x0),
141 fh1PtJetsLeadingRecInRan(0x0),
142 fh1NConstLeadingRecRan(0x0),
143 fh1PtJetsRecInRan(0x0),
144 fh1PtTracksGenIn(0x0),
a0d08b6d 145 fh1Nch(0x0),
aae7993b 146 fh1CentralityPhySel(0x0),
0341d0d8 147 fh1Centrality(0x0),
148 fh1CentralitySelect(0x0),
e2ca7519 149 fh1ZPhySel(0x0),
150 fh1Z(0x0),
151 fh1ZSelect(0x0),
dbf053f0 152 fh2NRecJetsPt(0x0),
153 fh2NRecTracksPt(0x0),
154 fh2NConstPt(0x0),
155 fh2NConstLeadingPt(0x0),
156 fh2JetPhiEta(0x0),
157 fh2LeadingJetPhiEta(0x0),
158 fh2JetEtaPt(0x0),
159 fh2LeadingJetEtaPt(0x0),
160 fh2TrackEtaPt(0x0),
161 fh2LeadingTrackEtaPt(0x0),
162 fh2JetsLeadingPhiEta(0x0),
163 fh2JetsLeadingPhiPt(0x0),
164 fh2TracksLeadingPhiEta(0x0),
165 fh2TracksLeadingPhiPt(0x0),
166 fh2TracksLeadingJetPhiPt(0x0),
167 fh2JetsLeadingPhiPtW(0x0),
168 fh2TracksLeadingPhiPtW(0x0),
169 fh2TracksLeadingJetPhiPtW(0x0),
170 fh2NRecJetsPtRan(0x0),
171 fh2NConstPtRan(0x0),
172 fh2NConstLeadingPtRan(0x0),
a0d08b6d 173 fh2PtNch(0x0),
174 fh2PtNchRan(0x0),
175 fh2PtNchN(0x0),
176 fh2PtNchNRan(0x0),
dbf053f0 177 fh2TracksLeadingJetPhiPtRan(0x0),
178 fh2TracksLeadingJetPhiPtWRan(0x0),
179 fHistList(0x0)
180{
82a35b1a 181 for(int i = 0;i<3;i++){
182 fh1BiARandomCones[i] = 0;
183 fh1BiARandomConesRan[i] = 0;
184 }
aae7993b 185 for(int i = 0;i<kMaxCent;i++){
186 fh2JetsLeadingPhiPtC[i] = 0;
187 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
188 fh2TracksLeadingJetPhiPtC[i] = 0;
189 fh2TracksLeadingJetPhiPtWC[i] = 0;
190 }
dbf053f0 191}
192
193AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
194 AliAnalysisTaskSE(name),
195 fAOD(0x0),
54424110 196 fAODExtension(0x0),
197 fRef(new TRefArray),
dbf053f0 198 fUseAODTrackInput(kFALSE),
199 fUseAODMCInput(kFALSE),
200 fUseGlobalSelection(kFALSE),
9fd9f55b 201 fUseBackgroundCalc(kFALSE),
dbf053f0 202 fFilterMask(0),
203 fTrackTypeRec(kTrackUndef),
204 fTrackTypeGen(kTrackUndef),
bd80a748 205 fNSkipLeadingRan(0),
c4856527 206 fNRandomCones(0),
dbf053f0 207 fAvgTrials(1),
208 fExternalWeight(1),
8a320ab9 209 fTrackEtaWindow(0.9),
dbf053f0 210 fRecEtaWindow(0.5),
54424110 211 fTrackPtCut(0.),
212 fJetOutputMinPt(1),
0341d0d8 213 fJetTriggerPtCut(0),
aae7993b 214 fCentCutUp(0),
215 fCentCutLo(0),
54424110 216 fNonStdBranch(""),
0341d0d8 217 fBackgroundBranch(""),
54424110 218 fNonStdFile(""),
dbf053f0 219 fRparam(1.0),
220 fAlgorithm(fastjet::kt_algorithm),
221 fStrategy(fastjet::Best),
222 fRecombScheme(fastjet::BIpt_scheme),
223 fAreaType(fastjet::active_area),
d6e66a82 224 fGhostArea(0.01),
225 fActiveAreaRepeats(1),
226 fGhostEtamax(1.5),
8f85b773 227 fTCAJetsOut(0x0),
228 fTCAJetsOutRan(0x0),
229 fTCARandomConesOut(0x0),
230 fTCARandomConesOutRan(0x0),
231 fAODJetBackgroundOut(0x0),
d7ca0e14 232 fRandom(0),
dbf053f0 233 fh1Xsec(0x0),
234 fh1Trials(0x0),
235 fh1PtHard(0x0),
236 fh1PtHardNoW(0x0),
237 fh1PtHardTrials(0x0),
238 fh1NJetsRec(0x0),
239 fh1NConstRec(0x0),
240 fh1NConstLeadingRec(0x0),
241 fh1PtJetsRecIn(0x0),
242 fh1PtJetsLeadingRecIn(0x0),
243 fh1PtJetConstRec(0x0),
244 fh1PtJetConstLeadingRec(0x0),
245 fh1PtTracksRecIn(0x0),
246 fh1PtTracksLeadingRecIn(0x0),
247 fh1NJetsRecRan(0x0),
248 fh1NConstRecRan(0x0),
249 fh1PtJetsLeadingRecInRan(0x0),
250 fh1NConstLeadingRecRan(0x0),
251 fh1PtJetsRecInRan(0x0),
252 fh1PtTracksGenIn(0x0),
a0d08b6d 253 fh1Nch(0x0),
aae7993b 254 fh1CentralityPhySel(0x0),
0341d0d8 255 fh1Centrality(0x0),
256 fh1CentralitySelect(0x0),
e2ca7519 257 fh1ZPhySel(0x0),
258 fh1Z(0x0),
259 fh1ZSelect(0x0),
dbf053f0 260 fh2NRecJetsPt(0x0),
261 fh2NRecTracksPt(0x0),
262 fh2NConstPt(0x0),
263 fh2NConstLeadingPt(0x0),
264 fh2JetPhiEta(0x0),
265 fh2LeadingJetPhiEta(0x0),
266 fh2JetEtaPt(0x0),
267 fh2LeadingJetEtaPt(0x0),
268 fh2TrackEtaPt(0x0),
269 fh2LeadingTrackEtaPt(0x0),
270 fh2JetsLeadingPhiEta(0x0),
271 fh2JetsLeadingPhiPt(0x0),
272 fh2TracksLeadingPhiEta(0x0),
273 fh2TracksLeadingPhiPt(0x0),
274 fh2TracksLeadingJetPhiPt(0x0),
275 fh2JetsLeadingPhiPtW(0x0),
276 fh2TracksLeadingPhiPtW(0x0),
277 fh2TracksLeadingJetPhiPtW(0x0),
278 fh2NRecJetsPtRan(0x0),
279 fh2NConstPtRan(0x0),
280 fh2NConstLeadingPtRan(0x0),
a0d08b6d 281 fh2PtNch(0x0),
282 fh2PtNchRan(0x0),
283 fh2PtNchN(0x0),
284 fh2PtNchNRan(0x0),
dbf053f0 285 fh2TracksLeadingJetPhiPtRan(0x0),
286 fh2TracksLeadingJetPhiPtWRan(0x0),
287 fHistList(0x0)
288{
82a35b1a 289 for(int i = 0;i<3;i++){
290 fh1BiARandomCones[i] = 0;
291 fh1BiARandomConesRan[i] = 0;
292 }
aae7993b 293 for(int i = 0;i<kMaxCent;i++){
294 fh2JetsLeadingPhiPtC[i] = 0;
295 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
296 fh2TracksLeadingJetPhiPtC[i] = 0;
297 fh2TracksLeadingJetPhiPtWC[i] = 0;
298 }
82a35b1a 299 DefineOutput(1, TList::Class());
dbf053f0 300}
301
302
303
304Bool_t AliAnalysisTaskJetCluster::Notify()
305{
306 //
307 // Implemented Notify() to read the cross sections
308 // and number of trials from pyxsec.root
309 //
310 return kTRUE;
311}
312
313void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
314{
315
316 //
317 // Create the output container
318 //
319
d7ca0e14 320 fRandom = new TRandom3(0);
54424110 321
322
dbf053f0 323 // Connect the AOD
324
325
326 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
327
d6e66a82 328
54424110 329
330 if(fNonStdBranch.Length()!=0)
331 {
332 // only create the output branch if we have a name
333 // Create a new branch for jets...
334 // -> cleared in the UserExec....
335 // here we can also have the case that the brnaches are written to a separate file
336
8f85b773 337 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
338 fTCAJetsOut->SetName(fNonStdBranch.Data());
339 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
d6e66a82 340
341
8f85b773 342 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
343 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
344 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
d6e66a82 345
9fd9f55b 346 if(fUseBackgroundCalc){
347 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
8f85b773 348 fAODJetBackgroundOut = new AliAODJetEventBackground();
349 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
350 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
9fd9f55b 351 }
c4856527 352 }
353 // create the branch for the random cones with the same R
354 TString cName = Form("%sRandomCone",fNonStdBranch.Data());
355
356 if(fNRandomCones>0){
d7ca0e14 357 if(!AODEvent()->FindListObject(cName.Data())){
8f85b773 358 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
359 fTCARandomConesOut->SetName(cName.Data());
360 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
d7ca0e14 361 }
d7ca0e14 362 // create the branch with the random for the random cones on the random event
363 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
364 if(!AODEvent()->FindListObject(cName.Data())){
8f85b773 365 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
366 fTCARandomConesOutRan->SetName(cName.Data());
367 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
d7ca0e14 368 }
9fd9f55b 369 }
c4856527 370
54424110 371 if(fNonStdFile.Length()!=0){
372 //
373 // case that we have an AOD extension we need to fetch the jets from the extended output
f840d64c 374 // we identify the extension aod event by looking for the branchname
54424110 375 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
8f85b773 376 // case that we have an AOD extension we need can fetch the background maybe from the extended output
377 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
54424110 378 }
379 }
380
381
dbf053f0 382 if(!fHistList)fHistList = new TList();
9c6bd749 383 fHistList->SetOwner();
8f85b773 384
dbf053f0 385 Bool_t oldStatus = TH1::AddDirectoryStatus();
386 TH1::AddDirectory(kFALSE);
387
388 //
389 // Histogram
390
391 const Int_t nBinPt = 200;
392 Double_t binLimitsPt[nBinPt+1];
393 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
394 if(iPt == 0){
395 binLimitsPt[iPt] = 0.0;
396 }
397 else {// 1.0
d7ca0e14 398 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
dbf053f0 399 }
400 }
401
402 const Int_t nBinPhi = 90;
403 Double_t binLimitsPhi[nBinPhi+1];
404 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
405 if(iPhi==0){
406 binLimitsPhi[iPhi] = -1.*TMath::Pi();
407 }
408 else{
409 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
410 }
411 }
412
413
414
415 const Int_t nBinEta = 40;
416 Double_t binLimitsEta[nBinEta+1];
417 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
418 if(iEta==0){
419 binLimitsEta[iEta] = -2.0;
420 }
421 else{
422 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
423 }
424 }
425
d7ca0e14 426 const int nChMax = 4000;
dbf053f0 427
428 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
429 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
430
431 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
432 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
433
434
435 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
436 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
437
438 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
439 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
440 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
441 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
442
443
444 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
445 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
446 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
447
448 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
449 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
450 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
451 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
452 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
453 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
8a320ab9 454 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
455 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
456 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
a0d08b6d 457 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
dbf053f0 458
0341d0d8 459 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
460 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
aae7993b 461 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
0341d0d8 462
e2ca7519 463 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
464 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
465 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
466
dbf053f0 467 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
468 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
469 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
470 //
471
472
473 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
474 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
475 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
476 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
477
a0d08b6d 478 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
479 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);
480 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);
481 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);
482
483
484
dbf053f0 485 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
486 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
487 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
488 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
489
490 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
491 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
492 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
493 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
494
495 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
496 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
497 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
498 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
499
500
501
502 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
503 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
504 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
505 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
506 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
507 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
508 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
509 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
510 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
511 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
512 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
513 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
514
515 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
516 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
517 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
518 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
519
520 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
521 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
522 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
523 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
524
525
c4856527 526 if(fNRandomCones>0&&fUseBackgroundCalc){
d7ca0e14 527 for(int i = 0;i<3;i++){
528 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
529 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
530 }
531 }
c4856527 532
aae7993b 533 for(int i = 0;i < kMaxCent;i++){
534 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
535 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
536 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
537 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
538 }
dbf053f0 539
540 const Int_t saveLevel = 3; // large save level more histos
541 if(saveLevel>0){
542 fHistList->Add(fh1Xsec);
543 fHistList->Add(fh1Trials);
544
545 fHistList->Add(fh1NJetsRec);
546 fHistList->Add(fh1NConstRec);
547 fHistList->Add(fh1NConstLeadingRec);
548 fHistList->Add(fh1PtJetsRecIn);
549 fHistList->Add(fh1PtJetsLeadingRecIn);
550 fHistList->Add(fh1PtTracksRecIn);
551 fHistList->Add(fh1PtTracksLeadingRecIn);
552 fHistList->Add(fh1PtJetConstRec);
553 fHistList->Add(fh1PtJetConstLeadingRec);
554 fHistList->Add(fh1NJetsRecRan);
555 fHistList->Add(fh1NConstRecRan);
556 fHistList->Add(fh1PtJetsLeadingRecInRan);
557 fHistList->Add(fh1NConstLeadingRecRan);
558 fHistList->Add(fh1PtJetsRecInRan);
a0d08b6d 559 fHistList->Add(fh1Nch);
0341d0d8 560 fHistList->Add(fh1Centrality);
561 fHistList->Add(fh1CentralitySelect);
aae7993b 562 fHistList->Add(fh1CentralityPhySel);
e2ca7519 563 fHistList->Add(fh1Z);
564 fHistList->Add(fh1ZSelect);
565 fHistList->Add(fh1ZPhySel);
8f85b773 566 if(fNRandomCones>0&&fUseBackgroundCalc){
d7ca0e14 567 for(int i = 0;i<3;i++){
568 fHistList->Add(fh1BiARandomCones[i]);
569 fHistList->Add(fh1BiARandomConesRan[i]);
570 }
82a35b1a 571 }
aae7993b 572 for(int i = 0;i < kMaxCent;i++){
573 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
574 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
575 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
576 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
577 }
578
dbf053f0 579 fHistList->Add(fh2NRecJetsPt);
580 fHistList->Add(fh2NRecTracksPt);
581 fHistList->Add(fh2NConstPt);
582 fHistList->Add(fh2NConstLeadingPt);
a0d08b6d 583 fHistList->Add(fh2PtNch);
584 fHistList->Add(fh2PtNchRan);
585 fHistList->Add(fh2PtNchN);
586 fHistList->Add(fh2PtNchNRan);
dbf053f0 587 fHistList->Add(fh2JetPhiEta);
588 fHistList->Add(fh2LeadingJetPhiEta);
589 fHistList->Add(fh2JetEtaPt);
590 fHistList->Add(fh2LeadingJetEtaPt);
591 fHistList->Add(fh2TrackEtaPt);
592 fHistList->Add(fh2LeadingTrackEtaPt);
593 fHistList->Add(fh2JetsLeadingPhiEta );
594 fHistList->Add(fh2JetsLeadingPhiPt);
595 fHistList->Add(fh2TracksLeadingPhiEta);
596 fHistList->Add(fh2TracksLeadingPhiPt);
597 fHistList->Add(fh2TracksLeadingJetPhiPt);
598 fHistList->Add(fh2JetsLeadingPhiPtW);
599 fHistList->Add(fh2TracksLeadingPhiPtW);
600 fHistList->Add(fh2TracksLeadingJetPhiPtW);
601 fHistList->Add(fh2NRecJetsPtRan);
602 fHistList->Add(fh2NConstPtRan);
603 fHistList->Add(fh2NConstLeadingPtRan);
604 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
605 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
85587ff2 606 }
dbf053f0 607
608 // =========== Switch on Sumw2 for all histos ===========
609 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
610 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
611 if (h1){
612 h1->Sumw2();
613 continue;
614 }
615 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
616 if(hn)hn->Sumw2();
617 }
618 TH1::AddDirectory(oldStatus);
619}
620
621void AliAnalysisTaskJetCluster::Init()
622{
623 //
624 // Initialization
625 //
626
dbf053f0 627 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
628
629}
630
631void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
632{
633
634 if(fUseGlobalSelection){
635 // no selection by the service task, we continue
636 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
637 PostData(1, fHistList);
638 return;
639 }
640
0341d0d8 641
642
54424110 643 // handle and reset the output jet branch
c4856527 644
8f85b773 645 if(fTCAJetsOut)fTCAJetsOut->Delete();
646 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
647 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
648 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
649 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
54424110 650
a18eedbb 651 AliAODJetEventBackground* externalBackground = 0;
0341d0d8 652 if(!externalBackground&&fBackgroundBranch.Length()){
653 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
8f85b773 654 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
0341d0d8 655 }
dbf053f0 656 //
657 // Execute analysis for current event
658 //
659 AliESDEvent *fESD = 0;
660 if(fUseAODTrackInput){
661 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
662 if(!fAOD){
663 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
664 return;
665 }
666 // fethc the header
667 }
668 else{
669 // assume that the AOD is in the general output...
670 fAOD = AODEvent();
671 if(!fAOD){
672 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
673 return;
674 }
675 if(fDebug>0){
676 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
677 }
678 }
679
680 Bool_t selectEvent = false;
54424110 681 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
dbf053f0 682
aae7993b 683 Float_t cent = 0;
e2ca7519 684 Float_t zVtx = 0;
aae7993b 685 Int_t cenClass = -1;
dbf053f0 686 if(fAOD){
687 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
54424110 688 TString vtxTitle(vtxAOD->GetTitle());
e2ca7519 689 zVtx = vtxAOD->GetZ();
690
aae7993b 691 cent = fAOD->GetHeader()->GetCentrality();
692 if(cent<10)cenClass = 0;
693 else if(cent<30)cenClass = 1;
694 else if(cent<50)cenClass = 2;
695 else if(cent<80)cenClass = 3;
e2ca7519 696 if(physicsSelection){
697 fh1CentralityPhySel->Fill(cent);
698 fh1ZPhySel->Fill(zVtx);
699 }
700
54424110 701
702 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
dbf053f0 703 Float_t yvtx = vtxAOD->GetY();
704 Float_t xvtx = vtxAOD->GetX();
705 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
0996d685 706 if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
e2ca7519 707 if(physicsSelection){
708 selectEvent = true;
709 }
dbf053f0 710 }
711 }
aae7993b 712 if(fCentCutUp>0){
713 if(cent<fCentCutLo||cent>fCentCutUp){
714 selectEvent = false;
715 }
716 }
717
dbf053f0 718 }
719 if(!selectEvent){
720 PostData(1, fHistList);
721 return;
722 }
aae7993b 723 fh1Centrality->Fill(cent);
e2ca7519 724 fh1Z->Fill(zVtx);
dbf053f0 725 fh1Trials->Fill("#sum{ntrials}",1);
726
727
728 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
729
730 // ==== General variables needed
731
732
733
734 // we simply fetch the tracks/mc particles as a list of AliVParticles
735
736 TList recParticles;
737 TList genParticles;
738
739 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
a0d08b6d 740 Float_t nCh = recParticles.GetEntries();
741 fh1Nch->Fill(nCh);
dbf053f0 742 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
743 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
744 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
745
746 // find the jets....
747
748 vector<fastjet::PseudoJet> inputParticlesRec;
749 vector<fastjet::PseudoJet> inputParticlesRecRan;
d7ca0e14 750
751 // Generate the random cones
752
753 AliAODJet vTmpRan(1,0,0,1);
dbf053f0 754 for(int i = 0; i < recParticles.GetEntries(); i++){
755 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
756 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
d7ca0e14 757 // we take total momentum here
dbf053f0 758 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
759 jInp.set_user_index(i);
760 inputParticlesRec.push_back(jInp);
761
762 // the randomized input changes eta and phi, but keeps the p_T
bd80a748 763 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
764 Double_t pT = vp->Pt();
8a320ab9 765 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
d7ca0e14 766 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
bd80a748 767
d7ca0e14 768 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
bd80a748 769 Double_t pZ = pT/TMath::Tan(theta);
770
771 Double_t pX = pT * TMath::Cos(phi);
772 Double_t pY = pT * TMath::Sin(phi);
773 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
774 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
d7ca0e14 775
bd80a748 776 jInpRan.set_user_index(i);
777 inputParticlesRecRan.push_back(jInpRan);
d7ca0e14 778 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
bd80a748 779 }
54424110 780
781 // fill the tref array, only needed when we write out jets
8f85b773 782 if(fTCAJetsOut){
54424110 783 if(i == 0){
784 fRef->Delete(); // make sure to delete before placement new...
785 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
786 }
787 fRef->Add(vp);
788 }
d7ca0e14 789 }// recparticles
d7ca0e14 790
dbf053f0 791 if(inputParticlesRec.size()==0){
792 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
793 PostData(1, fHistList);
794 return;
795 }
796
797 // run fast jet
4dd6159d 798 // employ setters for these...
799
d6e66a82 800
9fd9f55b 801 // now create the object that holds info about ghosts
802 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
803 // reduce CPU time...
804 fGhostArea = 0.5;
805 fActiveAreaRepeats = 0;
806 }
d7ca0e14 807
808 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
4dd6159d 809 fastjet::AreaType areaType = fastjet::active_area;
810 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
82a35b1a 811 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
dbf053f0 812 vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
4dd6159d 813 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
dbf053f0 814
d6e66a82 815 //range where to compute background
816 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
817 phiMin = 0;
818 phiMax = 2*TMath::Pi();
819 rapMax = fGhostEtamax - fRparam;
820 rapMin = - fGhostEtamax + fRparam;
821 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
822
823
824
dbf053f0 825 inclusiveJets = clustSeq.inclusive_jets();
826 sortedJets = sorted_by_pt(inclusiveJets);
827
dbf053f0 828 fh1NJetsRec->Fill(sortedJets.size());
829
830 // loop over all jets an fill information, first on is the leading jet
831
54424110 832 Int_t nRecOver = inclusiveJets.size();
833 Int_t nRec = inclusiveJets.size();
834 if(inclusiveJets.size()>0){
835 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
aae7993b 836 Double_t area = clustSeq.area(sortedJets[0]);
0341d0d8 837 leadingJet.SetEffArea(area,0);
dbf053f0 838 Float_t pt = leadingJet.Pt();
54424110 839 Int_t nAodOutJets = 0;
840 Int_t nAodOutTracks = 0;
841 AliAODJet *aodOutJet = 0;
dbf053f0 842
843 Int_t iCount = 0;
844 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
845 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
846 while(pt<ptCut&&iCount<nRec){
847 nRecOver--;
848 iCount++;
849 if(iCount<nRec){
850 pt = sortedJets[iCount].perp();
851 }
852 }
853 if(nRecOver<=0)break;
854 fh2NRecJetsPt->Fill(ptCut,nRecOver);
855 }
856 Float_t phi = leadingJet.Phi();
857 if(phi<0)phi+=TMath::Pi()*2.;
858 Float_t eta = leadingJet.Eta();
0341d0d8 859 Float_t pTback = 0;
860 if(externalBackground){
861 // carefull has to be filled in a task before
862 // todo, ReArrange to the botom
d73e1768 863 pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
0341d0d8 864 }
865 pt = leadingJet.Pt() - pTback;
dbf053f0 866 // correlation of leading jet with tracks
867 TIterator *recIter = recParticles.MakeIterator();
dbf053f0 868 recIter->Reset();
225f0094 869 AliVParticle *tmpRecTrack = 0;
dbf053f0 870 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
871 Float_t tmpPt = tmpRecTrack->Pt();
872 // correlation
873 Float_t tmpPhi = tmpRecTrack->Phi();
874 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
875 Float_t dPhi = phi - tmpPhi;
876 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
877 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 878 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
879 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
aae7993b 880 if(cenClass>=0){
881 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
882 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
883 }
884
54424110 885 }
886
887
d73e1768 888 TLorentzVector vecareab;
4096256e 889 for(int j = 0; j < nRec;j++){
890 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
891 aodOutJet = 0;
892 nAodOutTracks = 0;
893 Float_t tmpPt = tmpRec.Pt();
fb10c4b8 894
8f85b773 895 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
896 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
fb10c4b8 897 Double_t area1 = clustSeq.area(sortedJets[j]);
898 aodOutJet->SetEffArea(area1,0);
d73e1768 899 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
900 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
901 aodOutJet->SetVectorAreaCharged(&vecareab);
902
903
fb10c4b8 904 }
905
906
4096256e 907 Float_t tmpPtBack = 0;
908 if(externalBackground){
909 // carefull has to be filled in a task before
2b018f7a 910 // todo, ReArrange to the botom
4096256e 911 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
912 }
913 tmpPt = tmpPt - tmpPtBack;
914 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
915
916 fh1PtJetsRecIn->Fill(tmpPt);
917 // Fill Spectra with constituents
918 vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
fb10c4b8 919
4096256e 920 fh1NConstRec->Fill(constituents.size());
921 fh2PtNch->Fill(nCh,tmpPt);
922 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
923 fh2NConstPt->Fill(tmpPt,constituents.size());
924 // loop over constiutents and fill spectrum
4096256e 925
54424110 926 for(unsigned int ic = 0; ic < constituents.size();ic++){
927 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
928 fh1PtJetConstRec->Fill(part->Pt());
929 if(aodOutJet){
930 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
931 }
932 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
933 }
54424110 934
935 // correlation
936 Float_t tmpPhi = tmpRec.Phi();
937 Float_t tmpEta = tmpRec.Eta();
c4856527 938 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
54424110 939 if(j==0){
940 fh1PtJetsLeadingRecIn->Fill(tmpPt);
941 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
942 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
943 fh1NConstLeadingRec->Fill(constituents.size());
944 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
945 continue;
946 }
947 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
948 fh2JetEtaPt->Fill(tmpEta,tmpPt);
949 Float_t dPhi = phi - tmpPhi;
950 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
951 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
952 Float_t dEta = eta - tmpRec.Eta();
953 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
954 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
aae7993b 955 if(cenClass>=0){
956 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
957 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
958 }
54424110 959 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
8f85b773 960 }// loop over reconstructed jets
54424110 961 delete recIter;
8f85b773 962
963
964
965 // Add the random cones...
966 if(fNRandomCones>0&&fTCARandomConesOut){
967 // create a random jet within the acceptance
968 Double_t etaMax = 0.8 - fRparam;
969 Int_t nCone = 0;
970 Int_t nConeRan = 0;
971 Double_t pTC = 1; // small number
972 for(int ir = 0;ir < fNRandomCones;ir++){
973 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
974 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
975 // massless jet
976 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
977 Double_t pZC = pTC/TMath::Tan(thetaC);
978 Double_t pXC = pTC * TMath::Cos(phiC);
979 Double_t pYC = pTC * TMath::Sin(phiC);
980 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
981 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
982 bool skip = false;
983 for(int jj = 0; jj < TMath::Min(nRec,2);jj++){// test for overlap with leading jets
984 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
985 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
986 skip = true;
987 break;
988 }
989 }
990 // test for overlap with previous cones to avoid double counting
991 for(int iic = 0;iic<ir;iic++){
992 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
993 if(iicone){
994 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
995 skip = true;
996 break;
997 }
998 }
999 }
1000 if(skip)continue;
1001 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
e309f70f 1002 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nConeRan++]) AliAODJet(tmpRecC);
1003 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nCone++]) AliAODJet(tmpRecC);
8f85b773 1004 }// loop over random cones creation
1005
1006
1007 // loop over the reconstructed particles and add up the pT in the random cones
1008 // maybe better to loop over randomized particles not in the real jets...
1009 // but this by definition brings dow average energy in the whole event
1010 AliAODJet vTmpRanR(1,0,0,1);
1011 for(int i = 0; i < recParticles.GetEntries(); i++){
1012 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1013 if(fTCARandomConesOut){
1014 for(int ir = 0;ir < fNRandomCones;ir++){
1015 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1016 if(jC&&jC->DeltaR(vp)<fRparam){
1017 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1018 }
1019 }
1020 }// add up energy in cone
1021
1022 // the randomized input changes eta and phi, but keeps the p_T
1023 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1024 Double_t pTR = vp->Pt();
8a320ab9 1025 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
8f85b773 1026 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1027
1028 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1029 Double_t pZR = pTR/TMath::Tan(thetaR);
1030
1031 Double_t pXR = pTR * TMath::Cos(phiR);
1032 Double_t pYR = pTR * TMath::Sin(phiR);
1033 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1034 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1035 if(fTCARandomConesOutRan){
1036 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1037 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1038 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1039 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1040 }
1041 }
1042 }
1043 }
1044 }// loop over recparticles
1045
1046 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1047 if(fTCARandomConesOut){
1048 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1049 // rescale the momntum vectors for the random cones
1050
1051 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1052 if(rC){
1053 Double_t etaC = rC->Eta();
1054 Double_t phiC = rC->Phi();
1055 // massless jet, unit vector
1056 pTC = rC->ChargedBgEnergy();
1057 if(pTC<=0)pTC = 0.1; // for almost empty events
1058 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1059 Double_t pZC = pTC/TMath::Tan(thetaC);
1060 Double_t pXC = pTC * TMath::Cos(phiC);
1061 Double_t pYC = pTC * TMath::Sin(phiC);
1062 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1063 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1064 rC->SetBgEnergy(0,0);
1065 rC->SetEffArea(jetArea,0);
1066 }
1067 }
1068 }
1069 if(!fTCARandomConesOutRan){
1070 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1071 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1072 // same wit random
1073 if(rC){
1074 Double_t etaC = rC->Eta();
1075 Double_t phiC = rC->Phi();
1076 // massless jet, unit vector
1077 pTC = rC->ChargedBgEnergy();
1078 if(pTC<=0)pTC = 0.1;// for almost empty events
1079 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1080 Double_t pZC = pTC/TMath::Tan(thetaC);
1081 Double_t pXC = pTC * TMath::Cos(phiC);
1082 Double_t pYC = pTC * TMath::Sin(phiC);
1083 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1084 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1085 rC->SetBgEnergy(0,0);
1086 rC->SetEffArea(jetArea,0);
1087 }
1088 }
1089 }
1090 }// if(fNRandomCones
d7ca0e14 1091
1092 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
d6e66a82 1093
4096256e 1094
8f85b773 1095
1096
1097
1098 if(fAODJetBackgroundOut){
d6e66a82 1099 vector<fastjet::PseudoJet> jets2=sortedJets;
1100 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1101 Double_t bkg1=0;
1102 Double_t sigma1=0.;
1103 Double_t meanarea1=0.;
1104 Double_t bkg2=0;
1105 Double_t sigma2=0.;
1106 Double_t meanarea2=0.;
82a35b1a 1107
d73e1768 1108 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
8f85b773 1109 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
f974a156 1110
1111 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1112 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
82a35b1a 1113
4a33f620 1114 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
8f85b773 1115 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
f974a156 1116 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1117 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1118
d6e66a82 1119 }
54424110 1120 }
d7ca0e14 1121
82a35b1a 1122
1123
1124
dbf053f0 1125
82a35b1a 1126 // fill track information
1127 Int_t nTrackOver = recParticles.GetSize();
dbf053f0 1128 // do the same for tracks and jets
82a35b1a 1129
1130 if(nTrackOver>0){
dbf053f0 1131 TIterator *recIter = recParticles.MakeIterator();
1132 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1133 Float_t pt = tmpRec->Pt();
82a35b1a 1134
dbf053f0 1135 // Printf("Leading track p_t %3.3E",pt);
1136 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1137 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1138 while(pt<ptCut&&tmpRec){
1139 nTrackOver--;
1140 tmpRec = (AliVParticle*)(recIter->Next());
1141 if(tmpRec){
1142 pt = tmpRec->Pt();
1143 }
1144 }
1145 if(nTrackOver<=0)break;
1146 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1147 }
1148
1149 recIter->Reset();
1150 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1151 Float_t phi = leading->Phi();
1152 if(phi<0)phi+=TMath::Pi()*2.;
1153 Float_t eta = leading->Eta();
1154 pt = leading->Pt();
1155 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1156 Float_t tmpPt = tmpRec->Pt();
1157 Float_t tmpEta = tmpRec->Eta();
1158 fh1PtTracksRecIn->Fill(tmpPt);
1159 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1160 if(tmpRec==leading){
1161 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1162 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1163 continue;
1164 }
1165 // correlation
1166 Float_t tmpPhi = tmpRec->Phi();
1167
1168 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1169 Float_t dPhi = phi - tmpPhi;
1170 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1171 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1172 Float_t dEta = eta - tmpRec->Eta();
1173 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1174 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1175 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1176 }
1177 delete recIter;
1178 }
1179
1180 // find the random jets
1181 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
d6e66a82 1182 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
dbf053f0 1183
1184 inclusiveJetsRan = clustSeqRan.inclusive_jets();
1185 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1186
1187 // fill the jet information from random track
1188
1189 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1190
1191 // loop over all jets an fill information, first on is the leading jet
1192
1193 Int_t nRecOverRan = inclusiveJetsRan.size();
1194 Int_t nRecRan = inclusiveJetsRan.size();
d73e1768 1195
dbf053f0 1196 if(inclusiveJetsRan.size()>0){
1197 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1198 Float_t pt = leadingJet.Pt();
1199
1200 Int_t iCount = 0;
d73e1768 1201 TLorentzVector vecarearanb;
1202
dbf053f0 1203 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1204 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1205 while(pt<ptCut&&iCount<nRecRan){
1206 nRecOverRan--;
1207 iCount++;
1208 if(iCount<nRecRan){
1209 pt = sortedJetsRan[iCount].perp();
1210 }
1211 }
1212 if(nRecOverRan<=0)break;
1213 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1214 }
1215 Float_t phi = leadingJet.Phi();
1216 if(phi<0)phi+=TMath::Pi()*2.;
1217 pt = leadingJet.Pt();
1218
1219 // correlation of leading jet with random tracks
1220
1221 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1222 {
1223 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1224 // correlation
1225 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1226 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1227 Float_t dPhi = phi - tmpPhi;
1228 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1229 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1230 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1231 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1232 }
1233
d6e66a82 1234 Int_t nAodOutJetsRan = 0;
1235 AliAODJet *aodOutJetRan = 0;
dbf053f0 1236 for(int j = 0; j < nRecRan;j++){
1237 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1238 Float_t tmpPt = tmpRec.Pt();
1239 fh1PtJetsRecInRan->Fill(tmpPt);
1240 // Fill Spectra with constituents
1241 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1242 fh1NConstRecRan->Fill(constituents.size());
1243 fh2NConstPtRan->Fill(tmpPt,constituents.size());
a0d08b6d 1244 fh2PtNchRan->Fill(nCh,tmpPt);
1245 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
d6e66a82 1246
1247
8f85b773 1248 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1249 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
d6e66a82 1250 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1251
d73e1768 1252 aodOutJetRan->SetEffArea(arearan,0);
1253 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1254 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1255 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1256
1257 }
d6e66a82 1258
1259
1260
1261
1262 for(unsigned int ic = 0; ic < constituents.size();ic++){
1263 // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1264 // fh1PtJetConstRec->Fill(part->Pt());
1265 if(aodOutJetRan){
1266 aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
d7ca0e14 1267 }
1268 }
d6e66a82 1269
1270
1271
1272
dbf053f0 1273 // correlation
1274 Float_t tmpPhi = tmpRec.Phi();
1275 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1276
1277 if(j==0){
1278 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1279 fh1NConstLeadingRecRan->Fill(constituents.size());
1280 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1281 continue;
1282 }
1283 }
d6e66a82 1284
1285
8f85b773 1286 if(fAODJetBackgroundOut){
d6e66a82 1287 Double_t bkg3=0.;
1288 Double_t sigma3=0.;
1289 Double_t meanarea3=0.;
4a33f620 1290 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
8f85b773 1291 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
d7ca0e14 1292 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1293 /*
82a35b1a 1294 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
9c6bd749 1295 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
d7ca0e14 1296 */
d6e66a82 1297 }
1298
1299
1300
dbf053f0 1301 }
0341d0d8 1302
1303
1304 // do the event selection if activated
1305 if(fJetTriggerPtCut>0){
1306 bool select = false;
2b018f7a 1307 Float_t minPt = fJetTriggerPtCut;
1308 /*
0341d0d8 1309 // hard coded for now ...
2b018f7a 1310 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
0341d0d8 1311 if(cent<10)minPt = 50;
aae7993b 1312 else if(cent<30)minPt = 42;
0341d0d8 1313 else if(cent<50)minPt = 28;
1314 else if(cent<80)minPt = 18;
2b018f7a 1315 */
0341d0d8 1316 float rho = 0;
1317 if(externalBackground)rho = externalBackground->GetBackground(2);
8f85b773 1318 if(fTCAJetsOut){
1319 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1320 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
0341d0d8 1321 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1322 if(ptSub>=minPt){
1323 select = true;
1324 break;
1325 }
1326 }
1327 }
1328
1329 if(select){
1330 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1331 fh1CentralitySelect->Fill(cent);
e2ca7519 1332 fh1ZSelect->Fill(zVtx);
0341d0d8 1333 aodH->SetFillAOD(kTRUE);
1334 }
1335 }
dbf053f0 1336
1337 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1338 PostData(1, fHistList);
1339}
1340
1341void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1342{
1343// Terminate analysis
1344//
1345 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1346}
1347
1348
1349Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1350
1351 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1352
1353 Int_t iCount = 0;
f840d64c 1354 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1355 if(type!=kTrackAODextraonly) {
1356 AliAODEvent *aod = 0;
1357 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1358 else aod = AODEvent();
1359 if(!aod){
1360 return iCount;
1361 }
1362 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1363 AliAODTrack *tr = aod->GetTrack(it);
1364 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
8a320ab9 1365 if(TMath::Abs(tr->Eta())>fTrackEtaWindow)continue;
f840d64c 1366 if(tr->Pt()<fTrackPtCut)continue;
1367 list->Add(tr);
1368 iCount++;
1369 }
dbf053f0 1370 }
f840d64c 1371 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1372 AliAODEvent *aod = 0;
1373 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1374 else aod = AODEvent();
1375
1376 if(!aod){
1377 return iCount;
1378 }
1379 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
8ddcf605 1380 if(!aodExtraTracks)return iCount;
f840d64c 1381 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1382 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1383 if (!track) continue;
1384
1385 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
8ddcf605 1386 if(!trackAOD)continue;
f840d64c 1387 if(trackAOD->Pt()<fTrackPtCut) continue;
f840d64c 1388 list->Add(trackAOD);
1389 iCount++;
1390 }
dbf053f0 1391 }
1392 }
1393 else if (type == kTrackKineAll||type == kTrackKineCharged){
1394 AliMCEvent* mcEvent = MCEvent();
1395 if(!mcEvent)return iCount;
1396 // we want to have alivpartilces so use get track
1397 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1398 if(!mcEvent->IsPhysicalPrimary(it))continue;
1399 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1400 if(type == kTrackKineAll){
0341d0d8 1401 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1402 list->Add(part);
1403 iCount++;
1404 }
1405 else if(type == kTrackKineCharged){
1406 if(part->Particle()->GetPDG()->Charge()==0)continue;
0341d0d8 1407 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1408 list->Add(part);
1409 iCount++;
1410 }
1411 }
1412 }
1413 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1414 AliAODEvent *aod = 0;
959508f4 1415 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
dbf053f0 1416 else aod = AODEvent();
1417 if(!aod)return iCount;
1418 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1419 if(!tca)return iCount;
1420 for(int it = 0;it < tca->GetEntriesFast();++it){
225f0094 1421 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
dbf053f0 1422 if(!part->IsPhysicalPrimary())continue;
1423 if(type == kTrackAODMCAll){
0341d0d8 1424 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1425 list->Add(part);
1426 iCount++;
1427 }
1428 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1429 if(part->Charge()==0)continue;
0341d0d8 1430 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1431 if(kTrackAODMCCharged){
1432 list->Add(part);
1433 }
1434 else {
8a320ab9 1435 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
dbf053f0 1436 list->Add(part);
1437 }
1438 iCount++;
1439 }
1440 }
1441 }// AODMCparticle
1442 list->Sort();
1443 return iCount;
dbf053f0 1444}
1445
1446/*
1447Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1448 for(int i = 0; i < particles.GetEntries(); i++){
1449 AliVParticle *vp = (AliVParticle*)particles.At(i);
1450 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1451 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1452 jInp.set_user_index(i);
1453 inputParticles.push_back(jInp);
1454 }
1455
1456 return 0;
1457
1458}
1459*/