]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - JETAN/AliAnalysisTaskJetCluster.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Changes for reading Embedded tracks in the cluster task
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetCluster.cxx
... / ...
CommitLineData
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>
35#include <TList.h>
36#include <TLorentzVector.h>
37#include <TClonesArray.h>
38#include "TDatabasePDG.h"
39
40#include "AliAnalysisTaskJetCluster.h"
41#include "AliAnalysisManager.h"
42#include "AliJetFinder.h"
43#include "AliJetHeader.h"
44#include "AliJetReader.h"
45#include "AliESDEvent.h"
46#include "AliAODEvent.h"
47#include "AliAODHandler.h"
48#include "AliAODTrack.h"
49#include "AliAODJet.h"
50#include "AliAODMCParticle.h"
51#include "AliMCEventHandler.h"
52#include "AliMCEvent.h"
53#include "AliStack.h"
54#include "AliGenPythiaEventHeader.h"
55#include "AliJetKineReaderHeader.h"
56#include "AliGenCocktailEventHeader.h"
57#include "AliInputEventHandler.h"
58#include "AliAODJetEventBackground.h"
59
60#include "fastjet/PseudoJet.hh"
61#include "fastjet/ClusterSequenceArea.hh"
62#include "fastjet/AreaDefinition.hh"
63#include "fastjet/JetDefinition.hh"
64// get info on how fastjet was configured
65#include "fastjet/config.h"
66
67
68ClassImp(AliAnalysisTaskJetCluster)
69
70AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
71 delete fRef;
72 delete fRandom;
73}
74
75AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
76 fAOD(0x0),
77 fAODExtension(0x0),
78 fRef(new TRefArray),
79 fUseAODTrackInput(kFALSE),
80 fUseAODMCInput(kFALSE),
81 fUseGlobalSelection(kFALSE),
82 fUseBackgroundCalc(kFALSE),
83 fFilterMask(0),
84 fTrackTypeRec(kTrackUndef),
85 fTrackTypeGen(kTrackUndef),
86 fNSkipLeadingRan(0),
87 fNRandomCones(0),
88 fAvgTrials(1),
89 fExternalWeight(1),
90 fRecEtaWindow(0.5),
91 fTrackPtCut(0.),
92 fJetOutputMinPt(1),
93 fJetTriggerPtCut(0),
94 fCentCutUp(0),
95 fCentCutLo(0),
96 fNonStdBranch(""),
97 fBackgroundBranch(""),
98 fNonStdFile(""),
99 fRparam(1.0),
100 fAlgorithm(fastjet::kt_algorithm),
101 fStrategy(fastjet::Best),
102 fRecombScheme(fastjet::BIpt_scheme),
103 fAreaType(fastjet::active_area),
104 fGhostArea(0.01),
105 fActiveAreaRepeats(1),
106 fGhostEtamax(1.5),
107 fRandom(0),
108 fh1Xsec(0x0),
109 fh1Trials(0x0),
110 fh1PtHard(0x0),
111 fh1PtHardNoW(0x0),
112 fh1PtHardTrials(0x0),
113 fh1NJetsRec(0x0),
114 fh1NConstRec(0x0),
115 fh1NConstLeadingRec(0x0),
116 fh1PtJetsRecIn(0x0),
117 fh1PtJetsLeadingRecIn(0x0),
118 fh1PtJetConstRec(0x0),
119 fh1PtJetConstLeadingRec(0x0),
120 fh1PtTracksRecIn(0x0),
121 fh1PtTracksLeadingRecIn(0x0),
122 fh1NJetsRecRan(0x0),
123 fh1NConstRecRan(0x0),
124 fh1PtJetsLeadingRecInRan(0x0),
125 fh1NConstLeadingRecRan(0x0),
126 fh1PtJetsRecInRan(0x0),
127 fh1PtTracksGenIn(0x0),
128 fh1Nch(0x0),
129 fh1CentralityPhySel(0x0),
130 fh1Centrality(0x0),
131 fh1CentralitySelect(0x0),
132 fh1ZPhySel(0x0),
133 fh1Z(0x0),
134 fh1ZSelect(0x0),
135 fh2NRecJetsPt(0x0),
136 fh2NRecTracksPt(0x0),
137 fh2NConstPt(0x0),
138 fh2NConstLeadingPt(0x0),
139 fh2JetPhiEta(0x0),
140 fh2LeadingJetPhiEta(0x0),
141 fh2JetEtaPt(0x0),
142 fh2LeadingJetEtaPt(0x0),
143 fh2TrackEtaPt(0x0),
144 fh2LeadingTrackEtaPt(0x0),
145 fh2JetsLeadingPhiEta(0x0),
146 fh2JetsLeadingPhiPt(0x0),
147 fh2TracksLeadingPhiEta(0x0),
148 fh2TracksLeadingPhiPt(0x0),
149 fh2TracksLeadingJetPhiPt(0x0),
150 fh2JetsLeadingPhiPtW(0x0),
151 fh2TracksLeadingPhiPtW(0x0),
152 fh2TracksLeadingJetPhiPtW(0x0),
153 fh2NRecJetsPtRan(0x0),
154 fh2NConstPtRan(0x0),
155 fh2NConstLeadingPtRan(0x0),
156 fh2PtNch(0x0),
157 fh2PtNchRan(0x0),
158 fh2PtNchN(0x0),
159 fh2PtNchNRan(0x0),
160 fh2TracksLeadingJetPhiPtRan(0x0),
161 fh2TracksLeadingJetPhiPtWRan(0x0),
162 fHistList(0x0)
163{
164 for(int i = 0;i<3;i++){
165 fh1BiARandomCones[i] = 0;
166 fh1BiARandomConesRan[i] = 0;
167 }
168 for(int i = 0;i<kMaxCent;i++){
169 fh2JetsLeadingPhiPtC[i] = 0;
170 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
171 fh2TracksLeadingJetPhiPtC[i] = 0;
172 fh2TracksLeadingJetPhiPtWC[i] = 0;
173 }
174}
175
176AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
177 AliAnalysisTaskSE(name),
178 fAOD(0x0),
179 fAODExtension(0x0),
180 fRef(new TRefArray),
181 fUseAODTrackInput(kFALSE),
182 fUseAODMCInput(kFALSE),
183 fUseGlobalSelection(kFALSE),
184 fUseBackgroundCalc(kFALSE),
185 fFilterMask(0),
186 fTrackTypeRec(kTrackUndef),
187 fTrackTypeGen(kTrackUndef),
188 fNSkipLeadingRan(0),
189 fNRandomCones(0),
190 fAvgTrials(1),
191 fExternalWeight(1),
192 fRecEtaWindow(0.5),
193 fTrackPtCut(0.),
194 fJetOutputMinPt(1),
195 fJetTriggerPtCut(0),
196 fCentCutUp(0),
197 fCentCutLo(0),
198 fNonStdBranch(""),
199 fBackgroundBranch(""),
200 fNonStdFile(""),
201 fRparam(1.0),
202 fAlgorithm(fastjet::kt_algorithm),
203 fStrategy(fastjet::Best),
204 fRecombScheme(fastjet::BIpt_scheme),
205 fAreaType(fastjet::active_area),
206 fGhostArea(0.01),
207 fActiveAreaRepeats(1),
208 fGhostEtamax(1.5),
209 fRandom(0),
210 fh1Xsec(0x0),
211 fh1Trials(0x0),
212 fh1PtHard(0x0),
213 fh1PtHardNoW(0x0),
214 fh1PtHardTrials(0x0),
215 fh1NJetsRec(0x0),
216 fh1NConstRec(0x0),
217 fh1NConstLeadingRec(0x0),
218 fh1PtJetsRecIn(0x0),
219 fh1PtJetsLeadingRecIn(0x0),
220 fh1PtJetConstRec(0x0),
221 fh1PtJetConstLeadingRec(0x0),
222 fh1PtTracksRecIn(0x0),
223 fh1PtTracksLeadingRecIn(0x0),
224 fh1NJetsRecRan(0x0),
225 fh1NConstRecRan(0x0),
226 fh1PtJetsLeadingRecInRan(0x0),
227 fh1NConstLeadingRecRan(0x0),
228 fh1PtJetsRecInRan(0x0),
229 fh1PtTracksGenIn(0x0),
230 fh1Nch(0x0),
231 fh1CentralityPhySel(0x0),
232 fh1Centrality(0x0),
233 fh1CentralitySelect(0x0),
234 fh1ZPhySel(0x0),
235 fh1Z(0x0),
236 fh1ZSelect(0x0),
237 fh2NRecJetsPt(0x0),
238 fh2NRecTracksPt(0x0),
239 fh2NConstPt(0x0),
240 fh2NConstLeadingPt(0x0),
241 fh2JetPhiEta(0x0),
242 fh2LeadingJetPhiEta(0x0),
243 fh2JetEtaPt(0x0),
244 fh2LeadingJetEtaPt(0x0),
245 fh2TrackEtaPt(0x0),
246 fh2LeadingTrackEtaPt(0x0),
247 fh2JetsLeadingPhiEta(0x0),
248 fh2JetsLeadingPhiPt(0x0),
249 fh2TracksLeadingPhiEta(0x0),
250 fh2TracksLeadingPhiPt(0x0),
251 fh2TracksLeadingJetPhiPt(0x0),
252 fh2JetsLeadingPhiPtW(0x0),
253 fh2TracksLeadingPhiPtW(0x0),
254 fh2TracksLeadingJetPhiPtW(0x0),
255 fh2NRecJetsPtRan(0x0),
256 fh2NConstPtRan(0x0),
257 fh2NConstLeadingPtRan(0x0),
258 fh2PtNch(0x0),
259 fh2PtNchRan(0x0),
260 fh2PtNchN(0x0),
261 fh2PtNchNRan(0x0),
262 fh2TracksLeadingJetPhiPtRan(0x0),
263 fh2TracksLeadingJetPhiPtWRan(0x0),
264 fHistList(0x0)
265{
266 for(int i = 0;i<3;i++){
267 fh1BiARandomCones[i] = 0;
268 fh1BiARandomConesRan[i] = 0;
269 }
270 for(int i = 0;i<kMaxCent;i++){
271 fh2JetsLeadingPhiPtC[i] = 0;
272 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
273 fh2TracksLeadingJetPhiPtC[i] = 0;
274 fh2TracksLeadingJetPhiPtWC[i] = 0;
275 }
276 DefineOutput(1, TList::Class());
277}
278
279
280
281Bool_t AliAnalysisTaskJetCluster::Notify()
282{
283 //
284 // Implemented Notify() to read the cross sections
285 // and number of trials from pyxsec.root
286 //
287 return kTRUE;
288}
289
290void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
291{
292
293 //
294 // Create the output container
295 //
296
297 fRandom = new TRandom3(0);
298
299
300 // Connect the AOD
301
302
303 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
304
305
306
307 if(fNonStdBranch.Length()!=0)
308 {
309 // only create the output branch if we have a name
310 // Create a new branch for jets...
311 // -> cleared in the UserExec....
312 // here we can also have the case that the brnaches are written to a separate file
313
314 TClonesArray *tca = new TClonesArray("AliAODJet", 0);
315 tca->SetName(fNonStdBranch.Data());
316 AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
317
318
319 TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
320 tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
321 AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
322
323 if(fUseBackgroundCalc){
324 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
325 AliAODJetEventBackground* evBkg = new AliAODJetEventBackground();
326 evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
327 AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data());
328 }
329 }
330 // create the branch for the random cones with the same R
331 TString cName = Form("%sRandomCone",fNonStdBranch.Data());
332
333 if(fNRandomCones>0){
334 if(!AODEvent()->FindListObject(cName.Data())){
335 TClonesArray *tcaC = new TClonesArray("AliAODJet", 0);
336 tcaC->SetName(cName.Data());
337 AddAODBranch("TClonesArray",&tcaC,fNonStdFile.Data());
338 }
339 // create the branch with the random for the random cones on the random event
340 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
341 if(!AODEvent()->FindListObject(cName.Data())){
342 TClonesArray *tcaCran = new TClonesArray("AliAODJet", 0);
343 tcaCran->SetName(cName.Data());
344 AddAODBranch("TClonesArray",&tcaCran,fNonStdFile.Data());
345 }
346 }
347
348 if(fNonStdFile.Length()!=0){
349 //
350 // case that we have an AOD extension we need to fetch the jets from the extended output
351 // we identifay the extension aod event by looking for the branchname
352 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
353 TObjArray* extArray = aodH->GetExtensions();
354 if (extArray) {
355 TIter next(extArray);
356 while ((fAODExtension=(AliAODExtension*)next())){
357 TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data());
358 if(fDebug>10){
359 Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__);
360 fAODExtension->GetAOD()->Dump();
361 }
362 if(obj){
363 if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data());
364 break;
365 }
366 fAODExtension = 0;
367 }
368 }
369 }
370 }
371
372
373 if(!fHistList)fHistList = new TList();
374 fHistList->SetOwner();
375
376 Bool_t oldStatus = TH1::AddDirectoryStatus();
377 TH1::AddDirectory(kFALSE);
378
379 //
380 // Histogram
381
382 const Int_t nBinPt = 200;
383 Double_t binLimitsPt[nBinPt+1];
384 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
385 if(iPt == 0){
386 binLimitsPt[iPt] = 0.0;
387 }
388 else {// 1.0
389 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
390 }
391 }
392
393 const Int_t nBinPhi = 90;
394 Double_t binLimitsPhi[nBinPhi+1];
395 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
396 if(iPhi==0){
397 binLimitsPhi[iPhi] = -1.*TMath::Pi();
398 }
399 else{
400 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
401 }
402 }
403
404
405
406 const Int_t nBinEta = 40;
407 Double_t binLimitsEta[nBinEta+1];
408 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
409 if(iEta==0){
410 binLimitsEta[iEta] = -2.0;
411 }
412 else{
413 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
414 }
415 }
416
417 const int nChMax = 4000;
418
419 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
420 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
421
422 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
423 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
424
425
426 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
427 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
428
429 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
430 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
431 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
432 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
433
434
435 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
436 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
437 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
438
439 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
440 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
441 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
442 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
443 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
444 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
445 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
446 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
447 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
448 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
449
450 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
451 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
452 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
453
454 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
455 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
456 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
457
458 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
459 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
460 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
461 //
462
463
464 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
465 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
466 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
467 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
468
469 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
470 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);
471 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);
472 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);
473
474
475
476 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
477 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
478 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
479 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
480
481 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
482 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
483 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
484 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
485
486 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
487 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
488 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
489 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
490
491
492
493 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
494 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
495 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
496 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
497 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
498 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
499 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
500 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
501 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
502 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
503 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
504 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
505
506 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
507 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
508 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
509 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
510
511 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
512 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
513 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
514 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
515
516
517 if(fNRandomCones>0&&fUseBackgroundCalc){
518 for(int i = 0;i<3;i++){
519 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
520 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
521 }
522 }
523
524 for(int i = 0;i < kMaxCent;i++){
525 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
526 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
527 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
528 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
529 }
530
531 const Int_t saveLevel = 3; // large save level more histos
532 if(saveLevel>0){
533 fHistList->Add(fh1Xsec);
534 fHistList->Add(fh1Trials);
535
536 fHistList->Add(fh1NJetsRec);
537 fHistList->Add(fh1NConstRec);
538 fHistList->Add(fh1NConstLeadingRec);
539 fHistList->Add(fh1PtJetsRecIn);
540 fHistList->Add(fh1PtJetsLeadingRecIn);
541 fHistList->Add(fh1PtTracksRecIn);
542 fHistList->Add(fh1PtTracksLeadingRecIn);
543 fHistList->Add(fh1PtJetConstRec);
544 fHistList->Add(fh1PtJetConstLeadingRec);
545 fHistList->Add(fh1NJetsRecRan);
546 fHistList->Add(fh1NConstRecRan);
547 fHistList->Add(fh1PtJetsLeadingRecInRan);
548 fHistList->Add(fh1NConstLeadingRecRan);
549 fHistList->Add(fh1PtJetsRecInRan);
550 fHistList->Add(fh1Nch);
551 fHistList->Add(fh1Centrality);
552 fHistList->Add(fh1CentralitySelect);
553 fHistList->Add(fh1CentralityPhySel);
554 fHistList->Add(fh1Z);
555 fHistList->Add(fh1ZSelect);
556 fHistList->Add(fh1ZPhySel);
557 if(fNRandomCones&&fUseBackgroundCalc){
558 for(int i = 0;i<3;i++){
559 fHistList->Add(fh1BiARandomCones[i]);
560 fHistList->Add(fh1BiARandomConesRan[i]);
561 }
562 }
563 for(int i = 0;i < kMaxCent;i++){
564 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
565 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
566 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
567 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
568 }
569
570 fHistList->Add(fh2NRecJetsPt);
571 fHistList->Add(fh2NRecTracksPt);
572 fHistList->Add(fh2NConstPt);
573 fHistList->Add(fh2NConstLeadingPt);
574 fHistList->Add(fh2PtNch);
575 fHistList->Add(fh2PtNchRan);
576 fHistList->Add(fh2PtNchN);
577 fHistList->Add(fh2PtNchNRan);
578 fHistList->Add(fh2JetPhiEta);
579 fHistList->Add(fh2LeadingJetPhiEta);
580 fHistList->Add(fh2JetEtaPt);
581 fHistList->Add(fh2LeadingJetEtaPt);
582 fHistList->Add(fh2TrackEtaPt);
583 fHistList->Add(fh2LeadingTrackEtaPt);
584 fHistList->Add(fh2JetsLeadingPhiEta );
585 fHistList->Add(fh2JetsLeadingPhiPt);
586 fHistList->Add(fh2TracksLeadingPhiEta);
587 fHistList->Add(fh2TracksLeadingPhiPt);
588 fHistList->Add(fh2TracksLeadingJetPhiPt);
589 fHistList->Add(fh2JetsLeadingPhiPtW);
590 fHistList->Add(fh2TracksLeadingPhiPtW);
591 fHistList->Add(fh2TracksLeadingJetPhiPtW);
592 fHistList->Add(fh2NRecJetsPtRan);
593 fHistList->Add(fh2NConstPtRan);
594 fHistList->Add(fh2NConstLeadingPtRan);
595 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
596 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
597 }
598
599 // =========== Switch on Sumw2 for all histos ===========
600 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
601 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
602 if (h1){
603 h1->Sumw2();
604 continue;
605 }
606 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
607 if(hn)hn->Sumw2();
608 }
609 TH1::AddDirectory(oldStatus);
610}
611
612void AliAnalysisTaskJetCluster::Init()
613{
614 //
615 // Initialization
616 //
617
618 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
619
620}
621
622void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
623{
624
625 if(fUseGlobalSelection){
626 // no selection by the service task, we continue
627 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
628 PostData(1, fHistList);
629 return;
630 }
631
632
633
634 // handle and reset the output jet branch
635 // only need this once
636 TClonesArray* jarray = 0;
637 TClonesArray* jarrayran = 0;
638 TClonesArray* rConeArray = 0;
639 TClonesArray* rConeArrayRan = 0;
640 AliAODJetEventBackground* evBkg = 0;
641 if(fNonStdBranch.Length()!=0) {
642 if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
643 if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
644 if(jarray)jarray->Delete(); // this is our responsibility, clear before filling again
645 if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
646 if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
647 if(jarrayran)jarrayran->Delete(); // this is our responsibility, clear before filling again
648
649 if(fUseBackgroundCalc){
650 if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
651 if(!evBkg) evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
652 if(evBkg)evBkg->Reset();
653 }
654
655 if(fNRandomCones>0){
656 TString cName = Form("%sRandomCone",fNonStdBranch.Data());
657 if(AODEvent())rConeArray = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
658 if(!rConeArray)rConeArray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
659 if(rConeArray)rConeArray->Delete();
660
661 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
662 if(AODEvent())rConeArrayRan = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
663 if(!rConeArrayRan)rConeArrayRan = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
664 if(rConeArrayRan)rConeArrayRan->Delete();
665 }
666 }
667
668 AliAODJetEventBackground* externalBackground = 0;
669 if(!externalBackground&&fBackgroundBranch.Length()){
670 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
671 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
672 }
673 //
674 // Execute analysis for current event
675 //
676 AliESDEvent *fESD = 0;
677 if(fUseAODTrackInput){
678 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
679 if(!fAOD){
680 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
681 return;
682 }
683 // fethc the header
684 }
685 else{
686 // assume that the AOD is in the general output...
687 fAOD = AODEvent();
688 if(!fAOD){
689 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
690 return;
691 }
692 if(fDebug>0){
693 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
694 }
695 }
696
697 Bool_t selectEvent = false;
698 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
699
700 Float_t cent = 0;
701 Float_t zVtx = 0;
702 Int_t cenClass = -1;
703 if(fAOD){
704 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
705 TString vtxTitle(vtxAOD->GetTitle());
706 zVtx = vtxAOD->GetZ();
707
708 cent = fAOD->GetHeader()->GetCentrality();
709 if(cent<10)cenClass = 0;
710 else if(cent<30)cenClass = 1;
711 else if(cent<50)cenClass = 2;
712 else if(cent<80)cenClass = 3;
713 if(physicsSelection){
714 fh1CentralityPhySel->Fill(cent);
715 fh1ZPhySel->Fill(zVtx);
716 }
717
718
719 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
720 Float_t yvtx = vtxAOD->GetY();
721 Float_t xvtx = vtxAOD->GetX();
722 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
723 if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
724 if(physicsSelection){
725 selectEvent = true;
726 }
727 }
728 }
729 if(fCentCutUp>0){
730 if(cent<fCentCutLo||cent>fCentCutUp){
731 selectEvent = false;
732 }
733 }
734
735 }
736 if(!selectEvent){
737 PostData(1, fHistList);
738 return;
739 }
740 fh1Centrality->Fill(cent);
741 fh1Z->Fill(zVtx);
742 fh1Trials->Fill("#sum{ntrials}",1);
743
744
745 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
746
747 // ==== General variables needed
748
749
750
751 // we simply fetch the tracks/mc particles as a list of AliVParticles
752
753 TList recParticles;
754 TList genParticles;
755
756 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
757 Float_t nCh = recParticles.GetEntries();
758 fh1Nch->Fill(nCh);
759 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
760 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
761 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
762
763 // find the jets....
764
765 vector<fastjet::PseudoJet> inputParticlesRec;
766 vector<fastjet::PseudoJet> inputParticlesRecRan;
767
768 // Generate the random cones
769
770 AliAODJet vTmpRan(1,0,0,1);
771 for(int i = 0; i < recParticles.GetEntries(); i++){
772 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
773 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
774 // we take total momentum here
775 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
776 jInp.set_user_index(i);
777 inputParticlesRec.push_back(jInp);
778
779 // the randomized input changes eta and phi, but keeps the p_T
780 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
781 Double_t pT = vp->Pt();
782 Double_t eta = 1.8 * fRandom->Rndm() - 0.9;
783 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
784
785 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
786 Double_t pZ = pT/TMath::Tan(theta);
787
788 Double_t pX = pT * TMath::Cos(phi);
789 Double_t pY = pT * TMath::Sin(phi);
790 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
791 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
792
793 jInpRan.set_user_index(i);
794 inputParticlesRecRan.push_back(jInpRan);
795 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
796 }
797
798 // fill the tref array, only needed when we write out jets
799 if(jarray){
800 if(i == 0){
801 fRef->Delete(); // make sure to delete before placement new...
802 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
803 }
804 fRef->Add(vp);
805 }
806 }// recparticles
807
808 if(inputParticlesRec.size()==0){
809 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
810 PostData(1, fHistList);
811 return;
812 }
813
814 // run fast jet
815 // employ setters for these...
816
817
818 // now create the object that holds info about ghosts
819 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
820 // reduce CPU time...
821 fGhostArea = 0.5;
822 fActiveAreaRepeats = 0;
823 }
824
825 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
826 fastjet::AreaType areaType = fastjet::active_area;
827 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
828 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
829 vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
830 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
831
832 //range where to compute background
833 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
834 phiMin = 0;
835 phiMax = 2*TMath::Pi();
836 rapMax = fGhostEtamax - fRparam;
837 rapMin = - fGhostEtamax + fRparam;
838 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
839
840
841
842 inclusiveJets = clustSeq.inclusive_jets();
843 sortedJets = sorted_by_pt(inclusiveJets);
844
845 fh1NJetsRec->Fill(sortedJets.size());
846
847 // loop over all jets an fill information, first on is the leading jet
848
849 Int_t nRecOver = inclusiveJets.size();
850 Int_t nRec = inclusiveJets.size();
851 if(inclusiveJets.size()>0){
852 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
853 Double_t area = clustSeq.area(sortedJets[0]);
854 leadingJet.SetEffArea(area,0);
855 Float_t pt = leadingJet.Pt();
856 Int_t nAodOutJets = 0;
857 Int_t nAodOutTracks = 0;
858 AliAODJet *aodOutJet = 0;
859
860 Int_t iCount = 0;
861 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
862 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
863 while(pt<ptCut&&iCount<nRec){
864 nRecOver--;
865 iCount++;
866 if(iCount<nRec){
867 pt = sortedJets[iCount].perp();
868 }
869 }
870 if(nRecOver<=0)break;
871 fh2NRecJetsPt->Fill(ptCut,nRecOver);
872 }
873 Float_t phi = leadingJet.Phi();
874 if(phi<0)phi+=TMath::Pi()*2.;
875 Float_t eta = leadingJet.Eta();
876 Float_t pTback = 0;
877 if(externalBackground){
878 // carefull has to be filled in a task before
879 // todo, ReArrange to the botom
880 pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
881 }
882 pt = leadingJet.Pt() - pTback;
883 // correlation of leading jet with tracks
884 TIterator *recIter = recParticles.MakeIterator();
885 recIter->Reset();
886 AliVParticle *tmpRecTrack = 0;
887 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
888 Float_t tmpPt = tmpRecTrack->Pt();
889 // correlation
890 Float_t tmpPhi = tmpRecTrack->Phi();
891 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
892 Float_t dPhi = phi - tmpPhi;
893 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
894 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
895 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
896 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
897 if(cenClass>=0){
898 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
899 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
900 }
901
902 }
903
904
905 TLorentzVector vecareab;
906 for(int j = 0; j < nRec;j++){
907 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
908 aodOutJet = 0;
909 nAodOutTracks = 0;
910 Float_t tmpPt = tmpRec.Pt();
911
912 if(tmpPt>fJetOutputMinPt&&jarray){// cut on the non-background subtracted...
913 aodOutJet = new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
914 Double_t area1 = clustSeq.area(sortedJets[j]);
915 aodOutJet->SetEffArea(area1,0);
916 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
917 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
918 aodOutJet->SetVectorAreaCharged(&vecareab);
919
920
921 }
922
923
924 Float_t tmpPtBack = 0;
925 if(externalBackground){
926 // carefull has to be filled in a task before
927 // todo, ReArrange to the botom
928 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
929 }
930 tmpPt = tmpPt - tmpPtBack;
931 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
932
933 fh1PtJetsRecIn->Fill(tmpPt);
934 // Fill Spectra with constituents
935 vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
936
937 fh1NConstRec->Fill(constituents.size());
938 fh2PtNch->Fill(nCh,tmpPt);
939 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
940 fh2NConstPt->Fill(tmpPt,constituents.size());
941 // loop over constiutents and fill spectrum
942
943 // Add the jet information and the track references to the output AOD
944
945
946 if(fNRandomCones>0){
947 // create a random jet within the acceptance
948 Double_t etaMax = 0.8 - fRparam;
949 Int_t nCone = 0;
950 Int_t nConeRan = 0;
951 Double_t pTC = 1; // small number
952 for(int ir = 0;ir < fNRandomCones;ir++){
953 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
954 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
955 // massless jet
956 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
957 Double_t pZC = pTC/TMath::Tan(thetaC);
958 Double_t pXC = pTC * TMath::Cos(phiC);
959 Double_t pYC = pTC * TMath::Sin(phiC);
960 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
961 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
962 bool skip = false;
963 for(int jj = 0; jj < TMath::Min(nRec,2);jj++){// test for overlap with leading jets
964 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
965 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
966 skip = true;
967 break;
968 }
969 }
970 // test for overlap with previous cones to avoid double counting
971 for(int iic = 0;iic<ir;iic++){
972 AliAODJet *iicone = (AliAODJet*)rConeArray->At(iic);
973 if(iicone){
974 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
975 skip = true;
976 break;
977 }
978 }
979 }
980 if(skip)continue;
981 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
982 if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRecC);
983 if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRecC);
984 }// random cones
985
986
987 // loop over the reconstructed particles and add up the pT in the random cones
988 // maybe better to loop over randomized particles not in the real jets...
989 // but this by definition brings dow average energy in the whole event
990 AliAODJet vTmpRanR(1,0,0,1);
991 for(int i = 0; i < recParticles.GetEntries(); i++){
992 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
993 if(rConeArray){
994 for(int ir = 0;ir < fNRandomCones;ir++){
995 AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);
996 if(jC&&jC->DeltaR(vp)<fRparam){
997 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
998 }
999 }
1000 }// add up energy in cone
1001
1002 // the randomized input changes eta and phi, but keeps the p_T
1003 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1004 Double_t pTR = vp->Pt();
1005 Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
1006 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1007
1008 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1009 Double_t pZR = pTR/TMath::Tan(thetaR);
1010
1011 Double_t pXR = pTR * TMath::Cos(phiR);
1012 Double_t pYR = pTR * TMath::Sin(phiR);
1013 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1014 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1015 if(rConeArrayRan){
1016 for(int ir = 0;ir < fNRandomCones;ir++){
1017 AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);
1018 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1019 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1020 }
1021 }
1022 }
1023 }
1024 }// loop over recparticles
1025
1026 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1027 for(int ir = 0;ir < fNRandomCones;ir++){
1028 // rescale the momntum vectors for the random cones
1029 if(!rConeArray)continue;
1030 AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
1031 if(rC){
1032 Double_t etaC = rC->Eta();
1033 Double_t phiC = rC->Phi();
1034 // massless jet, unit vector
1035 pTC = rC->ChargedBgEnergy();
1036 if(pTC<=0)pTC = 0.1; // for almost empty events
1037 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1038 Double_t pZC = pTC/TMath::Tan(thetaC);
1039 Double_t pXC = pTC * TMath::Cos(phiC);
1040 Double_t pYC = pTC * TMath::Sin(phiC);
1041 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1042 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1043 rC->SetBgEnergy(0,0);
1044 rC->SetEffArea(jetArea,0);
1045 }
1046 rC = (AliAODJet*)rConeArrayRan->At(ir);
1047 // same wit random
1048 if(rC){
1049 Double_t etaC = rC->Eta();
1050 Double_t phiC = rC->Phi();
1051 // massless jet, unit vector
1052 pTC = rC->ChargedBgEnergy();
1053 if(pTC<=0)pTC = 0.1;// for almost empty events
1054 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1055 Double_t pZC = pTC/TMath::Tan(thetaC);
1056 Double_t pXC = pTC * TMath::Cos(phiC);
1057 Double_t pYC = pTC * TMath::Sin(phiC);
1058 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1059 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1060 rC->SetBgEnergy(0,0);
1061 rC->SetEffArea(jetArea,0);
1062 }
1063 }
1064 }// if(fUseBackgroundCalc
1065
1066 for(unsigned int ic = 0; ic < constituents.size();ic++){
1067 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1068 fh1PtJetConstRec->Fill(part->Pt());
1069 if(aodOutJet){
1070 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1071 }
1072 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1073 }
1074
1075 // correlation
1076 Float_t tmpPhi = tmpRec.Phi();
1077 Float_t tmpEta = tmpRec.Eta();
1078 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1079 if(j==0){
1080 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1081 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1082 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1083 fh1NConstLeadingRec->Fill(constituents.size());
1084 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1085 continue;
1086 }
1087 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1088 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1089 Float_t dPhi = phi - tmpPhi;
1090 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1091 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1092 Float_t dEta = eta - tmpRec.Eta();
1093 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1094 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1095 if(cenClass>=0){
1096 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1097 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1098 }
1099 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1100 }
1101 delete recIter;
1102
1103 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1104
1105
1106 if(evBkg){
1107 vector<fastjet::PseudoJet> jets2=sortedJets;
1108 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1109 Double_t bkg1=0;
1110 Double_t sigma1=0.;
1111 Double_t meanarea1=0.;
1112 Double_t bkg2=0;
1113 Double_t sigma2=0.;
1114 Double_t meanarea2=0.;
1115
1116 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1117 evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
1118
1119 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1120 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1121
1122 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1123 evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
1124 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1125 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1126
1127 }
1128 }
1129
1130
1131
1132
1133
1134 // fill track information
1135 Int_t nTrackOver = recParticles.GetSize();
1136 // do the same for tracks and jets
1137
1138 if(nTrackOver>0){
1139 TIterator *recIter = recParticles.MakeIterator();
1140 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1141 Float_t pt = tmpRec->Pt();
1142
1143 // Printf("Leading track p_t %3.3E",pt);
1144 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1145 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1146 while(pt<ptCut&&tmpRec){
1147 nTrackOver--;
1148 tmpRec = (AliVParticle*)(recIter->Next());
1149 if(tmpRec){
1150 pt = tmpRec->Pt();
1151 }
1152 }
1153 if(nTrackOver<=0)break;
1154 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1155 }
1156
1157 recIter->Reset();
1158 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1159 Float_t phi = leading->Phi();
1160 if(phi<0)phi+=TMath::Pi()*2.;
1161 Float_t eta = leading->Eta();
1162 pt = leading->Pt();
1163 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1164 Float_t tmpPt = tmpRec->Pt();
1165 Float_t tmpEta = tmpRec->Eta();
1166 fh1PtTracksRecIn->Fill(tmpPt);
1167 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1168 if(tmpRec==leading){
1169 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1170 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1171 continue;
1172 }
1173 // correlation
1174 Float_t tmpPhi = tmpRec->Phi();
1175
1176 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1177 Float_t dPhi = phi - tmpPhi;
1178 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1179 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1180 Float_t dEta = eta - tmpRec->Eta();
1181 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1182 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1183 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1184 }
1185 delete recIter;
1186 }
1187
1188 // find the random jets
1189 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
1190 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1191
1192 inclusiveJetsRan = clustSeqRan.inclusive_jets();
1193 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1194
1195 // fill the jet information from random track
1196
1197 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1198
1199 // loop over all jets an fill information, first on is the leading jet
1200
1201 Int_t nRecOverRan = inclusiveJetsRan.size();
1202 Int_t nRecRan = inclusiveJetsRan.size();
1203
1204 if(inclusiveJetsRan.size()>0){
1205 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1206 Float_t pt = leadingJet.Pt();
1207
1208 Int_t iCount = 0;
1209 TLorentzVector vecarearanb;
1210
1211 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1212 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1213 while(pt<ptCut&&iCount<nRecRan){
1214 nRecOverRan--;
1215 iCount++;
1216 if(iCount<nRecRan){
1217 pt = sortedJetsRan[iCount].perp();
1218 }
1219 }
1220 if(nRecOverRan<=0)break;
1221 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1222 }
1223 Float_t phi = leadingJet.Phi();
1224 if(phi<0)phi+=TMath::Pi()*2.;
1225 pt = leadingJet.Pt();
1226
1227 // correlation of leading jet with random tracks
1228
1229 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1230 {
1231 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1232 // correlation
1233 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1234 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1235 Float_t dPhi = phi - tmpPhi;
1236 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1237 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1238 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1239 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1240 }
1241
1242 Int_t nAodOutJetsRan = 0;
1243 AliAODJet *aodOutJetRan = 0;
1244 for(int j = 0; j < nRecRan;j++){
1245 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1246 Float_t tmpPt = tmpRec.Pt();
1247 fh1PtJetsRecInRan->Fill(tmpPt);
1248 // Fill Spectra with constituents
1249 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1250 fh1NConstRecRan->Fill(constituents.size());
1251 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1252 fh2PtNchRan->Fill(nCh,tmpPt);
1253 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1254
1255
1256 if(tmpPt>fJetOutputMinPt&&jarrayran){
1257 aodOutJetRan = new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1258 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1259
1260 aodOutJetRan->SetEffArea(arearan,0);
1261 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1262 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1263 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1264
1265 }
1266
1267
1268
1269
1270 for(unsigned int ic = 0; ic < constituents.size();ic++){
1271 // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1272 // fh1PtJetConstRec->Fill(part->Pt());
1273 if(aodOutJetRan){
1274 aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
1275 }
1276 }
1277
1278
1279
1280
1281 // correlation
1282 Float_t tmpPhi = tmpRec.Phi();
1283 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1284
1285 if(j==0){
1286 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1287 fh1NConstLeadingRecRan->Fill(constituents.size());
1288 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1289 continue;
1290 }
1291 }
1292
1293
1294 if(evBkg){
1295 Double_t bkg3=0.;
1296 Double_t sigma3=0.;
1297 Double_t meanarea3=0.;
1298 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1299 evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
1300 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1301 /*
1302 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1303 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1304 */
1305 }
1306
1307
1308
1309 }
1310
1311
1312 // do the event selection if activated
1313 if(fJetTriggerPtCut>0){
1314 bool select = false;
1315 Float_t minPt = fJetTriggerPtCut;
1316 /*
1317 // hard coded for now ...
1318 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1319 if(cent<10)minPt = 50;
1320 else if(cent<30)minPt = 42;
1321 else if(cent<50)minPt = 28;
1322 else if(cent<80)minPt = 18;
1323 */
1324 float rho = 0;
1325 if(externalBackground)rho = externalBackground->GetBackground(2);
1326 if(jarray){
1327 for(int i = 0;i < jarray->GetEntriesFast();i++){
1328 AliAODJet *jet = (AliAODJet*)jarray->At(i);
1329 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1330 if(ptSub>=minPt){
1331 select = true;
1332 break;
1333 }
1334 }
1335 }
1336
1337 if(select){
1338 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1339 fh1CentralitySelect->Fill(cent);
1340 fh1ZSelect->Fill(zVtx);
1341 aodH->SetFillAOD(kTRUE);
1342 }
1343 }
1344
1345 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1346 PostData(1, fHistList);
1347}
1348
1349void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1350{
1351// Terminate analysis
1352//
1353 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1354}
1355
1356
1357Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1358
1359 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1360
1361 Int_t iCount = 0;
1362 if(type==kTrackAOD){
1363 AliAODEvent *aod = 0;
1364 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1365 else aod = AODEvent();
1366 if(!aod){
1367 return iCount;
1368 }
1369 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1370 AliAODTrack *tr = aod->GetTrack(it);
1371 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1372 if(TMath::Abs(tr->Eta())>0.9)continue;
1373 if(tr->Pt()<fTrackPtCut)continue;
1374 list->Add(tr);
1375 iCount++;
1376 }
1377 }
1378 else if (type == kTrackKineAll||type == kTrackKineCharged){
1379 AliMCEvent* mcEvent = MCEvent();
1380 if(!mcEvent)return iCount;
1381 // we want to have alivpartilces so use get track
1382 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1383 if(!mcEvent->IsPhysicalPrimary(it))continue;
1384 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1385 if(type == kTrackKineAll){
1386 if(part->Pt()<fTrackPtCut)continue;
1387 list->Add(part);
1388 iCount++;
1389 }
1390 else if(type == kTrackKineCharged){
1391 if(part->Particle()->GetPDG()->Charge()==0)continue;
1392 if(part->Pt()<fTrackPtCut)continue;
1393 list->Add(part);
1394 iCount++;
1395 }
1396 }
1397 }
1398 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1399 AliAODEvent *aod = 0;
1400 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1401 else aod = AODEvent();
1402 if(!aod)return iCount;
1403 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1404 if(!tca)return iCount;
1405 for(int it = 0;it < tca->GetEntriesFast();++it){
1406 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1407 if(!part->IsPhysicalPrimary())continue;
1408 if(type == kTrackAODMCAll){
1409 if(part->Pt()<fTrackPtCut)continue;
1410 list->Add(part);
1411 iCount++;
1412 }
1413 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1414 if(part->Charge()==0)continue;
1415 if(part->Pt()<fTrackPtCut)continue;
1416 if(kTrackAODMCCharged){
1417 list->Add(part);
1418 }
1419 else {
1420 if(TMath::Abs(part->Eta())>0.9)continue;
1421 list->Add(part);
1422 }
1423 iCount++;
1424 }
1425 }
1426 }// AODMCparticle
1427 list->Sort();
1428 return iCount;
1429}
1430
1431/*
1432Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1433 for(int i = 0; i < particles.GetEntries(); i++){
1434 AliVParticle *vp = (AliVParticle*)particles.At(i);
1435 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1436 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1437 jInp.set_user_index(i);
1438 inputParticles.push_back(jInp);
1439 }
1440
1441 return 0;
1442
1443}
1444*/