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