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