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