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