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