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