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