Init the TList with refence objects when using the copy ctor, added protection agains...
[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
f840d64c 351 // we identify the extension aod event by looking for the branchname
54424110 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
d73e1768 880 pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
0341d0d8 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
d73e1768 905 TLorentzVector vecareab;
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);
d73e1768 916 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
917 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
918 aodOutJet->SetVectorAreaCharged(&vecareab);
919
920
fb10c4b8 921 }
922
923
4096256e 924 Float_t tmpPtBack = 0;
925 if(externalBackground){
926 // carefull has to be filled in a task before
2b018f7a 927 // todo, ReArrange to the botom
4096256e 928 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
929 }
930 tmpPt = tmpPt - tmpPtBack;
931 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
932
933 fh1PtJetsRecIn->Fill(tmpPt);
934 // Fill Spectra with constituents
935 vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
fb10c4b8 936
4096256e 937 fh1NConstRec->Fill(constituents.size());
938 fh2PtNch->Fill(nCh,tmpPt);
939 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
940 fh2NConstPt->Fill(tmpPt,constituents.size());
941 // loop over constiutents and fill spectrum
942
943 // Add the jet information and the track references to the output AOD
4096256e 944
54424110 945
c4856527 946 if(fNRandomCones>0){
4096256e 947 // create a random jet within the acceptance
c4856527 948 Double_t etaMax = 0.8 - fRparam;
4096256e 949 Int_t nCone = 0;
950 Int_t nConeRan = 0;
951 Double_t pTC = 1; // small number
952 for(int ir = 0;ir < fNRandomCones;ir++){
953 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
954 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
955 // massless jet
956 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
957 Double_t pZC = pTC/TMath::Tan(thetaC);
958 Double_t pXC = pTC * TMath::Cos(phiC);
959 Double_t pYC = pTC * TMath::Sin(phiC);
960 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
961 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
962 bool skip = false;
c4856527 963 for(int jj = 0; jj < TMath::Min(nRec,2);jj++){// test for overlap with leading jets
4096256e 964 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
965 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
966 skip = true;
967 break;
968 }
969 }
c4856527 970 // test for overlap with previous cones to avoid double counting
971 for(int iic = 0;iic<ir;iic++){
972 AliAODJet *iicone = (AliAODJet*)rConeArray->At(iic);
973 if(iicone){
974 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
975 skip = true;
976 break;
977 }
978 }
979 }
4096256e 980 if(skip)continue;
981 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
982 if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRecC);
983 if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRecC);
984 }// random cones
985
986
987 // loop over the reconstructed particles and add up the pT in the random cones
988 // maybe better to loop over randomized particles not in the real jets...
989 // but this by definition brings dow average energy in the whole event
990 AliAODJet vTmpRanR(1,0,0,1);
991 for(int i = 0; i < recParticles.GetEntries(); i++){
992 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
993 if(rConeArray){
994 for(int ir = 0;ir < fNRandomCones;ir++){
995 AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);
996 if(jC&&jC->DeltaR(vp)<fRparam){
997 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
998 }
999 }
1000 }// add up energy in cone
1001
1002 // the randomized input changes eta and phi, but keeps the p_T
1003 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1004 Double_t pTR = vp->Pt();
1005 Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
1006 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1007
1008 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1009 Double_t pZR = pTR/TMath::Tan(thetaR);
1010
1011 Double_t pXR = pTR * TMath::Cos(phiR);
1012 Double_t pYR = pTR * TMath::Sin(phiR);
1013 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1014 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1015 if(rConeArrayRan){
1016 for(int ir = 0;ir < fNRandomCones;ir++){
1017 AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);
1018 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1019 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1020 }
1021 }
1022 }
1023 }
1024 }// loop over recparticles
1025
1026 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1027 for(int ir = 0;ir < fNRandomCones;ir++){
1028 // rescale the momntum vectors for the random cones
1029 if(!rConeArray)continue;
1030 AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
1031 if(rC){
1032 Double_t etaC = rC->Eta();
1033 Double_t phiC = rC->Phi();
1034 // massless jet, unit vector
1035 pTC = rC->ChargedBgEnergy();
1036 if(pTC<=0)pTC = 0.1; // for almost empty events
1037 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1038 Double_t pZC = pTC/TMath::Tan(thetaC);
1039 Double_t pXC = pTC * TMath::Cos(phiC);
1040 Double_t pYC = pTC * TMath::Sin(phiC);
1041 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1042 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1043 rC->SetBgEnergy(0,0);
1044 rC->SetEffArea(jetArea,0);
1045 }
1046 rC = (AliAODJet*)rConeArrayRan->At(ir);
1047 // same wit random
1048 if(rC){
1049 Double_t etaC = rC->Eta();
1050 Double_t phiC = rC->Phi();
1051 // massless jet, unit vector
1052 pTC = rC->ChargedBgEnergy();
1053 if(pTC<=0)pTC = 0.1;// for almost empty events
1054 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1055 Double_t pZC = pTC/TMath::Tan(thetaC);
1056 Double_t pXC = pTC * TMath::Cos(phiC);
1057 Double_t pYC = pTC * TMath::Sin(phiC);
1058 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1059 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1060 rC->SetBgEnergy(0,0);
1061 rC->SetEffArea(jetArea,0);
1062 }
1063 }
1064 }// if(fUseBackgroundCalc
54424110 1065
1066 for(unsigned int ic = 0; ic < constituents.size();ic++){
1067 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1068 fh1PtJetConstRec->Fill(part->Pt());
1069 if(aodOutJet){
1070 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1071 }
1072 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1073 }
54424110 1074
1075 // correlation
1076 Float_t tmpPhi = tmpRec.Phi();
1077 Float_t tmpEta = tmpRec.Eta();
c4856527 1078 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
54424110 1079 if(j==0){
1080 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1081 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1082 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1083 fh1NConstLeadingRec->Fill(constituents.size());
1084 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1085 continue;
1086 }
1087 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1088 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1089 Float_t dPhi = phi - tmpPhi;
1090 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1091 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1092 Float_t dEta = eta - tmpRec.Eta();
1093 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1094 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
aae7993b 1095 if(cenClass>=0){
1096 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1097 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1098 }
54424110 1099 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
4096256e 1100 }
54424110 1101 delete recIter;
d7ca0e14 1102
1103 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
d6e66a82 1104
4096256e 1105
d6e66a82 1106 if(evBkg){
1107 vector<fastjet::PseudoJet> jets2=sortedJets;
1108 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1109 Double_t bkg1=0;
1110 Double_t sigma1=0.;
1111 Double_t meanarea1=0.;
1112 Double_t bkg2=0;
1113 Double_t sigma2=0.;
1114 Double_t meanarea2=0.;
82a35b1a 1115
d73e1768 1116 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
d6e66a82 1117 evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
f974a156 1118
1119 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1120 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
82a35b1a 1121
4a33f620 1122 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
d6e66a82 1123 evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
f974a156 1124 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1125 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1126
d6e66a82 1127 }
54424110 1128 }
d7ca0e14 1129
82a35b1a 1130
1131
1132
dbf053f0 1133
82a35b1a 1134 // fill track information
1135 Int_t nTrackOver = recParticles.GetSize();
dbf053f0 1136 // do the same for tracks and jets
82a35b1a 1137
1138 if(nTrackOver>0){
dbf053f0 1139 TIterator *recIter = recParticles.MakeIterator();
1140 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1141 Float_t pt = tmpRec->Pt();
82a35b1a 1142
dbf053f0 1143 // Printf("Leading track p_t %3.3E",pt);
1144 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1145 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1146 while(pt<ptCut&&tmpRec){
1147 nTrackOver--;
1148 tmpRec = (AliVParticle*)(recIter->Next());
1149 if(tmpRec){
1150 pt = tmpRec->Pt();
1151 }
1152 }
1153 if(nTrackOver<=0)break;
1154 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1155 }
1156
1157 recIter->Reset();
1158 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1159 Float_t phi = leading->Phi();
1160 if(phi<0)phi+=TMath::Pi()*2.;
1161 Float_t eta = leading->Eta();
1162 pt = leading->Pt();
1163 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1164 Float_t tmpPt = tmpRec->Pt();
1165 Float_t tmpEta = tmpRec->Eta();
1166 fh1PtTracksRecIn->Fill(tmpPt);
1167 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1168 if(tmpRec==leading){
1169 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1170 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1171 continue;
1172 }
1173 // correlation
1174 Float_t tmpPhi = tmpRec->Phi();
1175
1176 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1177 Float_t dPhi = phi - tmpPhi;
1178 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1179 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1180 Float_t dEta = eta - tmpRec->Eta();
1181 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1182 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1183 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1184 }
1185 delete recIter;
1186 }
1187
1188 // find the random jets
1189 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
d6e66a82 1190 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
dbf053f0 1191
1192 inclusiveJetsRan = clustSeqRan.inclusive_jets();
1193 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1194
1195 // fill the jet information from random track
1196
1197 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1198
1199 // loop over all jets an fill information, first on is the leading jet
1200
1201 Int_t nRecOverRan = inclusiveJetsRan.size();
1202 Int_t nRecRan = inclusiveJetsRan.size();
d73e1768 1203
dbf053f0 1204 if(inclusiveJetsRan.size()>0){
1205 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1206 Float_t pt = leadingJet.Pt();
1207
1208 Int_t iCount = 0;
d73e1768 1209 TLorentzVector vecarearanb;
1210
dbf053f0 1211 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1212 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1213 while(pt<ptCut&&iCount<nRecRan){
1214 nRecOverRan--;
1215 iCount++;
1216 if(iCount<nRecRan){
1217 pt = sortedJetsRan[iCount].perp();
1218 }
1219 }
1220 if(nRecOverRan<=0)break;
1221 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1222 }
1223 Float_t phi = leadingJet.Phi();
1224 if(phi<0)phi+=TMath::Pi()*2.;
1225 pt = leadingJet.Pt();
1226
1227 // correlation of leading jet with random tracks
1228
1229 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1230 {
1231 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1232 // correlation
1233 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1234 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1235 Float_t dPhi = phi - tmpPhi;
1236 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1237 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1238 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1239 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1240 }
1241
d6e66a82 1242 Int_t nAodOutJetsRan = 0;
1243 AliAODJet *aodOutJetRan = 0;
dbf053f0 1244 for(int j = 0; j < nRecRan;j++){
1245 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1246 Float_t tmpPt = tmpRec.Pt();
1247 fh1PtJetsRecInRan->Fill(tmpPt);
1248 // Fill Spectra with constituents
1249 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1250 fh1NConstRecRan->Fill(constituents.size());
1251 fh2NConstPtRan->Fill(tmpPt,constituents.size());
a0d08b6d 1252 fh2PtNchRan->Fill(nCh,tmpPt);
1253 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
d6e66a82 1254
1255
1256 if(tmpPt>fJetOutputMinPt&&jarrayran){
1257 aodOutJetRan = new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1258 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1259
d73e1768 1260 aodOutJetRan->SetEffArea(arearan,0);
1261 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1262 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1263 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1264
1265 }
d6e66a82 1266
1267
1268
1269
1270 for(unsigned int ic = 0; ic < constituents.size();ic++){
1271 // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1272 // fh1PtJetConstRec->Fill(part->Pt());
1273 if(aodOutJetRan){
1274 aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
d7ca0e14 1275 }
1276 }
d6e66a82 1277
1278
1279
1280
dbf053f0 1281 // correlation
1282 Float_t tmpPhi = tmpRec.Phi();
1283 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1284
1285 if(j==0){
1286 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1287 fh1NConstLeadingRecRan->Fill(constituents.size());
1288 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1289 continue;
1290 }
1291 }
d6e66a82 1292
1293
1294 if(evBkg){
1295 Double_t bkg3=0.;
1296 Double_t sigma3=0.;
1297 Double_t meanarea3=0.;
4a33f620 1298 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
d6e66a82 1299 evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
d7ca0e14 1300 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1301 /*
82a35b1a 1302 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
9c6bd749 1303 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
d7ca0e14 1304 */
d6e66a82 1305 }
1306
1307
1308
dbf053f0 1309 }
0341d0d8 1310
1311
1312 // do the event selection if activated
1313 if(fJetTriggerPtCut>0){
1314 bool select = false;
2b018f7a 1315 Float_t minPt = fJetTriggerPtCut;
1316 /*
0341d0d8 1317 // hard coded for now ...
2b018f7a 1318 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
0341d0d8 1319 if(cent<10)minPt = 50;
aae7993b 1320 else if(cent<30)minPt = 42;
0341d0d8 1321 else if(cent<50)minPt = 28;
1322 else if(cent<80)minPt = 18;
2b018f7a 1323 */
0341d0d8 1324 float rho = 0;
1325 if(externalBackground)rho = externalBackground->GetBackground(2);
2b018f7a 1326 if(jarray){
0341d0d8 1327 for(int i = 0;i < jarray->GetEntriesFast();i++){
1328 AliAODJet *jet = (AliAODJet*)jarray->At(i);
1329 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1330 if(ptSub>=minPt){
1331 select = true;
1332 break;
1333 }
1334 }
1335 }
1336
1337 if(select){
1338 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1339 fh1CentralitySelect->Fill(cent);
e2ca7519 1340 fh1ZSelect->Fill(zVtx);
0341d0d8 1341 aodH->SetFillAOD(kTRUE);
1342 }
1343 }
dbf053f0 1344
1345 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1346 PostData(1, fHistList);
1347}
1348
1349void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1350{
1351// Terminate analysis
1352//
1353 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1354}
1355
1356
1357Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1358
1359 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1360
1361 Int_t iCount = 0;
f840d64c 1362 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1363 if(type!=kTrackAODextraonly) {
1364 AliAODEvent *aod = 0;
1365 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1366 else aod = AODEvent();
1367 if(!aod){
1368 return iCount;
1369 }
1370 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1371 AliAODTrack *tr = aod->GetTrack(it);
1372 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1373 if(TMath::Abs(tr->Eta())>0.9)continue;
1374 if(tr->Pt()<fTrackPtCut)continue;
1375 list->Add(tr);
1376 iCount++;
1377 }
dbf053f0 1378 }
f840d64c 1379 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1380 AliAODEvent *aod = 0;
1381 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1382 else aod = AODEvent();
1383
1384 if(!aod){
1385 return iCount;
1386 }
1387 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1388 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1389 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1390 if (!track) continue;
1391
1392 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1393 if(trackAOD->Pt()<fTrackPtCut) continue;
1394
1395 list->Add(trackAOD);
1396 iCount++;
1397 }
dbf053f0 1398 }
1399 }
1400 else if (type == kTrackKineAll||type == kTrackKineCharged){
1401 AliMCEvent* mcEvent = MCEvent();
1402 if(!mcEvent)return iCount;
1403 // we want to have alivpartilces so use get track
1404 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1405 if(!mcEvent->IsPhysicalPrimary(it))continue;
1406 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1407 if(type == kTrackKineAll){
0341d0d8 1408 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1409 list->Add(part);
1410 iCount++;
1411 }
1412 else if(type == kTrackKineCharged){
1413 if(part->Particle()->GetPDG()->Charge()==0)continue;
0341d0d8 1414 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1415 list->Add(part);
1416 iCount++;
1417 }
1418 }
1419 }
1420 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1421 AliAODEvent *aod = 0;
959508f4 1422 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
dbf053f0 1423 else aod = AODEvent();
1424 if(!aod)return iCount;
1425 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1426 if(!tca)return iCount;
1427 for(int it = 0;it < tca->GetEntriesFast();++it){
225f0094 1428 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
dbf053f0 1429 if(!part->IsPhysicalPrimary())continue;
1430 if(type == kTrackAODMCAll){
0341d0d8 1431 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1432 list->Add(part);
1433 iCount++;
1434 }
1435 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1436 if(part->Charge()==0)continue;
0341d0d8 1437 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1438 if(kTrackAODMCCharged){
1439 list->Add(part);
1440 }
1441 else {
1442 if(TMath::Abs(part->Eta())>0.9)continue;
1443 list->Add(part);
1444 }
1445 iCount++;
1446 }
1447 }
1448 }// AODMCparticle
1449 list->Sort();
1450 return iCount;
dbf053f0 1451}
1452
1453/*
1454Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1455 for(int i = 0; i < particles.GetEntries(); i++){
1456 AliVParticle *vp = (AliVParticle*)particles.At(i);
1457 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1458 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1459 jInp.set_user_index(i);
1460 inputParticles.push_back(jInp);
1461 }
1462
1463 return 0;
1464
1465}
1466*/