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