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