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