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