PPR version of the JETAN code (M. Lopez Noriega)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Jan 2006 09:15:05 +0000 (09:15 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Jan 2006 09:15:05 +0000 (09:15 +0000)
36 files changed:
JETAN/AliJet.cxx
JETAN/AliJet.h
JETAN/AliJetAnalysis.cxx
JETAN/AliJetAnalysis.h
JETAN/AliJetControlPlots.cxx
JETAN/AliJetControlPlots.h
JETAN/AliJetESDReader.cxx
JETAN/AliJetESDReader.h
JETAN/AliJetESDReaderHeader.cxx
JETAN/AliJetESDReaderHeader.h
JETAN/AliJetFinder.cxx
JETAN/AliJetFinder.h
JETAN/AliJetHeader.cxx
JETAN/AliJetHeader.h
JETAN/AliJetKineReader.cxx
JETAN/AliJetKineReader.h
JETAN/AliJetKineReaderHeader.cxx
JETAN/AliJetKineReaderHeader.h
JETAN/AliJetProductionData.h
JETAN/AliJetReader.cxx
JETAN/AliJetReader.h
JETAN/AliJetReaderHeader.cxx
JETAN/AliJetReaderHeader.h
JETAN/AliLeading.cxx
JETAN/AliLeading.h
JETAN/AliPxconeJetFinder.cxx
JETAN/AliPxconeJetHeader.cxx
JETAN/AliPxconeJetHeader.h
JETAN/AliUA1JetFinder.cxx
JETAN/AliUA1JetHeader.cxx
JETAN/AliUA1JetHeader.h
JETAN/JetAnalysisLinkDef.h
JETAN/UA1Common.h
JETAN/libJETAN.pkg
JETAN/read_jets.C
JETAN/ua1_jet_finder.F

index bf6d177..450f5a0 100644 (file)
 
 #include "AliJet.h"
 ClassImp(AliJet)
-////////////////////////////////////////////////////////////////////////
+
 
 AliJet::AliJet() 
 {
-  //
   // Constructor
-  //
   fJets = new TClonesArray("TLorentzVector",1000);
   fNInput=0;
   fNJets=0;
+  fEtAvg = 0;
   fInJet=TArrayI();
+  fPtIn=TArrayF();
+  fEtaIn=TArrayF();
+  fPhiIn=TArrayF();
   fPtFromSignal=TArrayF();
   fMultiplicities=TArrayI();
+  fNCells=TArrayI();
 } 
 
 ////////////////////////////////////////////////////////////////////////
@@ -57,7 +58,7 @@ AliJet::~AliJet()
 
 Bool_t AliJet::OutOfRange(Int_t i, const char *s) const
 {
-  // checks if i is a valid index. s= name of calling method
+  // checks if i is a valid index. s = name of calling method
   if (i >= fNJets || i < 0) {
     cout << s << " Index " << i << " out of range" << endl;
     return kTRUE;
@@ -71,14 +72,13 @@ TLorentzVector* AliJet::GetJet(Int_t i)
 {
   // returns i-jet
   if (OutOfRange(i, "AliJet::GetJet:")) return 0;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv; 
 }
 
 ////////////////////////////////////////////////////////////////////////
 
-Int_t AliJet::GetMultiplicity(Int_t i)
+Int_t AliJet::GetMultiplicity(Int_t i) const
 {
   // gets multiplicity of i-jet
   if (OutOfRange(i, "AliJet::GetMultiplicity:")) return 0;
@@ -87,11 +87,19 @@ Int_t AliJet::GetMultiplicity(Int_t i)
 
 ////////////////////////////////////////////////////////////////////////
 
+Int_t AliJet::GetNCell(Int_t i) const
+{
+  // gets number of cell of i-jet
+  if (OutOfRange(i, "AliJet::GetNCell:")) return 0;
+  return fNCells[i];
+}
+
+////////////////////////////////////////////////////////////////////////
+
 Double_t AliJet::GetPx(Int_t i)
 {
 // Get Px component of jet i
   if (OutOfRange(i, "AliJet::GetPx:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv->Px();
 }
@@ -102,7 +110,6 @@ Double_t AliJet::GetPy(Int_t i)
 {
 // Get Py component of jet i
   if (OutOfRange(i, "AliJet::GetPy:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv->Py();
 }
@@ -113,7 +120,6 @@ Double_t AliJet::GetPz(Int_t i)
 {
 // Get Pz component of jet i
   if (OutOfRange(i, "AliJet::GetPz:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv->Pz();
 }
@@ -124,7 +130,6 @@ Double_t AliJet::GetP(Int_t i)
 {
 // Get momentum of jet i
   if (OutOfRange(i, "AliJet::GetP:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv->P();
 }
@@ -135,7 +140,6 @@ Double_t AliJet::GetE(Int_t i)
 {
 // Get energy of jet i
   if (OutOfRange(i, "AliJet::GetE:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv->E();
 }
@@ -146,7 +150,6 @@ Double_t AliJet::GetPt(Int_t i)
 {
 // Get transverse momentum of jet i
   if (OutOfRange(i, "AliJet::GetPt:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv->Pt();
 }
@@ -157,7 +160,6 @@ Double_t AliJet::GetEta(Int_t i)
 {
 // Get eta of jet i
   if (OutOfRange(i, "AliJet::GetEta:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv->Eta();
 }
@@ -168,7 +170,6 @@ Double_t AliJet::GetPhi(Int_t i)
 {
 // Get phi of jet i
   if (OutOfRange(i, "AliJet::GetPhi:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return ( (lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi());
 }
@@ -179,7 +180,6 @@ Double_t AliJet::GetTheta(Int_t i)
 {
 // Get theta of jet i
   if (OutOfRange(i, "AliJet::GetTheta:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv->Theta();
 }
@@ -190,7 +190,6 @@ Double_t AliJet::GetMass(Int_t i)
 {
 // Get invariant mass of jet i
   if (OutOfRange(i, "AliJet::GetMass:")) return -1e30;
-
   TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
   return lv->M();
 }
@@ -209,12 +208,37 @@ void AliJet::AddJet(Double_t px, Double_t py, Double_t pz, Double_t e)
 void AliJet::SetInJet(Int_t* j)
 {
   // set information of which input object belongs
-  // to each jet. filled in by AliJetFinder
+  // to each jet. If zero, object was not assigned to
+  // a jet, if n,positive, it was assiged to jet n
+  // if n, negative, it is within cone of jet n, but
+  // it did not passed the user cuts. filled in by AliJetFinder
+
   if (fNInput>0) fInJet.Set(fNInput, j);
 }
 
 ////////////////////////////////////////////////////////////////////////
 
+void AliJet::SetEtaIn(Float_t* r)
+{
+  if (fNInput>0) fEtaIn.Set(fNInput, r);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJet::SetPtIn(Float_t* pt)
+{
+  if (fNInput>0) fPtIn.Set(fNInput, pt);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJet::SetPhiIn(Float_t* x)
+{
+  if (fNInput>0) fPhiIn.Set(fNInput, x);
+}
+
+////////////////////////////////////////////////////////////////////////
+
 void AliJet::SetPtFromSignal(Float_t* p)
 {
   // set information of percentage of pt of jets
@@ -233,6 +257,13 @@ void AliJet::SetMultiplicities(Int_t* m)
 
 ////////////////////////////////////////////////////////////////////////
 
+void AliJet::SetNCells(Int_t* n)
+{
+  if (fNJets>0) fNCells.Set(fNJets, n);
+}
+
+////////////////////////////////////////////////////////////////////////
+
 void AliJet::ClearJets(Option_t *option)
 {
   // reset all values
@@ -241,6 +272,11 @@ void AliJet::ClearJets(Option_t *option)
   fNJets=0;
   fMultiplicities.Set(0);
   fInJet.Set(0); 
+  fPtFromSignal.Set(0);
+  fPhiIn.Set(0);
+  fEtaIn.Set(0);
+  fPtIn.Set(0);
+  fNCells.Set(0);
 }
 
 ////////////////////////////////////////////////////////////////////////
index d9b4389..c9c5da7 100644 (file)
@@ -31,10 +31,16 @@ class AliJet : public TObject
   TClonesArray* GetJets() const {return fJets;}
   TArrayI GetInJet() const {return fInJet;}
   TArrayI GetMultiplicities() const {return fMultiplicities;}
+  TArrayI GetNCells() const {return fNCells;}
   TArrayF GetPtFromSignal() const {return fPtFromSignal;}
+  TArrayF GetEtaIn() const { return fEtaIn; }
+  TArrayF GetPhiIn() const { return fPhiIn; }
+  TArrayF GetPtIn() const { return fPtIn; }
+  Double_t GetEtAvg() const { return fEtAvg; }
 
   TLorentzVector* GetJet(Int_t i);
-  Int_t GetMultiplicity(Int_t i);
+  Int_t GetMultiplicity(Int_t i) const;
+  Int_t GetNCell(Int_t i) const;
   Double_t GetPx(Int_t i);
   Double_t GetPy(Int_t i);
   Double_t GetPz(Int_t i);
@@ -49,10 +55,15 @@ class AliJet : public TObject
   // Setters
   void SetNinput(Int_t i) {fNInput = i;}
   void AddJet(Double_t px, Double_t py, Double_t pz, Double_t e);
-  void SetInJet(Int_t* j);
   void SetMultiplicities(Int_t* m);
+  void SetNCells(Int_t* n);
   void SetPtFromSignal(Float_t* p);
-
+  void SetEtaIn(Float_t* eta);
+  void SetPhiIn(Float_t* phi);
+  void SetPtIn(Float_t* pt);
+  void SetInJet(Int_t* idx);
+  void SetEtAvg(Double_t et) { fEtAvg = et; }
+  
   // others
   Bool_t OutOfRange(Int_t i, const char *s) const;
   void ClearJets(Option_t *option="");
@@ -62,11 +73,18 @@ class AliJet : public TObject
 
   Int_t fNInput;               // number of input objects
   Int_t fNJets;                // number of jets found
+  Double_t fEtAvg;             // average background et per cell
+
   TArrayI fInJet;              // i-input object belongs to k-jet 
   TArrayI fMultiplicities;     // Multiplicity of each jet
+  TArrayI fNCells;             // Number of cells in jet
   TArrayF fPtFromSignal;       // percentage of pt from signal
   TClonesArray* fJets;         // 4-momenta of jets
 
+  TArrayF fEtaIn;              // Arrays of input particles kine:Eta
+  TArrayF fPhiIn;              // Arrays of input particles kine:Phi
+  TArrayF fPtIn;               // Arrays of input particles kine:Pt
+  
   ClassDef(AliJet,1)
 };
  
index 2c7e0cc..0f9787d 100644 (file)
  
 //---------------------------------------------------------------------
 // JetAnalysis class 
-// Analyse Jets
-// Author: andreas.morsch@cern.ch
+// Analyse Jets (already found jets)
+// Authors: andreas.morsch@cern.ch, jgcn@mail.cern.ch
+//          mercedes.lopez.noriega@cern.ch
 //---------------------------------------------------------------------
  
+#include <Riostream.h>
 #include "AliJetAnalysis.h"
 ClassImp(AliJetAnalysis)
-////////////////////////////////////////////////////////////////////////
+  
+// root
+#include <Riostream.h>
 #include <TH1F.h>
+#include <TH2F.h>
 #include <TProfile.h>
-#include <TCanvas.h>
 #include <TFile.h>
 #include <TTree.h>
 #include <TStyle.h>
 #include <TSystem.h>
 #include <TLorentzVector.h>
-
+#include <TMath.h>
+// aliroot
 #include "AliJetProductionDataPDC2004.h"
 #include "AliJet.h"
+#include "AliJetKineReaderHeader.h"
 #include "AliJetESDReaderHeader.h"
 #include "AliUA1JetHeader.h"
 #include "AliLeading.h"
+  
+  
+AliJetAnalysis::AliJetAnalysis()
+{
+  // Constructor
+  fDirectory     = 0x0;   
+  fBkgdDirectory = 0x0;   
+  fFile          = "anajets.root";   
+  fEventMin      =  0;
+  fEventMax      = -1;
+  fRunMin        =  0;
+  fRunMax        = 11;
+  fminMult       =  0;
+  fPercentage    = -1.0;
+  fEfactor       = 1.0;
+  SetReaderHeader();
+  
+  // for Analize()
+  fPythia    = kFALSE;
+  fDoPart    = kTRUE;
+  fDoGenJ    = kTRUE;
+  fDoRecJ    = kTRUE;
+  fDoKine    = kTRUE;
+  fDoCorr    = kTRUE;
+  fDoCorr50  = kFALSE;
+  fDoShap    = kTRUE;
+  fDoFrag    = kTRUE;
+  fDoTrig    = kTRUE;
+  fDoJt      = kTRUE;
+  fDodNdxi   = kTRUE;
+  fDoBkgd    = kFALSE;
+
+  fPart      = 0;
+  fGenJ      = 0;
+  fRecJ      = 0;
+  fRecB      = 0;
+  fWeight    = 1.0;
+  fWdEdr     = 0.0;
+  fPartPtCut = 0.0;
+  fWJt       = 0.0;
+  fWdNdxi    = 0.0;
+
+  // for bkgd
+  fp0    = 0.0;
+  fPtJ   = 0.0;
+  fEJ    = 0.0;
+  fEtaJ  = 0.0;
+  fPhiJ  = 0.0;
+  fjv3X = 0.0;
+  fjv3Y = 0.0;
+  fjv3Z = 0.0;
+
+  // kine histos
+  fRKineEneH = 0;
+  fRKinePtH  = 0;
+  fRKinePhiH = 0;
+  fRKineEtaH = 0;
+  fGKineEneH = 0;
+  fGKinePtH  = 0;
+  fGKinePhiH = 0;
+  fGKineEtaH = 0;
+  fPKineEneH = 0;
+  fPKinePtH  = 0;
+  fPKinePhiH = 0;
+  fPKineEtaH = 0;
+  
+  // correlation histograms
+  fPGCorrEneH = 0;
+  fPGCorrPtH  = 0;
+  fPGCorrEtaH = 0;
+  fPGCorrPhiH = 0;
+  fPRCorrEneH = 0;
+  fPRCorrPtH  = 0;
+  fPRCorrEtaH = 0;
+  fPRCorrPhiH = 0;
+  fRGCorrEneH = 0;
+  fRGCorrPtH  = 0;
+  fRGCorrEtaH = 0;
+  fRGCorrPhiH = 0;
+  
+  // correlation histograms when one particle 
+  // has more than 50% of the energy of the jet
+  fPRCorr50EneH = 0;
+  fPRCorr50PtH  = 0;
+  fPRCorr50EtaH = 0;
+  fPRCorr50PhiH = 0;
+  fRGCorr50EneH = 0;
+  fRGCorr50PtH  = 0;
+  fRGCorr50EtaH = 0;
+  fRGCorr50PhiH = 0;
+
+  // shape histos
+  fWShapR    = 0.0;
+  fRShapSelH = 0;
+  fRShapRejH = 0;
+  fRShapAllH = 0;  
+  
+  // fragmentation function histos
+  fWFragR    = 0.0;
+  fRFragSelH = 0;
+  fRFragRejH = 0;
+  fRFragAllH = 0; 
+
+  // trigger bias histos
+  fGTriggerEneH  = 0;
+  fRTriggerEneH  = 0;
+  fGPTriggerEneH = 0;
+  fPTriggerEneH  = 0;
+  
+  // dE/dr histo and its dr
+  fdEdrH   = 0;
+  fdEdrB   = 0;
+  fPtEneH2 = 0;
+  fdEdrW   = 0;
+  fdrdEdr = 1.0;
+
+  // jt histo and its dr
+  fJtH   = 0;
+  fJtB   = 0;
+  fJtW   = 0;
+  fdrJt = 1.0;
+
+  // dN/dxi histo and its dr
+  fdNdxiH   = 0;
+  fdNdxiB   = 0;
+  fdNdxiW   = 0;
+  fPtEneH   = 0;
+  fdrdNdxi = 0.7;
+
+  // initialize weight for dE/dr histo
+  SetdEdrWeight();
+}
+
+AliJetAnalysis::~AliJetAnalysis()
+{
+  // Destructor
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// define histogrames 
+
+void AliJetAnalysis::DefineHistograms()
+{
+  // Define the histograms to be filled
+  if (fDoKine) DefineKineH();
+  if (fDoCorr) DefineCorrH();
+  if (fDoCorr50) DefineCorr50H();
+  if (fDoShap) DefineShapH();
+  if (fDoFrag) DefineFragH();
+  if (fDoTrig) DefineTrigH();
+  if (fDoJt) DefineJtH();
+  if (fDodNdxi) DefinedNdxiH();
+}
+
+void AliJetAnalysis::DefineKineH()
+{
+  // leading particle    
+  if (fDoPart) {
+    fPKineEneH = new TH1F("PKineEne","Energy of leading particle",50,0,200);
+    SetProperties(fPKineEneH,"Energy (GeV)","Entries");
+    fPKinePtH = new TH1F("PKinePt","Pt of leading particle",50,0,200);
+    SetProperties(fPKinePtH,"P_{T} (GeV)","Entries");
+    fPKinePhiH = new TH1F("PKinePhiH","Azimuthal angle of leading particle",
+                         90,0.,2.0*TMath::Pi());
+    SetProperties(fPKinePhiH,"#phi","Entries");
+    fPKineEtaH = new TH1F("PKineEtaH","Pseudorapidity of leading particle",
+                         40,-1.0,1.0);
+    SetProperties(fPKineEtaH,"#eta","Entries");
+  }
+  // leading generated jet
+  if (fDoGenJ) {
+    fGKineEneH = new TH1F("GKineEne","Energy of generated jet",50,0,200);
+    SetProperties(fGKineEneH,"Energy (GeV)","Entries");
+    fGKinePtH = new TH1F("GKinePt","Pt of generated jet",50,0,200);
+    SetProperties(fGKinePtH,"P_{T} (GeV)","Entries");
+    fGKinePhiH = new TH1F("GKinePhiH","Azimuthal angle of generated jet",
+                         90,0.,2.0*TMath::Pi());
+    SetProperties(fGKinePhiH,"#phi","Entries");
+    fGKineEtaH = new TH1F("GKineEtaH","Pseudorapidity of generated jet",
+                         40,-1.0,1.0);
+    SetProperties(fGKineEtaH,"#eta","Entries");
+  }
+  // leading reconstructed jet
+  if (fDoRecJ) {
+    fRKineEneH = new TH1F("RKineEne","Energy of reconstructed jet",50,0,200);
+    SetProperties(fRKineEneH,"Energy (GeV)","Entries");
+    fRKinePtH = new TH1F("RKinePt","Pt of reconstructed jet",50,0,200);
+    SetProperties(fRKinePtH,"P_{T} (GeV)","Entries");
+    fRKinePhiH = new TH1F("RKinePhiH","Azimuthal angle of reconstructed jet",
+                         90,0.,2.0*TMath::Pi());
+    SetProperties(fRKinePhiH,"#phi","Entries");
+    fRKineEtaH = new TH1F("RKineEtaH","Pseudorapidity of reconstructed jet",
+                         40,-1.0,1.0);
+    SetProperties(fRKineEtaH,"#eta","Entries");
+  } 
+}
 
-    AliJetAnalysis::AliJetAnalysis()
+void AliJetAnalysis::DefineCorrH()
 {
-    // Constructor
-    fDirectory    = 0x0;   
-    fEventMin     =  0;
-    fEventMax     = -1;
-    fRunMin       =  0;
-    fRunMax       = 11;
+  // correlation 
+  if (fDoPart && fDoGenJ) {
+    fPGCorrEneH = new TH2F("PGCorrEne","Energy correlation part-gen jet",
+                          40,0,200,40,0,200);
+    SetProperties(fPGCorrEneH,"Part Energy (GeV)","Gen Jet Energy (GeV)");
+    fPGCorrPtH = new TH2F("PGCorrPt","Pt correlation part-gen jet",
+                         40,0,200,40,0,200);
+    SetProperties(fPGCorrPtH,"Part P_{T} (GeV)","Gen Jet P_{T} (GeV)");
+    fPGCorrEtaH = new TH2F("PGCorrEta","Pseudorapidity correlation part-gen jet",
+                          40,-1.0,1.0,40,-1.0,1.0);
+    SetProperties(fPGCorrEtaH,"Part #eta","Gen Jet #eta");
+    fPGCorrPhiH = new TH2F("PGCorrPhi","Azimuthal angle correlation part-gen jet",
+                          90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+    SetProperties(fPGCorrPhiH,"Part #phi","Gen Jet #phi");
+  }  
+  if (fDoPart && fDoRecJ) {
+    fPRCorrEneH = new TH2F("PRCorrEne","Energy correlation part-rec jet",
+                          40,0,200,40,0,200);
+    SetProperties(fPRCorrEneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
+    fPRCorrPtH = new TH2F("PRCorrPt","Pt correlation part-rec jet",
+                         40,0,200,40,0,200);
+    SetProperties(fPRCorrPtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
+    fPRCorrEtaH = new TH2F("PRCorrEta","Pseudorapidity correlation part-rec jet",
+                          40,-1.0,1.0,40,-1.0,1.0);
+    SetProperties(fPRCorrEtaH,"part #eta","Rec Jet #eta");
+    fPRCorrPhiH = new TH2F("PRCorrPhi","Azimuthal angle correlation part-rec jet",
+                          90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+    SetProperties(fPRCorrPhiH,"Part #phi","Rec Jet #phi");
+  }  
+  if (fDoGenJ && fDoRecJ) {
+    fRGCorrEneH = new TH2F("RGCorrEne","Energy correlation rec jet-gen jet",
+                          40,0,200,40,0,200);
+    SetProperties(fRGCorrEneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
+    fRGCorrPtH = new TH2F("RGCorrPt","Pt correlation rec jet-gen jet",
+                         40,0,200,40,0,200);
+    SetProperties(fRGCorrPtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
+    fRGCorrEtaH = new TH2F("RGCorrEta","Pseudorapidity correlation rec jet-gen jet",
+                          40,-1.0,1.0,40,-1.0,1.0);
+    SetProperties(fRGCorrEtaH,"Rec Jet #eta","Gen Jet #eta");
+    fRGCorrPhiH = new TH2F("RGCorrPhi","Azimuthal angle correlation rec jet-gen jet",
+                          90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+    SetProperties(fRGCorrPhiH,"Rec Jet #phi","Gen Jet #phi");
+  }
 }
 
-void AliJetAnalysis::Analyse() 
+void AliJetAnalysis::DefineCorr50H()
 {
-    //
-    // Some histos
-    //
-    TH1F::AddDirectory(0);
-    TProfile::AddDirectory(0);
+  // correlation 
+  if (fDoPart && fDoRecJ) {
+    fPRCorr50EneH = new TH2F("PRCorr50Ene","Energy correlation part-rec jet",
+                            40,0,200,40,0,200);
+    SetProperties(fPRCorr50EneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
+    fPRCorr50PtH = new TH2F("PRCorr50Pt","Pt correlation part-rec jet",
+                           40,0,200,40,0,200);
+    SetProperties(fPRCorr50PtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
+    fPRCorr50EtaH = new TH2F("PRCorr50Eta","Pseudorapidity correlation part-rec jet",
+                            40,-1.0,1.0,40,-1.0,1.0);
+    SetProperties(fPRCorr50EtaH,"part #eta","Rec Jet #eta");
+    fPRCorr50PhiH = new TH2F("PRCorr50Phi","Azimuthal angle correlation part-rec jet",
+                            90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+    SetProperties(fPRCorr50PhiH,"Part #phi","Rec Jet #phi");
+  }
+  
+  if (fDoGenJ && fDoRecJ) {
+    fRGCorr50EneH = new TH2F("RGCorr50Ene","Energy correlation rec jet-gen jet",
+                            40,0,200,40,0,200);
+    SetProperties(fRGCorr50EneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
+    fRGCorr50PtH = new TH2F("RGCorr50Pt","Pt correlation rec jet-gen jet",
+                           40,0,200,40,0,200);
+    SetProperties(fRGCorr50PtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
+    fRGCorr50EtaH = new TH2F("RGCorr50Eta","Pseudorapidity correlation rec jet-gen jet",
+                            40,-1.0,1.0,40,-1.0,1.0);
+    SetProperties(fRGCorr50EtaH,"Rec Jet #eta","Gen Jet #eta");
+    fRGCorr50PhiH = new TH2F("RGCorr50Phi","Azimuthal angle correlation rec jet-gen jet",
+                            90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+    SetProperties(fRGCorr50PhiH,"Rec Jet #phi","Gen Jet #phi");
+  }
+}
+
+void AliJetAnalysis::DefineShapH()
+{
+  // leading reconstructed jet
+  if (fDoRecJ) {
+    fdEdrH = new TH2F("fdEdrH","dE/dr histo",20,0,1,40,0,200);
+    SetProperties(fdEdrH,"r","Rec Jet P_{T}");
+    fdEdrB = new TH2F("fdEdrB","dE/dr bkgdhisto",20,0,1,40,0,200);
+    SetProperties(fdEdrB,"r","Rec P_{T}");
+    fPtEneH2 = new TH2F("fPtEneH2","fPtEneH2",40,0,200,40,0,200);
+    SetProperties(fPtEneH2,"Rec Jet E","Rec Jet P_{T}");
+    fdEdrW = new TH1F("fdEdrW","weights for dE/dr",40,0,200);
+    SetProperties(fdEdrW,"Rec Jet P_{T}","weights");
     
-    TH1F* e0H    = new TH1F("e0H"  ,"Jet Energy (reconstructed)",      40,  0., 200.);
-    TH1F* e1H    = new TH1F("e1H" , "Jet Energy (generated)",          40,  0., 200.);
-    TH1F* e2H    = new TH1F("e2H" , "Jet Energy (generated, nrec = 0", 40,  0., 200.);
-    TH1F* e3H    = new TH1F("e3H" , "Jet Energy (leading)",            40,  0., 200.);
-    TH1F* e4H    = new TH1F("e4H" , "Jet Energy (reconstructed: 105 < Egen < 125", 40,  0., 200.);
+    fRShapSelH = new TH1F("RShapSel","Shape of generated jets (sel part)",20,0.,1.);
+    SetProperties(fRShapSelH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
+    fRShapRejH = new TH1F("RShapRej","Shape of generated jets (rej part)",20,0.,1.);
+    SetProperties(fRShapRejH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
+    fRShapAllH = new TH1F("RShapAll","Shape of generated jets (all part)",20,0.,1.);
+    SetProperties(fRShapAllH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
+  }
+}
 
-    TH1F* e5H    = new TH1F("e5H" , "Jet Energy (generated)", 40,  0., 200.);
-    TH1F* e6H    = new TH1F("e6H" , "Jet Energy (generated)", 40,  0., 200.);
-    TH1F* e7H    = new TH1F("e7H" , "Jet Energy (generated)", 40,  0., 200.);
-    TH1F* e8H    = new TH1F("e8H" , "Jet Energy (generated)", 40,  0., 200.);
+void AliJetAnalysis::DefineJtH()
+{
+  // Define the histogram for J_T
+  if (fDoRecJ) {
+    fJtH = new TH2F("fJtH","J_{T} histogram",80,0.,4.,40,0.,200.);
+    SetProperties(fJtH,"J_{T}","Rec Jet P_{T}");
+    fJtB = new TH2F("fJtB","J_{T} bkgd histogram",80,0.,4.,40,0.,200.);
+    SetProperties(fJtB,"J_{T}","Rec P_{T}");
+    fJtW = new TH1F("fJtW","J_{T} weight",40,0,200);
+    SetProperties(fJtW,"J_{T}W","weight");
+  }
+}
 
-    TProfile* r5H    = new TProfile("r5H" , "rec/generated", 20,  0., 200, 0., 1., "S");
-    TProfile* r6H    = new TProfile("r6H" , "rec/generated", 20,  0., 200, 0., 1., "S");
+void AliJetAnalysis::DefinedNdxiH()
+{
+  // Define the histogram for dN/dxi
+  if (fDoRecJ) {
+    fdNdxiH = new TH2F("fdNdxiH","dN/d#xi histo",200,0,10,40,0,200);
+    SetProperties(fdNdxiH,"xi","Rec Jet P_{T}");
+    fdNdxiB = new TH2F("fdNdxiB","dN/d#xi bkgd histo",200,0,10,40,0,200);
+    SetProperties(fdNdxiB,"xi","Rec P_{T}");
+    fdNdxiW = new TH1F("fdNdxiW","dN/d#xi histo",40,0,200);
+    SetProperties(fdNdxiW,"Rec Jet P_{T}","weights");
+    fPtEneH = new TH2F("fPtEneH","fPtEneH",40,0,200,40,0,200);
+    SetProperties(fPtEneH,"Rec Jet E","Rec Jet P_{T}");
+  }
+}
 
-    TProfile* r7H    = new TProfile("r7H" , "rec/generated", 20,  0., 200, 0., 1., "S");
-    TProfile* r8H    = new TProfile("r8H" , "rec/generated", 20,  0., 200, 0., 1., "S");
+void AliJetAnalysis::DefineFragH()
+{
+  // leading reconstructed jet
+  if (fDoRecJ) {
+    fRFragSelH = new TH1F("RFragSel","Frag Fun of reconstructed jets (sel part)",20,0.,10.);
+    SetProperties(fRFragSelH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+    fRFragRejH = new TH1F("RFragRej","Frag Fun of reconstructed jets (rej part)",20,0.,10.);
+    SetProperties(fRFragRejH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+    fRFragAllH = new TH1F("RFragAll","Frag Fun of reconstructed jets (all part)",20,0.,10.);
+    SetProperties(fRFragAllH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+  }
+}
 
+void AliJetAnalysis::DefineTrigH()
+{
+  // generated energy
+  fGTriggerEneH = new TProfile("GTriggerEne","Generated energy (trigger bias)", 
+                              20, 0., 200., 0., 1., "S");    
+  fGTriggerEneH->SetXTitle("E_{gen}");
+  fGTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
+  // reconstructed energy
+  fRTriggerEneH = new TProfile("RTriggerEne","Reconstructed energy (trigger bias)", 
+                              20, 0., 200., 0., 1., "S");  
+  fRTriggerEneH->SetXTitle("E_{rec}");
+  fRTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
+  // generated energy vs generated/leading
+  fGPTriggerEneH = new TProfile("GPTriggerEne","Generated energy (trigger bias)", 
+                               20, 0., 200., 0., 1., "S");    
+  fGPTriggerEneH->SetXTitle("E_{gen}");
+  fGPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
+  // leading particle energy
+  fPTriggerEneH = new TProfile("PTriggerEne","Leading particle energy (trigger bias)", 
+                              20, 0., 200., 0., 1., "S");  
+  fPTriggerEneH->SetXTitle("E_{L}/0.2");
+  fPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
 
-    TH1F* dr1H = new TH1F("dr1H", "delta R",  160, 0.,   2.);
-    TH1F* dr2H = new TH1F("dr2H", "delta R",  160, 0.,   2.);
-    TH1F* dr3H = new TH1F("dr4H", "delta R",  160, 0.,   2.);
+}
 
-    TH1F* etaH  = new TH1F("etaH",  "eta",  160, -2.,   2.);
-    TH1F* eta1H = new TH1F("eta1H", "eta",  160, -2.,   2.);
-    TH1F* eta2H = new TH1F("eta2H", "eta",  160, -2.,   2.);
+void AliJetAnalysis::SetProperties(TH1* h,const char* x, const char* y) const
+{
+  //Set properties of histos (x and y title and error propagation)
+  h->SetXTitle(x);
+  h->SetYTitle(y);
+  h->Sumw2();
+}
+////////////////////////////////////////////////////////////////////////
+// compute weights for dE/dr histogram
 
-    TH1F* phiH   = new TH1F("phiH",   "phi",  160, -3.,   3.);
-    TH1F* dphiH  = new TH1F("dphiH",  "phi",  160,  0.,   3.1415);
-    TH1F* phi1H  = new TH1F("phi1H", "phi",   160,  0.,   6.28);
-    TH1F* phi2H  = new TH1F("phi2H", "phi",   160,  0.,   6.28);
-    
+void AliJetAnalysis::SetdEdrWeight()
+{
+  // Due to the limited acceptance, each bin in the dE/dr has a different
+  // weight given by the ratio of the area of a disk dR to the area
+  // within the eta acceptance. The weight depends on the eta position 
+  // of the jet. Here a look up table for the weights is computed. It
+  // assumes |etaJet|<0.5 and |eta_lego|<0.9. It makes bins of 0.05
+  // units in eta. Note that this is fixed by the bin width chosen for
+  // the histogram --> this weights are tailored for the specific 
+  // histogram definition used here!
+  
+  // two dimensional table: first index, bin in eta of jet, second
+  // index bin of dE/dr histo
+
+  Int_t nBins = 20;
+  Float_t xMin = 0.0;
+  Float_t xMax = 1.0;
+  Float_t binSize = (xMax-xMin)/nBins;
 
-    TProfile* drP1    = new TProfile("drP1" , "Delta_R", 20,  0., 200, -1., 1., "S");
-    TProfile* drP2    = new TProfile("drP2" , "Delta_R", 20,  0., 200, -1., 1., "S");
+  Float_t da,ds,r1,r2,h1,h2,a1,a2,theta1,theta2,area1,area2;
+  Int_t ji;
+  for (Int_t i=0;i<(nBins/2);i++) {
+    // index of first histo bin needing a scale factor
+    ji=(nBins-2)-i;
+    for(Int_t j=0;j<nBins;j++) {
+      // area of ring.
+      da = TMath::Pi()*(binSize*binSize)*(1.0+2.0*j);
+      ds = 1.0;
+      if (j>=ji) { // compute missing area using segments of circle
+       r1=j*binSize;
+       r2=(j+1)*binSize;
+       h1=(j-ji)*binSize;
+       h2=(j+1-ji)*binSize;
+       a1=2.0*TMath::Sqrt(2.0*h1*r1-h1*h1);
+       a2=2.0*TMath::Sqrt(2.0*h2*r2-h2*h2);
+       theta1=2*TMath::ACos((r1-h1)/r1);
+       theta2=2*TMath::ACos((r2-h2)/r2);
+       area1=binSize*(r1*r1*theta1-a1*(r1-h1));
+       area2=binSize*(r2*r2*theta2-a2*(r2-h2));
+       ds = (da-(area2-area1))/da;
+      }
+      fWeightdEdr[i][j]=ds/da;
+    }
+  }
+}
+
+// get weight for dE/dr histogram
+Float_t AliJetAnalysis::GetdEdrWeight(Float_t etaJet, Float_t r)
+{
+  // Return the correponding weight for the dE/dr plot
+  Int_t nBins = 20;
+  Float_t xMin = 0.0;
+  Float_t xMax = 1.0;
+  Float_t binSize = (xMax-xMin)/nBins;
+
+  Float_t eta = TMath::Abs(etaJet);
+  if ((eta > 0.5) || (r > fdrdEdr)) return 0.0;
+  Int_t binJet = (Int_t) TMath::Floor(eta/binSize);
+  Int_t binR = (Int_t) TMath::Floor(r/binSize);
+  Float_t w = fWeightdEdr[binJet][binR];
+  return w;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// fill histograms 
+
+void AliJetAnalysis::FillHistograms()
+{
+  // fill histograms 
 
   // Run data 
-    AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004();
+  AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004();
+  
+  // Loop over runs
+  TFile* jFile = 0x0;
+  
+  for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++) {
+    // Open file
+    char fn[256];
+    sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data());
+    jFile = new TFile(fn);
+    printf("  Analyzing run: %d %s\n", iRun,fn);       
+    
+    // Get reader header and events to be looped over
+    AliJetReaderHeader *jReaderH = (AliJetReaderHeader*)(jFile->Get(fReaderHeader));
+    if (fEventMin == -1) fEventMin =  jReaderH->GetFirstEvent();
+    if (fEventMax == -1) {
+      fEventMax =  jReaderH->GetLastEvent();
+    } else {
+      fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent());
+    }
+    
+    AliUA1JetHeader *jH = (AliUA1JetHeader *) (jFile->Get("AliUA1JetHeader"));
     
+    // Calculate weight
+    fWeight = runData->GetWeight(iRun)/ ( (Float_t) (fEventMax - fEventMin + 1));
 
-    // Loop over runs
-    TFile* jFile = 0x0;
-  
-    for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++)
-    {
-
-       // Open file
-       char fn[20];
-       sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data());
-       
-       
-       jFile = new TFile(fn);
-
-       printf("  Analyzing run: %d %s\n", iRun,fn);    
-       // Get jet header and display parameters
-       AliUA1JetHeader* jHeader = 
-           (AliUA1JetHeader*) (jFile->Get("AliUA1JetHeader"));
-       // jHeader->PrintParameters();
-       
-       // Get reader header and events to be looped over
-       AliJetESDReaderHeader *jReaderH = 
-           (AliJetESDReaderHeader*)(jFile->Get("AliJetKineReaderHeader"));
-
-       if (fEventMin == -1) fEventMin =  jReaderH->GetFirstEvent();
-       if (fEventMax == -1) {
-           fEventMax =  jReaderH->GetLastEvent();
-       } else {
-           fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent());
+    // Loop over events
+    for (Int_t i = fEventMin; i < fEventMax; i++) {
+      if (i%100 == 0) printf("  Analyzing event %d / %d \n",i,fEventMax);
+      
+      // Get next tree with AliJet
+      char nameT[100];
+      sprintf(nameT, "TreeJ%d",i);
+      TTree *jetT =(TTree *)(jFile->Get(nameT));
+      if (fDoRecJ) jetT->SetBranchAddress("FoundJet",    &fRecJ);
+      if (fDoGenJ) jetT->SetBranchAddress("GenJet",      &fGenJ);
+      if (fDoPart) jetT->SetBranchAddress("LeadingPart", &fPart);
+      jetT->GetEntry(0);
+      
+      int nJets = fRecJ->GetNJets();  
+      
+      TArrayI inJet = fRecJ->GetInJet();
+      if(inJet.GetSize()>fminMult){     // removes events with low multiplicity
+       if (fDoKine) FillKineH();
+       if (fDoCorr) FillCorrH();
+       if (fDoCorr50) FillCorr50H();
+       if (fDoShap) FillShapH(jH->GetRadius());
+       if (fDoFrag) FillFragH();    
+       if (fDoTrig) FillTrigH();
+       if (fDoJt) FillJtH();
+       if (fDodNdxi) FilldNdxiH();
+      }
+      delete jetT;                       // jet should be deleted before creating a new one
+      if(inJet.GetSize()>fminMult){      // removes events with low multiplicity
+       if (fDoBkgd && nJets>0) FillBkgd(i,iRun);
+      }      
+    } // end loop over events in one file
+    if (jFile) jFile->Close();
+    delete jFile;
+  } // end loop over files
+}
+
+void AliJetAnalysis::FillKineH()
+{
+  // leading particle 
+  if (fDoPart && fPart->LeadingFound()){
+    fPKineEneH->Fill(fPart->GetE(),fWeight);
+    fPKinePtH->Fill(fPart->GetPt(),fWeight);
+    fPKinePhiH->Fill(fPart->GetPhi(),fWeight);
+    fPKineEtaH->Fill(fPart->GetEta(),fWeight);
+  } 
+  // leading generated jet
+  if (fDoGenJ && fGenJ->GetNJets()> 0){
+    fGKineEneH->Fill(fGenJ->GetE(0),fWeight);
+    fGKinePtH->Fill(fGenJ->GetPt(0),fWeight);
+    fGKinePhiH->Fill(fGenJ->GetPhi(0),fWeight);
+    fGKineEtaH->Fill(fGenJ->GetEta(0),fWeight);
+  } 
+  // leading reconstructed jet
+  if (fDoRecJ && fRecJ->GetNJets()> 0) {
+    TArrayF p=fRecJ->GetPtFromSignal();
+    if (p[0]>fPercentage) {
+      fRKineEneH->Fill(fRecJ->GetE(0)/fEfactor,fWeight);
+      fRKinePtH->Fill(fRecJ->GetPt(0),fWeight);
+      fRKinePhiH->Fill(fRecJ->GetPhi(0),fWeight);
+      fRKineEtaH->Fill(fRecJ->GetEta(0),fWeight);
+    }
+  }
+}
+
+void AliJetAnalysis::FillCorrH()
+{
+  // Fill correlation histograms
+  if (fDoPart && fPart->LeadingFound() && fDoGenJ && fGenJ->GetNJets()> 0) 
+    Correlation(fPart->GetLeading(),fGenJ->GetJet(0), 
+               fPGCorrEneH,fPGCorrPtH,fPGCorrEtaH,fPGCorrPhiH);
+  if (fDoPart && fPart->LeadingFound() && fDoRecJ && fRecJ->GetNJets()> 0) {
+    TArrayF p=fRecJ->GetPtFromSignal();
+    if (p[0]>fPercentage) 
+      Correlation(fPart->GetLeading(),fRecJ->GetJet(0), 
+                 fPRCorrEneH,fPRCorrPtH,fPRCorrEtaH,fPRCorrPhiH);
+  }
+  if (fDoGenJ && fGenJ->GetNJets()> 0  && fDoRecJ && fRecJ->GetNJets()> 0) {
+    TArrayF p=fRecJ->GetPtFromSignal();
+    if (p[0]>fPercentage) 
+      Correlation(fRecJ->GetJet(0), fGenJ->GetJet(0),
+                 fRGCorrEneH,fRGCorrPtH,fRGCorrEtaH,fRGCorrPhiH);
+  }
+}
+
+void AliJetAnalysis::Correlation(TLorentzVector *lv1,TLorentzVector *lv2,
+                                TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
+{
+  // Correlation histograms
+  h1->Fill(lv1->E(),lv2->E(),fWeight);
+  h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
+  h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
+  Float_t p1, p2;
+  p1 = ((lv1->Phi() < 0) ? (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
+  p2 = ((lv2->Phi() < 0) ? (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
+  h4->Fill(p1,p2,fWeight);
+}
+
+void AliJetAnalysis::FillCorr50H()
+{
+  // Fill correlation histograms when one particle has > 50% of jet energy
+  if (fDoRecJ && fRecJ->GetNJets()> 0) {
+    TArrayF p = fRecJ->GetPtFromSignal();
+    if (p[0]>fPercentage) {
+      if (fDoPart && fPart->LeadingFound()) 
+       Correlation50(fRecJ, fPart->GetLeading(),fRecJ->GetJet(0), 
+                     fPRCorr50EneH,fPRCorr50PtH,fPRCorr50EtaH,fPRCorr50PhiH);
+      if (fDoGenJ && fGenJ->GetNJets()> 0) 
+       Correlation50(fRecJ, fRecJ->GetJet(0), fGenJ->GetJet(0),
+                     fRGCorr50EneH,fRGCorr50PtH,fRGCorr50EtaH,fRGCorr50PhiH);
+    }
+  }
+}
+
+void AliJetAnalysis::Correlation50(AliJet *j,TLorentzVector *lv1,TLorentzVector *lv2,
+                                  TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
+{
+  // Correlation histograms when one particle has > 50% of jet energy
+    TArrayF ptin = j->GetPtIn();
+    TArrayF etain = j->GetEtaIn();
+    TArrayI inJet = j->GetInJet();
+    
+    Int_t flag = 0;
+    for(Int_t i=0;i<(inJet.GetSize());i++) {
+      if (inJet[i] == 1) {
+       Float_t t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+       Float_t x = (ptin[i]*TMath::Sqrt(1.+1./(t1*t1)))/(j->GetE(0));
+       if (x>0.5) flag++;
+      }
+    }
+    if (flag>1) cout << " Flag = " << flag << endl;
+    if (flag == 1) {
+      h1->Fill(lv1->E(),lv2->E(),fWeight);
+      h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
+      h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
+      Float_t p1, p2;
+      p1 = ((lv1->Phi() < 0) ? 
+           (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
+      p2 = ((lv2->Phi() < 0) ? 
+           (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
+      h4->Fill(p1,p2,fWeight);
+    }
+}
+
+void AliJetAnalysis::FillJtH()
+{
+  // Fill J_T histogram
+  if (fRecJ->GetNJets()> 0) {
+    fjv3X = 0.0; fjv3Y = 0.0; fjv3Z = 0.0;
+    TArrayF p=fRecJ->GetPtFromSignal();
+    if (p[0]>fPercentage) {
+      // initialize
+      const TVector3 kJv3 = fRecJ->GetJet(0)->Vect();
+      fjv3X =  kJv3.X(); fjv3Y =  kJv3.Y(); fjv3Z =  kJv3.Z();
+      TVector3 trk;
+      TArrayI inJet = fRecJ->GetInJet();
+      TArrayF etain = fRecJ->GetEtaIn();
+      TArrayF ptin = fRecJ->GetPtIn();
+      TArrayF phiin = fRecJ->GetPhiIn();
+      Float_t deta, dphi,jt, dr;
+      for(Int_t i=0;i<inJet.GetSize();i++) {
+       deta = etain[i] - fRecJ->GetEta(0);
+       dphi = phiin[i] - fRecJ->GetPhi(0);
+       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+       dr = TMath::Sqrt(deta * deta + dphi * dphi);
+       if ((dr<fdrJt) && (ptin[i] > fPartPtCut)) { 
+         trk.SetPtEtaPhi(ptin[i],etain[i],phiin[i]);
+         jt = TMath::Sin(trk.Angle(kJv3))*trk.Mag();
+         fJtH->Fill(jt,fRecJ->GetPt(0),fWeight);
        }
-       
-       
-       // Calculate weight
-       Float_t wgt = runData->GetWeight(iRun) / Float_t(fEventMax - fEventMin + 1);
-       Float_t ptmin, ptmax;
-       runData->GetPtHardLimits(iRun, ptmin, ptmax);
-       
-       
-       // Loop over events
-       AliJet *jets  = 0x0; 
-       AliJet *gjets = 0x0;
-       AliLeading *leading = 0x0;
-       Float_t egen, etag, econ, erec;
-       
-       
-       for (Int_t i = fEventMin; i < fEventMax; i++) {
-           printf("  Analyzing run: %d  Event %d / %d \n", 
-                  iRun, i, fEventMax);
-           
-           // Het next tree with AliJet
-           char nameT[100];
-           sprintf(nameT, "TreeJ%d",i);
-           TTree *jetT =(TTree *)(jFile->Get(nameT));
-           jetT->SetBranchAddress("FoundJet",    &jets);
-           jetT->SetBranchAddress("GenJet",      &gjets);
-           jetT->SetBranchAddress("LeadingPart", &leading);
-           jetT->GetEntry(0);
-           
-//
-//    Find the jet with the highest E_T within fiducial region
-//
-           Int_t njets = jets->GetNJets();
-           Int_t imax = -1;
-           Int_t jmax = -1;
-           Float_t emax = 0.;
-           
-           for (Int_t ij = 0; ij < njets; ij++) {
-               if (jets->GetPt(ij) > emax && 
-                   TMath::Abs(jets->GetEta(ij)) < 0.60) {
-                   emax = jets->GetPt(ij);
-                   jmax = imax;
-                   imax = ij;
-               }
-           }
-           
-           
-           if (imax == -1) {
-               Int_t ngen = gjets->GetNJets();
-               if(ngen>0) e2H->Fill(gjets->GetPt(0), wgt);
-           } else {
-//           printf("Reconstructed Jet %5d %13.3f %13.3f %13.3f\n", 
-//                    imax, emax, jets->GetEta(imax), jets->GetPhi(imax));
-               econ = jets->GetPt(imax);
-               erec = jets->GetPt(imax) / 0.65;
-               dr1H->Fill(jets->GetEta(imax) - gjets->GetEta(0));
-//
-//    Find the generated jet closest to the reconstructed 
-//
-               Float_t rmin;
-               Float_t etamin = 1.e6;
-               
-               Int_t   igen;
-               Float_t etaj = jets->GetEta(imax);
-               Float_t phij = jets->GetPhi(imax);
-               
-               Int_t ngen = gjets->GetNJets();
-               
-               if (ngen != 0) {
-                   rmin = 1.e6;
-                   igen = -1;
-                   for (Int_t j = 0; j < ngen; j++) {
-                       etag = gjets->GetEta(j);
-                       Float_t phig = gjets->GetPhi(j);
-                       Float_t deta = etag - etaj;
-                       Float_t dphi = TMath::Abs(phig - phij);
-                       if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
-                       Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
-                       if (r  < rmin) {
-                           rmin = r;
-                           etamin = deta;
-                           igen = j;
-                       }
-                   }
-                   
-                   egen = gjets->GetPt(igen);
-                   e1H->Fill(gjets->GetPt(igen), wgt);
-                   etag = gjets->GetEta(igen);
-                   Float_t phig = gjets->GetPhi(igen);
-                   Float_t dphi = phig - phij;
-                   
-//               if (econ < ptmax) {
-                   e0H->Fill(erec, wgt);
-//               } else {
-//                   e0H->Fill(erec, 6.7e-6);
-//               }
-                   
-                   
-                   
-                   if (egen > 20. && egen < 40.) {
-                       phiH->Fill((dphi));
-                       etaH->Fill(etag - etaj, wgt);
-                       phi1H->Fill(phig);
-                       phi2H->Fill(phij);                
-                       e4H->Fill(jets->GetPt(imax));
-                   }
-                   
-                   if (erec > 90. && erec < 110. && rmin < 0.1) {
-                       e5H->Fill(egen, wgt);
-                       dr2H->Fill(rmin);
-                       if (egen < 30.) {
-                           printf("Strange jet %6d %13.3f %13.3f %13.3f \n", 
-                                  imax, etaj, phij, erec);
-                           for (Int_t j = 0; j < ngen; j++) {
-                               printf("Generated %6d %13.3f %13.3f %13.3f\n", 
-                                      j, gjets->GetEta(j), 
-                                      gjets->GetPhi(j), gjets->GetPt(j));
-                               
-                           }
-                       }
-                   }
-                   
-                   if (rmin < 0.1) {
-                       r5H->Fill(egen, jets->GetPt(imax)/egen, wgt);
-                       r6H->Fill(jets->GetPt(imax) 
-                                 / 0.4, jets->GetPt(imax)/egen, wgt);
-                       e8H->Fill(erec, wgt);
-                   }
-                   
-                   if (rmin < 0.1) {
-                       drP1->Fill(egen, etamin, wgt);
-                   }
-               } // ngen !=0
-           } // has reconstructed jet
-           
-//
-//   Leading particle
-//
-           if (leading->LeadingFound()) {
-               Float_t etal = leading->GetLeading()->Eta();
-               Float_t phil = leading->GetLeading()->Phi();
-               Float_t el   = leading->GetLeading()->E();
-//       printf("Leading %f %f %f \n", etal, phil, el);
-         
-               e3H->Fill(el, wgt);
-         
-               Float_t rmin = 1.e6;
-               Float_t etamin = 1.e6;
-         
-               Int_t igen = -1;
-               Int_t ngen = gjets->GetNJets();
-               for (Int_t j = 0; j < ngen; j++) {
-                   etag = gjets->GetEta(j);
-                   Float_t phig = gjets->GetPhi(j);
-                   Float_t deta = etag-etal;
-                   Float_t dphi = TMath::Abs(phig - phil);
-                   if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
-                   
-                   Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
-             
-                   if (r  < rmin) {
-                       rmin = r;
-                       igen = j;
-                       etamin = deta;
-                   }
-               }
-               if (egen > 20. && egen < 40.) {
-                   dr3H->Fill(rmin);
-                   eta1H->Fill(etag-etal, wgt);
-                   
-               }
-         
-               if (el > 54. && el  < 66.) e6H->Fill(egen, wgt);
-               if (rmin < 0.3) {
-                   r7H->Fill(egen, el/egen, wgt);
-                   r8H->Fill(el / 0.2, el/egen, wgt);
-               }
-               
-               if (rmin < 0.1) {
-                   drP2->Fill(egen, etamin, wgt);
-               }
-           } // if leading particle
-           
-
-
-//
-// Generated Jet
-//
-           Int_t ngen = gjets->GetNJets();
-           emax = 0.;
-           imax = -1;
-           
-           for (Int_t j = 0; j < ngen; j++) {
-               if (gjets->GetPt(j) > 
-                   emax && TMath::Abs(gjets->GetEta(j)) < 0.5) {
-                   emax = gjets->GetPt(j);
-                   imax = j;
-               }
-           }
-           if (imax != -1) e7H->Fill(emax, wgt);
-           
-           
-           delete jetT;
-           
-       } // events
-       if (jFile) jFile->Close();
-       delete jFile;
-       
-    } // runs
+      }
+      fJtW->Fill(fRecJ->GetPt(0),fWeight);
+      fWJt+=fWeight;
+    }
+  }
+}
+
+void AliJetAnalysis::FilldNdxiH()
+{
+  // Fill dN/d#xi histograms
+  if (fRecJ->GetNJets()> 0) {
+    TArrayF p=fRecJ->GetPtFromSignal();
+    if (p[0]>fPercentage) {
+      fp0 = p[0];
+      TArrayI inJet = fRecJ->GetInJet();
+      TArrayF etain = fRecJ->GetEtaIn();
+      TArrayF ptin = fRecJ->GetPtIn();
+      TArrayF phiin = fRecJ->GetPhiIn();
+      Float_t xi,t1,ene,dr,deta,dphi;
+      for(Int_t i=0;i<inJet.GetSize();i++) {
+       deta = etain[i] - fRecJ->GetEta(0);
+       dphi = phiin[i] - fRecJ->GetPhi(0);
+       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+       dr = TMath::Sqrt(deta * deta + dphi * dphi);
+       if ((dr<fdrdNdxi) && (ptin[i] > fPartPtCut)) { 
+         t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+         ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
+         xi = (Float_t) TMath::Log((fRecJ->GetE(0)/fEfactor)/ene);
+         fdNdxiH->Fill(xi,fRecJ->GetPt(0),fWeight);
+       }  
+      } 
+      fdNdxiW->Fill(fRecJ->GetPt(0),fWeight);
+      fPtEneH->Fill(fRecJ->GetE(0)/fEfactor,fRecJ->GetPt(0),fWeight);
+      fWdNdxi+=fWeight;
+    }
+  }
+}
+
+void AliJetAnalysis::FillBkgd(Int_t eventN, Int_t runN)
+{
+  // Background calculated with hijing events (no pythia)  
+  if (fp0>fPercentage) {
+    fPtJ=0.,fEJ=0.,fEtaJ=0.,fPhiJ=0.;
+    // Background calculated with hijing events (no pythia)
+    AliJetProductionDataPDC2004* runDataB = new AliJetProductionDataPDC2004();
+    TFile* jFileB =0x0;;
+    char fnB[256];
     
-//  delete jFile;
-//  if (jFile) jFile->Close();  
-/*
-  TFile* f = new TFile("j.root", "recreate");
-  e0H->Write();
-  e1H->Write();
-  e2H->Write();
-  e3H->Write();
-  e4H->Write();
-  e7H->Write();
-  e8H->Write();
-  f->Close();
-*/
-
-    // Get Control Plots
-//  gStyle->SetOptStat(0);
+    sprintf(fnB, "%s/%s.root", fBkgdDirectory, (runDataB->GetRunTitle(runN)).Data());
+    jFileB = new TFile(fnB);
+
+    char nameB[100];
+    sprintf(nameB, "TreeJ%d",eventN);
+    TTree *jetB =(TTree *)(jFileB->Get(nameB));
+    jetB->SetBranchAddress("FoundJet",    &fRecB);
+    jetB->GetEntry(0);
     
-    TCanvas* c1 = new TCanvas("c1");
-  e0H->Draw();
-  e1H->SetLineColor(2);
-  e2H->SetLineColor(4);
-  e3H->SetLineColor(5);
-  e1H->Draw("same");
-  e3H->Draw("same");
-
-
-  TCanvas* c2 = new TCanvas("c2");
-//  dr1H->Draw();
-  dr2H->SetLineColor(2);
-  dr2H->Draw("");
-  
-  TCanvas* c3 = new TCanvas("c3");
-  dr2H->Draw();
-  dr3H->Draw("same");
+    TArrayI inJetB = fRecB->GetInJet();
+    TArrayF etainB = fRecB->GetEtaIn();
+    TArrayF ptinB = fRecB->GetPtIn();
+    TArrayF phiinB = fRecB->GetPhiIn();
+    fPtJ = fRecJ->GetPt(0);
+    fEJ = fRecJ->GetE(0);
+    fEtaJ = fRecJ->GetEta(0);
+    fPhiJ = fRecJ->GetPhi(0);
+    Float_t t1,ene,xi,detaB,dphiB,drB,jt,wB;
+    TVector3 trkB;
+    TVector3 jv3B;
+    jv3B.SetX(fjv3X); jv3B.SetY(fjv3Y); jv3B.SetZ(fjv3Z);
+    
+    for(Int_t k=0;k<inJetB.GetSize();k++){
+      if(ptinB[k] > fPartPtCut){
+       detaB = etainB[k] - fEtaJ;
+       dphiB = phiinB[k] - fPhiJ;
+       if (dphiB < -TMath::Pi()) dphiB= -dphiB - 2.0 * TMath::Pi();
+       if (dphiB > TMath::Pi()) dphiB = 2.0 * TMath::Pi() - dphiB;
+       drB = TMath::Sqrt(detaB * detaB + dphiB * dphiB);
+       t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etainB[k])));
+       ene = ptinB[k]*TMath::Sqrt(1.+1./(t1*t1));
+       trkB.SetPtEtaPhi(ptinB[k],etainB[k],phiinB[k]);
+       // --- dN/dxi
+       if (drB<fdrdNdxi) {
+         xi = (Float_t) TMath::Log(fEJ/ene);
+         fdNdxiB->Fill(xi,fPtJ,fWeight);
+       }
+       // --- Jt
+       if (drB<fdrJt) {
+         jt = TMath::Sin(trkB.Angle(jv3B))*(trkB.Mag());
+         fJtB->Fill(jt,fPtJ,fWeight);
+       }
+       // --- dE/dr
+       if (drB<fdrdEdr) {
+         wB = GetdEdrWeight(fEtaJ,drB)*fWeight*ene;
+         fdEdrB->Fill(drB,fPtJ,wB);
+       } 
+      }
+    }
+    delete jetB;
+    if (jFileB) jFileB->Close();
+    delete jFileB; 
+  }
+}
 
-  TCanvas* c4 = new TCanvas("c4");
-  e0H->Draw();
+void AliJetAnalysis::FillShapH(Float_t r)
+{
+  // Fill jet shape histograms
+  if (fDoRecJ && fRecJ->GetNJets()> 0) {
+    TArrayF p=fRecJ->GetPtFromSignal();
+    if (p[0]>fPercentage) {
+      Shape(fRecJ,fRShapSelH,fRShapRejH,fRShapAllH,fdEdrH,fPtEneH2,fdEdrW,r);
+      fWdEdr+=fWeight;
+      fWShapR+=fWeight;
+    }
+  }
+}
 
-  TCanvas* c5 = new TCanvas("c5");
-  etaH->SetXTitle("#eta_{rec} - #eta_{gen}");
-  
-  etaH->Draw();
-  eta1H->SetLineColor(2);
+void AliJetAnalysis::Shape(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha, 
+                          TH2F *hdedr, TH2F *hptene, TH1F* wdedr, Float_t radius)
+{
+  // initialize
+  TArrayI inJet = j->GetInJet();
+  TArrayF etain = j->GetEtaIn();
+  TArrayF ptin = j->GetPtIn();
+  TArrayF phiin = j->GetPhiIn();
   
-  eta1H->Draw("same");
+  // first compute dE/dr histo
+  Float_t etaj = j->GetEta(0);
+  Float_t ene,w,deta,dphi,dr,t1;
+  for(Int_t i=0;i<inJet.GetSize();i++) {
+    deta = etain[i] - j->GetEta(0);
+    dphi = phiin[i] - j->GetPhi(0);
+    if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+    if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+    dr = TMath::Sqrt(deta * deta + dphi * dphi);
+    if ((dr<fdrdEdr) && (ptin[i] > fPartPtCut)) {
+      t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+      ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
+      w = GetdEdrWeight(etaj,dr)*fWeight*ene;
+      hdedr->Fill(dr,j->GetPt(0),w);
+    }
+  }
+  hptene->Fill(fRecJ->GetE(0),fRecJ->GetPt(0),fWeight);
+  wdedr->Fill(j->GetPt(0),fWeight);
   
-  TCanvas* c5a = new TCanvas("c5a");
-  eta1H->Draw();
-
-  TCanvas* c5b = new TCanvas("c5b");
-  eta2H->Draw();
-
-  TCanvas* c6 = new TCanvas("c6");
-  e4H->Draw();
-  TCanvas* c7 = new TCanvas("c7");
-  phiH->Draw();
-
-  TCanvas* c7a = new TCanvas("c7a");
-  phi1H->Draw();
-  TCanvas* c7b = new TCanvas("c7b");
-  phi2H->Draw();
-
-  TCanvas* c8 = new TCanvas("c8");
-  e5H->SetXTitle("E_{gen} (GeV)");
-  e5H->Draw();
-  e6H->SetLineColor(2);
-  e6H->Draw("same");
+  // now compute shape histos
+  Int_t nBins = ha->GetNbinsX();
+  Float_t xMin = ha->GetXaxis()->GetXmin();
+  Float_t xMax = ha->GetXaxis()->GetXmax();
+  Float_t binSize,halfBin;
+  binSize = (xMax-xMin)/nBins;
+  halfBin = binSize/2;
+  Float_t rptS[20], rptR[20], rptA[20];
+  for(Int_t i=0;i<nBins;i++) rptS[i]=rptR[i]=rptA[i]=0.0;
+  // fill bins in r for leading jet
+  for(Int_t i=0;i<inJet.GetSize();i++) {
+    if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
+      deta = etain[i] - j->GetEta(0);
+      dphi = phiin[i] - j->GetPhi(0);
+      if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+      if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+      dr = TMath::Sqrt(deta * deta + dphi * dphi);
+      Float_t rR = dr/radius;
+      Int_t bin = (Int_t) TMath::Floor(rR/binSize);
+      rptA[bin]+=ptin[i]/(j->GetPt(0));
+      if (inJet[i] == 1) rptS[bin]+=ptin[i]/(j->GetPt(0));
+      if (fPythia && inJet[i] == -1) 
+       rptR[bin]+=ptin[i]/(j->GetPt(0));
+    }
+  }
   
-  TCanvas* c9 = new TCanvas("c9");
-  e6H->Draw();
+  // compute shape and fill histogram
+  Float_t ptS,ptR,ptA,r;
+  ptS=ptR=ptA=0.0;
+  for (Int_t i=0;i<nBins;i++) {
+    ptS+=rptS[i];
+    ptR+=rptR[i];
+    ptA+=rptA[i];
+    r=(i+1)*binSize-halfBin;
+    hs->Fill(r,ptS*fWeight);
+    if(fPythia) {
+      hr->Fill(r,ptR*fWeight);
+      ha->Fill(r,ptA*fWeight);
+    }
+  }
+}
+
+void AliJetAnalysis::FillFragH()
+{
+  // Fill fragmentation histogram
+  if (fDoRecJ && fRecJ->GetNJets()> 0) {
+    TArrayF p=fRecJ->GetPtFromSignal();
+    if (p[0]>fPercentage) {
+      FragFun(fRecJ,fRFragSelH,fRFragRejH,fRFragAllH);
+      fWFragR+=fWeight;
+    }
+  }
+}
 
+void AliJetAnalysis::FragFun(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha)
+{
+  // Calculate de fragmentation function
+  TArrayI inJet = j->GetInJet();
+  TArrayF etain = j->GetEtaIn();
+  TArrayF ptin = j->GetPtIn();
   
+  Float_t xi,t1,ene;
   
-  TCanvas* c10 = new TCanvas("c10");
+  for(Int_t i=0;i<(inJet.GetSize());i++) {
+    if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
+      t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+      ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
+      xi = (Float_t) TMath::Log((j->GetE(0))/ene);
+      if (fPythia) ha->Fill(xi,fWeight);
+      if (inJet[i] == 1) hs->Fill(xi,fWeight);
+      if (fPythia && inJet[i] == -1) hr->Fill(xi,fWeight);
+    }
+  }
+}
 
+void AliJetAnalysis:: FillTrigH()
+{
+  // Fill trigger bias histograms
+  if(fDoTrig && fGenJ->GetNJets()>0 && fRecJ->GetNJets()>0){   
+    TArrayF p=fRecJ->GetPtFromSignal();
+    if (p[0]>fPercentage) {
+      float genEne = fGenJ->GetE(0);
+      float recEne = fRecJ->GetE(0);
+      float eneRatio = (recEne/fEfactor)/genEne;
+      
+      fGTriggerEneH->Fill(genEne, eneRatio, fWeight);
+      fRTriggerEneH->Fill(recEne/fEfactor, eneRatio, fWeight);
+      
+      if (fPart->LeadingFound()){
+       float leaEne = fPart->GetE();
+       float eneRatio2 = leaEne/genEne;
+       fGPTriggerEneH->Fill(genEne, eneRatio2, fWeight);
+       fPTriggerEneH->Fill(leaEne/0.2, eneRatio2, fWeight);
+      }
+    }
+  }
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// Normalize histogrames 
+
+void AliJetAnalysis::NormHistograms()
+{
+  // normalize shape histograms
+  if (fDoShap) {
+    if (fDoRecJ && fWShapR>0.) {  // leading reconstructed jet
+      fRShapSelH->Scale(1.0/fWShapR);
+      fRShapRejH->Scale(1.0/fWShapR);
+      fRShapAllH->Scale(1.0/fWShapR);
+    }
+  }
   
-  r5H->SetMaximum(1.);
-  r5H->Draw();
-  r5H->SetXTitle("E_{gen} (GeV)");
-  r5H->SetYTitle("E_{leading} / E_{gen}");
-  r6H->SetLineColor(2);
-  r6H->Draw("same");
-
-  TCanvas* c11 = new TCanvas("c11");
-//
-// Leading particle rec/gen
-//
-  r7H->SetMaximum(1.);
-  
-  r7H->Draw();
-  r7H->SetXTitle("E_{gen} (GeV)");
-  r7H->SetYTitle("E_{leading} / E_{gen}");
+  // normalize fragmentation function histograms
+  if (fDoFrag) {
+    if (fDoRecJ && fWFragR>0.) {  // leading reconstructed jet
+      fRFragSelH->Scale(2.0/fWFragR);
+      fRFragRejH->Scale(2.0/fWFragR);
+      fRFragAllH->Scale(2.0/fWFragR);
+    }
+  }
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetAnalysis::PlotHistograms()
+{
+  // Plot histogramas (to be done...)
+  if (fDoKine) PlotKineH();
+  if (fDoCorr) PlotCorrH();
+  if (fDoCorr50) PlotCorr50H();
+  if (fDoShap) PlotShapH();
+  if (fDoFrag) PlotFragH();
+  if (fDoTrig) PlotTrigH();
+}
+void AliJetAnalysis::PlotKineH() const
+{
+    // missing    
+    if (fDoPart) ;
+    if (fDoGenJ) ;
+    if (fDoRecJ) ;
+}
+
+void AliJetAnalysis::PlotCorrH() const
+{
+    // missing    
+    if (fDoPart && fDoGenJ) ;
+    if (fDoPart && fDoRecJ) ; 
+    if (fDoGenJ && fDoRecJ) ; 
+}
+void AliJetAnalysis::PlotCorr50H() const
+{
+    // missing    
+    if (fDoPart && fDoGenJ) ;
+    if (fDoPart && fDoRecJ) ; 
+    if (fDoGenJ && fDoRecJ) ; 
+}
+
+void AliJetAnalysis::PlotShapH() const
+{
+    // missing    
+    if (fDoGenJ) ;
+    if (fDoRecJ) ;
+}
+
+void AliJetAnalysis::PlotFragH() const
+{
+    // missing    
+    if (fDoGenJ) ;
+    if (fDoRecJ) ;
+}
+
+void AliJetAnalysis::PlotTrigH()
+{
+    // missing    
+
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetAnalysis::SaveHistograms()
+{
+  // Save histograms
+  TFile *fOut = new TFile(fFile,"recreate");
+  fOut->cd();
+  if (fDoKine) SaveKineH();
+  if (fDoCorr) SaveCorrH();
+  if (fDoCorr50) SaveCorr50H();
+  if (fDoShap) SaveShapH();
+  if (fDoFrag) SaveFragH();
+  if (fDoTrig) SaveTrigH();
+  if (fDoJt) SaveJtH();
+  if (fDodNdxi) SavedNdxiH();
+  fOut->Close();
+}
+
+void AliJetAnalysis::SaveKineH()
+{
+  // Save kinematic histograms
+  if (fDoPart) {
+    fPKineEneH->Write();
+    fPKinePtH->Write();
+    fPKinePhiH->Write();
+    fPKineEtaH->Write();
+  }
   
-  r8H->SetLineColor(2);
-  r8H->Draw("same");
+  if (fDoGenJ) {
+    fGKineEneH->Write();
+    fGKinePtH->Write();
+    fGKinePhiH->Write();
+    fGKineEtaH->Write();
+  }
   
-  TCanvas* c12 = new TCanvas("c12");
-  drP1->SetXTitle("E_{gen} (GeV)");
-  drP1->SetYTitle("#eta_{rec} - #eta_{gen}");
-  drP1->Draw();
-
-  TCanvas* c13 = new TCanvas("c13");
-  drP2->SetXTitle("E_{gen} (GeV)");
-  drP2->SetYTitle("#eta_{leading} - #eta_{gen}");
+  if (fDoRecJ) {
+    fRKineEneH->Write();
+    fRKinePtH->Write();
+    fRKinePhiH->Write();
+    fRKineEtaH->Write();
+  }
+}
+
+void AliJetAnalysis::SaveCorrH()
+{
+  // Save correlation histograms
+  if (fDoPart && fDoGenJ) {
+    fPGCorrEneH->Write();
+    fPGCorrPtH->Write();
+    fPGCorrEtaH->Write();
+    fPGCorrPhiH->Write();
+  }
   
-  drP2->Draw();
+  if (fDoPart && fDoRecJ) {
+    fPRCorrEneH->Write();
+    fPRCorrPtH->Write();
+    fPRCorrEtaH->Write();
+    fPRCorrPhiH->Write();
+  }
   
-  TCanvas* c14 = new TCanvas("c14");
-  dphiH->Draw();
-
-/*
-  e1H->Write();
-  e2H->Write();
-  e3H->Write();
-  e4H->Write();
-*/
-} 
+  if (fDoGenJ && fDoRecJ) {
+    fRGCorrEneH->Write();
+    fRGCorrPtH->Write();
+    fRGCorrEtaH->Write();
+    fRGCorrPhiH->Write();
+  } 
+}
+
+void AliJetAnalysis::SaveCorr50H()
+{
+  // Save correlation histograms (special case)
+  if (fDoPart && fDoRecJ) {
+    fPRCorr50EneH->Write();
+    fPRCorr50PtH->Write();
+    fPRCorr50EtaH->Write();
+    fPRCorr50PhiH->Write();
+  }
+  if (fDoGenJ && fDoRecJ) {
+    fRGCorr50EneH->Write();
+    fRGCorr50PtH->Write();
+    fRGCorr50EtaH->Write();
+    fRGCorr50PhiH->Write();
+  } 
+}
+
+void AliJetAnalysis::SaveShapH()
+{
+  // Save jet shape histograms
+  if (fDoRecJ) {
+    fRShapSelH->Write();
+    fdEdrH->Write();
+    if(fDoBkgd) fdEdrB->Write();
+    fPtEneH2->Write();
+    fdEdrW->Write();
+    if (fPythia){
+      fRShapRejH->Write();
+      fRShapAllH->Write();  
+    }
+  }    
+}
+
+void AliJetAnalysis::SaveJtH()
+{
+  // Save J_T histograms
+  if (fDoRecJ) {
+    fJtH->Write();
+    if(fDoBkgd) fJtB->Write();
+    fJtW->Write();
+  }
+}
+
+void AliJetAnalysis::SavedNdxiH()
+{
+  // Save dN/d#xi histograms
+  if (fDoRecJ) {
+    fdNdxiH->Write();
+    if(fDoBkgd) fdNdxiB->Write();
+    fPtEneH->Write();
+    fdNdxiW->Write();
+  }
+}
+
+void AliJetAnalysis::SaveFragH()
+{
+  // Save fragmentation histograms
+  if (fDoRecJ) {
+    fRFragSelH->Write();
+    if(fPythia){
+      fRFragRejH->Write();
+      fRFragAllH->Write();  
+    }
+  }
+}
+
+void AliJetAnalysis::SaveTrigH()
+{
+  // Save trigger bias histograms
+  if(fDoTrig){
+    fGTriggerEneH->Write();
+    fRTriggerEneH->Write();
+    fGPTriggerEneH->Write();
+    fPTriggerEneH->Write();
+  }
+}
+
+////////////////////////////////////////////////////////////////////////
+// main Analysis function
+
+void AliJetAnalysis::Analyze()
+    
+{
+    // Kinematics for
+    //   leading particle
+    //   leading generated jet
+    //   leading reconstructed jet
+
+    // Correlations amd resolutions
+    //    a) correlations in energy, pt, phi, eta
+    //    b) resolutions in energy, pt, phi, eta, r 
+    //   leading particle and leading generated jet
+    //   leading particle and leading reconstructed jet
+    //   leading generated jet and leading reconstructed jet
+
+    // Fragmentation functions and Shapes
+    //    a) integrated over all pt
+    //    b) in 3 flavors:
+    //       b.1) only for user selected particles in jet
+    //       b.2) only for user rejected particles in jet
+    //       b.3) for all particles in jet
+
+    DefineHistograms();
+    FillHistograms();
+    NormHistograms();
+    PlotHistograms();
+    SaveHistograms();
+}
 
index cba1109..9c94dd4 100644 (file)
  
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
+
 //---------------------------------------------------------------------
 // JetAnalysis class 
-// Perform Jet Analysis
-// Author: amdreas.morsch@cern.ch
+// Perform Jet Analysis on already found jets
+// Author: andreas.morsch@cern.ch, jgcn@mail.cern.ch
+//         mercedes.lopez.noriega@cern.ch
 //---------------------------------------------------------------------
+
 #include <TObject.h> 
+class AliLeading;
+class AliJet;
+class TH1;
+class TH1F;
+class TH2F;
+class TProfile;
+class TLorentzVector;
+
 class AliJetAnalysis : public TObject
 {
  public:
  
     AliJetAnalysis();
-    ~AliJetAnalysis(){;}
+    virtual ~AliJetAnalysis();
 
-    void Analyse();
-    // Setter
-    void SetDirectory(char* directory) {fDirectory = directory;}
-    void SetEventRange(Int_t imin, Int_t imax) {fEventMin = imin; fEventMax = imax;}
-    void SetRunRange(Int_t imin, Int_t imax) {fRunMin = imin; fRunMax = imax;}
+    void Analyze();
+    // define histograms
+    void DefineHistograms();
+    void DefineKineH();
+    void DefineCorrH();
+    void DefineCorr50H();
+    void DefineShapH();
+    void DefineFragH();
+    void DefineTrigH();
+    void DefineJtH();
+    void DefinedNdxiH();
+    // fill histograms
+    void FillHistograms();
+    void FillKineH();
+    void FillCorrH();
+    void FillCorr50H();
+    void FillShapH(Float_t r);
+    void FillFragH();
+    void FillTrigH();
+    void FillJtH();
+    void FilldNdxiH();
+    void FillBkgd(Int_t eventN, Int_t runN);
+    // normalize histograms
+    void NormHistograms();
+    // plot histograms
+    void PlotHistograms();
+    void PlotKineH() const;
+    void PlotCorrH() const;
+    void PlotCorr50H() const;
+    void PlotShapH() const;
+    void PlotFragH() const;
+    void PlotTrigH();
+    // save histograms
+    void SaveHistograms();
+    void SaveKineH();
+    void SaveCorrH();
+    void SaveCorr50H();
+    void SaveShapH();
+    void SaveFragH();
+    void SaveTrigH();
+    void SaveJtH();
+    void SavedNdxiH();
+    // other functions
+    void Shape(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha, TH2F* hd, TH2F* hp, TH1F* wd, Float_t r);
+    void FragFun(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha);
+    void Correlation(TLorentzVector *lv1,TLorentzVector *lv2,TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4);
+    void Correlation50(AliJet *j,TLorentzVector *lv1,TLorentzVector *lv2,TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4);
+    // setters
+    void SetDirectory(char* directory) 
+      {fDirectory = directory;}                      // directory where file with jets is
+    void SetBkgdDirectory(char* directory) 
+      {fBkgdDirectory = directory;}                  // directory where file with background is
+    void SetOutputFile(char* file) {fFile = file;}   // file where plots will be saved
+    void SetPercentage(Float_t p) {fPercentage = p;} // minimum percentage of tracks coming from pythia (very aprox.)
+    void SetEventRange(Int_t imin, Int_t imax) 
+      {fEventMin = imin; fEventMax = imax;}           // first and last event
+    void SetRunRange(Int_t imin, Int_t imax) 
+      {fRunMin = imin; fRunMax = imax;}               // first and last run
+    void SetMinimumMult(Int_t m){fminMult = m;}       // minimum multiplicity cut  
+    void SetPythia(Bool_t f = kFALSE){fPythia = f;}   // If only pythia, to save everything...  
+    void SetDoJt(Bool_t f = kTRUE){fDoJt = f;}        // To get j_T distribution
+    void SetDodNdxi(Bool_t f = kTRUE){fDodNdxi = f;}  // To get #xi distribution
+    void SetDoBkgd(Bool_t f = kTRUE) {fDoBkgd = f;}   // To get the bkgd for j_T, xi and dEdr in a hijing event
+    void SetDoLeadPart(Bool_t f = kTRUE){fDoPart = f;}// To make plots for leading particle
+    void SetDoGenJet(Bool_t f = kTRUE){fDoGenJ = f;}  // To make plots for generated jets
+    void SetDoRecJet(Bool_t f = kTRUE){fDoRecJ = f;}  // To make plots for reconstructed jets
+    void SetDoKinematics(Bool_t f = kTRUE){fDoKine = f;}       // To make the kine plots
+    void SetDoCorrelations(Bool_t f = kTRUE){fDoCorr = f;}     // Correlation histograms 
+    void SetDoCorr50(Bool_t f = kFALSE){fDoCorr50 = f;}        // Correlation histograms when one particle has more than 50% E
+    void SetDoShape(Bool_t f = kTRUE){fDoShap = f;}            // Shape plots
+    void SetDoFragmentations(Bool_t f = kTRUE){fDoFrag = f;}   // Fragmentation
+    void SetDoTriggerBias(Bool_t f = kTRUE){fDoTrig = f;}      // Trigger bias plots
+    void SetDivideEnergy(Float_t Efactor){fEfactor = Efactor;} // Divides E of rec.jet by Efactor
+    void SetProperties(TH1* h,const char* x, const char* y) const;
+    void SetReaderHeader(char *s="AliJetKineReaderHeader"){fReaderHeader = s;}
+    void SetdEdrWeight();
+    void SetPartPtCut(Float_t c){fPartPtCut = c;}
+    void SetdrJt(Float_t r){fdrJt = r;}
+    void SetdrdNdxi(Float_t r){fdrdNdxi = r;}
+    void SetdrdEdr(Float_t r){fdrdEdr = r;}
+    // getters
+    Float_t GetdEdrWeight(Float_t eta, Float_t r);
+    
  private:
-    char*  fDirectory;   // Directory
-    Int_t  fEventMin;    // Minimum event number
-    Int_t  fEventMax;    // Maximum event number
-    Int_t  fRunMin;      // Minimum run number 
-    Int_t  fRunMax;      // Maximum run number
+    char*  fReaderHeader;    // Reader header
+    char*  fDirectory;       // Directory
+    char*  fBkgdDirectory;   // Directory for background
+    char*  fFile;            // Output file name
+    Int_t  fEventMin;        // Minimum event number
+    Int_t  fEventMax;        // Maximum event number
+    Int_t  fRunMin;          // Minimum run number 
+    Int_t  fRunMax;          // Maximum run number
+    Int_t  fminMult;         // Minimum multiplicity for events
+    Float_t fPercentage;     // percentage of pt from signal particles to accept a jet
+    Float_t fPartPtCut;      // cut in the pt of particles in histos
+    Float_t fdrJt;          // maximum dr for Jt plot
+    Float_t fdrdNdxi;       // maximum dr for dN/dxi plot
+    Float_t fdrdEdr;        // maximum dr for dE/dr plot
+    Float_t fEfactor;        // factor by which energy the reconstructed jet will be divided
+
+    Float_t fp0;    // percentage of tracks in reconstructed jet coming from pythia
+                    // so far calculated in aprox. way, it needs to be improved!
+    // for background from hijing events:
+    Float_t fPtJ;     // P_T of the pythia jet
+    Float_t fEJ;      // Energy of the pythia jet
+    Float_t fEtaJ;    // Eta of the pythia jet
+    Float_t fPhiJ;    // Phi of the pythia jet
+    Float_t fjv3X, fjv3Y, fjv3Z;     // x,y,z of the pythia jet
+
+    // user options    
+    Bool_t fPythia;      // if pythia events
+    Bool_t fDoPart;      // do analysis of leading particle
+    Bool_t fDoGenJ;      // do analysis of leading generated jet
+    Bool_t fDoRecJ;      // do analysis of leading rec jet
+    Bool_t fDoKine;      // do kinematic plots
+    Bool_t fDoCorr;      // do correlation plots
+    Bool_t fDoCorr50;    // do correlation plots when one track more than 50% of jet energy
+    Bool_t fDoShap;      // do shape plots
+    Bool_t fDoFrag;      // do fragmentation plots 
+    Bool_t fDoTrig;      // do trigger bias plots
+    Bool_t fDoJt;        // do jt histo
+    Bool_t fDodNdxi;     // do dN/dxi histo
+    Bool_t fDoBkgd;      // get dN/dxi bkgd using hijing tracks only
     
-       
+    // weights
+    Float_t fWeight;             // event weight
+    Float_t fWShapR;             // weighted number of jets 
+    Float_t fWFragR;             // weighted number of jets 
+    Float_t fWeightdEdr[10][20]; // weight for acceptance of dE/dr histo
+    Float_t fWdEdr;              // weighted number of events for dEdr histo
+    Float_t fWJt;                // weight for Jt
+    Float_t fWdNdxi;             // weight fro dNd#xi
+    
+    // leading hets and particles
+    AliLeading* fPart;   // pointer to leading particle
+    AliJet* fGenJ;       // pointer to leading generated jet
+    AliJet* fRecJ;       // pointer to leading reconstructed jet
+    AliJet* fRecB;       // pointer to leading reconstructed jet for background
+
+    // kine histos
+    TH1F *fRKineEneH;  // Reconstructed energy histo
+    TH1F *fRKinePtH;   // Reconstructed Pt histo
+    TH1F *fRKinePhiH;  // Reconstructed phi histo
+    TH1F *fRKineEtaH;  // Reconstructed eta histo
+    TH1F *fGKineEneH;  // Generated energy histo
+    TH1F *fGKinePtH;   // Generated Pt histo
+    TH1F *fGKinePhiH;  // Generated phi histo
+    TH1F *fGKineEtaH;  // Generated eta histo
+    TH1F *fPKineEneH;  // Pythia energy histo
+    TH1F *fPKinePtH;   // Pythia Pt histo
+    TH1F *fPKinePhiH;  // Pythia phi histo
+    TH1F *fPKineEtaH;  // Pythia eta histo
+
+    // correlation histograms
+    TH2F *fPGCorrEneH;  // Energy correlation part-gen jet
+    TH2F *fPGCorrPtH;   // Pt correlation part-gen jet
+    TH2F *fPGCorrEtaH;  // Pseudorapidity correlation part-gen jet
+    TH2F *fPGCorrPhiH;  // Azimuthal angle correlation part-gen jet
+    TH2F *fPRCorrEneH;  // Energy correlation part-rec jet
+    TH2F *fPRCorrPtH;   // Pt correlation part-rec jet
+    TH2F *fPRCorrEtaH;  // Pseudorapidity correlation part-rec jet
+    TH2F *fPRCorrPhiH;  // Azimuthal angle correlation part-rec jet
+    TH2F *fRGCorrEneH;  // Energy correlation rec jet-gen jet
+    TH2F *fRGCorrPtH;   // Pt correlation rec jet-gen jet
+    TH2F *fRGCorrEtaH;  // Pseudorapidity correlation rec jet-gen jet
+    TH2F *fRGCorrPhiH;  // Azimuthal angle correlation rec jet-gen jet
+   
+    // correlation histograms
+    TH2F *fPRCorr50EneH;  // Energy correlation part-rec jet
+    TH2F *fPRCorr50PtH;   // Pt correlation part-rec jet
+    TH2F *fPRCorr50EtaH;  // Pseudorapidity correlation part-rec jet
+    TH2F *fPRCorr50PhiH;  // Azimuthal angle correlation part-rec jet
+    TH2F *fRGCorr50EneH;  // Energy correlation rec jet-gen jet
+    TH2F *fRGCorr50PtH;   // Pt correlation rec jet-gen jet
+    TH2F *fRGCorr50EtaH;  // Pseudorapidity correlation rec jet-gen jet
+    TH2F *fRGCorr50PhiH;  // Azimuthal angle correlation rec jet-gen jet
+   
+    // fragmentation function and shape histos
+    TH1F *fRFragSelH;  // Frag Fun of reconstructed jets (sel part)
+    TH1F *fRFragRejH;  // Frag Fun of reconstructed jets (rej part)
+    TH1F *fRFragAllH;  // Frag Fun of reconstructed jets (all part)
+    TH1F *fRShapSelH;  // Shape of generated jets (sel part)
+    TH1F *fRShapRejH;  // Shape of generated jets (rej part)
+    TH1F *fRShapAllH;  // Shape of generated jets (all part)
+    
+    // trigger bias histos 
+    TProfile *fGTriggerEneH;  // Generated energy (trigger bias)
+    TProfile *fRTriggerEneH;  // Reconstructed energy (trigger bias)
+    TProfile *fGPTriggerEneH; // Generated energy (trigger bias)
+    TProfile *fPTriggerEneH;  // Leading particle energy (trigger bias)
+
+    // dE/dr histo
+    TH2F* fdEdrH;  // dE/dr histo
+    TH2F* fdEdrB;  // dE/dr bkgdhisto
+    TH2F* fPtEneH2;// fPtEneH2
+    TH1F* fdEdrW;  // weights for dE/dr
+
+    // Jt histo
+    TH2F* fJtH;  // J_{T} histogram
+    TH2F* fJtB;  // J_{T} bkgd histogram
+    TH1F* fJtW;  // J_{T} weight
+
+    // dN/dxi histo
+    TH2F* fdNdxiH;  // dN/d#xi histo
+    TH2F* fdNdxiB;  // dN/d#xi bkgd histo
+    TH1F* fdNdxiW;  // dN/d#xi weight histo
+    TH2F* fPtEneH;  // fPtEneH
+
     ClassDef(AliJetAnalysis,1)
 };
  
index 6e37ef5..a9c1030 100755 (executable)
@@ -22,9 +22,6 @@
 
 #include <TStyle.h>
 #include <TCanvas.h>
-#include <TLorentzVector.h>
-#include <TFile.h> 
-#include <TClonesArray.h>
 #include <TH1I.h>
 #include <TH1D.h>
 #include "AliJetReader.h"
@@ -37,17 +34,20 @@ ClassImp(AliJetControlPlots)
 
 AliJetControlPlots::AliJetControlPlots()
 {
-  //
   // Constructor
-  //
   fNJetT=0;
 
+  // general properties
   fNJetsH = new TH1I("fNJetsH","Number of Jets",12,0,11);
   SetProperties(fNJetsH,"Number of jets","Entries");
 
   fMultH = new TH1I("fMultH","Multiplicity of Jets",22,0,21);
   SetProperties(fMultH,"Multiplicity of jets","Entries");
 
+  fInJetH = new  TH1D("fInJetH","Percentage of particles in jets",50,0,1);
+  SetProperties(fInJetH,"percentage of particles in jets","Entries");
+
+  // kinematics
   fPtH = new TH1D("fPtH","Pt of Jets",50,0.,200.);
   SetProperties(fPtH,"P_{t} [GeV]","Entries");
 
@@ -57,85 +57,144 @@ AliJetControlPlots::AliJetControlPlots()
   fEneH = new TH1D("fEneH","Energy of Jets",50,0.,200.);
   SetProperties(fEneH,"Energy [GeV]","Entries");
 
-  fFragH = new TH1D("fFragH","Jet Fragmentation",20,0.,1.);
-  SetProperties(fFragH,"x=E_{i}/E_{JET}","1/N_{JET}dN_{ch}/dx");
-
-  fFragLnH = new TH1D("fFragLnH","Jet Fragmentation",20,0.,10);
-  SetProperties(fFragLnH,"#xi=ln(E_{JET}/E_{i}","1/N_{JET}dN_{ch}/d#xi");
-
   fPhiH = new TH1D("fPhiH","Azimuthal angle of Jets",
-                  60,-TMath::Pi(),TMath::Pi());
+                  60,0.,2.0*TMath::Pi());
   SetProperties(fPhiH,"#phi","Entries");
-  
-  fInJetH = new  TH1D("fInJetH","Percentage of particles in jets",
-                     50,0,1);
-  SetProperties(fInJetH,"percentage of particles in jets","Entries");
+
+  // fragmentation 
+  fFragH = new TH1D("fFragH","Leading Jet Fragmentation (selected tracks)",
+                   20,0.,1.);
+  SetProperties(fFragH,"x=E_{i}/E_{JET}","1/N_{JET}dN_{ch}/dx");
+
+  fFragrH = new TH1D("fFragrH","Leading Jet Fragmentation (rejected tracks)",
+                   20,0.,1.);
+  SetProperties(fFragrH,"x=E_{i}/E_{JET}","1/N_{JET}dN_{ch}/dx");
+
+  // fragmentation log
+  fFragLnH = new TH1D("fFragLnH","Leading Jet Fragmentation (selected tracks)",
+                     20,0.,10);
+  SetProperties(fFragLnH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+
+  fFragLnrH = new TH1D("fFragLnrH",
+                      "Leading Jet Fragmentation (rejected tracks)",
+                      20,0.,10);
+  SetProperties(fFragLnrH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+
+  // jet shape
+  fShapeH = new TH1D("fShapeH","Leading Jet Shape (selected tracks)",
+                    20,0.,1.);
+  SetProperties(fShapeH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R)");
+      
+  fShaperH = new TH1D("fShaperH","Leading Jet Shape (rejected tracks)",
+                    20,0.,1.);
+  SetProperties(fShaperH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R)");
 }
 
 ////////////////////////////////////////////////////////////////////////
 
 AliJetControlPlots::~AliJetControlPlots()
 {
-  //
   // Destructor
-  //
   delete fNJetsH;
   delete fMultH;
   delete fPtH;
   delete fEtaH;
   delete fEneH;
-  delete fFragH;
-  delete fFragLnH;
   delete fPhiH;
   delete fInJetH;
+  delete fFragH;
+  delete fFragLnH;
+  delete fFragrH;
+  delete fFragLnrH;
+  delete fShapeH;
+  delete fShaperH;
 }
 
 ////////////////////////////////////////////////////////////////////////
 
-void AliJetControlPlots::FillHistos(AliJet *j, AliJetReader *r)
+void AliJetControlPlots::FillHistos(AliJet *j)
 {
-// Fills the histograms
-
+  // Fills the histograms
   Int_t nj = j->GetNJets();
-  fNJetsH->Fill(nj);
-  fNJetT+=nj;
-
+  fNJetsH->Fill(nj,1);
+  if (nj == 0) return;
+  
   // kinematics, occupancy and multiplicities
   TArrayI mj = j->GetMultiplicities();
   Int_t mjTot=0;
   for (Int_t i=0;i<nj;i++) {
     mjTot+=mj[i];
-    fMultH->Fill(mj[i]);    
-    fPtH->Fill(j->GetPt(i));
-    fEtaH->Fill(j->GetEta(i));
-    fEneH->Fill(j->GetE(i));
-    fPhiH->Fill(j->GetPhi(i));
+    fMultH->Fill(mj[i],1);    
+    fPtH->Fill(j->GetPt(i),1);
+    fEtaH->Fill(j->GetEta(i),1);
+    fEneH->Fill(j->GetE(i),1);
+    fPhiH->Fill(j->GetPhi(i),1);
   }
-  if (nj>0) 
-    fInJetH->Fill(((Double_t) mjTot)/((Double_t) j->GetNinput()));
-
-  // fragmentation 
-  TClonesArray *lvArray = r->GetMomentumArray();
+  fInJetH->Fill(((Double_t) mjTot)/((Double_t) j->GetNinput()),1);
+  
+  // fragmentation of leading jet 
   TArrayI inJet = j->GetInJet();
-  Int_t nIn = j->GetNinput();
-  for(Int_t i=0;i<nIn;i++) {
-    if (inJet[i] == -1) continue;
-    TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
-    Double_t xe = (lv->E())/(j->GetE(inJet[i]));
-    fFragH->Fill(xe);
-    fFragLnH->Fill(TMath::Log(1.0/xe));
+  TArrayF etain = j->GetEtaIn();
+  TArrayF ptin = j->GetPtIn();
+  for(Int_t i=0;i<(inJet.GetSize());i++) {
+    Float_t t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+    Float_t ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
+    if (inJet[i] == 1) {
+      fFragH->Fill((Float_t) ene/(j->GetE(0)),1);
+      fFragLnH->Fill((Float_t) TMath::Log((j->GetE(0))/ene),1);
+    }
+    if (inJet[i] == -1) {
+      fFragrH->Fill((Float_t) ene/(j->GetE(0)),1);
+      fFragLnrH->Fill((Float_t) TMath::Log((j->GetE(0))/ene),1);
+    }
+  }
+  
+  // shape of leading jet 
+  // * CAREFUL: depends on number of bins and bin sizes
+  //   of shape histo. HARDWIRED at the moment *
+  //   Plot against dr NOT dr/R_jet!!!
+  TArrayF phiin = j->GetPhiIn();
+  Float_t rptS[20], rptR[20];
+  for(Int_t i=0;i<20;i++) rptS[i]=rptR[i]=0.0;
+  
+  for(Int_t i=0;i<inJet.GetSize();i++) {
+    if (inJet[i] == 1 || inJet[i] == -1) {
+      Float_t deta = etain[i] - j->GetEta(0);
+      Float_t dphi = phiin[i] - j->GetPhi(0);
+      if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+      if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+      Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+      if (dr>1) continue;
+      Int_t bin = (Int_t) TMath::Floor(dr/0.05);
+      if (inJet[i] == 1) rptS[bin]+=ptin[i]/(j->GetPt(0));
+      if (inJet[i] == -1) rptR[bin]+=ptin[i]/(j->GetPt(0));      
+    }
+  }
+  
+  Float_t ptS,ptR,r;
+  ptS=ptR=0.0;
+  for (Int_t i=0;i<20;i++) {
+    ptS+=rptS[i];
+    ptR+=rptR[i];
+    r=(i+1)*0.05-0.025;
+    fShapeH->Fill(r,ptS);
+    fShaperH->Fill(r,ptR);      
   }
+  fNJetT++;
 }
 
 ////////////////////////////////////////////////////////////////////////
 
 void AliJetControlPlots::Normalize()
 {
+// *CAREFUL: depends on histogram number of bins 
   if (fNJetT == 0) return;
-  fFragH->Sumw2();
   fFragH->Scale(20.0/((Double_t) fNJetT));
-  fFragLnH->Sumw2();
   fFragLnH->Scale(2.0/((Double_t) fNJetT));
+  fFragrH->Scale(20.0/((Double_t) fNJetT));
+  fFragLnrH->Scale(2.0/((Double_t) fNJetT));
+  fShapeH->Scale(1.0/((Double_t) fNJetT));
+  fShaperH->Scale(1.0/((Double_t) fNJetT));
 }
 
 ////////////////////////////////////////////////////////////////////////
@@ -146,8 +205,8 @@ void AliJetControlPlots::PlotHistos()
   gStyle->SetOptStat(kFALSE);
   gStyle->SetOptTitle(kFALSE);
   
-  TCanvas *c = new TCanvas("c","AliJetControlPlots",50,200,900,700);
-  c->Divide(3,3);
+  TCanvas *c = new TCanvas("c","AliJetControlPlots",50,200,1200,700);
+  c->Divide(4,3);
   c->cd(1);  fNJetsH->Draw("e1");
   c->cd(2);  fMultH->Draw("e1");
   c->cd(3);  fInJetH->Draw("e1");
@@ -155,21 +214,23 @@ void AliJetControlPlots::PlotHistos()
   c->cd(5);  fEtaH->Draw("e1");
   c->cd(6);  fPhiH->Draw("e1");
   c->cd(7);  fEneH->Draw("e1");
-  c->cd(8);  fFragH->Draw("e1");
-  c->cd(9);  fFragLnH->Draw("e1");
+  c->cd(8);  fFragLnH->Draw("e1");
+  c->cd(9);  fFragLnrH->Draw("e1");
+  c->cd(10);  fShapeH->Draw("e1");
+  c->cd(11);  fShaperH->Draw("e1");
 }
 
 ////////////////////////////////////////////////////////////////////////
 
 void AliJetControlPlots::SetProperties(TH1* h,const char* x, const char* y) const
 {
-//
 // Sets the histogram style properties
   h->SetMarkerStyle(20);
   h->SetMarkerSize(.5);
   h->SetMarkerColor(2);
   h->SetXTitle(x);
   h->SetYTitle(y);
+  h->Sumw2();
 }
 
 
index ff372db..d3b19de 100755 (executable)
@@ -31,31 +31,39 @@ class AliJetControlPlots : public TObject
   // getters
   TH1I *GetNJetsH() {return fNJetsH;}
   TH1I *GetMultH() {return fMultH;}
+  TH1D *GetPhiH() {return fPhiH;}
+  TH1D *GetFractionInJetH() {return fInJetH;}
   TH1D *GetEneH() {return fEneH;}
   TH1D *GetPtH() {return fPtH;}
   TH1D *GetEtaH() {return fEtaH;}
   TH1D *GetFragH() {return fFragH;}
   TH1D *GetFragLnH() {return fFragLnH;}
-  TH1D *GetPhiH() {return fPhiH;}
-  TH1D *GetFractionInJetH() {return fInJetH;}
+  TH1D *GetFragrH() {return fFragrH;}
+  TH1D *GetFragLnrH() {return fFragLnrH;}
+  TH1D *GetShapeH() {return fShapeH;}
+  TH1D *GetShaperH() {return fShaperH;}  
   
   // others
-  void FillHistos(AliJet *j,AliJetReader *r);
+  void FillHistos(AliJet *j);
   void PlotHistos();
   void SetProperties(TH1* h,const char* x, const char* y) const;
   void Normalize();
 
  protected:
-  TH1I *fNJetsH;    // distribution of number of jets
+  TH1I *fNJetsH;   // distribution of number of jets
   TH1I *fMultH;    // jet multiplicity
-  TH1D *fPtH;       // pt spectra
-  TH1D *fEtaH;      // eta distribution
-  TH1D *fEneH;      // energy distribution
-  TH1D *fFragH;      // jet fragmentation
-  TH1D *fFragLnH;      // jet fragmentation in ln scale
-  TH1D *fPhiH;      // phi distribution
-  TH1D *fInJetH;    // percentage of input particles in a jet
-  Int_t fNJetT;     // total number of jets for normalization
+  TH1D *fPtH;      // pt spectra
+  TH1D *fEtaH;     // eta distribution
+  TH1D *fEneH;     // energy distribution
+  TH1D *fFragH;    // leading jet fragmentation (selected part)
+  TH1D *fFragLnH;  // leading jet fragmentation in ln scale
+  TH1D *fFragrH;   // leading jet fragmentation (rejected part)
+  TH1D *fFragLnrH; // leading jet fragmentation in ln scale
+  TH1D *fShapeH;   // leading jet shape (selected part)
+  TH1D *fShaperH;  // leading jet shape (rejected part)  
+  TH1D *fPhiH;     // phi distribution
+  TH1D *fInJetH;   // percentage of input particles in a jet
+  Int_t fNJetT;    // total number of jets for normalization
 
   ClassDef(AliJetControlPlots,1)
 };
index 4a2ce34..a2988d2 100755 (executable)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
  
-//
 // Jet ESD Reader 
 // ESD reader for jet analysis
-// Author: Mercedes Lopez Noriega 
-// mercedes.lopez.noriega@cern.ch
-//
-
+// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
 
 #include <Riostream.h>
 #include <TSystem.h>
 #include <TLorentzVector.h>
+#include <TVector3.h>
 #include "AliJetESDReader.h"
 #include "AliJetESDReaderHeader.h"
 #include "AliESD.h"
@@ -34,7 +31,7 @@ ClassImp(AliJetESDReader)
 AliJetESDReader::AliJetESDReader()
 {
   // Constructor    
-
+  printf("\nIn reader constructor\n");
   fReaderHeader = 0x0;
   fMass = 0;
   fSign = 0;
@@ -45,14 +42,14 @@ AliJetESDReader::AliJetESDReader()
 AliJetESDReader::~AliJetESDReader()
 {
   // Destructor
-  //  delete fReaderHeader;
 }
 
 //____________________________________________________________________________
 
 void AliJetESDReader::OpenInputFiles()
-
 {
+  // Open input files
+  printf("\nOpening files\n");
   // chain for the ESDs
   fChain   = new TChain("esdTree");
   fChainMC = new TChain("mcStackTree");
@@ -69,8 +66,8 @@ void AliJetESDReader::OpenInputFiles()
       printf("Adding %s\n",name);
       char path[256];
       sprintf(path,"%s/%s",dirName,name);
-      fChain->AddFile(path,-1,"esdTree");
-      fChainMC->AddFile(path,-1,"mcStackTree");
+      fChain->AddFile(path,-1);
+      fChainMC->AddFile(path,-1);
      }
   }
 
@@ -95,50 +92,58 @@ void AliJetESDReader::OpenInputFiles()
 
 void AliJetESDReader::FillMomentumArray(Int_t event)
 {
-// Fills the momentum array
+  // Fill the momentum array for each track
   Int_t goodTrack = 0;
   Int_t nt = 0;
-  Float_t pt;
-  Double_t p[3]; // track 3 momentum vector
+  Float_t pt, eta;
+  TVector3 p3;
 
   // clear momentum array
   ClearArray();
+
   // get event from chain
   fChain->GetEntry(event);
   fChainMC->GetEntry(event);
+
   // get number of tracks in event (for the loop)
   nt = fESD->GetNumberOfTracks();
-  // tmporary storage of signal flags
-  Int_t* flag  = new Int_t[nt];
+ // temporary storage of signal and pt cut flag
+  Int_t* sflag  = new Int_t[nt];
+  Int_t* cflag  = new Int_t[nt];
 
   // get cuts set by user
-  Float_t ptMin = ((AliJetESDReaderHeader*) fReaderHeader)->GetPtCut();
+  Float_t ptMin = fReaderHeader->GetPtCut();
+  Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+  Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
 
   //loop over tracks
   for (Int_t it = 0; it < nt; it++) {
       AliESDtrack *track = fESD->GetTrack(it);
       UInt_t status = track->GetStatus();
-      if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
-
+      p3 = track->P3();
+      pt = p3.Pt();
+      if (((status & AliESDtrack::kITSrefit) == 0) ||
+         ((status & AliESDtrack::kTPCrefit) == 0)) continue;    // quality check
       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() 
-         && TMath::Abs(track->GetLabel()) > 10000)  continue;
-    
-      track->GetPxPyPz(p);
-      pt = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]); // pt of the track
-
-      if (pt < ptMin) continue; //check  cuts 
-
-      new ((*fMomentumArray)[goodTrack]) 
-         TLorentzVector(p[0], p[1], p[2],
-                        TMath::Sqrt(pt * pt +p[2] * p[2]));
-      flag[goodTrack]=0;
-      if (TMath::Abs(track->GetLabel()) < 10000) flag[goodTrack]=1;
+         && TMath::Abs(track->GetLabel()) > 10000)  continue;   // quality check
+      if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly() 
+         && TMath::Abs(track->GetLabel()) < 10000)  continue;   // quality check
+      eta = p3.Eta();
+      if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
+
+      new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
+      sflag[goodTrack]=0;
+      if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+      cflag[goodTrack]=0;
+      if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
       goodTrack++;
   }
   // set the signal flags
-  fSignalFlag.Set(goodTrack,flag);
+  fSignalFlag.Set(goodTrack,sflag);
+  fCutFlag.Set(goodTrack,cflag);
 
-  printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
+  //printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
 }
 
 
index 790c83a..e425f9d 100755 (executable)
@@ -4,9 +4,11 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
  
+//---------------------------------------------------------------------
 // Jet ESD Reader 
 // ESD reader for jet analysis
 // Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+//---------------------------------------------------------------------
 
 #include "AliJetReader.h"
 class AliJetESDReaderHeader;
@@ -19,8 +21,8 @@ class AliJetESDReader : public AliJetReader
   virtual ~AliJetESDReader();
 
   // Getters
-  Float_t GetTrackMass() const {return fMass;}     // returns mass of the track
-  Int_t   GetTrackSign() const {return fSign;}     // returns sign of the track
+  Float_t GetTrackMass() const {return fMass;}  // returns mass of the track
+  Int_t   GetTrackSign() const {return fSign;}  // returns sign of the track
 
   // Setters
   void FillMomentumArray(Int_t event); 
index 41fb5ae..bdd606d 100755 (executable)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-//
+
+//---------------------------------------------------------------------
 // Jet ESD Reader Header
 // Header for the ESD reader in the jet analysis
-// Author: Mercedes Lopez Noriega 
-// mercedes.lopez.noriega@cern.ch
-//
+// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+//---------------------------------------------------------------------
  
 #include "AliJetESDReaderHeader.h"
 
@@ -29,10 +28,8 @@ AliJetESDReaderHeader::AliJetESDReaderHeader():
  AliJetReaderHeader("AliJetESDReaderHeader") 
 {
   // Constructor
-  SetPtCut();
-  SetDCA();
-  SetTLength();
   SetReadSignalOnly(kFALSE);
+  SetReadBkgdOnly(kFALSE);
   
 }
 
index 56ac1e6..a0c36b8 100755 (executable)
@@ -18,25 +18,18 @@ class AliJetESDReaderHeader : public AliJetReaderHeader
   virtual ~AliJetESDReaderHeader();
   
   // Getters
-  Float_t GetPtCut()       const  {return fPtCut;}
-  Float_t GetDCA()         const  {return fDCA;} // not working so far..(always 0)
-  Float_t GetTLength()     const  {return fTLength;} // not working so far.. (always 0)
   Bool_t  ReadSignalOnly() const  {return fReadSignalOnly;}
-  
+  Bool_t  ReadBkgdOnly() const  {return fReadBkgdOnly;}
          
   // Setters
-  virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;}
-  virtual void SetDCA(Float_t dca = 3.0) {fDCA = dca;}
-  virtual void SetTLength(Float_t length = 0.0) {fTLength = length;}
   virtual void SetReadSignalOnly(Bool_t flag = kTRUE) {fReadSignalOnly = flag;}
+  virtual void SetReadBkgdOnly(Bool_t flag = kTRUE) {fReadBkgdOnly = flag;}
   
  protected:
   //parameters set by user
-  Float_t fPtCut;          // pt cut 
-  Float_t fDCA;            // dca cut
-  Float_t fTLength;        // track length cut
   Bool_t  fReadSignalOnly; // read particles from signal event only
-  ClassDef(AliJetESDReaderHeader,1);
+  Bool_t  fReadBkgdOnly;   // read particles from bkgd event only
+  ClassDef(AliJetESDReaderHeader,2);
 };
  
 #endif
index 7b866dd..37bfa56 100644 (file)
@@ -12,8 +12,7 @@
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-  
+   
 //---------------------------------------------------------------------
 // Jet finder base class
 // manages the search for jets 
 #include "AliLeading.h"
 #include "AliHeader.h"
 
-
 ClassImp(AliJetFinder)
 
-////////////////////////////////////////////////////////////////////////
 
 AliJetFinder::AliJetFinder()
 {
-  //
   // Constructor
-  //
   fOut     = 0x0;
   fJets    = new AliJet();
   fGenJets = new AliJet();
@@ -50,15 +45,11 @@ AliJetFinder::AliJetFinder()
   SetPlotMode(kFALSE);
 }
 
-
 ////////////////////////////////////////////////////////////////////////
 
 AliJetFinder::~AliJetFinder()
 {
-  //
   // destructor
-  //
-
   // here reset and delete jets
   fJets->ClearJets();
   delete fJets;
@@ -72,10 +63,8 @@ AliJetFinder::~AliJetFinder()
   delete fOut;
   // reset and delete control plots
   if (fPlots) delete fPlots;
-  // delete fLeading;
 }
 
-
 ////////////////////////////////////////////////////////////////////////
 
 void AliJetFinder::SetOutputFile(const char *name)
@@ -84,25 +73,22 @@ void AliJetFinder::SetOutputFile(const char *name)
   fOut = new TFile(name,"recreate");
 }
 
-
 ////////////////////////////////////////////////////////////////////////
 
 void AliJetFinder::PrintJets()
 {
-//
-// Print jet information
+  // Print jet information
   cout << " Jets found with jet algorithm:" << endl;
   fJets->PrintJets();
   cout << " Jets found by pythia:" << endl;
   fGenJets->PrintJets();
 }
 
-
 ////////////////////////////////////////////////////////////////////////
 
 void AliJetFinder::SetPlotMode(Bool_t b)
 {
-// Sets the plotting mode
+  // Sets the plotting mode
   fPlotMode=b;
   if (b && !fPlots) fPlots = new AliJetControlPlots(); 
 }
@@ -111,7 +97,7 @@ void AliJetFinder::SetPlotMode(Bool_t b)
 
 void AliJetFinder::WriteJetsToFile(Int_t i)
 {
-// Writes the jets to file
+  // Writes the jets to file
   fOut->cd();
   char hname[30];
   sprintf(hname,"TreeJ%d",i);
@@ -134,12 +120,11 @@ void AliJetFinder::WriteRHeaderToFile()
   rh->Write();
 }
 
-
 ////////////////////////////////////////////////////////////////////////
 
 void AliJetFinder::GetGenJets()
 {
-// Get the generated jet information from mc header
+  // Get the generated jet information from mc header
   AliHeader* alih = fReader->GetAliHeader(); 
   if (alih == 0) return;
   AliGenEventHeader * genh = alih->GenEventHeader();
@@ -175,7 +160,6 @@ void AliJetFinder::Run()
       WriteRHeaderToFile();
       WriteJHeaderToFile();
   }
-
   // loop over events
   Int_t nFirst,nLast;
   nFirst = fReader->GetReaderHeader()->GetFirstEvent();
@@ -187,7 +171,7 @@ void AliJetFinder::Run()
       GetGenJets();
       FindJets();
       if (fOut) WriteJetsToFile(i);
-      if (fPlots) fPlots->FillHistos(fJets,fReader);
+      if (fPlots) fPlots->FillHistos(fJets);
       fLeading->Reset();
       fGenJets->ClearJets();
       Reset();
@@ -203,4 +187,3 @@ void AliJetFinder::Run()
       fOut->Close();
   }
 }
-
index efefbcf..c3feee4 100755 (executable)
@@ -4,7 +4,6 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
  
 //---------------------------------------------------------------------
 // Jet finder base class
 // manages the search for jets 
@@ -19,7 +18,6 @@ class AliJetReader;
 class AliJetControlPlots;
 class AliLeading;
 
-
 class AliJetFinder : public TObject 
 {
  public:
@@ -28,9 +26,9 @@ class AliJetFinder : public TObject
   virtual ~AliJetFinder();
 
   // getters
-  virtual AliJet *GetJets()      {return fJets;}
-  virtual Bool_t GetPlotMode()   const {return fPlotMode;}
-  virtual TFile* GetOutputFile() { return fOut; }
+  virtual AliJet *GetJets() {return fJets;}
+  virtual Bool_t GetPlotMode() const {return fPlotMode;}
+  virtual TFile* GetOutputFile() {return fOut;}
   // setters
   virtual void SetPlotMode(Bool_t b);
   virtual void SetOutputFile(const char *name="jets.root");
@@ -51,7 +49,7 @@ class AliJetFinder : public TObject
  protected:
   Bool_t fPlotMode;              // do you want control plots?
   AliJet* fJets;                 // pointer to jet class
-  AliJet* fGenJets;               // pointer to generated jets
+  AliJet* fGenJets;              // pointer to generated jets
   AliLeading* fLeading;          // pointer to leading particle data 
   AliJetReader* fReader;         // pointer to reader
   AliJetControlPlots* fPlots;    // pointer to control plots
index 0c39bef..bd50406 100755 (executable)
@@ -32,6 +32,8 @@ AliJetHeader::AliJetHeader():
   //
   // Constructor
   //
+    fJetEtaMax = 0.5;
+    fJetEtaMin = -0.5;    
 }
  
 ////////////////////////////////////////////////////////////////////////
@@ -43,6 +45,8 @@ AliJetHeader::AliJetHeader(const char * name):
   //
   // Constructor
   //
+    fJetEtaMax = 0.5;
+    fJetEtaMin = -0.5;    
 }
 
 ////////////////////////////////////////////////////////////////////////
index d85de2f..55e45a5 100755 (executable)
@@ -23,16 +23,23 @@ class AliJetHeader : public TNamed
   virtual ~AliJetHeader() { }
 
   // Getters
-  virtual TString GetComment() {return fComment;} 
-   
+  virtual TString GetComment() const {return fComment;} 
+  virtual Float_t GetJetEtaMax() const {return fJetEtaMax;}
+  virtual Float_t GetJetEtaMin() const {return fJetEtaMin;}  
+  
   // Setters
-  virtual void SetComment(const char* com) {fComment=TString(com);} 
+  virtual void SetComment(const char* com) {fComment=TString(com);}
+  virtual void SetJetEtaMax(Float_t eta= 0.5) {fJetEtaMax = eta;}
+  virtual void SetJetEtaMin(Float_t eta= -0.5) {fJetEtaMin = eta;}  
+  
 
   // others
   
 protected:
   TString fComment; // a comment 
-
+  Float_t fJetEtaMax; // maximum eta for the jet
+  Float_t fJetEtaMin; // minimum eta for the jet
+  
   ClassDef(AliJetHeader,1)
 };
  
index 2b96d2c..363650c 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-// 
+//-------------------------------------------------------------------------
 // Jet kinematics reader 
 // MC reader for jet analysis
-// Author: Andreas Morsch 
-// andreas.morsch@cern.ch
-//
+// Author: Andreas Morsch (andreas.morsch@cern.ch)
+//-------------------------------------------------------------------------
 
 // From root ...
 #include <TClonesArray.h>
 #include <TPDGCode.h>
 #include <TParticle.h>
 #include <TParticlePDG.h>
-#include <TVector3.h>
+//#include <TVector3.h>
 #include <TLorentzVector.h>
 #include <TSystem.h>
 #include <TDatabasePDG.h>
@@ -59,166 +58,180 @@ AliJetKineReader::~AliJetKineReader()
 
 void AliJetKineReader::OpenInputFiles()
 {
-    // Opens the input file using the run loader
-    const char* dirName = fReaderHeader->GetDirectory();
-    char path[256];
-    sprintf(path, "%s/galice.root",dirName);
-    fRunLoader = AliRunLoader::Open(path);
-    fRunLoader->LoadKinematics();
-    fRunLoader->LoadHeader(); 
-    
-    Int_t nMax = fRunLoader->GetNumberOfEvents();
-    printf("\nTotal number of events = %d", nMax);
-    
+  // Opens the input file using the run loader
+  const char* dirName = fReaderHeader->GetDirectory();
+  char path[256];
+  sprintf(path, "%s/galice.root",dirName);
+  fRunLoader = AliRunLoader::Open(path);
+  fRunLoader->LoadKinematics();
+  fRunLoader->LoadHeader(); 
+  
+  // get number of events
+  Int_t nMax = fRunLoader->GetNumberOfEvents();
+  printf("\nTotal number of events = %d", nMax);
+  
   // set number of events in header
-    if (fReaderHeader->GetLastEvent() == -1)
-       fReaderHeader->SetLastEvent(nMax);
-    else {
-       Int_t nUsr = fReaderHeader->GetLastEvent();
-       fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
-    }
+  if (fReaderHeader->GetLastEvent() == -1)
+    fReaderHeader->SetLastEvent(nMax);
+  else {
+    Int_t nUsr = fReaderHeader->GetLastEvent();
+    fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
+  }
 }
 
 //____________________________________________________________________________
 
 void AliJetKineReader::FillMomentumArray(Int_t event)
 {
-//
-// Fill momentum array for event
-//
-    Int_t goodTrack = 0;
-    // Clear array
-    ClearArray();
-    // Get event from runloader
-    fRunLoader->GetEvent(event);
-    // Get the stack
-    AliStack* stack = fRunLoader->Stack();
-    // Number of primaries
-    Int_t nt = stack->GetNprimary();
-    // Get cuts set by user and header
-    Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
-    fAliHeader = fRunLoader->GetHeader();
-        
-    // Loop over particles
-    Int_t* flag  = new Int_t[nt];
-    for (Int_t it = 0; it < nt; it++) {
-       TParticle *part = stack->Particle(it);
-       Int_t   status  = part->GetStatusCode();
-       Int_t   pdg     = TMath::Abs(part->GetPdgCode());
-       Float_t pt      = part->Pt(); 
-       
-       // Skip non-final state particles, neutrinos and particles with pt < pt_min 
-       
-       if (
-           (status != 1)            
-           || (pdg == 12 || pdg == 14 || pdg == 16) 
-           || (pt < ptMin)             
-           ) continue; 
-
-
-       Float_t p       = part->P();
-       Float_t p0      = p;
-       Float_t eta     = part->Eta();
-       Float_t phi     = part->Phi();
-
-        // Fast simulation of TPC if requested
-       if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
-           // Charged particles only
-           Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
-           if (charge == 0)               continue;
-           // Simulate efficiency
-           if (!Efficiency(p0, eta, phi)) continue;
-           // Simulate resolution
-           p = SmearMomentum(4, p0);
-       } // Fast TPC
-       
-        // Fast simulation of EMCAL if requested
-       if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
-           Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
-           // Charged particles
-           if (charge != 0) {
-               // Simulate efficiency
-               if (!Efficiency(p0, eta, phi)) continue;
-               // Simulate resolution
-               p = SmearMomentum(4, p0);
-           } // charged
-           // Neutral particles
-           // Exclude K0L, n, nbar
-           if (pdg == kNeutron || pdg == kK0Long) continue;
-       } // Fast EMCAL
-       
-       fMass = part->GetCalcMass();
-       fPdgC = part->GetPdgCode();
-       // Fill momentum array
-       Float_t r  = p/p0;
-//             printf("smeared %13.3f %13.3f %13.3f\n", p0, p, r);
-       
-       Float_t px = r * part->Px(); 
-       Float_t py = r * part->Py(); 
-       Float_t pz = r * part->Pz();
-       Float_t m  = part->GetMass();
-       Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
-       
-       new ((*fMomentumArray)[goodTrack]) 
-           TLorentzVector(px, py, pz, e);
-       goodTrack++;
+  // Fill momentum array for event
+  Int_t goodTrack = 0;
+  // Clear array
+  ClearArray();
+  // Get event from runloader
+  fRunLoader->GetEvent(event);
+  // Get the stack
+  AliStack* stack = fRunLoader->Stack();
+  // Number of primaries
+  Int_t nt = stack->GetNprimary();
+
+  // Get cuts set by user and header
+  Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
+  Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+  Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
+  fAliHeader = fRunLoader->GetHeader();
+  
+  // temporary storage of signal and cut flags
+  Int_t* sflag  = new Int_t[nt];
+  Int_t* cflag  = new Int_t[nt];
+  
+  TLorentzVector p4;
+  // Loop over particles
+  for (Int_t it = 0; it < nt; it++) {
+    TParticle *part = stack->Particle(it);
+    Int_t   status  = part->GetStatusCode();
+    Int_t   pdg     = TMath::Abs(part->GetPdgCode());
+    Float_t pt      = part->Pt(); 
+    
+    // Skip non-final state particles, neutrinos 
+    // and particles with pt < pt_min 
+    if ((status != 1)            
+       || (pdg == 12 || pdg == 14 || pdg == 16)) continue; 
+    
+    Float_t p       = part->P();
+    Float_t p0      = p;
+    Float_t eta     = part->Eta();
+    Float_t phi     = part->Phi();
+    
+    // Fast simulation of TPC if requested
+    if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
+      // Charged particles only
+      Float_t charge = 
+       TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+      if (charge == 0)               continue;
+      // Simulate efficiency
+      if (!Efficiency(p0, eta, phi)) continue;
+      // Simulate resolution
+      p = SmearMomentum(4, p0);
+    } // End fast TPC
+    
+    // Fast simulation of EMCAL if requested
+    if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
+      Float_t charge = 
+       TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+      // Charged particles only
+      if (charge != 0){
+       // Simulate efficiency
+       if (!Efficiency(p0, eta, phi)) continue;
+       // Simulate resolution
+       p = SmearMomentum(4, p0);
+      } // end "if" charged particles
+      // Neutral particles (exclude K0L, n, nbar)
+      if (pdg == kNeutron || pdg == kK0Long) continue;
+    } // End fast EMCAL
+    
+    fMass = part->GetCalcMass();
+    fPdgC = part->GetPdgCode();
+
+    // Fill momentum array
+    Float_t r  = p/p0;
+    Float_t px = r * part->Px(); 
+    Float_t py = r * part->Py(); 
+    Float_t pz = r * part->Pz();
+    Float_t m  = part->GetMass();
+    Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
+    p4 = TLorentzVector(px, py, pz, e);
+    if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
+    new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
+    
+    // all tracks are signal ... ?
+    sflag[goodTrack]=1;
+    cflag[goodTrack]=0;
+    if (pt > ptMin) cflag[goodTrack]=1; // track surviving pt cut
+    goodTrack++;
   }
-    fSignalFlag.Set(goodTrack,flag);
-    printf("\nNumber of good tracks in event %d= %d \n",event,goodTrack);
+
+  // set the signal flags
+  fSignalFlag.Set(goodTrack,sflag);
+  fCutFlag.Set(goodTrack,cflag);
+  printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
+
 }
 
 
 Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
 {
-//
-//  The relative momentum error, i.e. 
-//  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
-//  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
-//  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
-//
-//  If we include information from TRD b will be by a factor 2/3 smaller.
-//
-//  ind = 1: high multiplicity
-//  ind = 2: low  multiplicity
-//  ind = 3: high multiplicity + TRD
-//  ind = 4: low  multiplicity + TRD
-
-    Float_t pSmeared;
-    Float_t a = 0.75;
-    Float_t b = 0.12;
-
-    if (ind == 1) b = 0.12;
-    if (ind == 2) b = 0.08;
-    if (ind == 3) b = 0.12;    
-    if (ind == 4) b = 0.08;    
-    
-    Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
-    pSmeared = p + gRandom->Gaus(0., sigma);
-    return pSmeared;
+  //  The relative momentum error, i.e. 
+  //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
+  //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
+  //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
+  //
+  //  If we include information from TRD, b will be a factor 2/3 smaller.
+  //
+  //  ind = 1: high multiplicity
+  //  ind = 2: low  multiplicity
+  //  ind = 3: high multiplicity + TRD
+  //  ind = 4: low  multiplicity + TRD
+  
+  Float_t pSmeared;
+  Float_t a = 0.75;
+  Float_t b = 0.12;
+  
+  if (ind == 1) b = 0.12;
+  if (ind == 2) b = 0.08;
+  if (ind == 3) b = 0.12;    
+  if (ind == 4) b = 0.08;    
+  
+  Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
+  pSmeared = p + gRandom->Gaus(0., sigma);
+  return pSmeared;
 }
+
+
 Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
 {
-//
-// Fast simulation of geometrical acceptance and tracking efficiency
-//
-//  Tracking
-    Float_t eff = 0.99;
-    if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
-// Geometry
-    if (p > 0.5) {
-       phi *= 180. / TMath::Pi();
-       // Sector number 0 - 17
-       Int_t isec  = Int_t(phi / 20.);
-       // Sector centre
-       Float_t phi0 = isec * 20. + 10.;
-       Float_t phir = TMath::Abs(phi-phi0);
-       // 2 deg of dead space
-       if (phir > 9.) eff = 0.;
-    }
-    if (gRandom->Rndm() > eff) {
-       return kFALSE;
-    } else {
-       return kTRUE;
-    }
+  // Fast simulation of geometrical acceptance and tracking efficiency
+  //  Tracking
+  Float_t eff = 0.99;
+  if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
+  // Geometry
+  if (p > 0.5) {
+    phi *= 180. / TMath::Pi();
+    // Sector number 0 - 17
+    Int_t isec  = Int_t(phi / 20.);
+    // Sector centre
+    Float_t phi0 = isec * 20. + 10.;
+    Float_t phir = TMath::Abs(phi-phi0);
+    // 2 deg of dead space
+    if (phir > 9.) eff = 0.;
+  }
+
+  if (gRandom->Rndm() > eff) {
+    return kFALSE;
+  } else {
+    return kTRUE;
+  }
+
 }
 
index ee26dbb..c507aca 100644 (file)
@@ -21,7 +21,6 @@ class AliJetKineReader : public AliJetReader
   // Getters
   Float_t GetParticleMass() const {return fMass;}        // returns mass of the Track
   Int_t   GetParticlePdgCode() const {return fPdgC;}     // returns Pdg code
-
   // Setters
   void FillMomentumArray(Int_t event);
   void OpenInputFiles();
index 6d2682f..e338e78 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
  
+//---------------------------------------------------------------------
 // Jet Kine Reader Header
 // Header for the MC kinematics reader in the jet analysis
 // Author: Andreas Morsch (andreas.morsch@cern.ch)
+//---------------------------------------------------------------------
 
  
 #include "AliJetKineReaderHeader.h"
@@ -28,8 +30,8 @@ AliJetKineReaderHeader::AliJetKineReaderHeader():
  AliJetReaderHeader("AliJetKineReaderHeader") 
 {
   // Constructor
-  SetPtCut();
   SetFastSimTPC(kFALSE);
+  SetFastSimEMCAL(kFALSE);
 }
 
 //____________________________________________________________________________
index 49ef2e1..e190359 100644 (file)
@@ -18,19 +18,19 @@ class AliJetKineReaderHeader : public AliJetReaderHeader
   virtual ~AliJetKineReaderHeader();
   
   // Setters
-  void SetFastSimTPC(Bool_t flag = kTRUE) {fFastSimTPC = flag;}
-  void SetFastSimEMCAL(Bool_t flag = kTRUE) {fFastSimEMCAL = flag;}
-  virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;}
+  void SetFastSimTPC(Bool_t flag = kTRUE) {fFastSimTPC = flag;} // if TPC fast simulation
+  void SetFastSimEMCAL(Bool_t flag = kTRUE) {fFastSimEMCAL = flag;} // if EMCAL fast simulation
+
   // Getter
   Bool_t  FastSimTPC() const  {return fFastSimTPC;}
   Bool_t  FastSimEMCAL() const  {return fFastSimEMCAL;}
-  Float_t GetPtCut()   const  {return fPtCut;}
+
          
  protected:
   //parameters set by user
   Bool_t fFastSimTPC;
   Bool_t fFastSimEMCAL;
-  Float_t fPtCut;
+
   ClassDef(AliJetKineReaderHeader,1);
 };
  
index ee58a84..eb5f24d 100644 (file)
@@ -23,7 +23,7 @@ class AliJetProductionData : public TObject
     void    GetPtHardLimits(Int_t bin, Float_t& ptmin, Float_t& ptmax);
     TString GetRunTitle(Int_t bin);
     Float_t GetWeight(Int_t bin);
- public:         
+
     Int_t     fNbins;         // Number of pt_hard bins used in the production
     Float_t*  fPtHardLimits;  //[fNbins+1]
     Float_t*  fWeights;       //[fNbins]
index 3835a8a..3ea3e46 100755 (executable)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
  
-//---------------------------------------------------------------------
+//------------------------------------------------------------------------
 // Jet reader base class
 // manages the reading of input for jet algorithms
 // Author: jgcn@mda.cinvestav.mx
-//---------------------------------------------------------------------
+//------------------------------------------------------------------------
 
+// root
 #include <TClonesArray.h>
-
+//AliRoot
 #include "AliJetReader.h"
 #include "AliJetReaderHeader.h"
 #include "AliESD.h"
@@ -40,6 +41,7 @@ AliJetReader::AliJetReader()
   fArrayMC = 0;
   fAliHeader = 0;
   fSignalFlag = TArrayI();
+  fCutFlag = TArrayI();
 }
 
 ////////////////////////////////////////////////////////////////////////
index cb8979d..1c71365 100755 (executable)
@@ -4,11 +4,9 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
  
-//---------------------------------------------------------------------
 // Jet reader base class
 // manages the reading of input for jet algorithms
 // Author: jgcn@mda.cinvestav.mx
-//---------------------------------------------------------------------
   
 #include <TObject.h>
 #include <TChain.h>
@@ -19,7 +17,6 @@ class AliJetReaderHeader;
 class AliESD;
 class AliHeader;
 
-
 class AliJetReader : public TObject 
 {
  public: 
@@ -30,8 +27,9 @@ class AliJetReader : public TObject
   virtual TClonesArray *GetMomentumArray() {return fMomentumArray;}
   virtual Int_t GetChainEntries() {return fChain->GetEntries();} 
   virtual AliJetReaderHeader* GetReaderHeader() { return fReaderHeader;}
-  virtual AliHeader* GetAliHeader() { return fAliHeader;}
+  virtual AliHeader* GetAliHeader() {return fAliHeader;}
   virtual Int_t GetSignalFlag(Int_t i) const {return fSignalFlag[i];}
+  virtual Int_t GetCutFlag(Int_t i) const {return fCutFlag[i];}
 
   // Setters
   virtual void FillMomentumArray(Int_t) {}
@@ -51,7 +49,8 @@ class AliJetReader : public TObject
   AliJetReaderHeader      *fReaderHeader;  // pointer to header
   AliHeader               *fAliHeader;     // pointer to event header
   TArrayI fSignalFlag;   // to flag if a particle comes from pythia or 
-                        // from the underlying event
+                         // from the underlying event
+  TArrayI fCutFlag;      // to flag if a particle passed the pt cut or not
 
   ClassDef(AliJetReader,1)
 };
index cf39a0a..269c67a 100755 (executable)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
  
-//---------------------------------------------------------------------
+//-------------------------------------------------------------------------
+//-------------------------------------------------------------------------
 // base class for Jet Reader Header 
 // Author: jgcn@mda.cinvestav.mx
-//---------------------------------------------------------------------
+//-------------------------------------------------------------------------
 
 #include "AliJetReaderHeader.h"
 
 ClassImp(AliJetReaderHeader)
@@ -29,13 +29,12 @@ AliJetReaderHeader::AliJetReaderHeader():
  TNamed("AliJetReaderHeader", "Jet Reader Header"),
  fComment("No comment"),fDir(""),fPattern("")
 {
-  //
   // Constructor
-  //
   fFirst = 0;
   fLast = -1;
-  fFiducialEtaMin = -0.5;
-  fFiducialEtaMax =  0.5;
+  fFiducialEtaMin = -0.9;
+  fFiducialEtaMax =  0.9;
+  SetPtCut();
 }
 
 ////////////////////////////////////////////////////////////////////////
@@ -44,13 +43,12 @@ AliJetReaderHeader::AliJetReaderHeader(const char * name):
  TNamed(name, "Jet Reader Header"),
  fComment("No comment"),fDir(""),fPattern("")
 {
-  //
   // Constructor
-  //
   fFirst = 0;
   fLast = -1;
-  fFiducialEtaMin = -0.5;
-  fFiducialEtaMax =  0.5;
+  fFiducialEtaMin = -0.9;
+  fFiducialEtaMax =  0.9;
+  SetPtCut();
 }
 
 ////////////////////////////////////////////////////////////////////////
@@ -58,7 +56,5 @@ AliJetReaderHeader::AliJetReaderHeader(const char * name):
 AliJetReaderHeader::~AliJetReaderHeader()
 {
   // destructor
-  //  delete fComment;
-  //delete fDir;
-  //delete fPattern;
+
 }
index fd9a780..7466b42 100755 (executable)
@@ -4,10 +4,10 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
  
-//---------------------------------------------------------------------
+//-------------------------------------------------------------------------
 // base class for Jet Reader Header 
 // Author: jgcn@mda.cinvestav.mx
-//---------------------------------------------------------------------
+//-------------------------------------------------------------------------
   
 #include <TNamed.h>
 #include <TString.h>
@@ -21,12 +21,14 @@ class AliJetReaderHeader : public TNamed
   virtual ~AliJetReaderHeader();
   
   // Getters
-  virtual TString GetComment() {return fComment;}
+  virtual TString GetComment() const {return fComment;}
   virtual const char* GetDirectory() {return fDir.Data();}
   virtual const char* GetPattern() {return fPattern.Data();}
   virtual Float_t     GetFiducialEtaMin() const {return fFiducialEtaMin;}
   virtual Float_t     GetFiducialEtaMax() const {return fFiducialEtaMax;}  
-  
+  virtual Float_t GetPtCut()       const  {return fPtCut;}
+  Float_t GetDCA()         const  {return fDCA;}       // not working so far..(always 0)
+  Float_t GetTLength()     const  {return fTLength;}   // not working so far.. (always 0)
   Int_t   GetNEvents() const {return fLast-fFirst;}
   Int_t   GetLastEvent() const {return fLast;}
   Int_t   GetFirstEvent() const {return fFirst;}
@@ -38,18 +40,25 @@ class AliJetReaderHeader : public TNamed
   virtual void SetFirstEvent(Int_t i=0) {fFirst=i;}
   virtual void SetLastEvent(Int_t i=-1) {fLast=i;}
   virtual void SetFiducialEta(Float_t etamin, Float_t etamax) 
-      { fFiducialEtaMin = etamin; fFiducialEtaMax = etamax;}  
+      { fFiducialEtaMin = etamin; fFiducialEtaMax = etamax;}
+  virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;}
+  virtual void SetDCA(Float_t dca = 0.0) {fDCA = dca;}
+  virtual void SetTLength(Float_t length = 0.0) {fTLength = length;}
  protected:
 
   Int_t fFirst;            // First and last events analyzed
   Int_t fLast;             // in current set of files
   Float_t fFiducialEtaMin; // Fiducial minimum eta
-  Float_t fFiducialEtaMax; // Fiducial maximum eta 
+  Float_t fFiducialEtaMax; // Fiducial maximum eta
+  Float_t fPtCut;          // pt cut
+  Float_t fDCA;            // dca cut
+  Float_t fTLength;        // track length cut
   TString fComment;        // a comment
   TString fDir;            // directory with input files
   TString fPattern;        // pattern to look for input files
   
-  ClassDef(AliJetReaderHeader,1);
+  ClassDef(AliJetReaderHeader,2);
 };
  
 #endif
index c0d2e89..f7f07ec 100644 (file)
@@ -34,9 +34,7 @@ ClassImp(AliLeading)
 
 AliLeading::AliLeading() 
 {
-  //
   // Constructor
-  //
   fNassoc =  0;
   fLeading = new TLorentzVector(0.,0.,0.,0.);
   fLow     = -TMath::Pi()/2.0;
@@ -49,66 +47,62 @@ AliLeading::AliLeading()
 
 AliLeading::~AliLeading()
 {
-  //
   // Destructor
-  //
   delete fLeading;
 }
 
 ////////////////////////////////////////////////////////////////////////
 
 void AliLeading::FindLeading(AliJetReader *reader)
-
 {
-  //
   // find leading particle in the array of lorentz vectors
   // lvArray and fill the correlation histogram
-  //
-
-    
+  
   AliJetReaderHeader* header = reader->GetReaderHeader();
-    
+  
   TClonesArray* lvArray = reader->GetMomentumArray();
   Int_t nIn = lvArray->GetEntries();
-  fNassoc = nIn-1;
-
-  if (fNassoc < 0) return;
-
+  
   // find max
   Double_t ptMax = 0.0;
   Int_t idxMax = -1;
   for (Int_t i = 0; i < nIn; i++){
-      TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
-      if (lv->Pt()   > ptMax                       && 
-         lv->Eta()  > header->GetFiducialEtaMin() &&
-         lv->Eta()  < header->GetFiducialEtaMax()) 
-      {
-         ptMax  = lv->Pt();
-         idxMax = i;
-      }
+    TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+    if ((reader->GetCutFlag(i) == 1)             &&
+       lv->Pt()   > ptMax                       && 
+       lv->Eta()  > header->GetFiducialEtaMin() &&
+       lv->Eta()  < header->GetFiducialEtaMax()){
+      ptMax  = lv->Pt();
+      idxMax = i;
+    }
   }
   
   if (idxMax == -1) {
-      fFound = kFALSE;
-      Reset();
-      return;
+    fFound = kFALSE;
+    Reset();
+    return;
   }
   
   // fill correlation array
   fLeading = (TLorentzVector*) lvArray->At(idxMax);
   fFound = kTRUE;
   
+  fNassoc = 0;  
   for (Int_t i = 0; i < nIn; i++) {
-      if (i == idxMax) continue;
-      TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+    if (i == idxMax) continue;
+    TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+    if ( (reader->GetCutFlag(i) == 1) &&
+        lv->Eta()  > header->GetFiducialEtaMin() &&
+        lv->Eta()  < header->GetFiducialEtaMax()) {
       Double_t dphi = fLeading->DeltaPhi(*lv);
       if (dphi < fLow) dphi = 2.0 * TMath::Pi() + dphi;
       // find bin and fill array
-      
       Int_t iBin = (Int_t) 
-         TMath::Floor((dphi - fLow)
-                      *((Double_t) fnBin) / (2.0 * TMath::Pi()));
+       TMath::Floor((dphi - fLow)
+                    *((Double_t) fnBin) / (2.0 * TMath::Pi()));
       fCorr.AddAt(fCorr.At(iBin)+1,iBin);
+      fNassoc++;
+    }
   }
 }
 
@@ -126,7 +120,6 @@ void AliLeading::Reset()
 ////////////////////////////////////////////////////////////////////////
 
 void AliLeading::PrintLeading()
-
 {
 // Print leading particle information
   if (fNassoc<0) {
@@ -141,3 +134,35 @@ void AliLeading::PrintLeading()
        << fLeading->Eta() << "," << fLeading->Phi() << ")" << endl;
   cout << "    " << fNassoc << " associated particles." << endl;
 }
+
+////////////////////////////////////////////////////////////////////////
+Double_t AliLeading::GetE()
+{
+  return fLeading->E();
+}
+////////////////////////////////////////////////////////////////////////
+Double_t AliLeading::GetPt()
+{
+  return fLeading->Pt();
+}
+////////////////////////////////////////////////////////////////////////
+Double_t AliLeading::GetEta()
+{
+  return fLeading->Eta();
+}
+////////////////////////////////////////////////////////////////////////
+Double_t AliLeading::GetPhi()
+{
+  // get phi of leading
+  return ( (fLeading->Phi() < 0) ? 
+          (fLeading->Phi()) + 2. * TMath::Pi() : 
+          fLeading->Phi());
+}
+                                                                                
index b01c5b4..cef8512 100644 (file)
@@ -4,7 +4,6 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
  
 //---------------------------------------------------------------------
 // Class to find and store the leading particle in event and
 // store its correlation to associated particles
@@ -17,7 +16,6 @@
 class TLorentzVector;
 class AliJetReader;
 
-
 class AliLeading : public TObject
 {
  public:
@@ -26,11 +24,15 @@ class AliLeading : public TObject
   ~AliLeading();
 
   // Getters
-  Int_t GetNassoc() const { return fNassoc;}
-  TLorentzVector* GetLeading() const { return fLeading;}
-  TArrayI GetCorr() const { return fCorr;}
-  Bool_t LeadingFound() const { return fFound;}
-  
+  Int_t GetNassoc() const {return fNassoc;}
+  TLorentzVector* GetLeading() const {return fLeading;}
+  TArrayI GetCorr() const {return fCorr;}
+  Bool_t LeadingFound() const {return fFound;}
+  Double_t GetE();
+  Double_t GetPt();
+  Double_t GetEta();
+  Double_t GetPhi();
+   
   // Setters
   // Others
   void FindLeading(AliJetReader* reader);
@@ -39,13 +41,13 @@ class AliLeading : public TObject
 
  protected:
 
-  Int_t fNassoc; // number of associated particles
+  Int_t fNassoc;            // number of associated particles
   TLorentzVector* fLeading; // leading particle
-  TArrayI fCorr; // array to store azimuthal correlation
-                // between leading and assoc particles
-  Int_t fnBin;  // number of bins in array
-  Double_t fLow; // value corresponding to lower bound of bin 0
-  Bool_t  fFound;
+  TArrayI fCorr;            // array to store azimuthal correlation
+                            // between leading and assoc particles
+  Int_t fnBin;              // number of bins in array
+  Double_t fLow;            // value corresponding to lower bound of bin 0
+  Bool_t  fFound;           // leading found
   
   ClassDef(AliLeading,1);
 };
index becf007..4e92d8e 100755 (executable)
@@ -23,6 +23,7 @@
 
 #include <Riostream.h>
 #include <TLorentzVector.h>
+#include <TFile.h>
 #include "AliPxconeJetFinder.h"
 #include "AliPxconeJetHeader.h"
 #include "AliJetReader.h"
index eb99684..32fb1f7 100755 (executable)
@@ -29,9 +29,7 @@ ClassImp(AliPxconeJetHeader)
 AliPxconeJetHeader::AliPxconeJetHeader():
   AliJetHeader("AliPxconeJetHeader")
 {
-  //
   // Constructor
-  //
   SetMode();
   SetRadius();
   SetMinPt();
@@ -41,13 +39,10 @@ AliPxconeJetHeader::AliPxconeJetHeader():
 
 ////////////////////////////////////////////////////////////////////////
  
-void AliPxconeJetHeader::PrintParameters()
+void AliPxconeJetHeader::PrintParameters() const
 
 {
-  //
   // prints out parameters of jet algorithm
-  //
-
   cout << " PXCONE jet algorithm " << endl;
   cout << "  Running mode: " << fMode << endl;
   cout << "  Cone size: " << fRadius << endl;
index 007b0e6..1e1acb5 100755 (executable)
@@ -34,7 +34,7 @@ class AliPxconeJetHeader : public AliJetHeader
   void SetOverlap(Double_t o=.75) {fOverlap=o;}
 
   // others
-  void PrintParameters();
+  void PrintParameters() const;
    
 protected:
   Int_t fMode;           // ee or pp mode
index 2e291de..c9c1e95 100755 (executable)
@@ -12,8 +12,6 @@
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-   
-  
  
 //---------------------------------------------------------------------
 // UA1 Jet finder 
@@ -24,6 +22,7 @@
 
 #include <Riostream.h>
 #include <TLorentzVector.h>
+#include <TFile.h>
 #include <TH2F.h>
 #include "AliUA1JetFinder.h"
 #include "AliUA1JetHeader.h"
@@ -40,9 +39,7 @@ ClassImp(AliUA1JetFinder)
 AliUA1JetFinder::AliUA1JetFinder()
 
 {
-  //
   // Constructor
-  //
   fHeader = 0x0;
   fLego   = 0x0;
 }
@@ -52,13 +49,7 @@ AliUA1JetFinder::AliUA1JetFinder()
 AliUA1JetFinder::~AliUA1JetFinder()
 
 {
-  //
   // destructor
-  //
-
-  // delete fLego;            
-  
-  // reset and delete header
 }
 
 ////////////////////////////////////////////////////////////////////////
@@ -90,7 +81,6 @@ void AliUA1JetFinder::FindJets()
 {
   // initialize event, look for jets, download jet info
 
-
   // initialization in 2 steps
   // 1) transform input to pt,eta,phi plus lego
   // 2) dump lego
@@ -101,16 +91,19 @@ void AliUA1JetFinder::FindJets()
   if (nIn == 0) return;
     
   // local arrays for input
+  Float_t* enT  = new Float_t[nIn];
   Float_t* ptT  = new Float_t[nIn];
   Float_t* etaT = new Float_t[nIn];
   Float_t* phiT = new Float_t[nIn];
 
   // load input vectors
   for (Int_t i = 0; i < nIn; i++){
-      TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
-      ptT[i]  = lv->Pt();
-      etaT[i] = lv->Eta();
-      phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
+    TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+    enT[i]  = lv->Energy();
+    ptT[i]  = lv->Pt();
+    etaT[i] = lv->Eta();
+    phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
+    if (fReader->GetCutFlag(i) == 1) 
       fLego->Fill(etaT[i], phiT[i], ptT[i]);
   }
   fJets->SetNinput(nIn);
@@ -125,8 +118,8 @@ void AliUA1JetFinder::FindJets()
   TAxis* xaxis = fLego->GetXaxis();
   TAxis* yaxis = fLego->GetYaxis();
   Float_t e = 0.0;
-  for (Int_t i = 1; i <= fHeader->GetNbinEta(); i++) {
-      for (Int_t j = 1; j <= fHeader->GetNbinPhi(); j++) {
+  for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
+      for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
          e = fLego->GetBinContent(i,j);
          if (e > 0.0) e -= fHeader->GetMinCellEt();
          if (e < 0.0) e = 0.;
@@ -140,7 +133,7 @@ void AliUA1JetFinder::FindJets()
   }
 
   // run the algo. Parameters from header
-  Int_t nTot      = (fHeader->GetNbinEta())*(fHeader->GetNbinPhi());
+  Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
   Float_t minmove = fHeader->GetMinMove();
   Float_t maxmove = fHeader->GetMaxMove();
   Int_t mode      = fHeader->GetMode();
@@ -149,60 +142,83 @@ void AliUA1JetFinder::FindJets()
   ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell, 
                 minmove, maxmove, mode, precbg, ierror);
 
-
   // sort jets
   Int_t * idx  = new Int_t[UA1JETS.njet];
   TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
 
   // download jet info.   
-  AliJetReaderHeader* rheader = fReader->GetReaderHeader();
   for(Int_t i = 0; i < UA1JETS.njet; i++) {
     // reject events outside acceptable eta range
-    if (((UA1JETS.etaj[1][idx[i]])> (rheader->GetFiducialEtaMax()))
-       || ((UA1JETS.etaj[1][idx[i]]) < (rheader->GetFiducialEtaMin())))
-      continue;
-    Float_t px, py,pz,en; // convert to 4-vector
-    px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
-    py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
-    pz = UA1JETS.etj[idx[i]] /
-      TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
-    en = TMath::Sqrt(px * px + py * py + pz * pz);
-    fJets->AddJet(px, py, pz, en);
+    if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
+       || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
+         continue;
+      Float_t px, py,pz,en; // convert to 4-vector
+      px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
+      py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
+      pz = UA1JETS.etj[idx[i]] /
+         TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
+      en = TMath::Sqrt(px * px + py * py + pz * pz);
+      fJets->AddJet(px, py, pz, en);
   }
   
   // find multiplicities and relationship jet-particle
   // find also percentage of pt from pythia
   Int_t* injet = new Int_t[nIn];
-  for (Int_t i = 0; i < nIn; i++) injet[i]= -1;
+  Int_t* sflag = new Int_t[nIn];
+  for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
   Int_t* mult  = new Int_t[fJets->GetNJets()];
+  Int_t* ncell  = new Int_t[fJets->GetNJets()];
   Float_t* percentage  = new Float_t[fJets->GetNJets()];
 
   for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
-    Float_t pt_sig = 0.0;
-    mult[i] = 0;
-    for (Int_t j = 0; j < nIn; j++) {
-      Float_t deta = etaT[j] - fJets->GetEta(i);
-      Float_t dphi = phiT[j] - fJets->GetPhi(i);
-      if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
-      if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
-      Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
-      if (dr < fHeader->GetRadius() && injet[j] == -1) {
-       if (fReader->GetSignalFlag(j) == 1) pt_sig+=ptT[j];
-       injet[j] = i;
-       mult[i]++;
+      Float_t pt_sig = 0.0;
+      mult[i] = 0;
+      ncell[i] = UA1JETS.ncellj[i];
+      for (Int_t j = 0; j < nIn; j++) {
+         Float_t deta = etaT[j] - fJets->GetEta(i);
+         Float_t dphi = phiT[j] - fJets->GetPhi(i);
+         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+         if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+         Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+         if (dr < fHeader->GetRadius() && injet[j] == 0) {
+             injet[j] = -(i+1);
+             if ((fReader->GetCutFlag(j)) == 1 &&
+                 (etaT[j] < fHeader->GetLegoEtaMax()) &&
+                 (etaT[j] > fHeader->GetLegoEtaMin())) {
+                 injet[j] = i+1;
+                 mult[i]++;
+                 if (fReader->GetSignalFlag(j) == 1) {
+                   pt_sig+=ptT[j];
+                   sflag[j]=1;
+                 }
+             }
+         }
       }
-    }
-    percentage[i] = pt_sig/((Double_t) fJets->GetPt(i));    
+      percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
+       ((Double_t) fJets->GetPt(i));    
   }
+  fJets->SetNCells(ncell);
   fJets->SetPtFromSignal(percentage);
   fJets->SetMultiplicities(mult);
   fJets->SetInJet(injet);
+  fJets->SetEtaIn(etaT);
+  fJets->SetPhiIn(phiT);
+  fJets->SetPtIn(ptT);
+  fJets->SetEtAvg(UA1JETS.etavg);
+  delete ncell;
+  delete enT;
+  delete ptT;
+  delete etaT;
+  delete phiT;
+  delete injet;
+  delete idx;
+  delete mult;
+  delete percentage;
 }
 
 ////////////////////////////////////////////////////////////////////////
 
 void AliUA1JetFinder::Reset()
-
 {
   fLego->Reset();
   fJets->ClearJets();
@@ -216,17 +232,16 @@ void AliUA1JetFinder::WriteJHeaderToFile()
   fHeader->Write();
 }
 
-
 ////////////////////////////////////////////////////////////////////////
 
 void AliUA1JetFinder::Init()
 {
   // initializes some variables
   Float_t dEta, dPhi;
-  dEta = (fHeader->GetEtaMax()-fHeader->GetEtaMin())
-    /((Float_t) fHeader->GetNbinEta());
-  dPhi = (fHeader->GetPhiMax()-fHeader->GetPhiMin())
-    /((Float_t) fHeader->GetNbinPhi());
+  dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
+    /((Float_t) fHeader->GetLegoNbinEta());
+  dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
+    /((Float_t) fHeader->GetLegoNbinPhi());
 
   UA1CELL.etaCellSize = dEta;
   UA1CELL.phiCellSize = dPhi;
@@ -240,6 +255,7 @@ void AliUA1JetFinder::Init()
   // book lego
   fLego = new 
     TH2F("legoH","eta-phi",
-        fHeader->GetNbinEta(), fHeader->GetEtaMin(), fHeader->GetEtaMax(),
-        fHeader->GetNbinPhi(), fHeader->GetPhiMin(), fHeader->GetPhiMax());
+        fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(), 
+        fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(), 
+        fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
 }
index 03566ca..6d4aed8 100755 (executable)
 #include "AliUA1JetHeader.h"
 ClassImp(AliUA1JetHeader)
  
 ////////////////////////////////////////////////////////////////////////
 
 AliUA1JetHeader::AliUA1JetHeader():
   AliJetHeader("AliUA1JetHeader")
 {
-  //
   // Constructor
-  //
   fConeRadius =  0.3;
   fEtSeed     =  3.0;
   fMinJetEt   = 10.0;
@@ -41,23 +38,21 @@ AliUA1JetHeader::AliUA1JetHeader():
   fMinMove    =  0.05;
   fMaxMove    =  0.15;
   fPrecBg     =  0.035;
-  fNbinEta    =  36;
-  fNbinPhi    = 124;
-  fPhiMin     =   0.;
-  fPhiMax     = 2. * TMath::Pi();
-  fEtaMin     = -0.9;
-  fEtaMax     =  0.9;
+  fLegoNbinEta    =  36;
+  fLegoNbinPhi    = 124;
+  fLegoPhiMin     =   0.;
+  fLegoPhiMax     = 2. * TMath::Pi();
+  fLegoEtaMin     = -0.9;
+  fLegoEtaMax     =  0.9;
+  fOnlySignal = kFALSE;
+  fOnlyBkgd = kFALSE;
 }
  
-
 ////////////////////////////////////////////////////////////////////////
  
 void AliUA1JetHeader::PrintParameters() const
-
 {
-  //
   // prints out parameters of jet algorithm
-  //
 
   cout << " UA1 jet algorithm " << endl;
   cout << " * Jet parameters: " << endl;
@@ -71,10 +66,10 @@ void AliUA1JetHeader::PrintParameters() const
   cout << "   Maximum allowed move: " << fMaxMove << endl;
   cout << "   Precision for background: " << fPrecBg << endl;
   cout << " * Lego parameters: " << endl;
-  cout << "   Number of bins in eta: " << fNbinEta<< endl;
-  cout << "   Number of bins in phi: " << fNbinPhi<< endl;
-  cout << "   Minimum azimuthal angle: " << fPhiMin<< endl;
-  cout << "   Maximum azimuthal angle: " << fPhiMax<< endl;
-  cout << "   Minimum rapidity angle: " << fEtaMin<< endl;
-  cout << "   Maximum rapidity angle: " << fEtaMax<< endl;
+  cout << "   Number of bins in eta: " << fLegoNbinEta<< endl;
+  cout << "   Number of bins in phi: " << fLegoNbinPhi<< endl;
+  cout << "   Minimum azimuthal angle: " << fLegoPhiMin<< endl;
+  cout << "   Maximum azimuthal angle: " << fLegoPhiMax<< endl;
+  cout << "   Minimum rapidity angle: " << fLegoEtaMin<< endl;
+  cout << "   Maximum rapidity angle: " << fLegoEtaMax<< endl;
 }
index 62eb5aa..f22c963 100755 (executable)
@@ -4,7 +4,6 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
  
 //---------------------------------------------------------------------
 // Jet header class for UA1 algorithm
 // Stores the parameters of the UA1 jet algorithm
@@ -29,12 +28,14 @@ class AliUA1JetHeader : public AliJetHeader
   Float_t GetMinMove()   const  {return fMinMove;}
   Float_t GetMaxMove()   const  {return fMaxMove;}
   Float_t GetPrecBg()    const  {return fPrecBg;}
-  Int_t   GetNbinEta()   const  {return fNbinEta;}
-  Int_t   GetNbinPhi()   const  {return fNbinPhi;}
-  Float_t GetEtaMin()    const  {return fEtaMin;}
-  Float_t GetEtaMax()    const  {return fEtaMax;}
-  Float_t GetPhiMin()    const  {return fPhiMin;}
-  Float_t GetPhiMax()    const  {return fPhiMax;}
+  Int_t   GetLegoNbinEta()   const  {return fLegoNbinEta;}
+  Int_t   GetLegoNbinPhi()   const  {return fLegoNbinPhi;}
+  Float_t GetLegoEtaMin()    const  {return fLegoEtaMin;}
+  Float_t GetLegoEtaMax()    const  {return fLegoEtaMax;}
+  Float_t GetLegoPhiMin()    const  {return fLegoPhiMin;}
+  Float_t GetLegoPhiMax()    const  {return fLegoPhiMax;}
+  Bool_t GetOnlySignal()     const {return fOnlySignal;}
+  Bool_t GetOnlyBkgd()     const {return fOnlyBkgd;}
 
   // Setters
   void SetRadius(Float_t f) {fConeRadius=f;}
@@ -45,12 +46,14 @@ class AliUA1JetHeader : public AliJetHeader
   void SetMinMove(Float_t f) {fMinMove=f;}
   void SetMaxMove(Float_t f) {fMaxMove=f;}
   void SetPrecBg(Float_t f) {fPrecBg=f;}
-  void SetNbinEta(Int_t f) {fNbinEta=f;}
-  void SetNbinPhi(Int_t f) {fNbinPhi=f;}
-  void SetEtaMin(Float_t f) {fEtaMin=f;}
-  void SetEtaMax(Float_t f) {fEtaMax=f;}
-  void SetPhiMin(Float_t f) {fPhiMin=f;}
-  void SetPhiMax(Float_t f) {fPhiMax=f;}
+  void SetLegoNbinEta(Int_t f) {fLegoNbinEta=f;}
+  void SetLegoNbinPhi(Int_t f) {fLegoNbinPhi=f;}
+  void SetLegoEtaMin(Float_t f) {fLegoEtaMin=f;}
+  void SetLegoEtaMax(Float_t f) {fLegoEtaMax=f;}
+  void SetLegoPhiMin(Float_t f) {fLegoPhiMin=f;}
+  void SetLegoPhiMax(Float_t f) {fLegoPhiMax=f;}
+  void SetOnlySignal(Bool_t b) {fOnlySignal= b;}
+  void SetOnlyBkgd(Bool_t b) {fOnlyBkgd= b;}
 
   // others
   void PrintParameters() const; 
@@ -68,12 +71,15 @@ protected:
   Float_t fMaxMove;         // max cone move
   Float_t fPrecBg;          // max value of change for BG (in %)
   // parameters for legos
-  Int_t   fNbinEta;         //! number of cells in eta
-  Int_t   fNbinPhi;         //! number of cells in phi
-  Float_t fEtaMin;          //! minimum eta  
-  Float_t fEtaMax;          //! maximum eta
-  Float_t fPhiMin;          //! minimun phi
-  Float_t fPhiMax;          //! maximum phi
+  Int_t   fLegoNbinEta;         //! number of cells in eta
+  Int_t   fLegoNbinPhi;         //! number of cells in phi
+  Float_t fLegoEtaMin;          //! minimum eta  
+  Float_t fLegoEtaMax;          //! maximum eta
+  Float_t fLegoPhiMin;          //! minimun phi
+  Float_t fLegoPhiMax;          //! maximum phi
+  // parameters for Jet search
+  Bool_t fOnlySignal; // use only signal
+  Bool_t fOnlyBkgd;   // use only background
 
   ClassDef(AliUA1JetHeader,1)
 };
index e575817..d8abac0 100644 (file)
 #pragma link C++ class AliUA1JetFinder+;
 #pragma link C++ class AliJetReaderHeader+;
 #pragma link C++ class AliJetESDReaderHeader+;
-#pragma link C++ class AliJetMCReaderHeader+;
 #pragma link C++ class AliJetReader+;
 #pragma link C++ class AliJetESDReader+;
-#pragma link C++ class AliJetMCReader+;
 #pragma link C++ class AliJetKineReader+;
 #pragma link C++ class AliJetKineReaderHeader+;
 #pragma link C++ class AliJetControlPlots+;
@@ -24,4 +22,5 @@
 #pragma link C++ class AliJetProductionData+;
 #pragma link C++ class AliJetProductionDataPDC2004+;
 #pragma link C++ class AliJetAnalysis+;
+#pragma link C++ class AliJetDistributions+;
 #endif
index e7f332d..4e1e611 100755 (executable)
@@ -14,6 +14,7 @@ extern "C" {
        Float_t etaj[2][100];
        Float_t phij[2][100];
        Int_t   ncellj[100];
+      Float_t etavg;
     } UA1JetsCommon;
 
 #define UA1JETS COMMON_BLOCK(UA1JETS,ua1jets)
index 80508af..6d0d1a0 100644 (file)
@@ -4,10 +4,10 @@ SRCS =        AliJet.cxx AliJetHeader.cxx AliPxconeJetHeader.cxx \
        AliJetFinder.cxx AliPxconeJetFinder.cxx AliJetReaderHeader.cxx \
        AliJetESDReaderHeader.cxx AliJetReader.cxx AliJetESDReader.cxx \
        AliJetControlPlots.cxx AliUA1JetHeader.cxx AliUA1JetFinder.cxx \
-       AliJetMCReaderHeader.cxx AliJetMCReader.cxx AliLeading.cxx \
+       AliLeading.cxx \
        AliJetKineReader.cxx AliJetKineReaderHeader.cxx \
        AliJetProductionData.cxx AliJetProductionDataPDC2004.cxx \
-       AliJetAnalysis.cxx
+       AliJetAnalysis.cxx AliJetDistributions.cxx
 
 FSRCS = pxcone.F ua1_jet_finder.F
 
index f7df5e7..bf11ac4 100644 (file)
@@ -13,7 +13,13 @@ void read_jets(const char* fn = "jets.root")
     TH1F* dr1H = new TH1F("dr1H", "delta R",  160., 0.,   2.);
     TH1F* dr2H = new TH1F("dr2H", "delta R",  160., 0.,   2.);
     TH1F* dr3H = new TH1F("dr4H", "delta R",  160., 0.,   2.);
-    TH1F* etaH = new TH1F("etaH", "eta",  160., -2.,   2.);
+    TH1F* etaH  = new TH1F("etaH",  "eta",  160., -2.,   2.);
+    TH1F* eta1H = new TH1F("eta1H", "eta",  160., -2.,   2.);
+    TH1F* eta2H = new TH1F("eta2H", "eta",  160., -2.,   2.);
+
+    TH1F* phiH  = new TH1F("phiH",  "phi",  160., -3.,   3.);
+    TH1F* phi1H = new TH1F("phi1H", "phi",  160.,  0.,   6.28);
+    TH1F* phi2H = new TH1F("phi2H", "phi",  160.,  0.,   6.28);
     
 
   // load jet library
@@ -97,10 +103,18 @@ void read_jets(const char* fn = "jets.root")
 
              Float_t egen = gjets->GetPt(igen);
              e1H->Fill(gjets->GetPt(igen));
-             
-             if (egen > 105. && egen < 125.) {
+             Float_t etag = gjets->GetEta(igen);
+             Float_t phig = gjets->GetPhi(igen);
+             Float_t dphi = phig - phij;
+
+             if (egen > 125. && egen < 150.) {
+                 phiH->Fill((dphi));
+                 etaH->Fill(etag - etaj);
+                 phi1H->Fill(phig);
+                 phi2H->Fill(phij);              
+                 eta1H->Fill(etag);
+                 eta2H->Fill(etaj);
                  e4H->Fill(emax);
-                 if (rmin > 0.3) etaH->Fill(etaj);
                  dr2H->Fill(rmin);
              }
          }
@@ -124,6 +138,7 @@ void read_jets(const char* fn = "jets.root")
        Float_t deta = etag-etal;
        Float_t dphi = TMath::Abs(phig - phil);
        if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
+
        Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
        
        if (r  < rmin) {
@@ -131,7 +146,7 @@ void read_jets(const char* fn = "jets.root")
                igen = j;
        }
     }
+    if (egen > 125. && egen < 150.) 
        dr3H->Fill(rmin);
 
 //    cout << " Generated Jets:" << endl;
@@ -156,8 +171,8 @@ void read_jets(const char* fn = "jets.root")
   dr2H->Draw("");
   
   TCanvas* c3 = new TCanvas("c3");
-  dr3H->Draw();
-  dr2H->Draw("same");
+  dr2H->Draw();
+  dr3H->Draw("same");
 
   TCanvas* c4 = new TCanvas("c4");
   eH->Draw();
@@ -165,6 +180,19 @@ void read_jets(const char* fn = "jets.root")
   TCanvas* c5 = new TCanvas("c5");
   etaH->Draw();
 
+  TCanvas* c5a = new TCanvas("c5a");
+  eta1H->Draw();
+
+  TCanvas* c5b = new TCanvas("c5b");
+  eta2H->Draw();
+
   TCanvas* c6 = new TCanvas("c6");
   e4H->Draw();
+  TCanvas* c7 = new TCanvas("c7");
+  phiH->Draw();
+
+  TCanvas* c7a = new TCanvas("c7a");
+  phi1H->Draw();
+  TCanvas* c7b = new TCanvas("c7b");
+  phi2H->Draw();
 }
index 1bf3c02..2b0a74f 100755 (executable)
@@ -38,10 +38,10 @@ c:>------------------------------------------------------------------
       integer kpri
       data    kpri /0/
       integer njet, ncellj
-      real    etj, etaj, phij
+      real    etj, etaj, phij, etavg
 *    Results
       COMMON /UA1JETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2), 
-     +     NCELLJ(100)
+     +     NCELLJ(100), ETAVG
 *    Cell Geometry
       COMMON /UA1CELL/ etaCellSize, phiCellSize
 *    Parameters
@@ -227,17 +227,17 @@ c*-sum up unused cells within required distance of given eta/phi
             
 c*-reject cluster below minimum Ej_min
 c* protection (am)
-            etas=eta+etas/ets
-            arg = 0.
-            if (ets .ne. 0.) then
-               if (abs(etas/ets) .lt. 23.719) then
-                  arg = ets * cosh(etas/ets)
-               else
-                  arg = 1.e10
-               endif
-            endif
+
+c            arg = 0.
+c            if (ets .ne. 0.) then
+c               if (abs(etas/ets) .lt. 23.719) then
+c                  arg = ets * cosh(etas/ets)
+c               else
+c                  arg = 1.e10
+c               endif
+c            endif
             
-            if(arg .lt. ej_min) then
+            if(ets .lt. ej_min) then
                do k=1,ncell
                   if(flag(k).le.0) flag(k)=0
                enddo
@@ -246,6 +246,7 @@ c*-eles, store flags and jet variables
                do k=1,ncell
                   if(flag(k).eq.-1) flag(k)=1
                enddo
+               etas=eta+etas/ets
                phi=phi+phis/ets
                do while(phi .ge. C_2PI)
                   phi=phi-C_2PI
@@ -260,6 +261,7 @@ c*-eles, store flags and jet variables
                etaj(njet,2)=etas
                phij(njet,2)=phi
                ncellj(njet)=nc
+               etavg = et_ave
             endif 
          endif
          i=i+1