]> 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
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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>
006b2a30 35#include <TF1.h>
dbf053f0 36#include <TList.h>
37#include <TLorentzVector.h>
38#include <TClonesArray.h>
39#include "TDatabasePDG.h"
8ad1ea82 40#include <TGrid.h>
dbf053f0 41
42#include "AliAnalysisTaskJetCluster.h"
43#include "AliAnalysisManager.h"
44#include "AliJetFinder.h"
45#include "AliJetHeader.h"
46#include "AliJetReader.h"
dbf053f0 47#include "AliESDEvent.h"
48#include "AliAODEvent.h"
49#include "AliAODHandler.h"
1240edf5 50#include "AliAODExtension.h"
dbf053f0 51#include "AliAODTrack.h"
52#include "AliAODJet.h"
53#include "AliAODMCParticle.h"
54#include "AliMCEventHandler.h"
55#include "AliMCEvent.h"
56#include "AliStack.h"
57#include "AliGenPythiaEventHeader.h"
58#include "AliJetKineReaderHeader.h"
59#include "AliGenCocktailEventHeader.h"
60#include "AliInputEventHandler.h"
d6e66a82 61#include "AliAODJetEventBackground.h"
dbf053f0 62
63#include "fastjet/PseudoJet.hh"
64#include "fastjet/ClusterSequenceArea.hh"
65#include "fastjet/AreaDefinition.hh"
66#include "fastjet/JetDefinition.hh"
67// get info on how fastjet was configured
68#include "fastjet/config.h"
69
393a3bd4 70using std::vector;
dbf053f0 71
72ClassImp(AliAnalysisTaskJetCluster)
73
54424110 74AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
f0659f11 75 //
76 // Destructor
77 //
006b2a30 78
54424110 79 delete fRef;
d7ca0e14 80 delete fRandom;
8f85b773 81
82 if(fTCAJetsOut)fTCAJetsOut->Delete();
83 delete fTCAJetsOut;
006b2a30 84
8f85b773 85 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
86 delete fTCAJetsOutRan;
006b2a30 87
8f85b773 88 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
89 delete fTCARandomConesOut;
006b2a30 90
8f85b773 91 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
92 delete fTCARandomConesOutRan;
006b2a30 93
8f85b773 94 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
95 delete fAODJetBackgroundOut;
54424110 96}
97
b2c1e955 98AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
99 AliAnalysisTaskSE(),
dbf053f0 100 fAOD(0x0),
54424110 101 fAODExtension(0x0),
102 fRef(new TRefArray),
dbf053f0 103 fUseAODTrackInput(kFALSE),
104 fUseAODMCInput(kFALSE),
9fd9f55b 105 fUseBackgroundCalc(kFALSE),
b2c1e955 106 fEventSelection(kFALSE),
d66558be 107 fRequireVZEROAC(kFALSE),
108 fRequireTZEROvtx(kFALSE),
1f7b59f6 109 fUseHFcuts(kFALSE),
dbf053f0 110 fFilterMask(0),
d7396b04 111 fFilterMaskBestPt(0),
b2c1e955 112 fFilterType(0),
bea24d2d 113 fJetTypes(kJet),
dbf053f0 114 fTrackTypeRec(kTrackUndef),
bd80a748 115 fTrackTypeGen(kTrackUndef),
116 fNSkipLeadingRan(0),
af399bc2 117 fNSkipLeadingCone(0),
c4856527 118 fNRandomCones(0),
bd80a748 119 fAvgTrials(1),
8a320ab9 120 fExternalWeight(1),
121 fTrackEtaWindow(0.9),
c18b64ea 122 fRequireITSRefit(0),
27c67edf 123 fApplySharedClusterCut(0),
dbf053f0 124 fRecEtaWindow(0.5),
54424110 125 fTrackPtCut(0.),
617932d0 126 fJetOutputMinPt(0.150),
65c43de8 127 fMaxTrackPtInJet(100.),
0341d0d8 128 fJetTriggerPtCut(0),
ca6d12f9 129 fVtxZCut(8),
130 fVtxR2Cut(1),
aae7993b 131 fCentCutUp(0),
132 fCentCutLo(0),
7bedcd31 133 fStoreRhoLeadingTrackCorr(kFALSE),
54424110 134 fNonStdBranch(""),
0341d0d8 135 fBackgroundBranch(""),
54424110 136 fNonStdFile(""),
006b2a30 137 fMomResH1(0x0),
138 fMomResH2(0x0),
139 fMomResH3(0x0),
140 fMomResH1Fit(0x0),
141 fMomResH2Fit(0x0),
142 fMomResH3Fit(0x0),
143 fhEffH1(0x0),
144 fhEffH2(0x0),
145 fhEffH3(0x0),
fa7c34ba 146 fUseTrPtResolutionSmearing(kFALSE),
006b2a30 147 fUseDiceEfficiency(kFALSE),
da8f7918 148 fDiceEfficiencyMinPt(-1.),
fa7c34ba 149 fUseTrPtResolutionFromOADB(kFALSE),
150 fUseTrEfficiencyFromOADB(kFALSE),
151 fPathTrPtResolution(""),
152 fPathTrEfficiency(""),
153 fChangeEfficiencyFraction(0.),
da8f7918 154 fEfficiencyFixed(1.),
dbf053f0 155 fRparam(1.0),
156 fAlgorithm(fastjet::kt_algorithm),
157 fStrategy(fastjet::Best),
158 fRecombScheme(fastjet::BIpt_scheme),
159 fAreaType(fastjet::active_area),
d6e66a82 160 fGhostArea(0.01),
161 fActiveAreaRepeats(1),
162 fGhostEtamax(1.5),
8f85b773 163 fTCAJetsOut(0x0),
164 fTCAJetsOutRan(0x0),
165 fTCARandomConesOut(0x0),
166 fTCARandomConesOutRan(0x0),
167 fAODJetBackgroundOut(0x0),
d7ca0e14 168 fRandom(0),
dbf053f0 169 fh1Xsec(0x0),
170 fh1Trials(0x0),
171 fh1PtHard(0x0),
172 fh1PtHardNoW(0x0),
173 fh1PtHardTrials(0x0),
174 fh1NJetsRec(0x0),
175 fh1NConstRec(0x0),
176 fh1NConstLeadingRec(0x0),
177 fh1PtJetsRecIn(0x0),
178 fh1PtJetsLeadingRecIn(0x0),
179 fh1PtJetConstRec(0x0),
180 fh1PtJetConstLeadingRec(0x0),
181 fh1PtTracksRecIn(0x0),
182 fh1PtTracksLeadingRecIn(0x0),
183 fh1NJetsRecRan(0x0),
184 fh1NConstRecRan(0x0),
185 fh1PtJetsLeadingRecInRan(0x0),
186 fh1NConstLeadingRecRan(0x0),
187 fh1PtJetsRecInRan(0x0),
188 fh1PtTracksGenIn(0x0),
a0d08b6d 189 fh1Nch(0x0),
aae7993b 190 fh1CentralityPhySel(0x0),
0341d0d8 191 fh1Centrality(0x0),
192 fh1CentralitySelect(0x0),
e2ca7519 193 fh1ZPhySel(0x0),
194 fh1Z(0x0),
195 fh1ZSelect(0x0),
dbf053f0 196 fh2NRecJetsPt(0x0),
197 fh2NRecTracksPt(0x0),
198 fh2NConstPt(0x0),
199 fh2NConstLeadingPt(0x0),
200 fh2JetPhiEta(0x0),
201 fh2LeadingJetPhiEta(0x0),
202 fh2JetEtaPt(0x0),
203 fh2LeadingJetEtaPt(0x0),
204 fh2TrackEtaPt(0x0),
205 fh2LeadingTrackEtaPt(0x0),
206 fh2JetsLeadingPhiEta(0x0),
207 fh2JetsLeadingPhiPt(0x0),
208 fh2TracksLeadingPhiEta(0x0),
209 fh2TracksLeadingPhiPt(0x0),
210 fh2TracksLeadingJetPhiPt(0x0),
211 fh2JetsLeadingPhiPtW(0x0),
212 fh2TracksLeadingPhiPtW(0x0),
213 fh2TracksLeadingJetPhiPtW(0x0),
214 fh2NRecJetsPtRan(0x0),
215 fh2NConstPtRan(0x0),
216 fh2NConstLeadingPtRan(0x0),
a0d08b6d 217 fh2PtNch(0x0),
218 fh2PtNchRan(0x0),
219 fh2PtNchN(0x0),
220 fh2PtNchNRan(0x0),
dbf053f0 221 fh2TracksLeadingJetPhiPtRan(0x0),
222 fh2TracksLeadingJetPhiPtWRan(0x0),
cbfccad4 223 fh3CentvsRhoLeadingTrackPt(0x0),
224 fh3CentvsSigmaLeadingTrackPt(0x0),
225 fh3MultvsRhoLeadingTrackPt(0x0),
226 fh3MultvsSigmaLeadingTrackPt(0x0),
7bedcd31 227 fh3CentvsRhoLeadingTrackPtQ1(0x0),
228 fh3CentvsRhoLeadingTrackPtQ2(0x0),
229 fh3CentvsRhoLeadingTrackPtQ3(0x0),
230 fh3CentvsRhoLeadingTrackPtQ4(0x0),
231 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
232 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
233 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
234 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
235 fh3MultvsRhoLeadingTrackPtQ1(0x0),
236 fh3MultvsRhoLeadingTrackPtQ2(0x0),
237 fh3MultvsRhoLeadingTrackPtQ3(0x0),
238 fh3MultvsRhoLeadingTrackPtQ4(0x0),
239 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
240 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
241 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
242 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
cbfccad4 243 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
244 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
245 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
246 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
006b2a30 247 fh2PtGenPtSmeared(0x0),
248 fp1Efficiency(0x0),
249 fp1PtResolution(0x0),
dbf053f0 250 fHistList(0x0)
251{
f0659f11 252 //
253 // Constructor
254 //
255
82a35b1a 256 for(int i = 0;i<3;i++){
257 fh1BiARandomCones[i] = 0;
258 fh1BiARandomConesRan[i] = 0;
259 }
aae7993b 260 for(int i = 0;i<kMaxCent;i++){
261 fh2JetsLeadingPhiPtC[i] = 0;
262 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
263 fh2TracksLeadingJetPhiPtC[i] = 0;
264 fh2TracksLeadingJetPhiPtWC[i] = 0;
265 }
dbf053f0 266}
267
268AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
269 AliAnalysisTaskSE(name),
270 fAOD(0x0),
54424110 271 fAODExtension(0x0),
272 fRef(new TRefArray),
dbf053f0 273 fUseAODTrackInput(kFALSE),
274 fUseAODMCInput(kFALSE),
9fd9f55b 275 fUseBackgroundCalc(kFALSE),
d66558be 276 fEventSelection(kFALSE),
277 fRequireVZEROAC(kFALSE),
278 fRequireTZEROvtx(kFALSE),
1f7b59f6 279 fUseHFcuts(kFALSE),
d66558be 280 fFilterMask(0),
d7396b04 281 fFilterMaskBestPt(0),
b2c1e955 282 fFilterType(0),
bea24d2d 283 fJetTypes(kJet),
dbf053f0 284 fTrackTypeRec(kTrackUndef),
285 fTrackTypeGen(kTrackUndef),
bd80a748 286 fNSkipLeadingRan(0),
af399bc2 287 fNSkipLeadingCone(0),
c4856527 288 fNRandomCones(0),
dbf053f0 289 fAvgTrials(1),
290 fExternalWeight(1),
8a320ab9 291 fTrackEtaWindow(0.9),
c18b64ea 292 fRequireITSRefit(0),
27c67edf 293 fApplySharedClusterCut(0),
dbf053f0 294 fRecEtaWindow(0.5),
54424110 295 fTrackPtCut(0.),
617932d0 296 fJetOutputMinPt(0.150),
65c43de8 297 fMaxTrackPtInJet(100.),
0341d0d8 298 fJetTriggerPtCut(0),
ca6d12f9 299 fVtxZCut(8),
300 fVtxR2Cut(1),
aae7993b 301 fCentCutUp(0),
302 fCentCutLo(0),
7bedcd31 303 fStoreRhoLeadingTrackCorr(kFALSE),
54424110 304 fNonStdBranch(""),
0341d0d8 305 fBackgroundBranch(""),
54424110 306 fNonStdFile(""),
006b2a30 307 fMomResH1(0x0),
308 fMomResH2(0x0),
309 fMomResH3(0x0),
310 fMomResH1Fit(0x0),
311 fMomResH2Fit(0x0),
312 fMomResH3Fit(0x0),
313 fhEffH1(0x0),
314 fhEffH2(0x0),
315 fhEffH3(0x0),
fa7c34ba 316 fUseTrPtResolutionSmearing(kFALSE),
006b2a30 317 fUseDiceEfficiency(kFALSE),
da8f7918 318 fDiceEfficiencyMinPt(-1.),
fa7c34ba 319 fUseTrPtResolutionFromOADB(kFALSE),
320 fUseTrEfficiencyFromOADB(kFALSE),
321 fPathTrPtResolution(""),
322 fPathTrEfficiency(""),
323 fChangeEfficiencyFraction(0.),
da8f7918 324 fEfficiencyFixed(1.),
dbf053f0 325 fRparam(1.0),
326 fAlgorithm(fastjet::kt_algorithm),
327 fStrategy(fastjet::Best),
328 fRecombScheme(fastjet::BIpt_scheme),
329 fAreaType(fastjet::active_area),
d6e66a82 330 fGhostArea(0.01),
331 fActiveAreaRepeats(1),
332 fGhostEtamax(1.5),
8f85b773 333 fTCAJetsOut(0x0),
334 fTCAJetsOutRan(0x0),
335 fTCARandomConesOut(0x0),
336 fTCARandomConesOutRan(0x0),
337 fAODJetBackgroundOut(0x0),
d7ca0e14 338 fRandom(0),
dbf053f0 339 fh1Xsec(0x0),
340 fh1Trials(0x0),
341 fh1PtHard(0x0),
342 fh1PtHardNoW(0x0),
343 fh1PtHardTrials(0x0),
344 fh1NJetsRec(0x0),
345 fh1NConstRec(0x0),
346 fh1NConstLeadingRec(0x0),
347 fh1PtJetsRecIn(0x0),
348 fh1PtJetsLeadingRecIn(0x0),
349 fh1PtJetConstRec(0x0),
350 fh1PtJetConstLeadingRec(0x0),
351 fh1PtTracksRecIn(0x0),
352 fh1PtTracksLeadingRecIn(0x0),
353 fh1NJetsRecRan(0x0),
354 fh1NConstRecRan(0x0),
355 fh1PtJetsLeadingRecInRan(0x0),
356 fh1NConstLeadingRecRan(0x0),
357 fh1PtJetsRecInRan(0x0),
358 fh1PtTracksGenIn(0x0),
a0d08b6d 359 fh1Nch(0x0),
aae7993b 360 fh1CentralityPhySel(0x0),
0341d0d8 361 fh1Centrality(0x0),
362 fh1CentralitySelect(0x0),
e2ca7519 363 fh1ZPhySel(0x0),
364 fh1Z(0x0),
365 fh1ZSelect(0x0),
dbf053f0 366 fh2NRecJetsPt(0x0),
367 fh2NRecTracksPt(0x0),
368 fh2NConstPt(0x0),
369 fh2NConstLeadingPt(0x0),
370 fh2JetPhiEta(0x0),
371 fh2LeadingJetPhiEta(0x0),
372 fh2JetEtaPt(0x0),
373 fh2LeadingJetEtaPt(0x0),
374 fh2TrackEtaPt(0x0),
375 fh2LeadingTrackEtaPt(0x0),
376 fh2JetsLeadingPhiEta(0x0),
377 fh2JetsLeadingPhiPt(0x0),
378 fh2TracksLeadingPhiEta(0x0),
379 fh2TracksLeadingPhiPt(0x0),
380 fh2TracksLeadingJetPhiPt(0x0),
381 fh2JetsLeadingPhiPtW(0x0),
382 fh2TracksLeadingPhiPtW(0x0),
383 fh2TracksLeadingJetPhiPtW(0x0),
384 fh2NRecJetsPtRan(0x0),
385 fh2NConstPtRan(0x0),
386 fh2NConstLeadingPtRan(0x0),
a0d08b6d 387 fh2PtNch(0x0),
388 fh2PtNchRan(0x0),
389 fh2PtNchN(0x0),
390 fh2PtNchNRan(0x0),
dbf053f0 391 fh2TracksLeadingJetPhiPtRan(0x0),
392 fh2TracksLeadingJetPhiPtWRan(0x0),
cbfccad4 393 fh3CentvsRhoLeadingTrackPt(0x0),
394 fh3CentvsSigmaLeadingTrackPt(0x0),
395 fh3MultvsRhoLeadingTrackPt(0x0),
396 fh3MultvsSigmaLeadingTrackPt(0x0),
7bedcd31 397 fh3CentvsRhoLeadingTrackPtQ1(0x0),
398 fh3CentvsRhoLeadingTrackPtQ2(0x0),
399 fh3CentvsRhoLeadingTrackPtQ3(0x0),
400 fh3CentvsRhoLeadingTrackPtQ4(0x0),
401 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
402 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
403 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
404 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
405 fh3MultvsRhoLeadingTrackPtQ1(0x0),
406 fh3MultvsRhoLeadingTrackPtQ2(0x0),
407 fh3MultvsRhoLeadingTrackPtQ3(0x0),
408 fh3MultvsRhoLeadingTrackPtQ4(0x0),
409 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
410 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
411 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
412 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
cbfccad4 413 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
414 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
415 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
416 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
006b2a30 417 fh2PtGenPtSmeared(0x0),
418 fp1Efficiency(0x0),
419 fp1PtResolution(0x0),
dbf053f0 420 fHistList(0x0)
421{
f0659f11 422 //
423 // named ctor
424 //
006b2a30 425
82a35b1a 426 for(int i = 0;i<3;i++){
427 fh1BiARandomCones[i] = 0;
428 fh1BiARandomConesRan[i] = 0;
429 }
aae7993b 430 for(int i = 0;i<kMaxCent;i++){
431 fh2JetsLeadingPhiPtC[i] = 0;
432 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
433 fh2TracksLeadingJetPhiPtC[i] = 0;
434 fh2TracksLeadingJetPhiPtWC[i] = 0;
435 }
8ad1ea82 436 DefineOutput(1, TList::Class());
dbf053f0 437}
438
439
440
441Bool_t AliAnalysisTaskJetCluster::Notify()
442{
443 //
444 // Implemented Notify() to read the cross sections
445 // and number of trials from pyxsec.root
446 //
447 return kTRUE;
448}
449
450void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
451{
452
453 //
454 // Create the output container
455 //
456
d7ca0e14 457 fRandom = new TRandom3(0);
54424110 458
459
dbf053f0 460 // Connect the AOD
461
462
463 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
464
d6e66a82 465
54424110 466
467 if(fNonStdBranch.Length()!=0)
468 {
469 // only create the output branch if we have a name
470 // Create a new branch for jets...
471 // -> cleared in the UserExec....
472 // here we can also have the case that the brnaches are written to a separate file
d508795b 473
bea24d2d 474 if(fJetTypes&kJet){
475 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
476 fTCAJetsOut->SetName(fNonStdBranch.Data());
477 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
478 }
006b2a30 479
bea24d2d 480 if(fJetTypes&kJetRan){
481 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
482 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
da8f7918 483 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
484 if( fEfficiencyFixed < 1.)
485 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dEffFixed%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
486 else
487 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
488 }
bea24d2d 489 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
490 }
491
9fd9f55b 492 if(fUseBackgroundCalc){
493 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
8f85b773 494 fAODJetBackgroundOut = new AliAODJetEventBackground();
495 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
da8f7918 496 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
497 if( fEfficiencyFixed < 1.)
498 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dEffFixed%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
499 else
500 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
501 }
8f85b773 502 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
9fd9f55b 503 }
c4856527 504 }
505 // create the branch for the random cones with the same R
af399bc2 506 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
da8f7918 507 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
508 if( fEfficiencyFixed < 1.)
509 cName = Form("%sDetector%d%dEffFixed%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.),fNSkipLeadingCone);
510 else
511 cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
512 }
c4856527 513 if(fNRandomCones>0){
bea24d2d 514 if(fJetTypes&kRC){
515 if(!AODEvent()->FindListObject(cName.Data())){
516 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
517 fTCARandomConesOut->SetName(cName.Data());
518 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
519 }
d7ca0e14 520 }
d7ca0e14 521 // create the branch with the random for the random cones on the random event
bea24d2d 522 if(fJetTypes&kRCRan){
523 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
524 if(!AODEvent()->FindListObject(cName.Data())){
525 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
526 fTCARandomConesOutRan->SetName(cName.Data());
527 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
528 }
d7ca0e14 529 }
9fd9f55b 530 }
c4856527 531
54424110 532 if(fNonStdFile.Length()!=0){
533 //
534 // case that we have an AOD extension we need to fetch the jets from the extended output
f840d64c 535 // we identify the extension aod event by looking for the branchname
54424110 536 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
8f85b773 537 // case that we have an AOD extension we need can fetch the background maybe from the extended output
538 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
54424110 539 }
540 }
541
542
dbf053f0 543 if(!fHistList)fHistList = new TList();
9c6bd749 544 fHistList->SetOwner();
750e22e8 545 PostData(1, fHistList); // post data in any case once
546
dbf053f0 547 Bool_t oldStatus = TH1::AddDirectoryStatus();
548 TH1::AddDirectory(kFALSE);
549
550 //
551 // Histogram
552
fe80c230 553 const Int_t nBinPt = 100;
dbf053f0 554 Double_t binLimitsPt[nBinPt+1];
555 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
556 if(iPt == 0){
557 binLimitsPt[iPt] = 0.0;
558 }
559 else {// 1.0
fe80c230 560 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
dbf053f0 561 }
562 }
563
564 const Int_t nBinPhi = 90;
565 Double_t binLimitsPhi[nBinPhi+1];
566 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
567 if(iPhi==0){
568 binLimitsPhi[iPhi] = -1.*TMath::Pi();
569 }
570 else{
571 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
572 }
573 }
574
575
576
577 const Int_t nBinEta = 40;
578 Double_t binLimitsEta[nBinEta+1];
579 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
580 if(iEta==0){
581 binLimitsEta[iEta] = -2.0;
582 }
583 else{
584 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
585 }
586 }
587
fe80c230 588 const int nChMax = 5000;
dbf053f0 589
590 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
591 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
592
593 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
594 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
595
596
597 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
598 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
599
600 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
601 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
602 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
603 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
604
605
606 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
607 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
608 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
609
610 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
611 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
612 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
613 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
614 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
615 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
8a320ab9 616 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
617 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
618 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
a0d08b6d 619 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
dbf053f0 620
0341d0d8 621 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
622 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
aae7993b 623 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
0341d0d8 624
e2ca7519 625 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
626 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
627 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
628
dbf053f0 629 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
630 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
631 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
632 //
633
634
635 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
636 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
637 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
638 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
639
a0d08b6d 640 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
641 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);
642 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);
643 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);
644
645
646
dbf053f0 647 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
648 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
649 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
650 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
651
652 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
653 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
654 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
655 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
656
657 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
8ad1ea82 658 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
dbf053f0 659 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
8ad1ea82 660 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
dbf053f0 661
662
663
664 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
8ad1ea82 665 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
dbf053f0 666 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
8ad1ea82 667 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
dbf053f0 668 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
669 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
670 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
8ad1ea82 671 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
dbf053f0 672 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
8ad1ea82 673 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
dbf053f0 674 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
8ad1ea82 675 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
dbf053f0 676
677 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
8ad1ea82 678 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
dbf053f0 679 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
8ad1ea82 680 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
dbf053f0 681
682 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
683 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
684 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
8ad1ea82 685 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
dbf053f0 686
7bedcd31 687 if(fStoreRhoLeadingTrackCorr) {
cbfccad4 688 fh3CentvsRhoLeadingTrackPt = new TH3F("fh3CentvsRhoLeadingTrackPt","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.,100,0.,100.);
689 fh3CentvsSigmaLeadingTrackPt = new TH3F("fh3CentvsSigmaLeadingTrackPt","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.,100,0.,100.);
690 fh3MultvsRhoLeadingTrackPt = new TH3F("fh3MultvsRhoLeadingTrackPt","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.,100,0.,100.);
691 fh3MultvsSigmaLeadingTrackPt = new TH3F("fh3MultvsSigmaLeadingTrackPt","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.,100,0.,100.);
7bedcd31 692
693
694 fh3CentvsRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
695 fh3CentvsRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
696 fh3CentvsRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
697 fh3CentvsRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
698
699 fh3CentvsSigmaLeadingTrackPtQ1 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ1","centrality vs background density Q1; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
700 fh3CentvsSigmaLeadingTrackPtQ2 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ2","centrality vs background density Q2; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
701 fh3CentvsSigmaLeadingTrackPtQ3 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ3","centrality vs background density Q3; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
702 fh3CentvsSigmaLeadingTrackPtQ4 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ4","centrality vs background density Q4; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
703
704 fh3MultvsRhoLeadingTrackPtQ1 = new TH3F("fh3MultvsRhoLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
705 fh3MultvsRhoLeadingTrackPtQ2 = new TH3F("fh3MultvsRhoLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
706 fh3MultvsRhoLeadingTrackPtQ3 = new TH3F("fh3MultvsRhoLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
707 fh3MultvsRhoLeadingTrackPtQ4 = new TH3F("fh3MultvsRhoLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
708
709 fh3MultvsSigmaLeadingTrackPtQ1 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
710 fh3MultvsSigmaLeadingTrackPtQ2 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
711 fh3MultvsSigmaLeadingTrackPtQ3 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
712 fh3MultvsSigmaLeadingTrackPtQ4 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
713
cbfccad4 714
715 fh3CentvsDeltaRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
716 fh3CentvsDeltaRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
717 fh3CentvsDeltaRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
718 fh3CentvsDeltaRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
719
720 fHistList->Add(fh3CentvsRhoLeadingTrackPt);
721 fHistList->Add(fh3CentvsSigmaLeadingTrackPt);
722 fHistList->Add(fh3MultvsRhoLeadingTrackPt);
723 fHistList->Add(fh3MultvsSigmaLeadingTrackPt);
7bedcd31 724
725 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ1);
726 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ2);
727 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ3);
728 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ4);
729
730 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ1);
731 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ2);
732 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ3);
733 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ4);
734
735 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ1);
736 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ2);
737 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ3);
738 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ4);
739
740 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ1);
741 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ2);
742 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3);
743 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4);
744
cbfccad4 745 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ1);
746 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ2);
747 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ3);
748 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ4);
749
7bedcd31 750 }
c0fdfe53 751
006b2a30 752 //Detector level effects histos
753 fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
754
755 fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
756 fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
757
758 fHistList->Add(fh2PtGenPtSmeared);
759 fHistList->Add(fp1Efficiency);
760 fHistList->Add(fp1PtResolution);
dbf053f0 761
c4856527 762 if(fNRandomCones>0&&fUseBackgroundCalc){
d7ca0e14 763 for(int i = 0;i<3;i++){
764 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
765 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
766 }
767 }
c4856527 768
aae7993b 769 for(int i = 0;i < kMaxCent;i++){
770 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
771 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
772 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
773 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
774 }
dbf053f0 775
776 const Int_t saveLevel = 3; // large save level more histos
777 if(saveLevel>0){
778 fHistList->Add(fh1Xsec);
779 fHistList->Add(fh1Trials);
780
781 fHistList->Add(fh1NJetsRec);
782 fHistList->Add(fh1NConstRec);
783 fHistList->Add(fh1NConstLeadingRec);
784 fHistList->Add(fh1PtJetsRecIn);
785 fHistList->Add(fh1PtJetsLeadingRecIn);
786 fHistList->Add(fh1PtTracksRecIn);
787 fHistList->Add(fh1PtTracksLeadingRecIn);
788 fHistList->Add(fh1PtJetConstRec);
789 fHistList->Add(fh1PtJetConstLeadingRec);
790 fHistList->Add(fh1NJetsRecRan);
791 fHistList->Add(fh1NConstRecRan);
792 fHistList->Add(fh1PtJetsLeadingRecInRan);
793 fHistList->Add(fh1NConstLeadingRecRan);
794 fHistList->Add(fh1PtJetsRecInRan);
a0d08b6d 795 fHistList->Add(fh1Nch);
0341d0d8 796 fHistList->Add(fh1Centrality);
797 fHistList->Add(fh1CentralitySelect);
aae7993b 798 fHistList->Add(fh1CentralityPhySel);
e2ca7519 799 fHistList->Add(fh1Z);
800 fHistList->Add(fh1ZSelect);
801 fHistList->Add(fh1ZPhySel);
8f85b773 802 if(fNRandomCones>0&&fUseBackgroundCalc){
d7ca0e14 803 for(int i = 0;i<3;i++){
804 fHistList->Add(fh1BiARandomCones[i]);
805 fHistList->Add(fh1BiARandomConesRan[i]);
806 }
82a35b1a 807 }
8ad1ea82 808 for(int i = 0;i < kMaxCent;i++){
809 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
810 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
811 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
812 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
813 }
aae7993b 814
dbf053f0 815 fHistList->Add(fh2NRecJetsPt);
816 fHistList->Add(fh2NRecTracksPt);
817 fHistList->Add(fh2NConstPt);
818 fHistList->Add(fh2NConstLeadingPt);
a0d08b6d 819 fHistList->Add(fh2PtNch);
820 fHistList->Add(fh2PtNchRan);
821 fHistList->Add(fh2PtNchN);
822 fHistList->Add(fh2PtNchNRan);
dbf053f0 823 fHistList->Add(fh2JetPhiEta);
824 fHistList->Add(fh2LeadingJetPhiEta);
825 fHistList->Add(fh2JetEtaPt);
826 fHistList->Add(fh2LeadingJetEtaPt);
827 fHistList->Add(fh2TrackEtaPt);
828 fHistList->Add(fh2LeadingTrackEtaPt);
829 fHistList->Add(fh2JetsLeadingPhiEta );
830 fHistList->Add(fh2JetsLeadingPhiPt);
831 fHistList->Add(fh2TracksLeadingPhiEta);
832 fHistList->Add(fh2TracksLeadingPhiPt);
833 fHistList->Add(fh2TracksLeadingJetPhiPt);
834 fHistList->Add(fh2JetsLeadingPhiPtW);
835 fHistList->Add(fh2TracksLeadingPhiPtW);
836 fHistList->Add(fh2TracksLeadingJetPhiPtW);
837 fHistList->Add(fh2NRecJetsPtRan);
838 fHistList->Add(fh2NConstPtRan);
839 fHistList->Add(fh2NConstLeadingPtRan);
840 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
841 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
85587ff2 842 }
dbf053f0 843
844 // =========== Switch on Sumw2 for all histos ===========
845 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
846 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
847 if (h1){
848 h1->Sumw2();
849 continue;
850 }
851 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
852 if(hn)hn->Sumw2();
853 }
854 TH1::AddDirectory(oldStatus);
855}
856
fa7c34ba 857void AliAnalysisTaskJetCluster::LocalInit()
dbf053f0 858{
859 //
860 // Initialization
861 //
862
dbf053f0 863 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
864
fa7c34ba 865 if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
866 if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB();
867
868
006b2a30 869 FitMomentumResolution();
870
dbf053f0 871}
872
873void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
874{
875
54424110 876 // handle and reset the output jet branch
c4856527 877
8f85b773 878 if(fTCAJetsOut)fTCAJetsOut->Delete();
879 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
880 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
881 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
882 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
54424110 883
a18eedbb 884 AliAODJetEventBackground* externalBackground = 0;
0341d0d8 885 if(!externalBackground&&fBackgroundBranch.Length()){
886 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
750e22e8 887 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
8f85b773 888 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
0341d0d8 889 }
dbf053f0 890 //
891 // Execute analysis for current event
892 //
893 AliESDEvent *fESD = 0;
894 if(fUseAODTrackInput){
895 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
896 if(!fAOD){
897 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
898 return;
899 }
006b2a30 900 // fetch the header
dbf053f0 901 }
902 else{
903 // assume that the AOD is in the general output...
904 fAOD = AODEvent();
905 if(!fAOD){
906 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
907 return;
908 }
909 if(fDebug>0){
910 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
911 }
912 }
006b2a30 913
914 //Check if information is provided detector level effects
fa7c34ba 915 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
b0b15bae 916 if( fEfficiencyFixed < 1. ) {
917 if (!fUseDiceEfficiency)
918 fUseDiceEfficiency = 1; // 1 is the default; 2 can be set by user
919 }
920 else {
921 if(!fhEffH1 || !fhEffH2 || !fhEffH3 ) fUseDiceEfficiency = kFALSE;
922 }
da8f7918 923
dbf053f0 924 Bool_t selectEvent = false;
54424110 925 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
dbf053f0 926
aae7993b 927 Float_t cent = 0;
e2ca7519 928 Float_t zVtx = 0;
aae7993b 929 Int_t cenClass = -1;
dbf053f0 930 if(fAOD){
931 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
54424110 932 TString vtxTitle(vtxAOD->GetTitle());
e2ca7519 933 zVtx = vtxAOD->GetZ();
934
46cc237f 935 cent = ((AliVAODHeader*)fAOD->GetHeader())->GetCentrality();
aae7993b 936 if(cent<10)cenClass = 0;
937 else if(cent<30)cenClass = 1;
938 else if(cent<50)cenClass = 2;
939 else if(cent<80)cenClass = 3;
e2ca7519 940 if(physicsSelection){
941 fh1CentralityPhySel->Fill(cent);
942 fh1ZPhySel->Fill(zVtx);
943 }
944
ca6d12f9 945 if(fEventSelection){
946 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
dbf053f0 947 Float_t yvtx = vtxAOD->GetY();
948 Float_t xvtx = vtxAOD->GetX();
949 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
ca6d12f9 950 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
e2ca7519 951 if(physicsSelection){
952 selectEvent = true;
953 }
dbf053f0 954 }
aae7993b 955 }
ca6d12f9 956 if(fCentCutUp>0){
957 if(cent<fCentCutLo||cent>fCentCutUp){
958 selectEvent = false;
959 }
960 }
961 }else{
962 selectEvent = true;
aae7993b 963 }
dbf053f0 964 }
d66558be 965
966
967 Bool_t T0 = false;
968 Bool_t V0 = false;
969 const AliAODVZERO *vzero = fAOD->GetVZEROData();
970 if(vzero){
971 if((vzero->GetTriggerChargeA()>0)&&(vzero->GetTriggerChargeC()>0)){
972 V0 = true;
973 }
974 }
ca6d12f9 975
d66558be 976 const AliAODTZERO *tzero = fAOD->GetTZEROData();
977 if(tzero){
978 if(TMath::Abs(tzero->GetT0VertexRaw())<100){
979 T0 = true;
980 }
981 }
982
983 if(fRequireVZEROAC&&fRequireTZEROvtx)selectEvent = selectEvent&&V0&&T0;
984 else if(fRequireTZEROvtx)selectEvent = selectEvent&&T0;
985 else if(fRequireVZEROAC)selectEvent = selectEvent&&V0;
986
ca6d12f9 987
dbf053f0 988 if(!selectEvent){
989 PostData(1, fHistList);
990 return;
991 }
aae7993b 992 fh1Centrality->Fill(cent);
e2ca7519 993 fh1Z->Fill(zVtx);
dbf053f0 994 fh1Trials->Fill("#sum{ntrials}",1);
995
996
997 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
998
999 // ==== General variables needed
1000
1001
1002
1003 // we simply fetch the tracks/mc particles as a list of AliVParticles
1004
1005 TList recParticles;
1006 TList genParticles;
1007
1008 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
a0d08b6d 1009 Float_t nCh = recParticles.GetEntries();
b0b15bae 1010 Float_t nGen=genParticles.GetEntries();
a0d08b6d 1011 fh1Nch->Fill(nCh);
dbf053f0 1012 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
1013 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
1014 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
1015
1016 // find the jets....
1017
1018 vector<fastjet::PseudoJet> inputParticlesRec;
1019 vector<fastjet::PseudoJet> inputParticlesRecRan;
d7ca0e14 1020
1021 // Generate the random cones
1022
1023 AliAODJet vTmpRan(1,0,0,1);
dbf053f0 1024 for(int i = 0; i < recParticles.GetEntries(); i++){
1025 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
006b2a30 1026
dbf053f0 1027 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
d7ca0e14 1028 // we take total momentum here
006b2a30 1029
fa7c34ba 1030 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
006b2a30 1031 //Add particles to fastjet in case we are not running toy model
1032 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1033 jInp.set_user_index(i);
1034 inputParticlesRec.push_back(jInp);
1035 }
1036 else if(fUseDiceEfficiency) {
1037
1038 // Dice to decide if particle is kept or not - toy model for efficiency
1039 //
da8f7918 1040 Double_t sumEff = 0.;
1041 Double_t pT = 0.;
006b2a30 1042 Double_t eff[3] = {0.};
006b2a30 1043 Int_t cat[3] = {0};
da8f7918 1044
1045 Double_t rnd = fRandom->Uniform(1.);
1046 if( fEfficiencyFixed<1. ) {
1047 sumEff = fEfficiencyFixed;
b0b15bae 1048 if (fUseDiceEfficiency == 2) {
1049 sumEff = (nCh - fEfficiencyFixed*nGen) / nCh; // rescale eff; fEfficiencyFixed is wrt to nGen, but dicing is fraction of nCh
1050 }
da8f7918 1051 } else {
1052
1053 pT = vp->Pt();
1054 Double_t pTtmp = pT;
1055 if(pT>10.) pTtmp = 10.;
1056 Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
1057 Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
1058 Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
1059
1060 //Sort efficiencies from large to small
1061 if(eff1>eff2 && eff1>eff3) {
1062 eff[0] = eff1;
1063 cat[0] = 1;
1064 if(eff2>eff3) {
1065 eff[1] = eff2;
1066 eff[2] = eff3;
1067 cat[1] = 2;
1068 cat[2] = 3;
1069 } else {
1070 eff[1] = eff3;
1071 eff[2] = eff2;
1072 cat[1] = 3;
1073 cat[2] = 2;
1074 }
006b2a30 1075 }
da8f7918 1076 else if(eff2>eff1 && eff2>eff3) {
1077 eff[0] = eff2;
1078 cat[0] = 2;
1079 if(eff1>eff3) {
1080 eff[1] = eff1;
1081 eff[2] = eff3;
1082 cat[1] = 1;
1083 cat[2] = 3;
1084 } else {
1085 eff[1] = eff3;
1086 eff[2] = eff1;
1087 cat[1] = 3;
1088 cat[2] = 1;
1089 }
006b2a30 1090 }
da8f7918 1091 else if(eff3>eff1 && eff3>eff2) {
1092 eff[0] = eff3;
1093 cat[0] = 3;
1094 if(eff1>eff2) {
1095 eff[1] = eff1;
1096 eff[2] = eff2;
1097 cat[1] = 1;
1098 cat[2] = 2;
1099 } else {
1100 eff[1] = eff2;
1101 eff[2] = eff1;
1102 cat[1] = 2;
1103 cat[2] = 1;
1104 }
006b2a30 1105 }
da8f7918 1106
1107 sumEff = eff[0]+eff[1]+eff[2];
006b2a30 1108 }
006b2a30 1109 fp1Efficiency->Fill(vp->Pt(),sumEff);
4bd66a92 1110 if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
006b2a30 1111
fa7c34ba 1112 if(fUseTrPtResolutionSmearing) {
006b2a30 1113 //Smear momentum of generated particle
1114 Double_t smear = 1.;
1115 //Select hybrid track category
1116 if(rnd<=eff[2])
1117 smear = GetMomentumSmearing(cat[2],pT);
1118 else if(rnd<=(eff[2]+eff[1]))
1119 smear = GetMomentumSmearing(cat[1],pT);
1120 else
1121 smear = GetMomentumSmearing(cat[0],pT);
1122
1123 fp1PtResolution->Fill(vp->Pt(),smear);
1124
1125 Double_t sigma = vp->Pt()*smear;
1126 Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
1127
1128 Double_t phi = vp->Phi();
1129 Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
1130 Double_t pX = pTrec * TMath::Cos(phi);
1131 Double_t pY = pTrec * TMath::Sin(phi);
1132 Double_t pZ = pTrec/TMath::Tan(theta);
1133 Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
1134
1135 fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
1136
1137 fastjet::PseudoJet jInp(pX,pY,pZ,p);
1138 jInp.set_user_index(i);
1139 inputParticlesRec.push_back(jInp);
1140
1141 }
1142 else {
1143 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1144 jInp.set_user_index(i);
1145 inputParticlesRec.push_back(jInp);
1146
1147 }
1148
1149 }
dbf053f0 1150
1151 // the randomized input changes eta and phi, but keeps the p_T
bd80a748 1152 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1153 Double_t pT = vp->Pt();
8a320ab9 1154 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
d7ca0e14 1155 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
bd80a748 1156
d7ca0e14 1157 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
bd80a748 1158 Double_t pZ = pT/TMath::Tan(theta);
1159
1160 Double_t pX = pT * TMath::Cos(phi);
1161 Double_t pY = pT * TMath::Sin(phi);
1162 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
1163 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
d7ca0e14 1164
bd80a748 1165 jInpRan.set_user_index(i);
1166 inputParticlesRecRan.push_back(jInpRan);
d7ca0e14 1167 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
bd80a748 1168 }
54424110 1169
1170 // fill the tref array, only needed when we write out jets
8f85b773 1171 if(fTCAJetsOut){
54424110 1172 if(i == 0){
1173 fRef->Delete(); // make sure to delete before placement new...
fa7c34ba 1174 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
006b2a30 1175 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
1176 }
54424110 1177 }
fa7c34ba 1178 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
54424110 1179 }
d7ca0e14 1180 }// recparticles
d7ca0e14 1181
dbf053f0 1182 if(inputParticlesRec.size()==0){
1183 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1184 PostData(1, fHistList);
1185 return;
1186 }
1187
1188 // run fast jet
4dd6159d 1189 // employ setters for these...
1190
d6e66a82 1191
9fd9f55b 1192 // now create the object that holds info about ghosts
bea24d2d 1193 /*
8ad1ea82 1194 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
9fd9f55b 1195 // reduce CPU time...
1196 fGhostArea = 0.5;
1197 fActiveAreaRepeats = 0;
8ad1ea82 1198 }
bea24d2d 1199 */
1200
1201 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
4dd6159d 1202 fastjet::AreaType areaType = fastjet::active_area;
1203 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
82a35b1a 1204 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
4dd6159d 1205 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
dbf053f0 1206
d6e66a82 1207 //range where to compute background
1208 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1209 phiMin = 0;
1210 phiMax = 2*TMath::Pi();
1211 rapMax = fGhostEtamax - fRparam;
1212 rapMin = - fGhostEtamax + fRparam;
1213 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1214
1215
fe80c230 1216 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1217 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
d6e66a82 1218
dbf053f0 1219
dbf053f0 1220 fh1NJetsRec->Fill(sortedJets.size());
1221
8ad1ea82 1222 // loop over all jets an fill information, first on is the leading jet
dbf053f0 1223
54424110 1224 Int_t nRecOver = inclusiveJets.size();
1225 Int_t nRec = inclusiveJets.size();
1226 if(inclusiveJets.size()>0){
1227 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
aae7993b 1228 Double_t area = clustSeq.area(sortedJets[0]);
0341d0d8 1229 leadingJet.SetEffArea(area,0);
dbf053f0 1230 Float_t pt = leadingJet.Pt();
54424110 1231 Int_t nAodOutJets = 0;
1232 Int_t nAodOutTracks = 0;
1233 AliAODJet *aodOutJet = 0;
dbf053f0 1234
1235 Int_t iCount = 0;
1236 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1237 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1238 while(pt<ptCut&&iCount<nRec){
1239 nRecOver--;
1240 iCount++;
1241 if(iCount<nRec){
1242 pt = sortedJets[iCount].perp();
1243 }
1244 }
1245 if(nRecOver<=0)break;
1246 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1247 }
1248 Float_t phi = leadingJet.Phi();
1249 if(phi<0)phi+=TMath::Pi()*2.;
1250 Float_t eta = leadingJet.Eta();
0341d0d8 1251 Float_t pTback = 0;
1252 if(externalBackground){
1253 // carefull has to be filled in a task before
006b2a30 1254 // todo, ReArrange to the botom
8ad1ea82 1255 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
0341d0d8 1256 }
1257 pt = leadingJet.Pt() - pTback;
dbf053f0 1258 // correlation of leading jet with tracks
1259 TIterator *recIter = recParticles.MakeIterator();
dbf053f0 1260 recIter->Reset();
225f0094 1261 AliVParticle *tmpRecTrack = 0;
dbf053f0 1262 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1263 Float_t tmpPt = tmpRecTrack->Pt();
1264 // correlation
1265 Float_t tmpPhi = tmpRecTrack->Phi();
1266 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1267 Float_t dPhi = phi - tmpPhi;
1268 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1269 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1270 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1271 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
aae7993b 1272 if(cenClass>=0){
1273 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1274 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1275 }
1276
54424110 1277 }
1278
1279
d73e1768 1280 TLorentzVector vecareab;
4096256e 1281 for(int j = 0; j < nRec;j++){
1282 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1283 aodOutJet = 0;
1284 nAodOutTracks = 0;
1285 Float_t tmpPt = tmpRec.Pt();
fb10c4b8 1286
8f85b773 1287 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1288 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
d508795b 1289 aodOutJet->GetRefTracks()->Clear();
fb10c4b8 1290 Double_t area1 = clustSeq.area(sortedJets[j]);
1291 aodOutJet->SetEffArea(area1,0);
d73e1768 1292 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1293 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1294 aodOutJet->SetVectorAreaCharged(&vecareab);
fb10c4b8 1295 }
1296
1297
4096256e 1298 Float_t tmpPtBack = 0;
1299 if(externalBackground){
1300 // carefull has to be filled in a task before
8ad1ea82 1301 // todo, ReArrange to the botom
4096256e 1302 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1303 }
1304 tmpPt = tmpPt - tmpPtBack;
1305 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1306
1307 fh1PtJetsRecIn->Fill(tmpPt);
ea950747 1308 // Fill Spectra with constituentsemacs
fe80c230 1309 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
fb10c4b8 1310
4096256e 1311 fh1NConstRec->Fill(constituents.size());
1312 fh2PtNch->Fill(nCh,tmpPt);
1313 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1314 fh2NConstPt->Fill(tmpPt,constituents.size());
1315 // loop over constiutents and fill spectrum
d508795b 1316
d7396b04 1317 AliVParticle *partLead = 0;
1318 Float_t ptLead = -1;
1319
d508795b 1320 for(unsigned int ic = 0; ic < constituents.size();ic++){
1321 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
006b2a30 1322 if(!part) continue;
d508795b 1323 fh1PtJetConstRec->Fill(part->Pt());
1324 if(aodOutJet){
fa7c34ba 1325 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
d7396b04 1326 if(part->Pt()>fMaxTrackPtInJet){
1327 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1328 }
1329 }
1330 if(part->Pt()>ptLead){
95808c2b 1331 ptLead = part->Pt();
d7396b04 1332 partLead = part;
d508795b 1333 }
1334 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1335 }
d7396b04 1336
1337 AliAODTrack *aodT = 0;
47373258 1338 if(partLead){
1339 if(aodOutJet){
1340 //set pT of leading constituent of jet
1341 aodOutJet->SetPtLeading(partLead->Pt());
1342 aodT = dynamic_cast<AliAODTrack*>(partLead);
1343 if(aodT){
1344 if(aodT->TestFilterBit(fFilterMaskBestPt)){
1345 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1346 }
d7396b04 1347 }
1348 }
1349 }
47373258 1350
8ad1ea82 1351 // correlation
1352 Float_t tmpPhi = tmpRec.Phi();
1353 Float_t tmpEta = tmpRec.Eta();
1354 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1355 if(j==0){
1356 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1357 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1358 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1359 fh1NConstLeadingRec->Fill(constituents.size());
1360 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1361 continue;
1362 }
1363 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1364 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1365 Float_t dPhi = phi - tmpPhi;
1366 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1367 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1368 Float_t dEta = eta - tmpRec.Eta();
1369 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1370 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1371 if(cenClass>=0){
1372 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1373 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1374 }
1375 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
8f85b773 1376 }// loop over reconstructed jets
8ad1ea82 1377 delete recIter;
1378
1379
1380
1381 // Add the random cones...
1382 if(fNRandomCones>0&&fTCARandomConesOut){
1383 // create a random jet within the acceptance
1384 Double_t etaMax = fTrackEtaWindow - fRparam;
1385 Int_t nCone = 0;
1386 Int_t nConeRan = 0;
1387 Double_t pTC = 1; // small number
1388 for(int ir = 0;ir < fNRandomCones;ir++){
1389 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1390 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1391 // massless jet
1392 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1393 Double_t pZC = pTC/TMath::Tan(thetaC);
1394 Double_t pXC = pTC * TMath::Cos(phiC);
1395 Double_t pYC = pTC * TMath::Sin(phiC);
1396 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1397 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1398 bool skip = false;
1399 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1400 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1401 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1402 skip = true;
1403 break;
1404 }
1405 }
1406 // test for overlap with previous cones to avoid double counting
1407 for(int iic = 0;iic<ir;iic++){
1408 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1409 if(iicone){
1410 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1411 skip = true;
1412 break;
1413 }
1414 }
1415 }
1416 if(skip)continue;
1417 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1418 tmpRecC.SetPtLeading(-1.);
1419 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1420 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1421 }// loop over random cones creation
8f85b773 1422
1423
8ad1ea82 1424 // loop over the reconstructed particles and add up the pT in the random cones
1425 // maybe better to loop over randomized particles not in the real jets...
1426 // but this by definition brings dow average energy in the whole event
1427 AliAODJet vTmpRanR(1,0,0,1);
1428 for(int i = 0; i < recParticles.GetEntries(); i++){
1429 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1430 if(fTCARandomConesOut){
1431 for(int ir = 0;ir < fNRandomCones;ir++){
1432 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1433 if(jC&&jC->DeltaR(vp)<fRparam){
1434 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1435 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1436 if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1437 }
1438 }
1439 }// add up energy in cone
1440
1441 // the randomized input changes eta and phi, but keeps the p_T
1442 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1443 Double_t pTR = vp->Pt();
1444 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1445 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
8f85b773 1446
8ad1ea82 1447 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1448 Double_t pZR = pTR/TMath::Tan(thetaR);
8f85b773 1449
8ad1ea82 1450 Double_t pXR = pTR * TMath::Cos(phiR);
1451 Double_t pYR = pTR * TMath::Sin(phiR);
1452 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1453 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1454 if(fTCARandomConesOutRan){
1455 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1456 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1457 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1458 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1459 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1460 if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1461 }
1462 }
1463 }
1464 }
1465 }// loop over recparticles
8f85b773 1466
8ad1ea82 1467 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1468 if(fTCARandomConesOut){
1469 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1470 // rescale the momentum vectors for the random cones
8f85b773 1471
8ad1ea82 1472 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1473 if(rC){
1474 Double_t etaC = rC->Eta();
1475 Double_t phiC = rC->Phi();
1476 // massless jet, unit vector
1477 pTC = rC->ChargedBgEnergy();
1478 if(pTC<=0)pTC = 0.001; // for almost empty events
1479 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1480 Double_t pZC = pTC/TMath::Tan(thetaC);
1481 Double_t pXC = pTC * TMath::Cos(phiC);
1482 Double_t pYC = pTC * TMath::Sin(phiC);
1483 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1484 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1485 rC->SetBgEnergy(0,0);
1486 rC->SetEffArea(jetArea,0);
1487 }
1488 }
1489 }
1490 if(fTCARandomConesOutRan){
1491 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1492 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1493 // same wit random
1494 if(rC){
1495 Double_t etaC = rC->Eta();
1496 Double_t phiC = rC->Phi();
1497 // massless jet, unit vector
1498 pTC = rC->ChargedBgEnergy();
1499 if(pTC<=0)pTC = 0.001;// for almost empty events
1500 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1501 Double_t pZC = pTC/TMath::Tan(thetaC);
1502 Double_t pXC = pTC * TMath::Cos(phiC);
1503 Double_t pYC = pTC * TMath::Sin(phiC);
1504 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1505 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1506 rC->SetBgEnergy(0,0);
1507 rC->SetEffArea(jetArea,0);
1508 }
1509 }
1510 }
1511 }// if(fNRandomCones
d7ca0e14 1512
8ad1ea82 1513 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
8f85b773 1514
8ad1ea82 1515 if(fAODJetBackgroundOut){
1516 vector<fastjet::PseudoJet> jets2=sortedJets;
1517 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
78912d73 1518
8ad1ea82 1519 Double_t bkg1=0;
1520 Double_t sigma1=0.;
1521 Double_t meanarea1=0.;
1522 Double_t bkg2=0;
1523 Double_t sigma2=0.;
1524 Double_t meanarea2=0.;
82a35b1a 1525
8ad1ea82 1526 Double_t bkg4=0;
1527 Double_t sigma4=0.;
1528 Double_t meanarea4=0.;
02151ccc 1529
8ad1ea82 1530 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1531 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
f974a156 1532
8ad1ea82 1533 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1534 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1535 clustSeq.get_median_rho_and_sigma(sortedJets, range, true, bkg4, sigma4, meanarea4, true);
1536 fAODJetBackgroundOut->SetBackground(3,bkg4,sigma4,meanarea4);
02151ccc 1537
8ad1ea82 1538 //////////////////////
82a35b1a 1539
8ad1ea82 1540 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1541 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1542 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1543 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
f974a156 1544
8ad1ea82 1545 }
78912d73 1546 }
7bedcd31 1547
78912d73 1548 if(fStoreRhoLeadingTrackCorr) {
1549 vector<fastjet::PseudoJet> jets3=sortedJets;
1550 if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2);
7bedcd31 1551
78912d73 1552 Double_t bkg2=0;
1553 Double_t sigma2=0.;
1554 Double_t meanarea2=0.;
c0fdfe53 1555
78912d73 1556 clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true);
1557 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1558
cbfccad4 1559 //Get leading track in event
1560 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1561
1562 fh3CentvsRhoLeadingTrackPt->Fill(cent,bkg2,leading->Pt());
1563 fh3CentvsSigmaLeadingTrackPt->Fill(cent,sigma2,leading->Pt());
1564 fh3MultvsRhoLeadingTrackPt->Fill(nCh,bkg2,leading->Pt());
1565 fh3MultvsSigmaLeadingTrackPt->Fill(nCh,sigma2,leading->Pt());
78912d73 1566
1567
1568 //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1569 //exclude 2 leading jets from event
1570 //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R (Near side to leading track)
1571 //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1572 //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R) (Away side to leading track)
1573 //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R) (Orthogonal to leading track)
1574
1575 //Assuming R=0.2 for background clusters
1576
1577 Double_t bkg2Q[4] = {0.};
1578 Double_t sigma2Q[4] = {0.};
1579 Double_t meanarea2Q[4] = {0.};
1580
1581 //Get phi of leading track in event
78912d73 1582 Float_t phiLeadingTrack = leading->Phi();
1583 Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1584
1585 //Quadrant1 - near side
1586 phiMin = phiLeadingTrack - phiStep;
1587 phiMax = phiLeadingTrack + phiStep;
1588 fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1589 clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1590
1591 //Quadrant2 - orthogonal
1592 Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1593 if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1594 phiMin = phiQ2 - phiStep;
1595 phiMax = phiQ2 + phiStep;
1596 fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1597 clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1598
1599 //Quadrant3 - away side
1600 Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1601 if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1602 phiMin = phiQ3 - phiStep;
1603 phiMax = phiQ3 + phiStep;
1604 fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1605 clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1606
1607 //Quadrant4 - othogonal
1608 Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1609 if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1610 phiMin = phiQ4 - phiStep;
1611 phiMax = phiQ4 + phiStep;
1612 fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1613 clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1614
1615 fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1616 fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1617 fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1618 fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1619
1620 fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1621 fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1622 fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1623 fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1624
1625 fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1626 fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1627 fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1628 fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1629
1630 fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1631 fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1632 fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1633 fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
c0fdfe53 1634
cbfccad4 1635 fh3CentvsDeltaRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0]-bkg2,leading->Pt());
1636 fh3CentvsDeltaRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1]-bkg2,leading->Pt());
1637 fh3CentvsDeltaRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2]-bkg2,leading->Pt());
1638 fh3CentvsDeltaRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3]-bkg2,leading->Pt());
1639
54424110 1640 }
d7ca0e14 1641
82a35b1a 1642
1643
1644
dbf053f0 1645
82a35b1a 1646 // fill track information
1647 Int_t nTrackOver = recParticles.GetSize();
dbf053f0 1648 // do the same for tracks and jets
82a35b1a 1649
1650 if(nTrackOver>0){
8ad1ea82 1651 TIterator *recIter = recParticles.MakeIterator();
1652 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1653 Float_t pt = tmpRec->Pt();
1654
1655 // Printf("Leading track p_t %3.3E",pt);
1656 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1657 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1658 while(pt<ptCut&&tmpRec){
1659 nTrackOver--;
1660 tmpRec = (AliVParticle*)(recIter->Next());
1661 if(tmpRec){
1662 pt = tmpRec->Pt();
1663 }
1664 }
1665 if(nTrackOver<=0)break;
1666 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1667 }
dbf053f0 1668
8ad1ea82 1669 recIter->Reset();
1670 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1671 Float_t phi = leading->Phi();
1672 if(phi<0)phi+=TMath::Pi()*2.;
1673 Float_t eta = leading->Eta();
1674 pt = leading->Pt();
1675 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1676 Float_t tmpPt = tmpRec->Pt();
1677 Float_t tmpEta = tmpRec->Eta();
1678 fh1PtTracksRecIn->Fill(tmpPt);
1679 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1680 if(tmpRec==leading){
1681 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1682 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1683 continue;
1684 }
dbf053f0 1685 // correlation
8ad1ea82 1686 Float_t tmpPhi = tmpRec->Phi();
dbf053f0 1687
8ad1ea82 1688 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1689 Float_t dPhi = phi - tmpPhi;
1690 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1691 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1692 Float_t dEta = eta - tmpRec->Eta();
1693 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1694 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1695 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1696 }
1697 delete recIter;
1698 }
1699
1700 // find the random jets
1701
1702 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
dbf053f0 1703
8ad1ea82 1704 // fill the jet information from random track
1705 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1706 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
dbf053f0 1707
1708 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1709
8ad1ea82 1710 // loop over all jets an fill information, first on is the leading jet
dbf053f0 1711
8ad1ea82 1712 Int_t nRecOverRan = inclusiveJetsRan.size();
1713 Int_t nRecRan = inclusiveJetsRan.size();
d73e1768 1714
8ad1ea82 1715 if(inclusiveJetsRan.size()>0){
1716 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1717 Float_t pt = leadingJet.Pt();
dbf053f0 1718
8ad1ea82 1719 Int_t iCount = 0;
1720 TLorentzVector vecarearanb;
d73e1768 1721
8ad1ea82 1722 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1723 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
dbf053f0 1724 while(pt<ptCut&&iCount<nRecRan){
1725 nRecOverRan--;
1726 iCount++;
1727 if(iCount<nRecRan){
1728 pt = sortedJetsRan[iCount].perp();
1729 }
1730 }
1731 if(nRecOverRan<=0)break;
1732 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1733 }
1734 Float_t phi = leadingJet.Phi();
1735 if(phi<0)phi+=TMath::Pi()*2.;
1736 pt = leadingJet.Pt();
1737
1738 // correlation of leading jet with random tracks
1739
1740 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1741 {
1742 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1743 // correlation
1744 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1745 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1746 Float_t dPhi = phi - tmpPhi;
1747 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1748 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1749 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1750 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1751 }
1752
d6e66a82 1753 Int_t nAodOutJetsRan = 0;
8ad1ea82 1754 AliAODJet *aodOutJetRan = 0;
dbf053f0 1755 for(int j = 0; j < nRecRan;j++){
1756 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1757 Float_t tmpPt = tmpRec.Pt();
1758 fh1PtJetsRecInRan->Fill(tmpPt);
1759 // Fill Spectra with constituents
fe80c230 1760 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
dbf053f0 1761 fh1NConstRecRan->Fill(constituents.size());
1762 fh2NConstPtRan->Fill(tmpPt,constituents.size());
a0d08b6d 1763 fh2PtNchRan->Fill(nCh,tmpPt);
1764 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
d6e66a82 1765
1766
8ad1ea82 1767 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1768 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1769 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1770 aodOutJetRan->GetRefTracks()->Clear();
1771 aodOutJetRan->SetEffArea(arearan,0);
1772 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1773 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1774 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
d73e1768 1775
8ad1ea82 1776 }
d6e66a82 1777
dbf053f0 1778 // correlation
1779 Float_t tmpPhi = tmpRec.Phi();
1780 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1781
1782 if(j==0){
1783 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1784 fh1NConstLeadingRecRan->Fill(constituents.size());
1785 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1786 continue;
1787 }
1788 }
d6e66a82 1789
1790
8f85b773 1791 if(fAODJetBackgroundOut){
8ad1ea82 1792 Double_t bkg3=0.;
1793 Double_t sigma3=0.;
1794 Double_t meanarea3=0.;
1795 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1796 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1797 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1798 /*
1799 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1800 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1801 */
d6e66a82 1802 }
1803
1804
1805
8ad1ea82 1806 }
1807
1808
1809 // do the event selection if activated
1810 if(fJetTriggerPtCut>0){
1811 bool select = false;
1812 Float_t minPt = fJetTriggerPtCut;
1813 /*
1814 // hard coded for now ...
1815 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1816 if(cent<10)minPt = 50;
1817 else if(cent<30)minPt = 42;
1818 else if(cent<50)minPt = 28;
1819 else if(cent<80)minPt = 18;
1820 */
1821 float rho = 0;
1822 if(externalBackground)rho = externalBackground->GetBackground(2);
1823 if(fTCAJetsOut){
1824 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1825 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1826 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1827 if(ptSub>=minPt){
1828 select = true;
1829 break;
1830 }
1831 }
1832 }
ea950747 1833
8ad1ea82 1834 if(select){
1835 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1836 fh1CentralitySelect->Fill(cent);
1837 fh1ZSelect->Fill(zVtx);
1838 aodH->SetFillAOD(kTRUE);
1839 }
1840 }
1841 if (fDebug > 2){
1842 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1843 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1844 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1845 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1846 }
1847 PostData(1, fHistList);
dbf053f0 1848}
1849
1850void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1851{
f0659f11 1852 //
1853 // Terminate analysis
1854 //
8ad1ea82 1855 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
006b2a30 1856
8ad1ea82 1857 if(fMomResH1Fit) delete fMomResH1Fit;
1858 if(fMomResH2Fit) delete fMomResH2Fit;
1859 if(fMomResH3Fit) delete fMomResH3Fit;
006b2a30 1860
dbf053f0 1861}
1862
1863
1864Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1865
f0659f11 1866 //
1867 // get list of tracks/particles for different types
1868 //
1869
79d8877b 1870 if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
dbf053f0 1871
1872 Int_t iCount = 0;
8ad1ea82 1873 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly || type==kTrackAODMCextra || type==kTrackAODMCextraonly){
1874
1875 if(type!=kTrackAODextraonly && type!=kTrackAODMCextraonly) {
f840d64c 1876 AliAODEvent *aod = 0;
1877 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1878 else aod = AODEvent();
1879 if(!aod){
617932d0 1880 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
f840d64c 1881 return iCount;
1882 }
d7396b04 1883
f840d64c 1884 for(int it = 0;it < aod->GetNumberOfTracks();++it){
c47debfd 1885 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
1886 if(!tr) AliFatal("Not a standard AOD");
b2c1e955 1887 Bool_t bGood = false;
1888 if(fFilterType == 0)bGood = true;
8a0f6112 1889 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1890 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
b2c1e955 1891 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
617932d0 1892 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1893 continue;
1894 }
1f7b59f6
S		 1895
1896 // heavy flavor jets
1897 if(fFilterMask==528 && fUseHFcuts){
1898 Double_t ntpcClus = tr->GetTPCNcls();
1899 Double_t trPt=tr->Pt();
1900 TFormula NTPCClsCut("f1NClustersTPCLinearPtDep","70.+30./20.*x");
1901
1902 if (trPt <= 20. && (ntpcClus < NTPCClsCut.Eval(trPt))) continue;
1903 else if (trPt > 20. && ntpcClus < 100) continue;
1904
3bd35123 1905 if (AvoidDoubleCountingHF(aod,tr)) continue;
1f7b59f6 1906 }
3bd35123 1907 //
1f7b59f6 1908
f662dc87 1909 if(fRequireITSRefit){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
27c67edf
ML		 1910 if (fApplySharedClusterCut) {
1911 Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
1912 if (frac > 0.4) continue;
1913 }
617932d0 1914 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1915 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1916 continue;
1917 }
1918 if(tr->Pt()<fTrackPtCut){
1919 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1920 continue;
1921 }
1922 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
f840d64c 1923 list->Add(tr);
1924 iCount++;
1925 }
dbf053f0 1926 }
f840d64c 1927 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1928 AliAODEvent *aod = 0;
1929 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1930 else aod = AODEvent();
1931
1932 if(!aod){
1933 return iCount;
1934 }
1935 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
8ddcf605 1936 if(!aodExtraTracks)return iCount;
f840d64c 1937 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1938 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1939 if (!track) continue;
1940
1941 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
8ddcf605 1942 if(!trackAOD)continue;
b2c1e955 1943 Bool_t bGood = false;
1944 if(fFilterType == 0)bGood = true;
8a0f6112 1945 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1946 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
b2c1e955 1947 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
f662dc87 1948 if(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
9c7b62e0 1949 if (fApplySharedClusterCut) {
1950 Double_t frac = Double_t(trackAOD->GetTPCnclsS()) /Double_t(trackAOD->GetTPCncls());
1951 if (frac > 0.4) continue;
1952 }
1953
1954
8ad1ea82 1955 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1956 if(trackAOD->Pt()<fTrackPtCut) continue;
79d8877b 1957 if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
8ad1ea82 1958 list->Add(trackAOD);
1959 iCount++;
1960 }
1961 }
1962
1963 if(type==kTrackAODMCextra || type==kTrackAODMCextraonly) { //embed generator level particles
1964 AliAODEvent *aod = 0;
1965 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1966 else aod = AODEvent();
1967 if(!aod){
1968 return iCount;
1969 }
1970 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraMCparticles"));
1971 if(!aodExtraTracks)return iCount;
1972 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1973 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
a998fc5a 1974 AliAODMCParticle *partmc = dynamic_cast<AliAODMCParticle*> ((*aodExtraTracks)[it]);
79d8877b 1975 if (!track) {
66ca3666 1976 if(fDebug>10) printf("track %d does not exist\n",it);
79d8877b 1977 continue;
1978 }
8ad1ea82 1979
7e4edf34 1980 if(!partmc) continue;
1981 if(!partmc->IsPhysicalPrimary())continue;
a998fc5a 1982
1983 if (track->Pt()<fTrackPtCut) {
7e4edf34 1984 if(fDebug>10) printf("track %d has too low pt %.2f\n",it,track->Pt());
a998fc5a 1985 continue;
1986 }
1987
7e4edf34 1988 /*
1989 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*>((*aodExtraTracks)[it]);//(track);
a998fc5a 1990
79d8877b 1991 if(!trackAOD) {
66ca3666 1992 if(fDebug>10) printf("trackAOD %d does not exist\n",it);
79d8877b 1993 continue;
1994 }
1995
1996 trackAOD->SetFlags(AliESDtrack::kEmbedded);
79d8877b 1997 trackAOD->SetFilterMap(fFilterMask);
7e4edf34 1998 */
1999 if(fDebug>10) printf("pt extra track %.2f \n", track->Pt());
79d8877b 2000
7e4edf34 2001 if(TMath::Abs(track->Eta())>fTrackEtaWindow) continue;
2002 if(track->Pt()<fTrackPtCut) continue;
2003 list->Add(track);
2004
f840d64c 2005 iCount++;
2006 }
dbf053f0 2007 }
8ad1ea82 2008
dbf053f0 2009 }
2010 else if (type == kTrackKineAll||type == kTrackKineCharged){
2011 AliMCEvent* mcEvent = MCEvent();
2012 if(!mcEvent)return iCount;
2013 // we want to have alivpartilces so use get track
2014 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
2015 if(!mcEvent->IsPhysicalPrimary(it))continue;
2016 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
2017 if(type == kTrackKineAll){
0341d0d8 2018 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 2019 list->Add(part);
2020 iCount++;
2021 }
2022 else if(type == kTrackKineCharged){
2023 if(part->Particle()->GetPDG()->Charge()==0)continue;
0341d0d8 2024 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 2025 list->Add(part);
2026 iCount++;
2027 }
2028 }
2029 }
2030 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
2031 AliAODEvent *aod = 0;
959508f4 2032 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
dbf053f0 2033 else aod = AODEvent();
2034 if(!aod)return iCount;
2035 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2036 if(!tca)return iCount;
2037 for(int it = 0;it < tca->GetEntriesFast();++it){
225f0094 2038 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
dbf053f0 2039 if(!part->IsPhysicalPrimary())continue;
2040 if(type == kTrackAODMCAll){
0341d0d8 2041 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 2042 list->Add(part);
2043 iCount++;
2044 }
2045 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
2046 if(part->Charge()==0)continue;
0341d0d8 2047 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 2048 if(kTrackAODMCCharged){
2049 list->Add(part);
2050 }
2051 else {
8a320ab9 2052 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
dbf053f0 2053 list->Add(part);
2054 }
2055 iCount++;
2056 }
2057 }
2058 }// AODMCparticle
d8589aad 2059 else if (type == kTrackAODMCHF){
2060
2061 AliAODEvent *aod = 0;
2062 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2063 else aod = AODEvent();
2064 if(!aod)return iCount;
2065 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2066 if(!tca)return iCount;
2067 for(int it = 0;it < tca->GetEntriesFast();++it){
2068 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
2069 if(!part) continue;
2070 int partpdg= part->PdgCode();
2071 if(!part->IsPhysicalPrimary() && !IsBMeson(partpdg) && !IsDMeson(partpdg) )continue;
2072
2073 if (IsBMeson(partpdg) || IsDMeson(partpdg)) {
2074 iCount+= AddDaughters( list , part,tca);
2075 }
2076 else {
2077
2078 if(part->Pt()<fTrackPtCut) continue;
2079 if(TMath::Abs(part->Eta())>fTrackEtaWindow) continue;
2080 if(part->Charge()==0) continue;
2081
2082 if((part->Pt()>=fTrackPtCut) && (TMath::Abs(part->Eta())<=fTrackEtaWindow) && (part->Charge()!=0))list->Add(part);
2083 iCount++;
2084 }
2085 }
2086 }
2087
dbf053f0 2088 list->Sort();
2089 return iCount;
dbf053f0 2090}
2091
d8589aad 2092Int_t AliAnalysisTaskJetCluster::AddDaughters(TList * list, AliAODMCParticle *part, TClonesArray * tca){
2093 Int_t count=0;
2094 Int_t nDaugthers = part->GetNDaughters();
2095 for(Int_t i=0;i< nDaugthers;++i){
2096 AliAODMCParticle *partdaughter = (AliAODMCParticle*)(tca->At(i));
2097 if(!partdaughter) continue;
2098 if(partdaughter->Pt()<fTrackPtCut)continue;
2099 if(TMath::Abs(partdaughter->Eta())>fTrackEtaWindow)continue;
2100 if(partdaughter->Charge()==0)continue;
2101
2102 if(!IsDMeson(partdaughter->PdgCode()) && !IsBMeson(partdaughter->PdgCode()) ){
2103 list->Add(partdaughter);
2104 count++;
2105 }
2106 else count+=AddDaughters(list,part,tca);
2107 }
2108return count;
2109}
3bd35123 2110
2111
2112Bool_t AliAnalysisTaskJetCluster::AvoidDoubleCountingHF(AliAODEvent *aod, AliAODTrack *tr1){
1f7b59f6 2113
3bd35123 2114 if(!(tr1->TestFilterBit(BIT(9)))) return kFALSE;
d8589aad 2115
3bd35123 2116 Int_t idtr1 = tr1->GetID();
d8589aad 2117
3bd35123 2118 for(int jt = 0;jt < aod->GetNumberOfTracks();++jt){
1f7b59f6 2119
c47debfd 2120 AliAODTrack *tr2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(jt));
2121 if(!tr2) AliFatal("Not a standard AOD");
3bd35123 2122 Int_t idtr2 = tr2->GetID();
2123
2124 if (!(tr2->TestFilterBit(BIT(4)))) continue;
2125 if (idtr1==(idtr2+1)*-1.) return kTRUE;
1f7b59f6 2126
3bd35123 2127 }
2128 return kFALSE;
1f7b59f6 2129}
fa7c34ba 2130void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
2131
8ad1ea82 2132 if (!gGrid) {
2133 AliInfo("Trying to connect to AliEn ...");
2134 TGrid::Connect("alien://");
2135 }
2136
575dee18 2137 TFile *f = TFile::Open(fPathTrPtResolution.Data());
fa7c34ba 2138
0e5bbbbc 2139 if(!f)return;
2140
fa7c34ba 2141 TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
2142 TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
2143 TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
2144
2145 SetSmearResolution(kTRUE);
2146 SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
575dee18 2147
fa7c34ba 2148
2149}
2150
2151void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
2152
8ad1ea82 2153 if (!gGrid) {
2154 AliInfo("Trying to connect to AliEn ...");
2155 TGrid::Connect("alien://");
2156 }
2157
575dee18 2158 TFile *f = TFile::Open(fPathTrEfficiency.Data());
0e5bbbbc 2159 if(!f)return;
fa7c34ba 2160
2161 TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
2162 TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
2163 TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
2164
2165 SetDiceEfficiency(kTRUE);
2166
2167 if(fChangeEfficiencyFraction>0.) {
2168
2169 TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
2170
2171 for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
2172 Double_t content = hEffPosGlobSt->GetBinContent(bin);
2173 hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
2174 }
2175
2176 SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2177
2178 }
2179 else
2180 SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2181
fa7c34ba 2182}
2183
006b2a30 2184void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
2185
2186 //
2187 // set mom res profiles
2188 //
2189
575dee18 2190 if(fMomResH1) delete fMomResH1;
2191 if(fMomResH2) delete fMomResH2;
2192 if(fMomResH3) delete fMomResH3;
2193
2194 fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
2195 fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
2196 fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
2197
006b2a30 2198}
2199
2200void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
2201 //
2202 // set tracking efficiency histos
2203 //
2204
2205 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
2206 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
2207 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
2208}
2209
2210Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
2211
2212 //
2213 // Get smearing on generated momentum
2214 //
2215
2216 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
2217
2218 TProfile *fMomRes = 0x0;
2219 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
2220 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
2221 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
2222
2223 if(!fMomRes) {
2224 return 0.;
2225 }
2226
2227
2228 Double_t smear = 0.;
2229
2230 if(pt>20.) {
2231 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
2232 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
2233 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
2234 }
2235 else {
2236
2237 Int_t bin = fMomRes->FindBin(pt);
2238
2239 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
2240
2241 }
2242
2243 if(fMomRes) delete fMomRes;
2244
2245 return smear;
2246}
2247
2248void AliAnalysisTaskJetCluster::FitMomentumResolution() {
2249 //
2250 // Fit linear function on momentum resolution at high pT
2251 //
2252
2253 if(!fMomResH1Fit && fMomResH1) {
2254 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
2255 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
2256 fMomResH1Fit ->SetRange(5.,100.);
2257 }
2258
2259 if(!fMomResH2Fit && fMomResH2) {
2260 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2261 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2262 fMomResH2Fit ->SetRange(5.,100.);
2263 }
2264
2265 if(!fMomResH3Fit && fMomResH3) {
2266 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2267 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2268 fMomResH3Fit ->SetRange(5.,100.);
2269 }
2270
2271}
2272
dbf053f0 2273/*
8ad1ea82 2274 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
dbf053f0 2275 for(int i = 0; i < particles.GetEntries(); i++){
8ad1ea82 2276 AliVParticle *vp = (AliVParticle*)particles.At(i);
2277 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2278 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2279 jInp.set_user_index(i);
2280 inputParticles.push_back(jInp);
dbf053f0 2281 }
2282
2283 return 0;
2284
8ad1ea82 2285 }
dbf053f0 2286*/
d8589aad 2287
2288
2289bool AliAnalysisTaskJetCluster::IsBMeson(int pc){
2290 int bPdG[] = {511,521,10511,10521,513,523,10513,10523,20513,20523,20513,20523,515,525,531,
2291 10531,533,10533,20533,535,541,10541,543,10543,20543,545,551,10551,100551,
2292 110551,200551,210551,553,10553,20553,30553,100553,110553,120553,130553,200553,210553,220553,
2293 300553,9000533,9010553,555,10555,20555,100555,110555,120555,200555,557,100557};
2294 for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true;
2295return false;
2296}
2297bool AliAnalysisTaskJetCluster::IsDMeson(int pc){
2298 int bPdG[] = {411,421,10411,10421,413,423,10413,10423,20431,20423,415,
2299 425,431,10431,433,10433,20433,435,441,10441,100441,443,10443,20443,
2300 100443,30443,9000443,9010443,9020443,445,100445};
2301 for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true;
2302return false;
2303}
d43a2b0c 2304