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