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