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