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