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