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