TPC event task updated
[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)
d6e66a82 1431
4096256e 1432
8f85b773 1433
1434
1435
1436 if(fAODJetBackgroundOut){
d6e66a82 1437 vector<fastjet::PseudoJet> jets2=sortedJets;
1438 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1439 Double_t bkg1=0;
1440 Double_t sigma1=0.;
1441 Double_t meanarea1=0.;
1442 Double_t bkg2=0;
1443 Double_t sigma2=0.;
1444 Double_t meanarea2=0.;
82a35b1a 1445
d73e1768 1446 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
8f85b773 1447 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
f974a156 1448
1449 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1450 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
82a35b1a 1451
4a33f620 1452 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
8f85b773 1453 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
f974a156 1454 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1455 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1456
7bedcd31 1457 if(fStoreRhoLeadingTrackCorr) {
1458 fh2CentvsRho->Fill(cent,bkg2);
1459 fh2CentvsSigma->Fill(cent,sigma2);
1460 fh2MultvsRho->Fill(nCh,bkg2);
1461 fh2MultvsSigma->Fill(nCh,sigma2);
1462
1463
1464 //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1465 //exclude 2 leading jets from event
1466 //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R (Near side to leading track)
1467 //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1468 //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R) (Away side to leading track)
1469 //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R) (Orthogonal to leading track)
1470
1471 //Assuming R=0.2 for background clusters
1472
1473 Double_t bkg2Q[4] = {0.};
1474 Double_t sigma2Q[4] = {0.};
1475 Double_t meanarea2Q[4] = {0.};
1476
1477 //Get phi of leading track in event
1478 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1479 Float_t phiLeadingTrack = leading->Phi();
1480 Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1481
1482 //Quadrant1 - near side
1483 phiMin = phiLeadingTrack - phiStep;
1484 phiMax = phiLeadingTrack + phiStep;
1485 fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1486 clustSeq.get_median_rho_and_sigma(jets2, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1487
1488 //Quadrant2 - orthogonal
1489 Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1490 if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1491 phiMin = phiQ2 - phiStep;
1492 phiMax = phiQ2 + phiStep;
1493 fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1494 clustSeq.get_median_rho_and_sigma(jets2, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1495
1496 //Quadrant3 - away side
1497 Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1498 if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1499 phiMin = phiQ3 - phiStep;
1500 phiMax = phiQ3 + phiStep;
1501 fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1502 clustSeq.get_median_rho_and_sigma(jets2, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1503
1504 //Quadrant4 - othogonal
1505 Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1506 if(phiQ3>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1507 phiMin = phiQ4 - phiStep;
1508 phiMax = phiQ4 + phiStep;
1509 fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1510 clustSeq.get_median_rho_and_sigma(jets2, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1511
1512 fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1513 fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1514 fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1515 fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1516
1517 fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1518 fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1519 fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1520 fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1521
1522 fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1523 fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1524 fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1525 fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1526
1527 fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1528 fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1529 fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1530 fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
1531
1532 }
c0fdfe53 1533
c0fdfe53 1534
d6e66a82 1535 }
54424110 1536 }
d7ca0e14 1537
82a35b1a 1538
1539
1540
dbf053f0 1541
82a35b1a 1542 // fill track information
1543 Int_t nTrackOver = recParticles.GetSize();
dbf053f0 1544 // do the same for tracks and jets
82a35b1a 1545
1546 if(nTrackOver>0){
dbf053f0 1547 TIterator *recIter = recParticles.MakeIterator();
1548 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1549 Float_t pt = tmpRec->Pt();
82a35b1a 1550
dbf053f0 1551 // Printf("Leading track p_t %3.3E",pt);
1552 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1553 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1554 while(pt<ptCut&&tmpRec){
1555 nTrackOver--;
1556 tmpRec = (AliVParticle*)(recIter->Next());
1557 if(tmpRec){
1558 pt = tmpRec->Pt();
1559 }
1560 }
1561 if(nTrackOver<=0)break;
1562 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1563 }
1564
1565 recIter->Reset();
1566 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1567 Float_t phi = leading->Phi();
1568 if(phi<0)phi+=TMath::Pi()*2.;
1569 Float_t eta = leading->Eta();
1570 pt = leading->Pt();
1571 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1572 Float_t tmpPt = tmpRec->Pt();
1573 Float_t tmpEta = tmpRec->Eta();
1574 fh1PtTracksRecIn->Fill(tmpPt);
1575 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1576 if(tmpRec==leading){
1577 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1578 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1579 continue;
1580 }
1581 // correlation
1582 Float_t tmpPhi = tmpRec->Phi();
1583
1584 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1585 Float_t dPhi = phi - tmpPhi;
1586 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1587 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1588 Float_t dEta = eta - tmpRec->Eta();
1589 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1590 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1591 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1592 }
1593 delete recIter;
1594 }
1595
1596 // find the random jets
fe80c230 1597
d6e66a82 1598 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
dbf053f0 1599
dbf053f0 1600 // fill the jet information from random track
fe80c230 1601 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1602 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
dbf053f0 1603
1604 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1605
1606 // loop over all jets an fill information, first on is the leading jet
1607
1608 Int_t nRecOverRan = inclusiveJetsRan.size();
1609 Int_t nRecRan = inclusiveJetsRan.size();
d73e1768 1610
dbf053f0 1611 if(inclusiveJetsRan.size()>0){
1612 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1613 Float_t pt = leadingJet.Pt();
1614
1615 Int_t iCount = 0;
d73e1768 1616 TLorentzVector vecarearanb;
1617
dbf053f0 1618 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1619 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1620 while(pt<ptCut&&iCount<nRecRan){
1621 nRecOverRan--;
1622 iCount++;
1623 if(iCount<nRecRan){
1624 pt = sortedJetsRan[iCount].perp();
1625 }
1626 }
1627 if(nRecOverRan<=0)break;
1628 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1629 }
1630 Float_t phi = leadingJet.Phi();
1631 if(phi<0)phi+=TMath::Pi()*2.;
1632 pt = leadingJet.Pt();
1633
1634 // correlation of leading jet with random tracks
1635
1636 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1637 {
1638 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1639 // correlation
1640 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1641 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1642 Float_t dPhi = phi - tmpPhi;
1643 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1644 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1645 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1646 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1647 }
1648
d6e66a82 1649 Int_t nAodOutJetsRan = 0;
1650 AliAODJet *aodOutJetRan = 0;
dbf053f0 1651 for(int j = 0; j < nRecRan;j++){
1652 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1653 Float_t tmpPt = tmpRec.Pt();
1654 fh1PtJetsRecInRan->Fill(tmpPt);
1655 // Fill Spectra with constituents
fe80c230 1656 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
dbf053f0 1657 fh1NConstRecRan->Fill(constituents.size());
1658 fh2NConstPtRan->Fill(tmpPt,constituents.size());
a0d08b6d 1659 fh2PtNchRan->Fill(nCh,tmpPt);
1660 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
d6e66a82 1661
1662
8f85b773 1663 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1664 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
d6e66a82 1665 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
d508795b 1666 aodOutJetRan->GetRefTracks()->Clear();
d73e1768 1667 aodOutJetRan->SetEffArea(arearan,0);
1668 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1669 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1670 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1671
1672 }
d6e66a82 1673
dbf053f0 1674 // correlation
1675 Float_t tmpPhi = tmpRec.Phi();
1676 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1677
1678 if(j==0){
1679 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1680 fh1NConstLeadingRecRan->Fill(constituents.size());
1681 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1682 continue;
1683 }
1684 }
d6e66a82 1685
1686
8f85b773 1687 if(fAODJetBackgroundOut){
d6e66a82 1688 Double_t bkg3=0.;
1689 Double_t sigma3=0.;
1690 Double_t meanarea3=0.;
4a33f620 1691 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
8f85b773 1692 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
d7ca0e14 1693 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1694 /*
82a35b1a 1695 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
9c6bd749 1696 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
d7ca0e14 1697 */
d6e66a82 1698 }
1699
1700
1701
dbf053f0 1702 }
0341d0d8 1703
1704
1705 // do the event selection if activated
1706 if(fJetTriggerPtCut>0){
1707 bool select = false;
2b018f7a 1708 Float_t minPt = fJetTriggerPtCut;
1709 /*
0341d0d8 1710 // hard coded for now ...
2b018f7a 1711 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
0341d0d8 1712 if(cent<10)minPt = 50;
aae7993b 1713 else if(cent<30)minPt = 42;
0341d0d8 1714 else if(cent<50)minPt = 28;
1715 else if(cent<80)minPt = 18;
2b018f7a 1716 */
0341d0d8 1717 float rho = 0;
1718 if(externalBackground)rho = externalBackground->GetBackground(2);
8f85b773 1719 if(fTCAJetsOut){
1720 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1721 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
0341d0d8 1722 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1723 if(ptSub>=minPt){
1724 select = true;
1725 break;
1726 }
1727 }
1728 }
ea950747 1729
0341d0d8 1730 if(select){
1731 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1732 fh1CentralitySelect->Fill(cent);
e2ca7519 1733 fh1ZSelect->Fill(zVtx);
0341d0d8 1734 aodH->SetFillAOD(kTRUE);
1735 }
1736 }
ea950747 1737 if (fDebug > 2){
1738 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1739 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1740 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1741 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1742 }
dbf053f0 1743 PostData(1, fHistList);
1744}
1745
1746void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1747{
f0659f11 1748 //
1749 // Terminate analysis
1750 //
006b2a30 1751 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1752
1753 if(fMomResH1Fit) delete fMomResH1Fit;
1754 if(fMomResH2Fit) delete fMomResH2Fit;
1755 if(fMomResH3Fit) delete fMomResH3Fit;
1756
dbf053f0 1757}
1758
1759
1760Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1761
f0659f11 1762 //
1763 // get list of tracks/particles for different types
1764 //
1765
dbf053f0 1766 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1767
1768 Int_t iCount = 0;
f840d64c 1769 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1770 if(type!=kTrackAODextraonly) {
1771 AliAODEvent *aod = 0;
1772 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1773 else aod = AODEvent();
1774 if(!aod){
617932d0 1775 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
f840d64c 1776 return iCount;
1777 }
d7396b04 1778
f840d64c 1779 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1780 AliAODTrack *tr = aod->GetTrack(it);
b2c1e955 1781 Bool_t bGood = false;
1782 if(fFilterType == 0)bGood = true;
8a0f6112 1783 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1784 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
b2c1e955 1785 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
617932d0 1786 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1787 continue;
1788 }
1789 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1790 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1791 continue;
1792 }
1793 if(tr->Pt()<fTrackPtCut){
1794 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1795 continue;
1796 }
1797 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
f840d64c 1798 list->Add(tr);
1799 iCount++;
1800 }
dbf053f0 1801 }
f840d64c 1802 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1803 AliAODEvent *aod = 0;
1804 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1805 else aod = AODEvent();
1806
1807 if(!aod){
1808 return iCount;
1809 }
1810 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
8ddcf605 1811 if(!aodExtraTracks)return iCount;
f840d64c 1812 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1813 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1814 if (!track) continue;
1815
1816 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
8ddcf605 1817 if(!trackAOD)continue;
b2c1e955 1818 Bool_t bGood = false;
1819 if(fFilterType == 0)bGood = true;
8a0f6112 1820 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1821 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
b2c1e955 1822 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
ca6d12f9 1823 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
f840d64c 1824 if(trackAOD->Pt()<fTrackPtCut) continue;
f840d64c 1825 list->Add(trackAOD);
1826 iCount++;
1827 }
dbf053f0 1828 }
1829 }
1830 else if (type == kTrackKineAll||type == kTrackKineCharged){
1831 AliMCEvent* mcEvent = MCEvent();
1832 if(!mcEvent)return iCount;
1833 // we want to have alivpartilces so use get track
1834 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1835 if(!mcEvent->IsPhysicalPrimary(it))continue;
1836 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1837 if(type == kTrackKineAll){
0341d0d8 1838 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1839 list->Add(part);
1840 iCount++;
1841 }
1842 else if(type == kTrackKineCharged){
1843 if(part->Particle()->GetPDG()->Charge()==0)continue;
0341d0d8 1844 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1845 list->Add(part);
1846 iCount++;
1847 }
1848 }
1849 }
1850 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1851 AliAODEvent *aod = 0;
959508f4 1852 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
dbf053f0 1853 else aod = AODEvent();
1854 if(!aod)return iCount;
1855 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1856 if(!tca)return iCount;
1857 for(int it = 0;it < tca->GetEntriesFast();++it){
225f0094 1858 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
dbf053f0 1859 if(!part->IsPhysicalPrimary())continue;
1860 if(type == kTrackAODMCAll){
0341d0d8 1861 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1862 list->Add(part);
1863 iCount++;
1864 }
1865 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1866 if(part->Charge()==0)continue;
0341d0d8 1867 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1868 if(kTrackAODMCCharged){
1869 list->Add(part);
1870 }
1871 else {
8a320ab9 1872 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
dbf053f0 1873 list->Add(part);
1874 }
1875 iCount++;
1876 }
1877 }
1878 }// AODMCparticle
1879 list->Sort();
1880 return iCount;
dbf053f0 1881}
1882
fa7c34ba 1883void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
1884
575dee18 1885 TFile *f = TFile::Open(fPathTrPtResolution.Data());
fa7c34ba 1886
0e5bbbbc 1887 if(!f)return;
1888
fa7c34ba 1889 TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
1890 TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
1891 TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
1892
1893 SetSmearResolution(kTRUE);
1894 SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
575dee18 1895
fa7c34ba 1896
1897}
1898
1899void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
1900
575dee18 1901 TFile *f = TFile::Open(fPathTrEfficiency.Data());
0e5bbbbc 1902 if(!f)return;
fa7c34ba 1903
1904 TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
1905 TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
1906 TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
1907
1908 SetDiceEfficiency(kTRUE);
1909
1910 if(fChangeEfficiencyFraction>0.) {
1911
1912 TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
1913
1914 for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
1915 Double_t content = hEffPosGlobSt->GetBinContent(bin);
1916 hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
1917 }
1918
1919 SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1920
1921 }
1922 else
1923 SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1924
fa7c34ba 1925}
1926
006b2a30 1927void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1928
1929 //
1930 // set mom res profiles
1931 //
1932
575dee18 1933 if(fMomResH1) delete fMomResH1;
1934 if(fMomResH2) delete fMomResH2;
1935 if(fMomResH3) delete fMomResH3;
1936
1937 fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
1938 fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
1939 fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
1940
006b2a30 1941}
1942
1943void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1944 //
1945 // set tracking efficiency histos
1946 //
1947
1948 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1949 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1950 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1951}
1952
1953Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1954
1955 //
1956 // Get smearing on generated momentum
1957 //
1958
1959 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1960
1961 TProfile *fMomRes = 0x0;
1962 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1963 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1964 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1965
1966 if(!fMomRes) {
1967 return 0.;
1968 }
1969
1970
1971 Double_t smear = 0.;
1972
1973 if(pt>20.) {
1974 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1975 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1976 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1977 }
1978 else {
1979
1980 Int_t bin = fMomRes->FindBin(pt);
1981
1982 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1983
1984 }
1985
1986 if(fMomRes) delete fMomRes;
1987
1988 return smear;
1989}
1990
1991void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1992 //
1993 // Fit linear function on momentum resolution at high pT
1994 //
1995
1996 if(!fMomResH1Fit && fMomResH1) {
1997 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1998 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1999 fMomResH1Fit ->SetRange(5.,100.);
2000 }
2001
2002 if(!fMomResH2Fit && fMomResH2) {
2003 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2004 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2005 fMomResH2Fit ->SetRange(5.,100.);
2006 }
2007
2008 if(!fMomResH3Fit && fMomResH3) {
2009 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2010 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2011 fMomResH3Fit ->SetRange(5.,100.);
2012 }
2013
2014}
2015
dbf053f0 2016/*
2017Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
2018 for(int i = 0; i < particles.GetEntries(); i++){
2019 AliVParticle *vp = (AliVParticle*)particles.At(i);
2020 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2021 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2022 jInp.set_user_index(i);
2023 inputParticles.push_back(jInp);
2024 }
2025
2026 return 0;
2027
2028}
2029*/