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