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