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