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