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