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