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