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