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