]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/DEV/AliAnalysisTaskJetCluster.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Compatibility with ROOT trunk
[u/mrichter/AliRoot.git] / JETAN / DEV / AliAnalysisTaskJetCluster.cxx
CommitLineData
d89b8229 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>
24#include <TRandom3.h>
25#include <TSystem.h>
26#include <TInterpreter.h>
27#include <TChain.h>
28#include <TRefArray.h>
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>
d89b8229 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"
46#include "AliESDEvent.h"
47#include "AliAODEvent.h"
48#include "AliAODHandler.h"
49#include "AliAODExtension.h"
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"
60#include "AliAODJetEventBackground.h"
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;
95808c2b 70
d89b8229 71ClassImp(AliAnalysisTaskJetCluster)
72
73AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
f0659f11 74 //
75 // Destructor
76 //
006b2a30 77
d89b8229 78 delete fRef;
79 delete fRandom;
80
81 if(fTCAJetsOut)fTCAJetsOut->Delete();
82 delete fTCAJetsOut;
006b2a30 83
d89b8229 84 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
85 delete fTCAJetsOutRan;
006b2a30 86
d89b8229 87 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
88 delete fTCARandomConesOut;
006b2a30 89
d89b8229 90 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
91 delete fTCARandomConesOutRan;
006b2a30 92
d89b8229 93 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
94 delete fAODJetBackgroundOut;
95}
96
97AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
98 AliAnalysisTaskSE(),
99 fAOD(0x0),
100 fAODExtension(0x0),
101 fRef(new TRefArray),
102 fUseAODTrackInput(kFALSE),
103 fUseAODMCInput(kFALSE),
104 fUseBackgroundCalc(kFALSE),
105 fEventSelection(kFALSE),
106 fFilterMask(0),
d7396b04 107 fFilterMaskBestPt(0),
d89b8229 108 fFilterType(0),
109 fJetTypes(kJet),
110 fTrackTypeRec(kTrackUndef),
111 fTrackTypeGen(kTrackUndef),
112 fNSkipLeadingRan(0),
113 fNSkipLeadingCone(0),
114 fNRandomCones(0),
115 fAvgTrials(1),
116 fExternalWeight(1),
117 fTrackEtaWindow(0.9),
118 fRecEtaWindow(0.5),
119 fTrackPtCut(0.),
120 fJetOutputMinPt(0.150),
121 fMaxTrackPtInJet(100.),
122 fJetTriggerPtCut(0),
123 fVtxZCut(8),
124 fVtxR2Cut(1),
125 fCentCutUp(0),
126 fCentCutLo(0),
127 fNonStdBranch(""),
128 fBackgroundBranch(""),
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),
d89b8229 141 fRparam(1.0),
142 fAlgorithm(fastjet::kt_algorithm),
143 fStrategy(fastjet::Best),
144 fRecombScheme(fastjet::BIpt_scheme),
145 fAreaType(fastjet::active_area),
146 fGhostArea(0.01),
147 fActiveAreaRepeats(1),
148 fGhostEtamax(1.5),
149 fTCAJetsOut(0x0),
150 fTCAJetsOutRan(0x0),
151 fTCARandomConesOut(0x0),
152 fTCARandomConesOutRan(0x0),
153 fAODJetBackgroundOut(0x0),
154 fRandom(0),
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),
175 fh1Nch(0x0),
176 fh1CentralityPhySel(0x0),
177 fh1Centrality(0x0),
178 fh1CentralitySelect(0x0),
179 fh1ZPhySel(0x0),
180 fh1Z(0x0),
181 fh1ZSelect(0x0),
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),
203 fh2PtNch(0x0),
204 fh2PtNchRan(0x0),
205 fh2PtNchN(0x0),
206 fh2PtNchNRan(0x0),
207 fh2TracksLeadingJetPhiPtRan(0x0),
208 fh2TracksLeadingJetPhiPtWRan(0x0),
006b2a30 209 fh2PtGenPtSmeared(0x0),
210 fp1Efficiency(0x0),
211 fp1PtResolution(0x0),
d89b8229 212 fHistList(0x0)
213{
f0659f11 214 //
215 // Constructor
216 //
217
d89b8229 218 for(int i = 0;i<3;i++){
219 fh1BiARandomCones[i] = 0;
220 fh1BiARandomConesRan[i] = 0;
221 }
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 }
228}
229
230AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
231 AliAnalysisTaskSE(name),
232 fAOD(0x0),
233 fAODExtension(0x0),
234 fRef(new TRefArray),
235 fUseAODTrackInput(kFALSE),
236 fUseAODMCInput(kFALSE),
237 fUseBackgroundCalc(kFALSE),
d7396b04 238 fEventSelection(kFALSE), fFilterMask(0),
239 fFilterMaskBestPt(0),
d89b8229 240 fFilterType(0),
241 fJetTypes(kJet),
242 fTrackTypeRec(kTrackUndef),
243 fTrackTypeGen(kTrackUndef),
244 fNSkipLeadingRan(0),
245 fNSkipLeadingCone(0),
246 fNRandomCones(0),
247 fAvgTrials(1),
248 fExternalWeight(1),
249 fTrackEtaWindow(0.9),
250 fRecEtaWindow(0.5),
251 fTrackPtCut(0.),
252 fJetOutputMinPt(0.150),
253 fMaxTrackPtInJet(100.),
254 fJetTriggerPtCut(0),
255 fVtxZCut(8),
256 fVtxR2Cut(1),
257 fCentCutUp(0),
258 fCentCutLo(0),
259 fNonStdBranch(""),
260 fBackgroundBranch(""),
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),
d89b8229 273 fRparam(1.0),
274 fAlgorithm(fastjet::kt_algorithm),
275 fStrategy(fastjet::Best),
276 fRecombScheme(fastjet::BIpt_scheme),
277 fAreaType(fastjet::active_area),
278 fGhostArea(0.01),
279 fActiveAreaRepeats(1),
280 fGhostEtamax(1.5),
281 fTCAJetsOut(0x0),
282 fTCAJetsOutRan(0x0),
283 fTCARandomConesOut(0x0),
284 fTCARandomConesOutRan(0x0),
285 fAODJetBackgroundOut(0x0),
286 fRandom(0),
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),
307 fh1Nch(0x0),
308 fh1CentralityPhySel(0x0),
309 fh1Centrality(0x0),
310 fh1CentralitySelect(0x0),
311 fh1ZPhySel(0x0),
312 fh1Z(0x0),
313 fh1ZSelect(0x0),
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),
335 fh2PtNch(0x0),
336 fh2PtNchRan(0x0),
337 fh2PtNchN(0x0),
338 fh2PtNchNRan(0x0),
339 fh2TracksLeadingJetPhiPtRan(0x0),
340 fh2TracksLeadingJetPhiPtWRan(0x0),
006b2a30 341 fh2PtGenPtSmeared(0x0),
342 fp1Efficiency(0x0),
343 fp1PtResolution(0x0),
d89b8229 344 fHistList(0x0)
345{
f0659f11 346 //
347 // named ctor
348 //
006b2a30 349
d89b8229 350 for(int i = 0;i<3;i++){
351 fh1BiARandomCones[i] = 0;
352 fh1BiARandomConesRan[i] = 0;
353 }
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 }
360 DefineOutput(1, TList::Class());
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
381 fRandom = new TRandom3(0);
382
383
384 // Connect the AOD
385
386
387 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
388
389
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
397
398 if(fJetTypes&kJet){
399 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
400 fTCAJetsOut->SetName(fNonStdBranch.Data());
401 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
402 }
006b2a30 403
d89b8229 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));
d89b8229 409 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
410 }
411
412 if(fUseBackgroundCalc){
413 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
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
d89b8229 419 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
420 }
421 }
422 // create the branch for the random cones with the same R
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
d89b8229 427 if(fNRandomCones>0){
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 }
434 }
435 // create the branch with the random for the random cones on the random event
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 }
443 }
444 }
445
446 if(fNonStdFile.Length()!=0){
447 //
448 // case that we have an AOD extension we need to fetch the jets from the extended output
449 // we identify the extension aod event by looking for the branchname
450 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
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);
453 }
454 }
455
006b2a30 456 // FitMomentumResolution();
457
d89b8229 458
459 if(!fHistList)fHistList = new TList();
460 fHistList->SetOwner();
461 PostData(1, fHistList); // post data in any case once
462
463 Bool_t oldStatus = TH1::AddDirectoryStatus();
464 TH1::AddDirectory(kFALSE);
465
466 //
467 // Histogram
468
469 const Int_t nBinPt = 100;
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
476 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
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
504 const int nChMax = 5000;
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);
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);
535 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
536
537 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
538 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
539 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
540
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
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
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
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);
d89b8229 612
613 if(fNRandomCones>0&&fUseBackgroundCalc){
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 }
619
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 }
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);
646 fHistList->Add(fh1Nch);
647 fHistList->Add(fh1Centrality);
648 fHistList->Add(fh1CentralitySelect);
649 fHistList->Add(fh1CentralityPhySel);
650 fHistList->Add(fh1Z);
651 fHistList->Add(fh1ZSelect);
652 fHistList->Add(fh1ZPhySel);
653 if(fNRandomCones>0&&fUseBackgroundCalc){
654 for(int i = 0;i<3;i++){
655 fHistList->Add(fh1BiARandomCones[i]);
656 fHistList->Add(fh1BiARandomConesRan[i]);
657 }
658 }
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
666 fHistList->Add(fh2NRecJetsPt);
667 fHistList->Add(fh2NRecTracksPt);
668 fHistList->Add(fh2NConstPt);
669 fHistList->Add(fh2NConstLeadingPt);
670 fHistList->Add(fh2PtNch);
671 fHistList->Add(fh2PtNchRan);
672 fHistList->Add(fh2PtNchN);
673 fHistList->Add(fh2PtNchNRan);
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);
693 }
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
714 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
715
006b2a30 716 FitMomentumResolution();
717
d89b8229 718}
719
720void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
721{
722
723 // handle and reset the output jet branch
724
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();
730
731 AliAODJetEventBackground* externalBackground = 0;
732 if(!externalBackground&&fBackgroundBranch.Length()){
733 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
734 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
735 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
736 }
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
d89b8229 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;
d89b8229 764
765 Bool_t selectEvent = false;
766 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
767
768 Float_t cent = 0;
769 Float_t zVtx = 0;
770 Int_t cenClass = -1;
771 if(fAOD){
772 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
773 TString vtxTitle(vtxAOD->GetTitle());
774 zVtx = vtxAOD->GetZ();
775
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;
781 if(physicsSelection){
782 fh1CentralityPhySel->Fill(cent);
783 fh1ZPhySel->Fill(zVtx);
784 }
785
786 if(fEventSelection){
787 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
788 Float_t yvtx = vtxAOD->GetY();
789 Float_t xvtx = vtxAOD->GetX();
790 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
791 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
792 if(physicsSelection){
793 selectEvent = true;
794 }
795 }
796 }
797 if(fCentCutUp>0){
798 if(cent<fCentCutLo||cent>fCentCutUp){
799 selectEvent = false;
800 }
801 }
802 }else{
803 selectEvent = true;
804 }
805 }
806
807
808 if(!selectEvent){
809 PostData(1, fHistList);
810 return;
811 }
812 fh1Centrality->Fill(cent);
813 fh1Z->Fill(zVtx);
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);
829 Float_t nCh = recParticles.GetEntries();
830 fh1Nch->Fill(nCh);
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;
839
840 // Generate the random cones
841
842 AliAODJet vTmpRan(1,0,0,1);
843 for(int i = 0; i < recParticles.GetEntries(); i++){
844 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
006b2a30 845
d89b8229 846 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
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 }
d89b8229 957
958 // the randomized input changes eta and phi, but keeps the p_T
959 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
960 Double_t pT = vp->Pt();
961 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
962 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
963
964 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
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);
971
972 jInpRan.set_user_index(i);
973 inputParticlesRecRan.push_back(jInpRan);
974 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
975 }
976
977 // fill the tref array, only needed when we write out jets
978 if(fTCAJetsOut){
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 }
d89b8229 984 }
006b2a30 985 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
d89b8229 986 }
987 }// recparticles
988
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
996 // employ setters for these...
997
998
999 // now create the object that holds info about ghosts
1000 /*
1001 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1002 // reduce CPU time...
1003 fGhostArea = 0.5;
1004 fActiveAreaRepeats = 0;
1005 }
1006 */
1007
1008 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1009 fastjet::AreaType areaType = fastjet::active_area;
1010 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1011 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1012 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1013
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
1023 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1024 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1025
1026
1027 fh1NJetsRec->Fill(sortedJets.size());
1028
1029 // loop over all jets an fill information, first on is the leading jet
1030
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());
1035 Double_t area = clustSeq.area(sortedJets[0]);
1036 leadingJet.SetEffArea(area,0);
1037 Float_t pt = leadingJet.Pt();
1038 Int_t nAodOutJets = 0;
1039 Int_t nAodOutTracks = 0;
1040 AliAODJet *aodOutJet = 0;
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();
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();
d89b8229 1063 }
1064 pt = leadingJet.Pt() - pTback;
1065 // correlation of leading jet with tracks
1066 TIterator *recIter = recParticles.MakeIterator();
1067 recIter->Reset();
1068 AliVParticle *tmpRecTrack = 0;
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();
1077 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1078 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1079 if(cenClass>=0){
1080 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1081 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1082 }
1083
1084 }
1085
1086
1087 TLorentzVector vecareab;
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();
1093
1094 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1095 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1096 aodOutJet->GetRefTracks()->Clear();
1097 Double_t area1 = clustSeq.area(sortedJets[j]);
1098 aodOutJet->SetEffArea(area1,0);
1099 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1100 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1101 aodOutJet->SetVectorAreaCharged(&vecareab);
1102 }
1103
1104
1105 Float_t tmpPtBack = 0;
1106 if(externalBackground){
1107 // carefull has to be filled in a task before
1108 // todo, ReArrange to the botom
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);
1115 // Fill Spectra with constituentsemacs
1116 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1117
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
1123
d7396b04 1124 AliVParticle *partLead = 0;
1125 Float_t ptLead = -1;
1126
d89b8229 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;
d89b8229 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;
d89b8229 1140 }
1141 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1142 }
ae1bdbd2 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 }
d89b8229 1155
1156 // correlation
1157 Float_t tmpPhi = tmpRec.Phi();
1158 Float_t tmpEta = tmpRec.Eta();
1159 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
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);
1176 if(cenClass>=0){
1177 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1178 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1179 }
1180 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1181 }// loop over reconstructed jets
1182 delete recIter;
1183
1184
1185
1186 // Add the random cones...
1187 if(fNRandomCones>0&&fTCARandomConesOut){
1188 // create a random jet within the acceptance
1189 Double_t etaMax = fTrackEtaWindow - fRparam;
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;
1204 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
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.);
d89b8229 1224 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1225 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
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){
1239 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1240 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
95808c2b 1241 if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
d89b8229 1242 }
1243 }
1244 }// add up energy in cone
95808c2b 1245
d89b8229 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();
1249 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
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){
1263 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1264 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
95808c2b 1265 if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
d89b8229 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
d89b8229 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();
1283 if(pTC<=0)pTC = 0.001; // for almost empty events
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 }
1295 if(fTCARandomConesOutRan){
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();
1304 if(pTC<=0)pTC = 0.001;// for almost empty events
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
1317
1318 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1319
1320
1321
1322
1323
1324 if(fAODJetBackgroundOut){
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.;
1333
1334 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1335 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1336
1337 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1338 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1339
1340 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1341 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1342 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1343 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1344
1345 }
1346 }
1347
1348
1349
1350
1351
1352 // fill track information
1353 Int_t nTrackOver = recParticles.GetSize();
1354 // do the same for tracks and jets
1355
1356 if(nTrackOver>0){
1357 TIterator *recIter = recParticles.MakeIterator();
1358 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1359 Float_t pt = tmpRec->Pt();
1360
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
1407
1408 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1409
1410 // fill the jet information from random track
1411 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1412 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
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();
1420
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;
1426 TLorentzVector vecarearanb;
1427
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();
1455 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1456 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1457 }
1458
1459 Int_t nAodOutJetsRan = 0;
1460 AliAODJet *aodOutJetRan = 0;
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
1466 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1467 fh1NConstRecRan->Fill(constituents.size());
1468 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1469 fh2PtNchRan->Fill(nCh,tmpPt);
1470 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1471
1472
1473 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1474 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1475 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1476 aodOutJetRan->GetRefTracks()->Clear();
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 }
1483
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 }
1495
1496
1497 if(fAODJetBackgroundOut){
1498 Double_t bkg3=0.;
1499 Double_t sigma3=0.;
1500 Double_t meanarea3=0.;
1501 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1502 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1503 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1504 /*
1505 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1506 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1507 */
1508 }
1509
1510
1511
1512 }
1513
1514
1515 // do the event selection if activated
1516 if(fJetTriggerPtCut>0){
1517 bool select = false;
1518 Float_t minPt = fJetTriggerPtCut;
1519 /*
1520 // hard coded for now ...
1521 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1522 if(cent<10)minPt = 50;
1523 else if(cent<30)minPt = 42;
1524 else if(cent<50)minPt = 28;
1525 else if(cent<80)minPt = 18;
1526 */
1527 float rho = 0;
1528 if(externalBackground)rho = externalBackground->GetBackground(2);
1529 if(fTCAJetsOut){
1530 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1531 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1532 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1533 if(ptSub>=minPt){
1534 select = true;
1535 break;
1536 }
1537 }
1538 }
1539
1540 if(select){
1541 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1542 fh1CentralitySelect->Fill(cent);
1543 fh1ZSelect->Fill(zVtx);
1544 aodH->SetFillAOD(kTRUE);
1545 }
1546 }
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 }
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
d89b8229 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
d89b8229 1576 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1577
1578 Int_t iCount = 0;
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){
1585 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1586 return iCount;
1587 }
d7396b04 1588
d89b8229 1589 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1590 AliAODTrack *tr = aod->GetTrack(it);
1591 Bool_t bGood = false;
1592 if(fFilterType == 0)bGood = true;
1593 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1594 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1595 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
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());
1608 list->Add(tr);
1609 iCount++;
1610 }
1611 }
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"));
1621 if(!aodExtraTracks)return iCount;
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);
1627 if(!trackAOD)continue;
1628 Bool_t bGood = false;
1629 if(fFilterType == 0)bGood = true;
1630 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1631 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1632 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1633 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1634 if(trackAOD->Pt()<fTrackPtCut) continue;
1635 list->Add(trackAOD);
1636 iCount++;
1637 }
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){
1648 if(part->Pt()<fTrackPtCut)continue;
1649 list->Add(part);
1650 iCount++;
1651 }
1652 else if(type == kTrackKineCharged){
1653 if(part->Particle()->GetPDG()->Charge()==0)continue;
1654 if(part->Pt()<fTrackPtCut)continue;
1655 list->Add(part);
1656 iCount++;
1657 }
1658 }
1659 }
1660 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1661 AliAODEvent *aod = 0;
1662 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
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){
1668 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1669 if(!part->IsPhysicalPrimary())continue;
1670 if(type == kTrackAODMCAll){
1671 if(part->Pt()<fTrackPtCut)continue;
1672 list->Add(part);
1673 iCount++;
1674 }
1675 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1676 if(part->Charge()==0)continue;
1677 if(part->Pt()<fTrackPtCut)continue;
1678 if(kTrackAODMCCharged){
1679 list->Add(part);
1680 }
1681 else {
1682 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1683 list->Add(part);
1684 }
1685 iCount++;
1686 }
1687 }
1688 }// AODMCparticle
1689 list->Sort();
1690 return iCount;
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
d89b8229 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*/