]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliJetAnalysis.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- configure adapted to the new directory structure of the HOMER module in PubSub
[u/mrichter/AliRoot.git] / JETAN / AliJetAnalysis.cxx
CommitLineData
62567a80 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//---------------------------------------------------------------------
17// JetAnalysis class
83a444b1 18// Analyse Jets (already found jets)
19// Authors: andreas.morsch@cern.ch, jgcn@mail.cern.ch
20// mercedes.lopez.noriega@cern.ch
62567a80 21//---------------------------------------------------------------------
22
83a444b1 23#include <Riostream.h>
62567a80 24#include "AliJetAnalysis.h"
25ClassImp(AliJetAnalysis)
83a444b1 26
27// root
28#include <Riostream.h>
62567a80 29#include <TH1F.h>
83a444b1 30#include <TH2F.h>
62567a80 31#include <TProfile.h>
62567a80 32#include <TFile.h>
33#include <TTree.h>
34#include <TStyle.h>
35#include <TSystem.h>
36#include <TLorentzVector.h>
83a444b1 37#include <TMath.h>
38// aliroot
62567a80 39#include "AliJetProductionDataPDC2004.h"
40#include "AliJet.h"
83a444b1 41#include "AliJetKineReaderHeader.h"
62567a80 42#include "AliJetESDReaderHeader.h"
43#include "AliUA1JetHeader.h"
44#include "AliLeading.h"
83a444b1 45
46
1b7d5d7e 47AliJetAnalysis::AliJetAnalysis():
48 fReaderHeader(0x0),
49 fDirectory(0x0),
50 fBkgdDirectory(0x0),
51 fFile(0x0),
52 fEventMin(0),
53 fEventMax(-1),
54 fRunMin(0),
55 fRunMax(11),
56 fminMult(0),
57 fPercentage(-1.0),
58 fPartPtCut(0.0),
59 fdrJt(1.0),
60 fdrdNdxi(0.7),
61 fdrdEdr(1.0),
62 fEfactor(1.0),
63 fp0(0.0),
64 fPtJ(0.0),
65 fEJ(0.0),
66 fEtaJ(0.0),
67 fPhiJ(0.0),
68 fjv3X(0.0),
69 fjv3Y(0.0),
70 fjv3Z(0.0),
71 fPythia(kFALSE),
72 fDoPart(kTRUE),
73 fDoGenJ(kTRUE),
74 fDoRecJ(kTRUE),
75 fDoKine(kTRUE),
76 fDoCorr(kTRUE),
77 fDoCorr50(kFALSE),
78 fDoShap(kTRUE),
79 fDoFrag(kTRUE),
80 fDoTrig(kTRUE),
81 fDoJt(kTRUE),
82 fDodNdxi(kTRUE),
83 fDoBkgd(kTRUE),
84 fWeight(1.0),
85 fWShapR(0.0),
86 fWFragR(0.0),
87 fWdEdr(0.0),
88 fWJt(0.0),
89 fWdNdxi(0.0),
90 fPart(0),
91 fGenJ(0),
92 fRecJ(0),
93 fRecB(0),
94 fRKineEneH(0),
95 fRKinePtH(0),
96 fRKinePhiH(0),
97 fRKineEtaH(0),
98 fGKineEneH(0),
99 fGKinePtH(0),
100 fGKinePhiH(0),
101 fGKineEtaH(0),
102 fPKineEneH(0),
103 fPKinePtH(0),
104 fPKinePhiH(0),
105 fPKineEtaH(0),
106 fPGCorrEneH(0),
107 fPGCorrPtH(0),
108 fPGCorrEtaH(0),
109 fPGCorrPhiH(0),
110 fPRCorrEneH(0),
111 fPRCorrPtH(0),
112 fPRCorrEtaH(0),
113 fPRCorrPhiH(0),
114 fRGCorrEneH(0),
115 fRGCorrPtH(0),
116 fRGCorrEtaH(0),
117 fRGCorrPhiH(0),
118 fPRCorr50EneH(0),
119 fPRCorr50PtH(0),
120 fPRCorr50EtaH(0),
121 fPRCorr50PhiH(0),
122 fRGCorr50EneH(0),
123 fRGCorr50PtH(0),
124 fRGCorr50EtaH(0),
125 fRGCorr50PhiH(0),
126 fRFragSelH(0),
127 fRFragRejH(0),
128 fRFragAllH(0),
129 fRShapSelH(0),
130 fRShapRejH(0),
131 fRShapAllH(0),
132 fGTriggerEneH(0),
133 fRTriggerEneH(0),
134 fGPTriggerEneH(0),
135 fPTriggerEneH(0),
136 fdEdrH(0),
137 fdEdrB(0),
138 fPtEneH2(0),
139 fdEdrW(0),
140 fJtH(0),
141 fJtB(0),
142 fJtW(0),
143 fdNdxiH(0),
144 fdNdxiB(0),
145 fdNdxiW(0),
146 fPtEneH(0)
83a444b1 147{
1b7d5d7e 148 // Default constructor
149 fFile = "anaJets.root";
83a444b1 150
83a444b1 151 // initialize weight for dE/dr histo
152 SetdEdrWeight();
153}
154
155AliJetAnalysis::~AliJetAnalysis()
156{
157 // Destructor
158}
159
160
161////////////////////////////////////////////////////////////////////////
162// define histogrames
163
164void AliJetAnalysis::DefineHistograms()
165{
166 // Define the histograms to be filled
167 if (fDoKine) DefineKineH();
168 if (fDoCorr) DefineCorrH();
169 if (fDoCorr50) DefineCorr50H();
170 if (fDoShap) DefineShapH();
171 if (fDoFrag) DefineFragH();
172 if (fDoTrig) DefineTrigH();
173 if (fDoJt) DefineJtH();
174 if (fDodNdxi) DefinedNdxiH();
175}
176
177void AliJetAnalysis::DefineKineH()
178{
179 // leading particle
180 if (fDoPart) {
181 fPKineEneH = new TH1F("PKineEne","Energy of leading particle",50,0,200);
182 SetProperties(fPKineEneH,"Energy (GeV)","Entries");
183 fPKinePtH = new TH1F("PKinePt","Pt of leading particle",50,0,200);
184 SetProperties(fPKinePtH,"P_{T} (GeV)","Entries");
185 fPKinePhiH = new TH1F("PKinePhiH","Azimuthal angle of leading particle",
186 90,0.,2.0*TMath::Pi());
187 SetProperties(fPKinePhiH,"#phi","Entries");
188 fPKineEtaH = new TH1F("PKineEtaH","Pseudorapidity of leading particle",
189 40,-1.0,1.0);
190 SetProperties(fPKineEtaH,"#eta","Entries");
191 }
192 // leading generated jet
193 if (fDoGenJ) {
194 fGKineEneH = new TH1F("GKineEne","Energy of generated jet",50,0,200);
195 SetProperties(fGKineEneH,"Energy (GeV)","Entries");
196 fGKinePtH = new TH1F("GKinePt","Pt of generated jet",50,0,200);
197 SetProperties(fGKinePtH,"P_{T} (GeV)","Entries");
198 fGKinePhiH = new TH1F("GKinePhiH","Azimuthal angle of generated jet",
199 90,0.,2.0*TMath::Pi());
200 SetProperties(fGKinePhiH,"#phi","Entries");
201 fGKineEtaH = new TH1F("GKineEtaH","Pseudorapidity of generated jet",
202 40,-1.0,1.0);
203 SetProperties(fGKineEtaH,"#eta","Entries");
204 }
205 // leading reconstructed jet
206 if (fDoRecJ) {
207 fRKineEneH = new TH1F("RKineEne","Energy of reconstructed jet",50,0,200);
208 SetProperties(fRKineEneH,"Energy (GeV)","Entries");
209 fRKinePtH = new TH1F("RKinePt","Pt of reconstructed jet",50,0,200);
210 SetProperties(fRKinePtH,"P_{T} (GeV)","Entries");
211 fRKinePhiH = new TH1F("RKinePhiH","Azimuthal angle of reconstructed jet",
212 90,0.,2.0*TMath::Pi());
213 SetProperties(fRKinePhiH,"#phi","Entries");
214 fRKineEtaH = new TH1F("RKineEtaH","Pseudorapidity of reconstructed jet",
215 40,-1.0,1.0);
216 SetProperties(fRKineEtaH,"#eta","Entries");
217 }
218}
62567a80 219
83a444b1 220void AliJetAnalysis::DefineCorrH()
62567a80 221{
83a444b1 222 // correlation
223 if (fDoPart && fDoGenJ) {
224 fPGCorrEneH = new TH2F("PGCorrEne","Energy correlation part-gen jet",
225 40,0,200,40,0,200);
226 SetProperties(fPGCorrEneH,"Part Energy (GeV)","Gen Jet Energy (GeV)");
227 fPGCorrPtH = new TH2F("PGCorrPt","Pt correlation part-gen jet",
228 40,0,200,40,0,200);
229 SetProperties(fPGCorrPtH,"Part P_{T} (GeV)","Gen Jet P_{T} (GeV)");
230 fPGCorrEtaH = new TH2F("PGCorrEta","Pseudorapidity correlation part-gen jet",
231 40,-1.0,1.0,40,-1.0,1.0);
232 SetProperties(fPGCorrEtaH,"Part #eta","Gen Jet #eta");
233 fPGCorrPhiH = new TH2F("PGCorrPhi","Azimuthal angle correlation part-gen jet",
234 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
235 SetProperties(fPGCorrPhiH,"Part #phi","Gen Jet #phi");
236 }
237 if (fDoPart && fDoRecJ) {
238 fPRCorrEneH = new TH2F("PRCorrEne","Energy correlation part-rec jet",
239 40,0,200,40,0,200);
240 SetProperties(fPRCorrEneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
241 fPRCorrPtH = new TH2F("PRCorrPt","Pt correlation part-rec jet",
242 40,0,200,40,0,200);
243 SetProperties(fPRCorrPtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
244 fPRCorrEtaH = new TH2F("PRCorrEta","Pseudorapidity correlation part-rec jet",
245 40,-1.0,1.0,40,-1.0,1.0);
246 SetProperties(fPRCorrEtaH,"part #eta","Rec Jet #eta");
247 fPRCorrPhiH = new TH2F("PRCorrPhi","Azimuthal angle correlation part-rec jet",
248 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
249 SetProperties(fPRCorrPhiH,"Part #phi","Rec Jet #phi");
250 }
251 if (fDoGenJ && fDoRecJ) {
252 fRGCorrEneH = new TH2F("RGCorrEne","Energy correlation rec jet-gen jet",
253 40,0,200,40,0,200);
254 SetProperties(fRGCorrEneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
255 fRGCorrPtH = new TH2F("RGCorrPt","Pt correlation rec jet-gen jet",
256 40,0,200,40,0,200);
257 SetProperties(fRGCorrPtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
258 fRGCorrEtaH = new TH2F("RGCorrEta","Pseudorapidity correlation rec jet-gen jet",
259 40,-1.0,1.0,40,-1.0,1.0);
260 SetProperties(fRGCorrEtaH,"Rec Jet #eta","Gen Jet #eta");
261 fRGCorrPhiH = new TH2F("RGCorrPhi","Azimuthal angle correlation rec jet-gen jet",
262 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
263 SetProperties(fRGCorrPhiH,"Rec Jet #phi","Gen Jet #phi");
264 }
62567a80 265}
266
83a444b1 267void AliJetAnalysis::DefineCorr50H()
62567a80 268{
83a444b1 269 // correlation
270 if (fDoPart && fDoRecJ) {
271 fPRCorr50EneH = new TH2F("PRCorr50Ene","Energy correlation part-rec jet",
272 40,0,200,40,0,200);
273 SetProperties(fPRCorr50EneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
274 fPRCorr50PtH = new TH2F("PRCorr50Pt","Pt correlation part-rec jet",
275 40,0,200,40,0,200);
276 SetProperties(fPRCorr50PtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
277 fPRCorr50EtaH = new TH2F("PRCorr50Eta","Pseudorapidity correlation part-rec jet",
278 40,-1.0,1.0,40,-1.0,1.0);
279 SetProperties(fPRCorr50EtaH,"part #eta","Rec Jet #eta");
280 fPRCorr50PhiH = new TH2F("PRCorr50Phi","Azimuthal angle correlation part-rec jet",
281 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
282 SetProperties(fPRCorr50PhiH,"Part #phi","Rec Jet #phi");
283 }
284
285 if (fDoGenJ && fDoRecJ) {
286 fRGCorr50EneH = new TH2F("RGCorr50Ene","Energy correlation rec jet-gen jet",
287 40,0,200,40,0,200);
288 SetProperties(fRGCorr50EneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
289 fRGCorr50PtH = new TH2F("RGCorr50Pt","Pt correlation rec jet-gen jet",
290 40,0,200,40,0,200);
291 SetProperties(fRGCorr50PtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
292 fRGCorr50EtaH = new TH2F("RGCorr50Eta","Pseudorapidity correlation rec jet-gen jet",
293 40,-1.0,1.0,40,-1.0,1.0);
294 SetProperties(fRGCorr50EtaH,"Rec Jet #eta","Gen Jet #eta");
295 fRGCorr50PhiH = new TH2F("RGCorr50Phi","Azimuthal angle correlation rec jet-gen jet",
296 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
297 SetProperties(fRGCorr50PhiH,"Rec Jet #phi","Gen Jet #phi");
298 }
299}
300
301void AliJetAnalysis::DefineShapH()
302{
303 // leading reconstructed jet
304 if (fDoRecJ) {
305 fdEdrH = new TH2F("fdEdrH","dE/dr histo",20,0,1,40,0,200);
306 SetProperties(fdEdrH,"r","Rec Jet P_{T}");
307 fdEdrB = new TH2F("fdEdrB","dE/dr bkgdhisto",20,0,1,40,0,200);
308 SetProperties(fdEdrB,"r","Rec P_{T}");
309 fPtEneH2 = new TH2F("fPtEneH2","fPtEneH2",40,0,200,40,0,200);
310 SetProperties(fPtEneH2,"Rec Jet E","Rec Jet P_{T}");
311 fdEdrW = new TH1F("fdEdrW","weights for dE/dr",40,0,200);
312 SetProperties(fdEdrW,"Rec Jet P_{T}","weights");
62567a80 313
83a444b1 314 fRShapSelH = new TH1F("RShapSel","Shape of generated jets (sel part)",20,0.,1.);
315 SetProperties(fRShapSelH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
316 fRShapRejH = new TH1F("RShapRej","Shape of generated jets (rej part)",20,0.,1.);
317 SetProperties(fRShapRejH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
318 fRShapAllH = new TH1F("RShapAll","Shape of generated jets (all part)",20,0.,1.);
319 SetProperties(fRShapAllH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
320 }
321}
62567a80 322
83a444b1 323void AliJetAnalysis::DefineJtH()
324{
325 // Define the histogram for J_T
326 if (fDoRecJ) {
327 fJtH = new TH2F("fJtH","J_{T} histogram",80,0.,4.,40,0.,200.);
328 SetProperties(fJtH,"J_{T}","Rec Jet P_{T}");
329 fJtB = new TH2F("fJtB","J_{T} bkgd histogram",80,0.,4.,40,0.,200.);
330 SetProperties(fJtB,"J_{T}","Rec P_{T}");
331 fJtW = new TH1F("fJtW","J_{T} weight",40,0,200);
332 SetProperties(fJtW,"J_{T}W","weight");
333 }
334}
62567a80 335
83a444b1 336void AliJetAnalysis::DefinedNdxiH()
337{
338 // Define the histogram for dN/dxi
339 if (fDoRecJ) {
340 fdNdxiH = new TH2F("fdNdxiH","dN/d#xi histo",200,0,10,40,0,200);
341 SetProperties(fdNdxiH,"xi","Rec Jet P_{T}");
342 fdNdxiB = new TH2F("fdNdxiB","dN/d#xi bkgd histo",200,0,10,40,0,200);
343 SetProperties(fdNdxiB,"xi","Rec P_{T}");
344 fdNdxiW = new TH1F("fdNdxiW","dN/d#xi histo",40,0,200);
345 SetProperties(fdNdxiW,"Rec Jet P_{T}","weights");
346 fPtEneH = new TH2F("fPtEneH","fPtEneH",40,0,200,40,0,200);
347 SetProperties(fPtEneH,"Rec Jet E","Rec Jet P_{T}");
348 }
349}
62567a80 350
83a444b1 351void AliJetAnalysis::DefineFragH()
352{
353 // leading reconstructed jet
354 if (fDoRecJ) {
355 fRFragSelH = new TH1F("RFragSel","Frag Fun of reconstructed jets (sel part)",20,0.,10.);
356 SetProperties(fRFragSelH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
357 fRFragRejH = new TH1F("RFragRej","Frag Fun of reconstructed jets (rej part)",20,0.,10.);
358 SetProperties(fRFragRejH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
359 fRFragAllH = new TH1F("RFragAll","Frag Fun of reconstructed jets (all part)",20,0.,10.);
360 SetProperties(fRFragAllH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
361 }
362}
62567a80 363
83a444b1 364void AliJetAnalysis::DefineTrigH()
365{
366 // generated energy
367 fGTriggerEneH = new TProfile("GTriggerEne","Generated energy (trigger bias)",
368 20, 0., 200., 0., 1., "S");
369 fGTriggerEneH->SetXTitle("E_{gen}");
370 fGTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
371 // reconstructed energy
372 fRTriggerEneH = new TProfile("RTriggerEne","Reconstructed energy (trigger bias)",
373 20, 0., 200., 0., 1., "S");
374 fRTriggerEneH->SetXTitle("E_{rec}");
375 fRTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
376 // generated energy vs generated/leading
377 fGPTriggerEneH = new TProfile("GPTriggerEne","Generated energy (trigger bias)",
378 20, 0., 200., 0., 1., "S");
379 fGPTriggerEneH->SetXTitle("E_{gen}");
380 fGPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
381 // leading particle energy
382 fPTriggerEneH = new TProfile("PTriggerEne","Leading particle energy (trigger bias)",
383 20, 0., 200., 0., 1., "S");
384 fPTriggerEneH->SetXTitle("E_{L}/0.2");
385 fPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
62567a80 386
83a444b1 387}
62567a80 388
83a444b1 389void AliJetAnalysis::SetProperties(TH1* h,const char* x, const char* y) const
390{
391 //Set properties of histos (x and y title and error propagation)
392 h->SetXTitle(x);
393 h->SetYTitle(y);
394 h->Sumw2();
395}
396
397////////////////////////////////////////////////////////////////////////
398// compute weights for dE/dr histogram
62567a80 399
83a444b1 400void AliJetAnalysis::SetdEdrWeight()
401{
402 // Due to the limited acceptance, each bin in the dE/dr has a different
403 // weight given by the ratio of the area of a disk dR to the area
404 // within the eta acceptance. The weight depends on the eta position
405 // of the jet. Here a look up table for the weights is computed. It
406 // assumes |etaJet|<0.5 and |eta_lego|<0.9. It makes bins of 0.05
407 // units in eta. Note that this is fixed by the bin width chosen for
408 // the histogram --> this weights are tailored for the specific
409 // histogram definition used here!
410
411 // two dimensional table: first index, bin in eta of jet, second
412 // index bin of dE/dr histo
413
414 Int_t nBins = 20;
415 Float_t xMin = 0.0;
416 Float_t xMax = 1.0;
417 Float_t binSize = (xMax-xMin)/nBins;
62567a80 418
83a444b1 419 Float_t da,ds,r1,r2,h1,h2,a1,a2,theta1,theta2,area1,area2;
420 Int_t ji;
421 for (Int_t i=0;i<(nBins/2);i++) {
422 // index of first histo bin needing a scale factor
423 ji=(nBins-2)-i;
424 for(Int_t j=0;j<nBins;j++) {
425 // area of ring.
426 da = TMath::Pi()*(binSize*binSize)*(1.0+2.0*j);
427 ds = 1.0;
428 if (j>=ji) { // compute missing area using segments of circle
429 r1=j*binSize;
430 r2=(j+1)*binSize;
431 h1=(j-ji)*binSize;
432 h2=(j+1-ji)*binSize;
433 a1=2.0*TMath::Sqrt(2.0*h1*r1-h1*h1);
434 a2=2.0*TMath::Sqrt(2.0*h2*r2-h2*h2);
435 theta1=2*TMath::ACos((r1-h1)/r1);
436 theta2=2*TMath::ACos((r2-h2)/r2);
437 area1=binSize*(r1*r1*theta1-a1*(r1-h1));
438 area2=binSize*(r2*r2*theta2-a2*(r2-h2));
439 ds = (da-(area2-area1))/da;
440 }
441 fWeightdEdr[i][j]=ds/da;
442 }
443 }
444}
445
446// get weight for dE/dr histogram
447Float_t AliJetAnalysis::GetdEdrWeight(Float_t etaJet, Float_t r)
448{
449 // Return the correponding weight for the dE/dr plot
450 Int_t nBins = 20;
451 Float_t xMin = 0.0;
452 Float_t xMax = 1.0;
453 Float_t binSize = (xMax-xMin)/nBins;
454
455 Float_t eta = TMath::Abs(etaJet);
456 if ((eta > 0.5) || (r > fdrdEdr)) return 0.0;
457 Int_t binJet = (Int_t) TMath::Floor(eta/binSize);
458 Int_t binR = (Int_t) TMath::Floor(r/binSize);
459 Float_t w = fWeightdEdr[binJet][binR];
460 return w;
461}
462
463
464////////////////////////////////////////////////////////////////////////
465// fill histograms
466
467void AliJetAnalysis::FillHistograms()
468{
469 // fill histograms
62567a80 470
471 // Run data
83a444b1 472 AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004();
473
474 // Loop over runs
475 TFile* jFile = 0x0;
476
477 for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++) {
478 // Open file
479 char fn[256];
480 sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data());
481 jFile = new TFile(fn);
482 printf(" Analyzing run: %d %s\n", iRun,fn);
483
484 // Get reader header and events to be looped over
485 AliJetReaderHeader *jReaderH = (AliJetReaderHeader*)(jFile->Get(fReaderHeader));
486 if (fEventMin == -1) fEventMin = jReaderH->GetFirstEvent();
487 if (fEventMax == -1) {
488 fEventMax = jReaderH->GetLastEvent();
489 } else {
490 fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent());
491 }
492
493 AliUA1JetHeader *jH = (AliUA1JetHeader *) (jFile->Get("AliUA1JetHeader"));
62567a80 494
83a444b1 495 // Calculate weight
496 fWeight = runData->GetWeight(iRun)/ ( (Float_t) (fEventMax - fEventMin + 1));
62567a80 497
83a444b1 498 // Loop over events
499 for (Int_t i = fEventMin; i < fEventMax; i++) {
500 if (i%100 == 0) printf(" Analyzing event %d / %d \n",i,fEventMax);
501
502 // Get next tree with AliJet
503 char nameT[100];
504 sprintf(nameT, "TreeJ%d",i);
505 TTree *jetT =(TTree *)(jFile->Get(nameT));
506 if (fDoRecJ) jetT->SetBranchAddress("FoundJet", &fRecJ);
507 if (fDoGenJ) jetT->SetBranchAddress("GenJet", &fGenJ);
508 if (fDoPart) jetT->SetBranchAddress("LeadingPart", &fPart);
509 jetT->GetEntry(0);
510
511 int nJets = fRecJ->GetNJets();
512
513 TArrayI inJet = fRecJ->GetInJet();
514 if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
515 if (fDoKine) FillKineH();
516 if (fDoCorr) FillCorrH();
517 if (fDoCorr50) FillCorr50H();
518 if (fDoShap) FillShapH(jH->GetRadius());
519 if (fDoFrag) FillFragH();
520 if (fDoTrig) FillTrigH();
521 if (fDoJt) FillJtH();
522 if (fDodNdxi) FilldNdxiH();
523 }
524 delete jetT; // jet should be deleted before creating a new one
525 if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
526 if (fDoBkgd && nJets>0) FillBkgd(i,iRun);
527 }
528 } // end loop over events in one file
529 if (jFile) jFile->Close();
530 delete jFile;
531 } // end loop over files
532}
533
534void AliJetAnalysis::FillKineH()
535{
536 // leading particle
537 if (fDoPart && fPart->LeadingFound()){
538 fPKineEneH->Fill(fPart->GetE(),fWeight);
539 fPKinePtH->Fill(fPart->GetPt(),fWeight);
540 fPKinePhiH->Fill(fPart->GetPhi(),fWeight);
541 fPKineEtaH->Fill(fPart->GetEta(),fWeight);
542 }
543 // leading generated jet
544 if (fDoGenJ && fGenJ->GetNJets()> 0){
545 fGKineEneH->Fill(fGenJ->GetE(0),fWeight);
546 fGKinePtH->Fill(fGenJ->GetPt(0),fWeight);
547 fGKinePhiH->Fill(fGenJ->GetPhi(0),fWeight);
548 fGKineEtaH->Fill(fGenJ->GetEta(0),fWeight);
549 }
550 // leading reconstructed jet
551 if (fDoRecJ && fRecJ->GetNJets()> 0) {
552 TArrayF p=fRecJ->GetPtFromSignal();
553 if (p[0]>fPercentage) {
554 fRKineEneH->Fill(fRecJ->GetE(0)/fEfactor,fWeight);
555 fRKinePtH->Fill(fRecJ->GetPt(0),fWeight);
556 fRKinePhiH->Fill(fRecJ->GetPhi(0),fWeight);
557 fRKineEtaH->Fill(fRecJ->GetEta(0),fWeight);
558 }
559 }
560}
561
562void AliJetAnalysis::FillCorrH()
563{
564 // Fill correlation histograms
565 if (fDoPart && fPart->LeadingFound() && fDoGenJ && fGenJ->GetNJets()> 0)
566 Correlation(fPart->GetLeading(),fGenJ->GetJet(0),
567 fPGCorrEneH,fPGCorrPtH,fPGCorrEtaH,fPGCorrPhiH);
568 if (fDoPart && fPart->LeadingFound() && fDoRecJ && fRecJ->GetNJets()> 0) {
569 TArrayF p=fRecJ->GetPtFromSignal();
570 if (p[0]>fPercentage)
571 Correlation(fPart->GetLeading(),fRecJ->GetJet(0),
572 fPRCorrEneH,fPRCorrPtH,fPRCorrEtaH,fPRCorrPhiH);
573 }
574 if (fDoGenJ && fGenJ->GetNJets()> 0 && fDoRecJ && fRecJ->GetNJets()> 0) {
575 TArrayF p=fRecJ->GetPtFromSignal();
576 if (p[0]>fPercentage)
577 Correlation(fRecJ->GetJet(0), fGenJ->GetJet(0),
578 fRGCorrEneH,fRGCorrPtH,fRGCorrEtaH,fRGCorrPhiH);
579 }
580}
581
582void AliJetAnalysis::Correlation(TLorentzVector *lv1,TLorentzVector *lv2,
583 TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
584{
585 // Correlation histograms
586 h1->Fill(lv1->E(),lv2->E(),fWeight);
587 h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
588 h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
589 Float_t p1, p2;
590 p1 = ((lv1->Phi() < 0) ? (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
591 p2 = ((lv2->Phi() < 0) ? (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
592 h4->Fill(p1,p2,fWeight);
593}
594
595void AliJetAnalysis::FillCorr50H()
596{
597 // Fill correlation histograms when one particle has > 50% of jet energy
598 if (fDoRecJ && fRecJ->GetNJets()> 0) {
599 TArrayF p = fRecJ->GetPtFromSignal();
600 if (p[0]>fPercentage) {
601 if (fDoPart && fPart->LeadingFound())
602 Correlation50(fRecJ, fPart->GetLeading(),fRecJ->GetJet(0),
603 fPRCorr50EneH,fPRCorr50PtH,fPRCorr50EtaH,fPRCorr50PhiH);
604 if (fDoGenJ && fGenJ->GetNJets()> 0)
605 Correlation50(fRecJ, fRecJ->GetJet(0), fGenJ->GetJet(0),
606 fRGCorr50EneH,fRGCorr50PtH,fRGCorr50EtaH,fRGCorr50PhiH);
607 }
608 }
609}
610
611void AliJetAnalysis::Correlation50(AliJet *j,TLorentzVector *lv1,TLorentzVector *lv2,
612 TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
613{
614 // Correlation histograms when one particle has > 50% of jet energy
615 TArrayF ptin = j->GetPtIn();
616 TArrayF etain = j->GetEtaIn();
617 TArrayI inJet = j->GetInJet();
618
619 Int_t flag = 0;
620 for(Int_t i=0;i<(inJet.GetSize());i++) {
621 if (inJet[i] == 1) {
622 Float_t t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
623 Float_t x = (ptin[i]*TMath::Sqrt(1.+1./(t1*t1)))/(j->GetE(0));
624 if (x>0.5) flag++;
625 }
626 }
627 if (flag>1) cout << " Flag = " << flag << endl;
628 if (flag == 1) {
629 h1->Fill(lv1->E(),lv2->E(),fWeight);
630 h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
631 h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
632 Float_t p1, p2;
633 p1 = ((lv1->Phi() < 0) ?
634 (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
635 p2 = ((lv2->Phi() < 0) ?
636 (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
637 h4->Fill(p1,p2,fWeight);
638 }
639}
640
641void AliJetAnalysis::FillJtH()
642{
643 // Fill J_T histogram
644 if (fRecJ->GetNJets()> 0) {
645 fjv3X = 0.0; fjv3Y = 0.0; fjv3Z = 0.0;
646 TArrayF p=fRecJ->GetPtFromSignal();
647 if (p[0]>fPercentage) {
648 // initialize
649 const TVector3 kJv3 = fRecJ->GetJet(0)->Vect();
650 fjv3X = kJv3.X(); fjv3Y = kJv3.Y(); fjv3Z = kJv3.Z();
651 TVector3 trk;
652 TArrayI inJet = fRecJ->GetInJet();
653 TArrayF etain = fRecJ->GetEtaIn();
654 TArrayF ptin = fRecJ->GetPtIn();
655 TArrayF phiin = fRecJ->GetPhiIn();
656 Float_t deta, dphi,jt, dr;
657 for(Int_t i=0;i<inJet.GetSize();i++) {
658 deta = etain[i] - fRecJ->GetEta(0);
659 dphi = phiin[i] - fRecJ->GetPhi(0);
660 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
661 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
662 dr = TMath::Sqrt(deta * deta + dphi * dphi);
663 if ((dr<fdrJt) && (ptin[i] > fPartPtCut)) {
664 trk.SetPtEtaPhi(ptin[i],etain[i],phiin[i]);
665 jt = TMath::Sin(trk.Angle(kJv3))*trk.Mag();
666 fJtH->Fill(jt,fRecJ->GetPt(0),fWeight);
62567a80 667 }
83a444b1 668 }
669 fJtW->Fill(fRecJ->GetPt(0),fWeight);
670 fWJt+=fWeight;
671 }
672 }
673}
674
675void AliJetAnalysis::FilldNdxiH()
676{
677 // Fill dN/d#xi histograms
678 if (fRecJ->GetNJets()> 0) {
679 TArrayF p=fRecJ->GetPtFromSignal();
680 if (p[0]>fPercentage) {
681 fp0 = p[0];
682 TArrayI inJet = fRecJ->GetInJet();
683 TArrayF etain = fRecJ->GetEtaIn();
684 TArrayF ptin = fRecJ->GetPtIn();
685 TArrayF phiin = fRecJ->GetPhiIn();
686 Float_t xi,t1,ene,dr,deta,dphi;
687 for(Int_t i=0;i<inJet.GetSize();i++) {
688 deta = etain[i] - fRecJ->GetEta(0);
689 dphi = phiin[i] - fRecJ->GetPhi(0);
690 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
691 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
692 dr = TMath::Sqrt(deta * deta + dphi * dphi);
693 if ((dr<fdrdNdxi) && (ptin[i] > fPartPtCut)) {
694 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
695 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
696 xi = (Float_t) TMath::Log((fRecJ->GetE(0)/fEfactor)/ene);
697 fdNdxiH->Fill(xi,fRecJ->GetPt(0),fWeight);
698 }
699 }
700 fdNdxiW->Fill(fRecJ->GetPt(0),fWeight);
701 fPtEneH->Fill(fRecJ->GetE(0)/fEfactor,fRecJ->GetPt(0),fWeight);
702 fWdNdxi+=fWeight;
703 }
704 }
705}
706
707void AliJetAnalysis::FillBkgd(Int_t eventN, Int_t runN)
708{
709 // Background calculated with hijing events (no pythia)
710 if (fp0>fPercentage) {
711 fPtJ=0.,fEJ=0.,fEtaJ=0.,fPhiJ=0.;
712 // Background calculated with hijing events (no pythia)
713 AliJetProductionDataPDC2004* runDataB = new AliJetProductionDataPDC2004();
714 TFile* jFileB =0x0;;
715 char fnB[256];
62567a80 716
83a444b1 717 sprintf(fnB, "%s/%s.root", fBkgdDirectory, (runDataB->GetRunTitle(runN)).Data());
718 jFileB = new TFile(fnB);
719
720 char nameB[100];
721 sprintf(nameB, "TreeJ%d",eventN);
722 TTree *jetB =(TTree *)(jFileB->Get(nameB));
723 jetB->SetBranchAddress("FoundJet", &fRecB);
724 jetB->GetEntry(0);
62567a80 725
83a444b1 726 TArrayI inJetB = fRecB->GetInJet();
727 TArrayF etainB = fRecB->GetEtaIn();
728 TArrayF ptinB = fRecB->GetPtIn();
729 TArrayF phiinB = fRecB->GetPhiIn();
730 fPtJ = fRecJ->GetPt(0);
731 fEJ = fRecJ->GetE(0);
732 fEtaJ = fRecJ->GetEta(0);
733 fPhiJ = fRecJ->GetPhi(0);
734 Float_t t1,ene,xi,detaB,dphiB,drB,jt,wB;
735 TVector3 trkB;
736 TVector3 jv3B;
737 jv3B.SetX(fjv3X); jv3B.SetY(fjv3Y); jv3B.SetZ(fjv3Z);
738
739 for(Int_t k=0;k<inJetB.GetSize();k++){
740 if(ptinB[k] > fPartPtCut){
741 detaB = etainB[k] - fEtaJ;
742 dphiB = phiinB[k] - fPhiJ;
743 if (dphiB < -TMath::Pi()) dphiB= -dphiB - 2.0 * TMath::Pi();
744 if (dphiB > TMath::Pi()) dphiB = 2.0 * TMath::Pi() - dphiB;
745 drB = TMath::Sqrt(detaB * detaB + dphiB * dphiB);
746 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etainB[k])));
747 ene = ptinB[k]*TMath::Sqrt(1.+1./(t1*t1));
748 trkB.SetPtEtaPhi(ptinB[k],etainB[k],phiinB[k]);
749 // --- dN/dxi
750 if (drB<fdrdNdxi) {
751 xi = (Float_t) TMath::Log(fEJ/ene);
752 fdNdxiB->Fill(xi,fPtJ,fWeight);
753 }
754 // --- Jt
755 if (drB<fdrJt) {
756 jt = TMath::Sin(trkB.Angle(jv3B))*(trkB.Mag());
757 fJtB->Fill(jt,fPtJ,fWeight);
758 }
759 // --- dE/dr
760 if (drB<fdrdEdr) {
761 wB = GetdEdrWeight(fEtaJ,drB)*fWeight*ene;
762 fdEdrB->Fill(drB,fPtJ,wB);
763 }
764 }
765 }
766 delete jetB;
767 if (jFileB) jFileB->Close();
768 delete jFileB;
769 }
770}
62567a80 771
83a444b1 772void AliJetAnalysis::FillShapH(Float_t r)
773{
774 // Fill jet shape histograms
775 if (fDoRecJ && fRecJ->GetNJets()> 0) {
776 TArrayF p=fRecJ->GetPtFromSignal();
777 if (p[0]>fPercentage) {
778 Shape(fRecJ,fRShapSelH,fRShapRejH,fRShapAllH,fdEdrH,fPtEneH2,fdEdrW,r);
779 fWdEdr+=fWeight;
780 fWShapR+=fWeight;
781 }
782 }
783}
62567a80 784
83a444b1 785void AliJetAnalysis::Shape(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha,
786 TH2F *hdedr, TH2F *hptene, TH1F* wdedr, Float_t radius)
787{
788 // initialize
789 TArrayI inJet = j->GetInJet();
790 TArrayF etain = j->GetEtaIn();
791 TArrayF ptin = j->GetPtIn();
792 TArrayF phiin = j->GetPhiIn();
62567a80 793
83a444b1 794 // first compute dE/dr histo
795 Float_t etaj = j->GetEta(0);
796 Float_t ene,w,deta,dphi,dr,t1;
797 for(Int_t i=0;i<inJet.GetSize();i++) {
798 deta = etain[i] - j->GetEta(0);
799 dphi = phiin[i] - j->GetPhi(0);
800 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
801 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
802 dr = TMath::Sqrt(deta * deta + dphi * dphi);
803 if ((dr<fdrdEdr) && (ptin[i] > fPartPtCut)) {
804 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
805 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
806 w = GetdEdrWeight(etaj,dr)*fWeight*ene;
807 hdedr->Fill(dr,j->GetPt(0),w);
808 }
809 }
810 hptene->Fill(fRecJ->GetE(0),fRecJ->GetPt(0),fWeight);
811 wdedr->Fill(j->GetPt(0),fWeight);
62567a80 812
83a444b1 813 // now compute shape histos
814 Int_t nBins = ha->GetNbinsX();
815 Float_t xMin = ha->GetXaxis()->GetXmin();
816 Float_t xMax = ha->GetXaxis()->GetXmax();
817 Float_t binSize,halfBin;
818 binSize = (xMax-xMin)/nBins;
819 halfBin = binSize/2;
820 Float_t rptS[20], rptR[20], rptA[20];
821 for(Int_t i=0;i<nBins;i++) rptS[i]=rptR[i]=rptA[i]=0.0;
822 // fill bins in r for leading jet
823 for(Int_t i=0;i<inJet.GetSize();i++) {
824 if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
825 deta = etain[i] - j->GetEta(0);
826 dphi = phiin[i] - j->GetPhi(0);
827 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
828 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
829 dr = TMath::Sqrt(deta * deta + dphi * dphi);
830 Float_t rR = dr/radius;
831 Int_t bin = (Int_t) TMath::Floor(rR/binSize);
832 rptA[bin]+=ptin[i]/(j->GetPt(0));
833 if (inJet[i] == 1) rptS[bin]+=ptin[i]/(j->GetPt(0));
834 if (fPythia && inJet[i] == -1)
835 rptR[bin]+=ptin[i]/(j->GetPt(0));
836 }
837 }
62567a80 838
83a444b1 839 // compute shape and fill histogram
840 Float_t ptS,ptR,ptA,r;
841 ptS=ptR=ptA=0.0;
842 for (Int_t i=0;i<nBins;i++) {
843 ptS+=rptS[i];
844 ptR+=rptR[i];
845 ptA+=rptA[i];
846 r=(i+1)*binSize-halfBin;
847 hs->Fill(r,ptS*fWeight);
848 if(fPythia) {
849 hr->Fill(r,ptR*fWeight);
850 ha->Fill(r,ptA*fWeight);
851 }
852 }
853}
854
855void AliJetAnalysis::FillFragH()
856{
857 // Fill fragmentation histogram
858 if (fDoRecJ && fRecJ->GetNJets()> 0) {
859 TArrayF p=fRecJ->GetPtFromSignal();
860 if (p[0]>fPercentage) {
861 FragFun(fRecJ,fRFragSelH,fRFragRejH,fRFragAllH);
862 fWFragR+=fWeight;
863 }
864 }
865}
62567a80 866
83a444b1 867void AliJetAnalysis::FragFun(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha)
868{
869 // Calculate de fragmentation function
870 TArrayI inJet = j->GetInJet();
871 TArrayF etain = j->GetEtaIn();
872 TArrayF ptin = j->GetPtIn();
62567a80 873
83a444b1 874 Float_t xi,t1,ene;
62567a80 875
83a444b1 876 for(Int_t i=0;i<(inJet.GetSize());i++) {
877 if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
878 t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
879 ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
880 xi = (Float_t) TMath::Log((j->GetE(0))/ene);
881 if (fPythia) ha->Fill(xi,fWeight);
882 if (inJet[i] == 1) hs->Fill(xi,fWeight);
883 if (fPythia && inJet[i] == -1) hr->Fill(xi,fWeight);
884 }
885 }
886}
62567a80 887
83a444b1 888void AliJetAnalysis:: FillTrigH()
889{
890 // Fill trigger bias histograms
891 if(fDoTrig && fGenJ->GetNJets()>0 && fRecJ->GetNJets()>0){
892 TArrayF p=fRecJ->GetPtFromSignal();
893 if (p[0]>fPercentage) {
894 float genEne = fGenJ->GetE(0);
895 float recEne = fRecJ->GetE(0);
896 float eneRatio = (recEne/fEfactor)/genEne;
897
898 fGTriggerEneH->Fill(genEne, eneRatio, fWeight);
899 fRTriggerEneH->Fill(recEne/fEfactor, eneRatio, fWeight);
900
901 if (fPart->LeadingFound()){
902 float leaEne = fPart->GetE();
903 float eneRatio2 = leaEne/genEne;
904 fGPTriggerEneH->Fill(genEne, eneRatio2, fWeight);
905 fPTriggerEneH->Fill(leaEne/0.2, eneRatio2, fWeight);
906 }
907 }
908 }
909}
910
911
912////////////////////////////////////////////////////////////////////////
913// Normalize histogrames
914
915void AliJetAnalysis::NormHistograms()
916{
917 // normalize shape histograms
918 if (fDoShap) {
919 if (fDoRecJ && fWShapR>0.) { // leading reconstructed jet
920 fRShapSelH->Scale(1.0/fWShapR);
921 fRShapRejH->Scale(1.0/fWShapR);
922 fRShapAllH->Scale(1.0/fWShapR);
923 }
924 }
62567a80 925
83a444b1 926 // normalize fragmentation function histograms
927 if (fDoFrag) {
928 if (fDoRecJ && fWFragR>0.) { // leading reconstructed jet
929 fRFragSelH->Scale(2.0/fWFragR);
930 fRFragRejH->Scale(2.0/fWFragR);
931 fRFragAllH->Scale(2.0/fWFragR);
932 }
933 }
934}
935
936////////////////////////////////////////////////////////////////////////
937
938void AliJetAnalysis::PlotHistograms()
939{
940 // Plot histogramas (to be done...)
941 if (fDoKine) PlotKineH();
942 if (fDoCorr) PlotCorrH();
943 if (fDoCorr50) PlotCorr50H();
944 if (fDoShap) PlotShapH();
945 if (fDoFrag) PlotFragH();
946 if (fDoTrig) PlotTrigH();
947}
948
949void AliJetAnalysis::PlotKineH() const
950{
951 // missing
952 if (fDoPart) ;
953 if (fDoGenJ) ;
954 if (fDoRecJ) ;
955}
956
957void AliJetAnalysis::PlotCorrH() const
958{
959 // missing
960 if (fDoPart && fDoGenJ) ;
961 if (fDoPart && fDoRecJ) ;
962 if (fDoGenJ && fDoRecJ) ;
963}
964void AliJetAnalysis::PlotCorr50H() const
965{
966 // missing
967 if (fDoPart && fDoGenJ) ;
968 if (fDoPart && fDoRecJ) ;
969 if (fDoGenJ && fDoRecJ) ;
970}
971
972void AliJetAnalysis::PlotShapH() const
973{
974 // missing
975 if (fDoGenJ) ;
976 if (fDoRecJ) ;
977}
978
979void AliJetAnalysis::PlotFragH() const
980{
981 // missing
982 if (fDoGenJ) ;
983 if (fDoRecJ) ;
984}
985
986void AliJetAnalysis::PlotTrigH()
987{
988 // missing
989
990}
991
992////////////////////////////////////////////////////////////////////////
993
994void AliJetAnalysis::SaveHistograms()
995{
996 // Save histograms
997 TFile *fOut = new TFile(fFile,"recreate");
998 fOut->cd();
999 if (fDoKine) SaveKineH();
1000 if (fDoCorr) SaveCorrH();
1001 if (fDoCorr50) SaveCorr50H();
1002 if (fDoShap) SaveShapH();
1003 if (fDoFrag) SaveFragH();
1004 if (fDoTrig) SaveTrigH();
1005 if (fDoJt) SaveJtH();
1006 if (fDodNdxi) SavedNdxiH();
1007 fOut->Close();
1008}
1009
1010void AliJetAnalysis::SaveKineH()
1011{
1012 // Save kinematic histograms
1013 if (fDoPart) {
1014 fPKineEneH->Write();
1015 fPKinePtH->Write();
1016 fPKinePhiH->Write();
1017 fPKineEtaH->Write();
1018 }
62567a80 1019
83a444b1 1020 if (fDoGenJ) {
1021 fGKineEneH->Write();
1022 fGKinePtH->Write();
1023 fGKinePhiH->Write();
1024 fGKineEtaH->Write();
1025 }
62567a80 1026
83a444b1 1027 if (fDoRecJ) {
1028 fRKineEneH->Write();
1029 fRKinePtH->Write();
1030 fRKinePhiH->Write();
1031 fRKineEtaH->Write();
1032 }
1033}
1034
1035void AliJetAnalysis::SaveCorrH()
1036{
1037 // Save correlation histograms
1038 if (fDoPart && fDoGenJ) {
1039 fPGCorrEneH->Write();
1040 fPGCorrPtH->Write();
1041 fPGCorrEtaH->Write();
1042 fPGCorrPhiH->Write();
1043 }
62567a80 1044
83a444b1 1045 if (fDoPart && fDoRecJ) {
1046 fPRCorrEneH->Write();
1047 fPRCorrPtH->Write();
1048 fPRCorrEtaH->Write();
1049 fPRCorrPhiH->Write();
1050 }
62567a80 1051
83a444b1 1052 if (fDoGenJ && fDoRecJ) {
1053 fRGCorrEneH->Write();
1054 fRGCorrPtH->Write();
1055 fRGCorrEtaH->Write();
1056 fRGCorrPhiH->Write();
1057 }
1058}
1059
1060void AliJetAnalysis::SaveCorr50H()
1061{
1062 // Save correlation histograms (special case)
1063 if (fDoPart && fDoRecJ) {
1064 fPRCorr50EneH->Write();
1065 fPRCorr50PtH->Write();
1066 fPRCorr50EtaH->Write();
1067 fPRCorr50PhiH->Write();
1068 }
1069 if (fDoGenJ && fDoRecJ) {
1070 fRGCorr50EneH->Write();
1071 fRGCorr50PtH->Write();
1072 fRGCorr50EtaH->Write();
1073 fRGCorr50PhiH->Write();
1074 }
1075}
1076
1077void AliJetAnalysis::SaveShapH()
1078{
1079 // Save jet shape histograms
1080 if (fDoRecJ) {
1081 fRShapSelH->Write();
1082 fdEdrH->Write();
1083 if(fDoBkgd) fdEdrB->Write();
1084 fPtEneH2->Write();
1085 fdEdrW->Write();
1086 if (fPythia){
1087 fRShapRejH->Write();
1088 fRShapAllH->Write();
1089 }
1090 }
1091}
1092
1093void AliJetAnalysis::SaveJtH()
1094{
1095 // Save J_T histograms
1096 if (fDoRecJ) {
1097 fJtH->Write();
1098 if(fDoBkgd) fJtB->Write();
1099 fJtW->Write();
1100 }
1101}
1102
1103void AliJetAnalysis::SavedNdxiH()
1104{
1105 // Save dN/d#xi histograms
1106 if (fDoRecJ) {
1107 fdNdxiH->Write();
1108 if(fDoBkgd) fdNdxiB->Write();
1109 fPtEneH->Write();
1110 fdNdxiW->Write();
1111 }
1112}
1113
1114void AliJetAnalysis::SaveFragH()
1115{
1116 // Save fragmentation histograms
1117 if (fDoRecJ) {
1118 fRFragSelH->Write();
1119 if(fPythia){
1120 fRFragRejH->Write();
1121 fRFragAllH->Write();
1122 }
1123 }
1124}
1125
1126void AliJetAnalysis::SaveTrigH()
1127{
1128 // Save trigger bias histograms
1129 if(fDoTrig){
1130 fGTriggerEneH->Write();
1131 fRTriggerEneH->Write();
1132 fGPTriggerEneH->Write();
1133 fPTriggerEneH->Write();
1134 }
1135}
1136
1137////////////////////////////////////////////////////////////////////////
1138// main Analysis function
1139
1140void AliJetAnalysis::Analyze()
1141
1142{
1143 // Kinematics for
1144 // leading particle
1145 // leading generated jet
1146 // leading reconstructed jet
1147
1148 // Correlations amd resolutions
1149 // a) correlations in energy, pt, phi, eta
1150 // b) resolutions in energy, pt, phi, eta, r
1151 // leading particle and leading generated jet
1152 // leading particle and leading reconstructed jet
1153 // leading generated jet and leading reconstructed jet
1154
1155 // Fragmentation functions and Shapes
1156 // a) integrated over all pt
1157 // b) in 3 flavors:
1158 // b.1) only for user selected particles in jet
1159 // b.2) only for user rejected particles in jet
1160 // b.3) for all particles in jet
1161
1162 DefineHistograms();
1163 FillHistograms();
1164 NormHistograms();
1165 PlotHistograms();
1166 SaveHistograms();
1167}
62567a80 1168