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