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