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