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