]> 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
Correct the setting of geometry matrices from ESDs
[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
d7ca0e14 953 // float areaRandomCone = rRandomCone2 *TMath::Pi();
4a33f620 954 clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
d6e66a82 955 evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
d7ca0e14 956 /*
957 fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
9c6bd749 958 fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
82a35b1a 959
4a33f620 960 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
d6e66a82 961 evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
82a35b1a 962 fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
9c6bd749 963 fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
d7ca0e14 964 */
d6e66a82 965 }
54424110 966 }
d7ca0e14 967
82a35b1a 968
969
970
dbf053f0 971
82a35b1a 972 // fill track information
973 Int_t nTrackOver = recParticles.GetSize();
dbf053f0 974 // do the same for tracks and jets
82a35b1a 975
976 if(nTrackOver>0){
dbf053f0 977 TIterator *recIter = recParticles.MakeIterator();
978 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
979 Float_t pt = tmpRec->Pt();
82a35b1a 980
dbf053f0 981 // Printf("Leading track p_t %3.3E",pt);
982 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
983 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
984 while(pt<ptCut&&tmpRec){
985 nTrackOver--;
986 tmpRec = (AliVParticle*)(recIter->Next());
987 if(tmpRec){
988 pt = tmpRec->Pt();
989 }
990 }
991 if(nTrackOver<=0)break;
992 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
993 }
994
995 recIter->Reset();
996 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
997 Float_t phi = leading->Phi();
998 if(phi<0)phi+=TMath::Pi()*2.;
999 Float_t eta = leading->Eta();
1000 pt = leading->Pt();
1001 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1002 Float_t tmpPt = tmpRec->Pt();
1003 Float_t tmpEta = tmpRec->Eta();
1004 fh1PtTracksRecIn->Fill(tmpPt);
1005 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1006 if(tmpRec==leading){
1007 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1008 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1009 continue;
1010 }
1011 // correlation
1012 Float_t tmpPhi = tmpRec->Phi();
1013
1014 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1015 Float_t dPhi = phi - tmpPhi;
1016 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1017 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1018 Float_t dEta = eta - tmpRec->Eta();
1019 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1020 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1021 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1022 }
1023 delete recIter;
1024 }
1025
1026 // find the random jets
1027 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
d6e66a82 1028 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
dbf053f0 1029
1030 inclusiveJetsRan = clustSeqRan.inclusive_jets();
1031 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1032
1033 // fill the jet information from random track
1034
1035 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1036
1037 // loop over all jets an fill information, first on is the leading jet
1038
1039 Int_t nRecOverRan = inclusiveJetsRan.size();
1040 Int_t nRecRan = inclusiveJetsRan.size();
1041 if(inclusiveJetsRan.size()>0){
1042 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1043 Float_t pt = leadingJet.Pt();
1044
1045 Int_t iCount = 0;
1046 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1047 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1048 while(pt<ptCut&&iCount<nRecRan){
1049 nRecOverRan--;
1050 iCount++;
1051 if(iCount<nRecRan){
1052 pt = sortedJetsRan[iCount].perp();
1053 }
1054 }
1055 if(nRecOverRan<=0)break;
1056 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1057 }
1058 Float_t phi = leadingJet.Phi();
1059 if(phi<0)phi+=TMath::Pi()*2.;
1060 pt = leadingJet.Pt();
1061
1062 // correlation of leading jet with random tracks
1063
1064 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1065 {
1066 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1067 // correlation
1068 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1069 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1070 Float_t dPhi = phi - tmpPhi;
1071 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1072 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1073 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1074 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1075 }
1076
d6e66a82 1077 Int_t nAodOutJetsRan = 0;
1078 AliAODJet *aodOutJetRan = 0;
dbf053f0 1079 for(int j = 0; j < nRecRan;j++){
1080 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1081 Float_t tmpPt = tmpRec.Pt();
1082 fh1PtJetsRecInRan->Fill(tmpPt);
1083 // Fill Spectra with constituents
1084 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1085 fh1NConstRecRan->Fill(constituents.size());
1086 fh2NConstPtRan->Fill(tmpPt,constituents.size());
a0d08b6d 1087 fh2PtNchRan->Fill(nCh,tmpPt);
1088 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
d6e66a82 1089
1090
1091 if(tmpPt>fJetOutputMinPt&&jarrayran){
1092 aodOutJetRan = new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1093 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1094
1095 aodOutJetRan->SetEffArea(arearan,0); }
1096
1097
1098
1099
1100 for(unsigned int ic = 0; ic < constituents.size();ic++){
1101 // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1102 // fh1PtJetConstRec->Fill(part->Pt());
1103 if(aodOutJetRan){
1104 aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
d7ca0e14 1105 }
1106 }
d6e66a82 1107
1108
1109
1110
dbf053f0 1111 // correlation
1112 Float_t tmpPhi = tmpRec.Phi();
1113 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1114
1115 if(j==0){
1116 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1117 fh1NConstLeadingRecRan->Fill(constituents.size());
1118 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1119 continue;
1120 }
1121 }
d6e66a82 1122
1123
1124 if(evBkg){
1125 Double_t bkg3=0.;
1126 Double_t sigma3=0.;
1127 Double_t meanarea3=0.;
4a33f620 1128 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
d6e66a82 1129 evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
d7ca0e14 1130 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1131 /*
82a35b1a 1132 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
9c6bd749 1133 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
d7ca0e14 1134 */
d6e66a82 1135 }
1136
1137
1138
dbf053f0 1139 }
1140
1141
1142 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1143 PostData(1, fHistList);
1144}
1145
1146void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1147{
1148// Terminate analysis
1149//
1150 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1151}
1152
1153
1154Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1155
1156 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1157
1158 Int_t iCount = 0;
1159 if(type==kTrackAOD){
1160 AliAODEvent *aod = 0;
1161 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1162 else aod = AODEvent();
1163 if(!aod){
1164 return iCount;
1165 }
1166 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1167 AliAODTrack *tr = aod->GetTrack(it);
1168 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1169 if(TMath::Abs(tr->Eta())>0.9)continue;
1170 // if(tr->Pt()<0.3)continue;
1171 list->Add(tr);
1172 iCount++;
1173 }
1174 }
1175 else if (type == kTrackKineAll||type == kTrackKineCharged){
1176 AliMCEvent* mcEvent = MCEvent();
1177 if(!mcEvent)return iCount;
1178 // we want to have alivpartilces so use get track
1179 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1180 if(!mcEvent->IsPhysicalPrimary(it))continue;
1181 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1182 if(type == kTrackKineAll){
1183 list->Add(part);
1184 iCount++;
1185 }
1186 else if(type == kTrackKineCharged){
1187 if(part->Particle()->GetPDG()->Charge()==0)continue;
1188 list->Add(part);
1189 iCount++;
1190 }
1191 }
1192 }
1193 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1194 AliAODEvent *aod = 0;
959508f4 1195 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
dbf053f0 1196 else aod = AODEvent();
1197 if(!aod)return iCount;
1198 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1199 if(!tca)return iCount;
1200 for(int it = 0;it < tca->GetEntriesFast();++it){
225f0094 1201 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
dbf053f0 1202 if(!part->IsPhysicalPrimary())continue;
1203 if(type == kTrackAODMCAll){
1204 list->Add(part);
1205 iCount++;
1206 }
1207 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1208 if(part->Charge()==0)continue;
1209 if(kTrackAODMCCharged){
1210 list->Add(part);
1211 }
1212 else {
1213 if(TMath::Abs(part->Eta())>0.9)continue;
1214 list->Add(part);
1215 }
1216 iCount++;
1217 }
1218 }
1219 }// AODMCparticle
1220 list->Sort();
1221 return iCount;
1222
1223}
1224
1225/*
1226Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1227 for(int i = 0; i < particles.GetEntries(); i++){
1228 AliVParticle *vp = (AliVParticle*)particles.At(i);
1229 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1230 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1231 jInp.set_user_index(i);
1232 inputParticles.push_back(jInp);
1233 }
1234
1235 return 0;
1236
1237}
1238*/