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