Fixed memory leak when using random cones
[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
1000 AliAODJet *rCone = 0;
1001 AliAODJet *rConeRan = 0;
1002 if(fTCARandomConesOut)rConeRan = new ((*fTCARandomConesOut)[nConeRan++]) AliAODJet(tmpRecC);
1003 if(fTCARandomConesOutRan)rCone = new ((*fTCARandomConesOutRan)[nCone++]) AliAODJet(tmpRecC);
1004 }// loop over random cones creation
1005
1006
1007 // loop over the reconstructed particles and add up the pT in the random cones
1008 // maybe better to loop over randomized particles not in the real jets...
1009 // but this by definition brings dow average energy in the whole event
1010 AliAODJet vTmpRanR(1,0,0,1);
1011 for(int i = 0; i < recParticles.GetEntries(); i++){
1012 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1013 if(fTCARandomConesOut){
1014 for(int ir = 0;ir < fNRandomCones;ir++){
1015 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1016 if(jC&&jC->DeltaR(vp)<fRparam){
1017 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1018 }
1019 }
1020 }// add up energy in cone
1021
1022 // the randomized input changes eta and phi, but keeps the p_T
1023 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1024 Double_t pTR = vp->Pt();
1025 Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
1026 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1027
1028 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1029 Double_t pZR = pTR/TMath::Tan(thetaR);
1030
1031 Double_t pXR = pTR * TMath::Cos(phiR);
1032 Double_t pYR = pTR * TMath::Sin(phiR);
1033 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1034 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1035 if(fTCARandomConesOutRan){
1036 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1037 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1038 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1039 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1040 }
1041 }
1042 }
1043 }
1044 }// loop over recparticles
1045
1046 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1047 if(fTCARandomConesOut){
1048 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1049 // rescale the momntum vectors for the random cones
1050
1051 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1052 if(rC){
1053 Double_t etaC = rC->Eta();
1054 Double_t phiC = rC->Phi();
1055 // massless jet, unit vector
1056 pTC = rC->ChargedBgEnergy();
1057 if(pTC<=0)pTC = 0.1; // for almost empty events
1058 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1059 Double_t pZC = pTC/TMath::Tan(thetaC);
1060 Double_t pXC = pTC * TMath::Cos(phiC);
1061 Double_t pYC = pTC * TMath::Sin(phiC);
1062 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1063 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1064 rC->SetBgEnergy(0,0);
1065 rC->SetEffArea(jetArea,0);
1066 }
1067 }
1068 }
1069 if(!fTCARandomConesOutRan){
1070 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1071 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1072 // same wit random
1073 if(rC){
1074 Double_t etaC = rC->Eta();
1075 Double_t phiC = rC->Phi();
1076 // massless jet, unit vector
1077 pTC = rC->ChargedBgEnergy();
1078 if(pTC<=0)pTC = 0.1;// for almost empty events
1079 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1080 Double_t pZC = pTC/TMath::Tan(thetaC);
1081 Double_t pXC = pTC * TMath::Cos(phiC);
1082 Double_t pYC = pTC * TMath::Sin(phiC);
1083 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1084 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1085 rC->SetBgEnergy(0,0);
1086 rC->SetEffArea(jetArea,0);
1087 }
1088 }
1089 }
1090 }// if(fNRandomCones
d7ca0e14 1091
1092 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
d6e66a82 1093
4096256e 1094
8f85b773 1095
1096
1097
1098 if(fAODJetBackgroundOut){
d6e66a82 1099 vector<fastjet::PseudoJet> jets2=sortedJets;
1100 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1101 Double_t bkg1=0;
1102 Double_t sigma1=0.;
1103 Double_t meanarea1=0.;
1104 Double_t bkg2=0;
1105 Double_t sigma2=0.;
1106 Double_t meanarea2=0.;
82a35b1a 1107
d73e1768 1108 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
8f85b773 1109 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
f974a156 1110
1111 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1112 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
82a35b1a 1113
4a33f620 1114 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
8f85b773 1115 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
f974a156 1116 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1117 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1118
d6e66a82 1119 }
54424110 1120 }
d7ca0e14 1121
82a35b1a 1122
1123
1124
dbf053f0 1125
82a35b1a 1126 // fill track information
1127 Int_t nTrackOver = recParticles.GetSize();
dbf053f0 1128 // do the same for tracks and jets
82a35b1a 1129
1130 if(nTrackOver>0){
dbf053f0 1131 TIterator *recIter = recParticles.MakeIterator();
1132 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1133 Float_t pt = tmpRec->Pt();
82a35b1a 1134
dbf053f0 1135 // Printf("Leading track p_t %3.3E",pt);
1136 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1137 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1138 while(pt<ptCut&&tmpRec){
1139 nTrackOver--;
1140 tmpRec = (AliVParticle*)(recIter->Next());
1141 if(tmpRec){
1142 pt = tmpRec->Pt();
1143 }
1144 }
1145 if(nTrackOver<=0)break;
1146 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1147 }
1148
1149 recIter->Reset();
1150 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1151 Float_t phi = leading->Phi();
1152 if(phi<0)phi+=TMath::Pi()*2.;
1153 Float_t eta = leading->Eta();
1154 pt = leading->Pt();
1155 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1156 Float_t tmpPt = tmpRec->Pt();
1157 Float_t tmpEta = tmpRec->Eta();
1158 fh1PtTracksRecIn->Fill(tmpPt);
1159 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1160 if(tmpRec==leading){
1161 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1162 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1163 continue;
1164 }
1165 // correlation
1166 Float_t tmpPhi = tmpRec->Phi();
1167
1168 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1169 Float_t dPhi = phi - tmpPhi;
1170 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1171 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1172 Float_t dEta = eta - tmpRec->Eta();
1173 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1174 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1175 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1176 }
1177 delete recIter;
1178 }
1179
1180 // find the random jets
1181 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
d6e66a82 1182 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
dbf053f0 1183
1184 inclusiveJetsRan = clustSeqRan.inclusive_jets();
1185 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1186
1187 // fill the jet information from random track
1188
1189 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1190
1191 // loop over all jets an fill information, first on is the leading jet
1192
1193 Int_t nRecOverRan = inclusiveJetsRan.size();
1194 Int_t nRecRan = inclusiveJetsRan.size();
d73e1768 1195
dbf053f0 1196 if(inclusiveJetsRan.size()>0){
1197 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1198 Float_t pt = leadingJet.Pt();
1199
1200 Int_t iCount = 0;
d73e1768 1201 TLorentzVector vecarearanb;
1202
dbf053f0 1203 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1204 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1205 while(pt<ptCut&&iCount<nRecRan){
1206 nRecOverRan--;
1207 iCount++;
1208 if(iCount<nRecRan){
1209 pt = sortedJetsRan[iCount].perp();
1210 }
1211 }
1212 if(nRecOverRan<=0)break;
1213 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1214 }
1215 Float_t phi = leadingJet.Phi();
1216 if(phi<0)phi+=TMath::Pi()*2.;
1217 pt = leadingJet.Pt();
1218
1219 // correlation of leading jet with random tracks
1220
1221 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1222 {
1223 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1224 // correlation
1225 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1226 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1227 Float_t dPhi = phi - tmpPhi;
1228 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1229 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1230 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1231 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1232 }
1233
d6e66a82 1234 Int_t nAodOutJetsRan = 0;
1235 AliAODJet *aodOutJetRan = 0;
dbf053f0 1236 for(int j = 0; j < nRecRan;j++){
1237 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1238 Float_t tmpPt = tmpRec.Pt();
1239 fh1PtJetsRecInRan->Fill(tmpPt);
1240 // Fill Spectra with constituents
1241 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1242 fh1NConstRecRan->Fill(constituents.size());
1243 fh2NConstPtRan->Fill(tmpPt,constituents.size());
a0d08b6d 1244 fh2PtNchRan->Fill(nCh,tmpPt);
1245 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
d6e66a82 1246
1247
8f85b773 1248 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1249 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
d6e66a82 1250 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1251
d73e1768 1252 aodOutJetRan->SetEffArea(arearan,0);
1253 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1254 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1255 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1256
1257 }
d6e66a82 1258
1259
1260
1261
1262 for(unsigned int ic = 0; ic < constituents.size();ic++){
1263 // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1264 // fh1PtJetConstRec->Fill(part->Pt());
1265 if(aodOutJetRan){
1266 aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
d7ca0e14 1267 }
1268 }
d6e66a82 1269
1270
1271
1272
dbf053f0 1273 // correlation
1274 Float_t tmpPhi = tmpRec.Phi();
1275 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1276
1277 if(j==0){
1278 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1279 fh1NConstLeadingRecRan->Fill(constituents.size());
1280 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1281 continue;
1282 }
1283 }
d6e66a82 1284
1285
8f85b773 1286 if(fAODJetBackgroundOut){
d6e66a82 1287 Double_t bkg3=0.;
1288 Double_t sigma3=0.;
1289 Double_t meanarea3=0.;
4a33f620 1290 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
8f85b773 1291 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
d7ca0e14 1292 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1293 /*
82a35b1a 1294 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
9c6bd749 1295 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
d7ca0e14 1296 */
d6e66a82 1297 }
1298
1299
1300
dbf053f0 1301 }
0341d0d8 1302
1303
1304 // do the event selection if activated
1305 if(fJetTriggerPtCut>0){
1306 bool select = false;
2b018f7a 1307 Float_t minPt = fJetTriggerPtCut;
1308 /*
0341d0d8 1309 // hard coded for now ...
2b018f7a 1310 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
0341d0d8 1311 if(cent<10)minPt = 50;
aae7993b 1312 else if(cent<30)minPt = 42;
0341d0d8 1313 else if(cent<50)minPt = 28;
1314 else if(cent<80)minPt = 18;
2b018f7a 1315 */
0341d0d8 1316 float rho = 0;
1317 if(externalBackground)rho = externalBackground->GetBackground(2);
8f85b773 1318 if(fTCAJetsOut){
1319 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1320 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
0341d0d8 1321 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1322 if(ptSub>=minPt){
1323 select = true;
1324 break;
1325 }
1326 }
1327 }
1328
1329 if(select){
1330 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1331 fh1CentralitySelect->Fill(cent);
e2ca7519 1332 fh1ZSelect->Fill(zVtx);
0341d0d8 1333 aodH->SetFillAOD(kTRUE);
1334 }
1335 }
dbf053f0 1336
1337 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1338 PostData(1, fHistList);
1339}
1340
1341void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1342{
1343// Terminate analysis
1344//
1345 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1346}
1347
1348
1349Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1350
1351 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1352
1353 Int_t iCount = 0;
f840d64c 1354 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1355 if(type!=kTrackAODextraonly) {
1356 AliAODEvent *aod = 0;
1357 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1358 else aod = AODEvent();
1359 if(!aod){
1360 return iCount;
1361 }
1362 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1363 AliAODTrack *tr = aod->GetTrack(it);
1364 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1365 if(TMath::Abs(tr->Eta())>0.9)continue;
1366 if(tr->Pt()<fTrackPtCut)continue;
1367 list->Add(tr);
1368 iCount++;
1369 }
dbf053f0 1370 }
f840d64c 1371 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1372 AliAODEvent *aod = 0;
1373 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1374 else aod = AODEvent();
1375
1376 if(!aod){
1377 return iCount;
1378 }
1379 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
8ddcf605 1380 if(!aodExtraTracks)return iCount;
f840d64c 1381 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1382 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1383 if (!track) continue;
1384
1385 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
8ddcf605 1386 if(!trackAOD)continue;
f840d64c 1387 if(trackAOD->Pt()<fTrackPtCut) continue;
f840d64c 1388 list->Add(trackAOD);
1389 iCount++;
1390 }
dbf053f0 1391 }
1392 }
1393 else if (type == kTrackKineAll||type == kTrackKineCharged){
1394 AliMCEvent* mcEvent = MCEvent();
1395 if(!mcEvent)return iCount;
1396 // we want to have alivpartilces so use get track
1397 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1398 if(!mcEvent->IsPhysicalPrimary(it))continue;
1399 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1400 if(type == kTrackKineAll){
0341d0d8 1401 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1402 list->Add(part);
1403 iCount++;
1404 }
1405 else if(type == kTrackKineCharged){
1406 if(part->Particle()->GetPDG()->Charge()==0)continue;
0341d0d8 1407 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1408 list->Add(part);
1409 iCount++;
1410 }
1411 }
1412 }
1413 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1414 AliAODEvent *aod = 0;
959508f4 1415 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
dbf053f0 1416 else aod = AODEvent();
1417 if(!aod)return iCount;
1418 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1419 if(!tca)return iCount;
1420 for(int it = 0;it < tca->GetEntriesFast();++it){
225f0094 1421 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
dbf053f0 1422 if(!part->IsPhysicalPrimary())continue;
1423 if(type == kTrackAODMCAll){
0341d0d8 1424 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1425 list->Add(part);
1426 iCount++;
1427 }
1428 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1429 if(part->Charge()==0)continue;
0341d0d8 1430 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1431 if(kTrackAODMCCharged){
1432 list->Add(part);
1433 }
1434 else {
1435 if(TMath::Abs(part->Eta())>0.9)continue;
1436 list->Add(part);
1437 }
1438 iCount++;
1439 }
1440 }
1441 }// AODMCparticle
1442 list->Sort();
1443 return iCount;
dbf053f0 1444}
1445
1446/*
1447Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1448 for(int i = 0; i < particles.GetEntries(); i++){
1449 AliVParticle *vp = (AliVParticle*)particles.At(i);
1450 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1451 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1452 jInp.set_user_index(i);
1453 inputParticles.push_back(jInp);
1454 }
1455
1456 return 0;
1457
1458}
1459*/