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