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