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