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