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