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