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