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