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