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