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