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