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