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