First commit of new jet reconstruction and analysis code to be used for the
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Jun 2005 15:20:15 +0000 (15:20 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Jun 2005 15:20:15 +0000 (15:20 +0000)
ALICE PPR.

43 files changed:
JETAN/AliJet.cxx [new file with mode: 0644]
JETAN/AliJet.h [new file with mode: 0644]
JETAN/AliJetControlPlots.cxx [new file with mode: 0755]
JETAN/AliJetControlPlots.h [new file with mode: 0755]
JETAN/AliJetESDReader.cxx [new file with mode: 0755]
JETAN/AliJetESDReader.h [new file with mode: 0755]
JETAN/AliJetESDReaderHeader.cxx [new file with mode: 0755]
JETAN/AliJetESDReaderHeader.h [new file with mode: 0755]
JETAN/AliJetFinder.cxx [new file with mode: 0644]
JETAN/AliJetFinder.h [new file with mode: 0755]
JETAN/AliJetHeader.cxx [new file with mode: 0755]
JETAN/AliJetHeader.h [new file with mode: 0755]
JETAN/AliJetKineReader.cxx [new file with mode: 0644]
JETAN/AliJetKineReader.h [new file with mode: 0644]
JETAN/AliJetKineReaderHeader.cxx [new file with mode: 0644]
JETAN/AliJetKineReaderHeader.h [new file with mode: 0644]
JETAN/AliJetMCReader.cxx [new file with mode: 0644]
JETAN/AliJetMCReader.h [new file with mode: 0644]
JETAN/AliJetMCReaderHeader.cxx [new file with mode: 0644]
JETAN/AliJetMCReaderHeader.h [new file with mode: 0644]
JETAN/AliJetReader.cxx [new file with mode: 0755]
JETAN/AliJetReader.h [new file with mode: 0755]
JETAN/AliJetReaderHeader.cxx [new file with mode: 0755]
JETAN/AliJetReaderHeader.h [new file with mode: 0755]
JETAN/AliLeading.cxx [new file with mode: 0644]
JETAN/AliLeading.h [new file with mode: 0644]
JETAN/AliPxconeJetFinder.cxx [new file with mode: 0755]
JETAN/AliPxconeJetFinder.h [new file with mode: 0755]
JETAN/AliPxconeJetHeader.cxx [new file with mode: 0755]
JETAN/AliPxconeJetHeader.h [new file with mode: 0755]
JETAN/AliUA1JetFinder.cxx [new file with mode: 0755]
JETAN/AliUA1JetFinder.h [new file with mode: 0755]
JETAN/AliUA1JetHeader.cxx [new file with mode: 0755]
JETAN/AliUA1JetHeader.h [new file with mode: 0755]
JETAN/JetAnalysisLinkDef.h [new file with mode: 0644]
JETAN/UA1Common.h [new file with mode: 0755]
JETAN/libJETAN.pkg [new file with mode: 0644]
JETAN/pxcone.F [new file with mode: 0755]
JETAN/read_jets.C [new file with mode: 0644]
JETAN/read_jets_kine.C [new file with mode: 0644]
JETAN/testpxcone.C [new file with mode: 0755]
JETAN/testua1.C [new file with mode: 0755]
JETAN/ua1_jet_finder.F [new file with mode: 0755]

diff --git a/JETAN/AliJet.cxx b/JETAN/AliJet.cxx
new file mode 100644 (file)
index 0000000..7a6cd63
--- /dev/null
@@ -0,0 +1,258 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//---------------------------------------------------------------------
+// Jet class 
+// Stores the output of a jet algorithm
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+#include <Riostream.h>
+#include <TClonesArray.h>
+#include <TLorentzVector.h>
+
+#include "AliJet.h"
+ClassImp(AliJet)
+////////////////////////////////////////////////////////////////////////
+
+AliJet::AliJet() 
+{
+  //
+  // Constructor
+  //
+  fJets = new TClonesArray("TLorentzVector",1000);
+  fNInput=0;
+  fNJets=0;
+  fInJet=TArrayI();
+  fMultiplicities=TArrayI();
+} 
+
+////////////////////////////////////////////////////////////////////////
+
+AliJet::~AliJet()
+{
+  // destructor
+  if (fJets) {
+    fJets->Delete();
+    delete fJets;
+  }
+}
+
+////////////////////////////////////////////////////////////////////////
+
+Bool_t AliJet::OutOfRange(Int_t i, const char *s) const
+{
+  // 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;
+  }
+  return kFALSE;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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)
+{
+  // gets multiplicity of i-jet
+  if (OutOfRange(i, "AliJet::GetMultiplicity:")) return 0;
+  return fMultiplicities[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();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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());
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJet::AddJet(Double_t px, Double_t py, Double_t pz, Double_t e)
+{
+// Add new jet to the list
+  new ((*fJets)[fNJets++]) TLorentzVector(px,py,pz,e);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJet::SetInJet(Int_t* j)
+{
+  // set information of which input object belongs
+  // to each jet. filled in by AliJetFinder
+  if (fNInput>0) fInJet.Set(fNInput, j);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJet::SetMultiplicities(Int_t* m)
+{
+  // set information of jet multiplicities
+  // filled in by AliJetFinder
+  if (fNJets>0) fMultiplicities.Set(fNJets, m);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJet::ClearJets(Option_t *option)
+{
+  // reset all values
+  fJets->Clear(option);
+  fNInput=0;
+  fNJets=0;
+  fMultiplicities.Set(0);
+  fInJet.Set(0); 
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJet::PrintJets()
+{
+// Print jet information
+  if (fNJets == 0) {
+    cout << " AliJet::PrintJets: There are no jets in this event " << endl;
+    return;
+  }
+  cout << " AliJet::PrintJets: There are " << fNJets
+       << " jets in this event" << endl;
+  for(Int_t i=0;i<fNJets;i++) {
+    cout << "   Jet " << i << " (px,py,pz,en)=(" << GetPx(i)
+        << "," << GetPy(i)
+        << "," << GetPz(i)
+        << "," << GetE(i) 
+        << ")" << endl;
+    cout << "         (pt,eta,phi)=(" << GetPt(i)
+        << "," << GetEta(i)
+        << "," << GetPhi(i) << ")" << endl;
+    cout << "         # of tracks =" << GetMultiplicity(i) << endl;
+  }
+}
diff --git a/JETAN/AliJet.h b/JETAN/AliJet.h
new file mode 100644 (file)
index 0000000..c4b6dac
--- /dev/null
@@ -0,0 +1,70 @@
+#ifndef ALIJET_H
+#define ALIJET_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Jet class 
+// Stores the output of a jet algorithm
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+#include <TObject.h>
+#include <TArrayI.h>
+
+
+class TClonesArray;
+class TLorentzVector;
+class AliJet : public TObject
+{
+ public:
+  AliJet();
+  ~AliJet();
+
+  // Getters
+  Int_t GetNinput() const { return fNInput; }
+  Int_t GetNJets() const {return fNJets;}
+  TClonesArray* GetJets() const {return fJets;}
+  TArrayI GetInJet() const {return fInJet;}
+  TArrayI GetMultiplicities() const {return fMultiplicities;}
+
+  TLorentzVector* GetJet(Int_t i);
+  Int_t GetMultiplicity(Int_t i);
+  Double_t GetPx(Int_t i);
+  Double_t GetPy(Int_t i);
+  Double_t GetPz(Int_t i);
+  Double_t GetP(Int_t i);
+  Double_t GetE(Int_t i);
+  Double_t GetPt(Int_t i);
+  Double_t GetEta(Int_t i);
+  Double_t GetPhi(Int_t i);
+  Double_t GetTheta(Int_t i);
+  Double_t GetMass(Int_t i);
+   
+  // 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);
+
+  // others
+  Bool_t OutOfRange(Int_t i, const char *s) const;
+  void ClearJets(Option_t *option="");
+  void PrintJets();
+
+ protected:
+
+  Int_t fNInput;               // number of input objects
+  Int_t fNJets;                // number of jets found
+  TArrayI fInJet;              // i-input object belongs to k-jet 
+  TArrayI fMultiplicities;     // Multiplicity of each jet
+  TClonesArray* fJets;         // 4-momenta of jets
+
+  ClassDef(AliJet,1)
+};
+#endif
diff --git a/JETAN/AliJetControlPlots.cxx b/JETAN/AliJetControlPlots.cxx
new file mode 100755 (executable)
index 0000000..6e37ef5
--- /dev/null
@@ -0,0 +1,175 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+  
+//---------------------------------------------------------------------
+// Jet Control Plots class 
+// manages histograms with control plots of jet searching
+// Stores the output of a jet algorithm
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TLorentzVector.h>
+#include <TFile.h> 
+#include <TClonesArray.h>
+#include <TH1I.h>
+#include <TH1D.h>
+#include "AliJetReader.h"
+#include "AliJet.h"
+
+#include "AliJetControlPlots.h"
+ClassImp(AliJetControlPlots)
+  
+////////////////////////////////////////////////////////////////////////
+
+AliJetControlPlots::AliJetControlPlots()
+{
+  //
+  // Constructor
+  //
+  fNJetT=0;
+
+  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");
+
+  fPtH = new TH1D("fPtH","Pt of Jets",50,0.,200.);
+  SetProperties(fPtH,"P_{t} [GeV]","Entries");
+
+  fEtaH = new TH1D("fEtaH","Pseudorapidity of Jets",30,-1.5,1.5);
+  SetProperties(fEtaH,"#eta","Entries");
+
+  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());
+  SetProperties(fPhiH,"#phi","Entries");
+  
+  fInJetH = new  TH1D("fInJetH","Percentage of particles in jets",
+                     50,0,1);
+  SetProperties(fInJetH,"percentage of particles in jets","Entries");
+}
+
+////////////////////////////////////////////////////////////////////////
+
+AliJetControlPlots::~AliJetControlPlots()
+{
+  //
+  // Destructor
+  //
+  delete fNJetsH;
+  delete fMultH;
+  delete fPtH;
+  delete fEtaH;
+  delete fEneH;
+  delete fFragH;
+  delete fFragLnH;
+  delete fPhiH;
+  delete fInJetH;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetControlPlots::FillHistos(AliJet *j, AliJetReader *r)
+{
+// Fills the histograms
+
+  Int_t nj = j->GetNJets();
+  fNJetsH->Fill(nj);
+  fNJetT+=nj;
+
+  // 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));
+  }
+  if (nj>0) 
+    fInJetH->Fill(((Double_t) mjTot)/((Double_t) j->GetNinput()));
+
+  // fragmentation 
+  TClonesArray *lvArray = r->GetMomentumArray();
+  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));
+  }
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetControlPlots::Normalize()
+{
+  if (fNJetT == 0) return;
+  fFragH->Sumw2();
+  fFragH->Scale(20.0/((Double_t) fNJetT));
+  fFragLnH->Sumw2();
+  fFragLnH->Scale(2.0/((Double_t) fNJetT));
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetControlPlots::PlotHistos()
+{
+  // style
+  gStyle->SetOptStat(kFALSE);
+  gStyle->SetOptTitle(kFALSE);
+  
+  TCanvas *c = new TCanvas("c","AliJetControlPlots",50,200,900,700);
+  c->Divide(3,3);
+  c->cd(1);  fNJetsH->Draw("e1");
+  c->cd(2);  fMultH->Draw("e1");
+  c->cd(3);  fInJetH->Draw("e1");
+  c->cd(4);  fPtH->Draw("e1");
+  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");
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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);
+}
+
+
diff --git a/JETAN/AliJetControlPlots.h b/JETAN/AliJetControlPlots.h
new file mode 100755 (executable)
index 0000000..ff372db
--- /dev/null
@@ -0,0 +1,64 @@
+#ifndef ALIJETCONTROLPLOTS_H
+#define ALIJETCONTROLPLOTS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Jet Control Plots class 
+// manages histograms with control plots of jet searching
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+
+#include <TObject.h>
+
+class TFile;
+class TClonesArray;
+class TH1I;
+class TH1D;
+class AliJetReader;
+class AliJet;
+
+class AliJetControlPlots : public TObject
+{
+ public:
+  AliJetControlPlots();
+  ~AliJetControlPlots();
+
+  // setter
+  // getters
+  TH1I *GetNJetsH() {return fNJetsH;}
+  TH1I *GetMultH() {return fMultH;}
+  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;}
+  
+  // others
+  void FillHistos(AliJet *j,AliJetReader *r);
+  void PlotHistos();
+  void SetProperties(TH1* h,const char* x, const char* y) const;
+  void Normalize();
+
+ protected:
+  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
+
+  ClassDef(AliJetControlPlots,1)
+};
+#endif
+
diff --git a/JETAN/AliJetESDReader.cxx b/JETAN/AliJetESDReader.cxx
new file mode 100755 (executable)
index 0000000..956bb06
--- /dev/null
@@ -0,0 +1,136 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * 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
+//
+
+
+#include <Riostream.h>
+#include <TSystem.h>
+#include <TLorentzVector.h>
+#include "AliJetESDReader.h"
+#include "AliJetESDReaderHeader.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+
+ClassImp(AliJetESDReader)
+
+AliJetESDReader::AliJetESDReader()
+{
+  // Constructor    
+
+  fReaderHeader = 0x0;
+  fMass = 0;
+  fSign = 0;
+}
+
+//____________________________________________________________________________
+
+AliJetESDReader::~AliJetESDReader()
+{
+  // Destructor
+  delete fReaderHeader;
+}
+
+//____________________________________________________________________________
+
+void AliJetESDReader::OpenInputFiles()
+
+{
+  // chain for the ESDs
+  fChain   = new TChain("esdTree");
+  fChainMC = new TChain("mcStackTree");
+
+  // get directory and pattern name from the header
+  const char* dirName=fReaderHeader->GetDirectory();
+  const char* pattern=fReaderHeader->GetPattern();
+    
+  // Add files matching patters to the chain
+  void *dir  = gSystem->OpenDirectory(dirName);
+  const char *name = 0x0;
+  while ((name = gSystem->GetDirEntry(dir))){
+    if (strstr(name,pattern)){
+      printf("Adding %s\n",name);
+      char path[256];
+      sprintf(path,"%s/%s",dirName,name);
+      fChain->AddFile(path,-1,"esdTree");
+      fChainMC->AddFile(path,-1,"mcStackTree");
+     }
+  }
+
+  gSystem ->FreeDirectory(dir);
+  fChain  ->SetBranchAddress("ESD",    &fESD);
+  fChainMC->SetBranchAddress("Header", &fAliHeader);
+  fChainMC->SetBranchAddress("Stack",  &fArrayMC);
+
+  int nMax = fChain->GetEntries(); 
+  printf("\nTotal number of events in chain= %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));
+  }
+}
+
+//____________________________________________________________________________
+
+void AliJetESDReader::FillMomentumArray(Int_t event)
+{
+// Fills the momentum array
+  Int_t goodTrack = 0;
+  Int_t nt = 0;
+  Float_t pt;
+  Double_t p[3]; // track 3 momentum vector
+
+  // 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();
+  // get cuts set by user
+  Float_t ptMin = ((AliJetESDReaderHeader*) fReaderHeader)->GetPtCut();
+
+  //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
+
+      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]));
+      goodTrack++;
+  }
+
+  printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
+}
+
+
diff --git a/JETAN/AliJetESDReader.h b/JETAN/AliJetESDReader.h
new file mode 100755 (executable)
index 0000000..790c83a
--- /dev/null
@@ -0,0 +1,36 @@
+#ifndef ALIJETESDREADER_H
+#define ALIJETESDREADER_H
+/* 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;
+
+
+class AliJetESDReader : public AliJetReader
+{
+ public: 
+  AliJetESDReader();
+  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
+
+  // Setters
+  void FillMomentumArray(Int_t event); 
+  void OpenInputFiles();
+   
+ protected:
+  Float_t fMass;    // Particle mass
+  Int_t   fSign;    // Particle sign
+
+  ClassDef(AliJetESDReader,1)
+};
+#endif
diff --git a/JETAN/AliJetESDReaderHeader.cxx b/JETAN/AliJetESDReaderHeader.cxx
new file mode 100755 (executable)
index 0000000..41fb5ae
--- /dev/null
@@ -0,0 +1,43 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * 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
+//
+#include "AliJetESDReaderHeader.h"
+
+ClassImp(AliJetESDReaderHeader)
+
+//____________________________________________________________________________
+AliJetESDReaderHeader::AliJetESDReaderHeader():
+ AliJetReaderHeader("AliJetESDReaderHeader") 
+{
+  // Constructor
+  SetPtCut();
+  SetDCA();
+  SetTLength();
+  SetReadSignalOnly(kFALSE);
+  
+}
+
+//____________________________________________________________________________
+AliJetESDReaderHeader::~AliJetESDReaderHeader()
+{
+  // destructor
+}
diff --git a/JETAN/AliJetESDReaderHeader.h b/JETAN/AliJetESDReaderHeader.h
new file mode 100755 (executable)
index 0000000..8b8493e
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALIJETESDREADERHEADER_H
+#define ALIJETESDREADERHEADER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+// Jet ESD Reader Header
+// Header for the ESD reader in the jet analysis
+// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+
+#include "AliJetReaderHeader.h"
+class AliJetESDReaderHeader : public AliJetReaderHeader
+{
+
+ public:
+  AliJetESDReaderHeader();
+  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;}
+  
+         
+  // Setters
+  virtual void SetPtCut(const Float_t par = 2.0) {fPtCut = par;}
+  virtual void SetDCA(const Float_t dca = 3.0) {fDCA = dca;}
+  virtual void SetTLength(const Float_t length = 0.0) {fTLength = length;}
+  virtual void SetReadSignalOnly(Bool_t flag = kTRUE) {fReadSignalOnly = 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);
+};
+#endif
diff --git a/JETAN/AliJetFinder.cxx b/JETAN/AliJetFinder.cxx
new file mode 100644 (file)
index 0000000..4682dcb
--- /dev/null
@@ -0,0 +1,206 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * 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 
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TFile.h>
+#include "AliGenPythiaEventHeader.h"
+#include "AliJetFinder.h"
+#include "AliJet.h"
+#include "AliJetReader.h"
+#include "AliJetReaderHeader.h"
+#include "AliJetControlPlots.h"
+#include "AliLeading.h"
+#include "AliHeader.h"
+
+
+ClassImp(AliJetFinder)
+
+////////////////////////////////////////////////////////////////////////
+
+AliJetFinder::AliJetFinder()
+{
+  //
+  // Constructor
+  //
+  fOut = 0;
+  fJets = new AliJet();
+  fGenJets = new AliJet();
+  fLeading = new AliLeading();
+  fReader = 0;
+  fPlots = 0;
+  SetPlotMode(kFALSE);
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+AliJetFinder::~AliJetFinder()
+{
+  //
+  // destructor
+  //
+
+  // here reset and delete jets
+  fJets->ClearJets();
+  delete fJets;
+  fGenJets->ClearJets();
+  delete fGenJets;
+  // close file
+  if (fOut) {
+    fOut->Close();
+    fOut->Delete();
+  }
+  delete fOut;
+  // reset and delete control plots
+  if (fPlots) delete fPlots;
+  delete fLeading;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetFinder::SetOutputFile(const char *name)
+{
+  // opens output file 
+  fOut = new TFile(name,"recreate");
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetFinder::PrintJets()
+{
+//
+// 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
+  fPlotMode=b;
+  if (b && !fPlots) fPlots = new AliJetControlPlots(); 
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetFinder::WriteJetsToFile(Int_t i)
+{
+// Writes the jets to file
+  fOut->cd();
+  char hname[30];
+  sprintf(hname,"TreeJ%d",i);
+  TTree* jetT = new TTree(hname,"AliJet");
+  jetT->Branch("FoundJet",&fJets,1000);
+  jetT->Branch("GenJet",&fGenJets,1000);
+  jetT->Branch("LeadingPart",&fLeading,1000);
+  jetT->Fill();
+  jetT->Write(hname);
+  delete jetT;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetFinder::WriteRHeaderToFile()
+{
+  // write reader header
+  fOut->cd();
+  AliJetReaderHeader *rh = fReader->GetReaderHeader();
+  rh->Write();
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetFinder::GetGenJets()
+{
+// Get the generated jet information from mc header
+  fGenJets->SetNinput(1);
+  AliHeader* alih = fReader->GetAliHeader(); 
+  if (alih == 0) return;
+  AliGenEventHeader * genh = alih->GenEventHeader();
+  if (genh == 0) return;
+  Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets(); 
+  Int_t* m = new Int_t[nj];
+  Int_t* k = new Int_t[nj];
+  for (Int_t i=0; i< nj; i++) {
+    Float_t p[4];
+    ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
+    fGenJets->AddJet(p[0],p[1],p[2],p[3]);
+    m[i]=1;
+    k[i]=i;
+  }
+  fGenJets->SetMultiplicities(m);
+  fGenJets->SetInJet(k);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetFinder::Run()
+{
+  // do some initialization
+  Init();
+
+  // connect files
+  fReader->OpenInputFiles();
+
+  // write headers
+  if (fOut) {
+      fOut->cd();
+      WriteRHeaderToFile();
+      WriteJHeaderToFile();
+  }
+
+  // loop over events
+  Int_t nFirst,nLast;
+  nFirst = fReader->GetReaderHeader()->GetFirstEvent();
+  nLast = fReader->GetReaderHeader()->GetLastEvent();
+  // loop over events
+  for (Int_t i=nFirst;i<nLast;i++) {
+      fReader->FillMomentumArray(i);
+      fLeading->FindLeading(fReader);
+      GetGenJets();
+      FindJets();
+      if (fOut) WriteJetsToFile(i);
+      if (fPlots) fPlots->FillHistos(fJets,fReader);
+      fLeading->Reset();
+      fGenJets->ClearJets();
+      Reset();
+  } 
+  // write out
+  if (fPlots) {
+      fPlots->Normalize();
+      fPlots->PlotHistos();
+  }
+  if (fOut) {
+      fOut->cd();
+      fPlots->Write();
+      fOut->Close();
+  }
+}
+
diff --git a/JETAN/AliJetFinder.h b/JETAN/AliJetFinder.h
new file mode 100755 (executable)
index 0000000..efefbcf
--- /dev/null
@@ -0,0 +1,63 @@
+#ifndef ALIJETFINDER_H
+#define ALIJETFINDER_H
+/* 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 
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include <TObject.h>
+
+class TFile;
+class AliJet;
+class AliJetReader;
+class AliJetControlPlots;
+class AliLeading;
+
+
+class AliJetFinder : public TObject 
+{
+ public:
+
+  AliJetFinder();
+  virtual ~AliJetFinder();
+
+  // getters
+  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");
+  virtual void SetJetReader(AliJetReader* r) {fReader=r;}
+
+  // others
+  virtual void PrintJets();
+  virtual void Run();
+  virtual void WriteJetsToFile(Int_t i);
+  virtual void WriteRHeaderToFile();
+  // the following have to be implemented for each specific finder
+  virtual void Init() { }
+  virtual void Reset() { }
+  virtual void FindJets() { }
+  virtual void WriteJHeaderToFile() { }
+  virtual void GetGenJets();
+
+ protected:
+  Bool_t fPlotMode;              // do you want control plots?
+  AliJet* fJets;                 // pointer to jet class
+  AliJet* fGenJets;               // pointer to generated jets
+  AliLeading* fLeading;          // pointer to leading particle data 
+  AliJetReader* fReader;         // pointer to reader
+  AliJetControlPlots* fPlots;    // pointer to control plots
+  TFile* fOut;                   // output file
+
+  ClassDef(AliJetFinder,1)
+};
+
+#endif
diff --git a/JETAN/AliJetHeader.cxx b/JETAN/AliJetHeader.cxx
new file mode 100755 (executable)
index 0000000..0c39bef
--- /dev/null
@@ -0,0 +1,50 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//---------------------------------------------------------------------
+// Jet header base class 
+// Stores a comment which describes the jet analysis
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+#include "AliJetHeader.h"
+ClassImp(AliJetHeader)
+////////////////////////////////////////////////////////////////////////
+
+AliJetHeader::AliJetHeader():
+  TNamed("AliJetHeader", "Jet Header"),
+  fComment("No comment")
+{
+  //
+  // Constructor
+  //
+}
+////////////////////////////////////////////////////////////////////////
+
+AliJetHeader::AliJetHeader(const char * name):
+  TNamed(name, "Jet Header"),
+  fComment("No comment")
+{
+  //
+  // Constructor
+  //
+}
+
+////////////////////////////////////////////////////////////////////////
+
diff --git a/JETAN/AliJetHeader.h b/JETAN/AliJetHeader.h
new file mode 100755 (executable)
index 0000000..d85de2f
--- /dev/null
@@ -0,0 +1,39 @@
+#ifndef ALIJETHEADER_H
+#define ALIJETHEADER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Jet header base class 
+// Stores a comment which describes the jet analysis
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+#include <TNamed.h>
+#include <TString.h>
+class AliJetHeader : public TNamed
+{
+ public:
+  AliJetHeader(const char* name);
+  AliJetHeader();
+  virtual ~AliJetHeader() { }
+
+  // Getters
+  virtual TString GetComment() {return fComment;} 
+   
+  // Setters
+  virtual void SetComment(const char* com) {fComment=TString(com);} 
+
+  // others
+  
+protected:
+  TString fComment; // a comment 
+
+  ClassDef(AliJetHeader,1)
+};
+#endif
diff --git a/JETAN/AliJetKineReader.cxx b/JETAN/AliJetKineReader.cxx
new file mode 100644 (file)
index 0000000..a09767e
--- /dev/null
@@ -0,0 +1,132 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+// 
+// Jet kinematics reader 
+// MC reader for jet analysis
+// 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 <TLorentzVector.h>
+#include <TSystem.h>
+#include <TDatabasePDG.h>
+#include <TRandom.h>
+// From AliRoot ...
+#include "AliJetKineReaderHeader.h"
+#include "AliJetKineReader.h"
+#include "AliRunLoader.h"
+#include "AliStack.h"
+
+ClassImp(AliJetKineReader)
+
+AliJetKineReader::AliJetKineReader()
+{
+  // Constructor
+  fReaderHeader = 0x0;
+  fRunLoader    = 0x0;
+  fMass         = 0;
+  fPdgC         = 0;
+}
+
+//____________________________________________________________________________
+
+AliJetKineReader::~AliJetKineReader()
+{
+  // Destructor
+  delete fReaderHeader;
+}
+
+//____________________________________________________________________________
+
+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);
+    
+  // 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));
+    }
+}
+
+//____________________________________________________________________________
+
+void AliJetKineReader::FillMomentumArray(Int_t event)
+{
+//
+// Fill momentum array for event
+//
+    Int_t goodTrack = 0;
+    // Clear array
+    ClearArray();
+
+    fRunLoader->GetEvent(event);
+    AliStack* stack = fRunLoader->Stack();
+    Int_t nt = stack->GetNprimary();
+    // Get cuts set by user
+    Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
+    fAliHeader = fRunLoader->GetHeader();
+        
+    // 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(); 
+       
+       if (
+           (status != 1)            
+           || (pdg == 12 || pdg == 14) 
+           || (pt < ptMin)             
+           ) continue; 
+
+       if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
+           Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+             if (
+                 (charge == 0)  
+                 || (gRandom->Rndm() < 0.1) 
+                 ) continue;
+       }
+       
+       fMass = part->GetCalcMass();
+       fPdgC = part->GetPdgCode();
+       // Fill momentum array
+       
+       new ((*fMomentumArray)[goodTrack]) 
+           TLorentzVector(part->Px(),part->Py(),part->Pz(), part->Energy());
+       goodTrack++;
+  }
+    printf("\nNumber of good tracks %d \n", goodTrack);
+}
+
+
diff --git a/JETAN/AliJetKineReader.h b/JETAN/AliJetKineReader.h
new file mode 100644 (file)
index 0000000..366868f
--- /dev/null
@@ -0,0 +1,38 @@
+#ifndef ALIJETKINEREADER_H
+#define ALIJETKINEREADER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+// Jet Kine Reader 
+// MC Kinematics reader for jet analysis
+// Author: Andreas Morsch (andreas.morsch@cern.ch)
+
+#include "AliJetReader.h"
+
+class AliRunLoader;
+
+class AliJetKineReader : public AliJetReader
+{
+ public: 
+  AliJetKineReader();
+  virtual ~AliJetKineReader();
+
+  // 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();
+
+ protected:
+  AliRunLoader           *fRunLoader; // Pointer to the run loader
+  
+  Float_t fMass;                      // Particle mass
+  Int_t   fPdgC;                      // Pdg code
+  ClassDef(AliJetKineReader,1)
+};
+#endif
diff --git a/JETAN/AliJetKineReaderHeader.cxx b/JETAN/AliJetKineReaderHeader.cxx
new file mode 100644 (file)
index 0000000..6d2682f
--- /dev/null
@@ -0,0 +1,39 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * 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"
+
+ClassImp(AliJetKineReaderHeader)
+
+//____________________________________________________________________________
+
+AliJetKineReaderHeader::AliJetKineReaderHeader():
+ AliJetReaderHeader("AliJetKineReaderHeader") 
+{
+  // Constructor
+  SetPtCut();
+  SetFastSimTPC(kFALSE);
+}
+
+//____________________________________________________________________________
+AliJetKineReaderHeader::~AliJetKineReaderHeader()
+{
+  // destructor
+}
diff --git a/JETAN/AliJetKineReaderHeader.h b/JETAN/AliJetKineReaderHeader.h
new file mode 100644 (file)
index 0000000..69caeaa
--- /dev/null
@@ -0,0 +1,34 @@
+#ifndef ALIJETKINEREADERHEADER_H
+#define ALIJETKINEREADERHEADER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+// Jet Kinematics Reader Header
+// Header for the Kine reader in the jet analysis
+// Author: Andreas Morsch (andreas.morsch@cern.ch)
+
+#include "AliJetReaderHeader.h"
+class AliJetKineReaderHeader : public AliJetReaderHeader
+{
+
+ public:
+  AliJetKineReaderHeader();
+  virtual ~AliJetKineReaderHeader();
+  
+  // Setters
+  void SetFastSimTPC(Bool_t flag = kTRUE) {fFastSimTPC = flag;}
+  virtual void SetPtCut(const Float_t par = 2.0) {fPtCut = par;}
+  // Getter
+  Bool_t  FastSimTPC() const  {return fFastSimTPC;}
+  Float_t GetPtCut()   const  {return fPtCut;}
+         
+ protected:
+  //parameters set by user
+  Bool_t fFastSimTPC;
+  Float_t fPtCut;
+  ClassDef(AliJetKineReaderHeader,1);
+};
+#endif
diff --git a/JETAN/AliJetMCReader.cxx b/JETAN/AliJetMCReader.cxx
new file mode 100644 (file)
index 0000000..173170c
--- /dev/null
@@ -0,0 +1,97 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+// Jet MC Reader 
+// MC reader for jet analysis
+// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+
+// From root ...
+#include <TClonesArray.h>
+#include <TPDGCode.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
+#include <TVector3.h>
+#include <TLorentzVector.h>
+#include <TSystem.h>
+// From AliRoot ...
+#include "AliJetMCReader.h"
+#include "AliJetMCReaderHeader.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+
+ClassImp(AliJetMCReader)
+
+AliJetMCReader::AliJetMCReader()
+{
+  // Constructor
+  fReaderHeader = 0x0;
+  fMass = 0;
+  fPdgC = 0;
+}
+
+//____________________________________________________________________________
+
+AliJetMCReader::~AliJetMCReader()
+{
+  // Destructor
+  delete fReaderHeader;
+}
+
+//____________________________________________________________________________
+
+
+void AliJetMCReader::FillMomentumArray(Int_t event)
+{
+// Fill momentum array
+  TClonesArray &arrayMC = *fArrayMC;
+  Int_t goodTrack = 0;
+  Int_t nt = 0;
+  Float_t pt, e;
+  TVector3 p;
+
+  // clear array
+  ClearArray();
+  // get event from chains
+  fChain->GetEntry(event);
+  fChainMC->GetEntry(event);
+  // get number of tracks in event (for the loop)
+  nt = fESD->GetNumberOfTracks();
+
+  // get cuts set by user
+  Double_t ptMin = ((AliJetMCReaderHeader*) fReaderHeader)->GetPtCut();
+
+  //loop over particles
+  for (Int_t it = 0; it < nt; it++) {
+    AliESDtrack *track = fESD->GetTrack(it); //track
+    UInt_t status = track->GetStatus();
+    if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
+    //    track->GetImpactParameters(dca,z);
+    // if (dca > dcaMax) continue; // check track is reasonable 
+    Int_t label = TMath::Abs(track->GetLabel());
+    TParticle *part = (TParticle*)arrayMC[label]; //particle
+    pt = part->Pt(); // pt of the particle
+    if (pt < ptMin) continue; //check  cuts 
+    p = part->P();
+    e = part->Energy();
+    fMass = part->GetCalcMass();
+    fPdgC = part->GetPdgCode();
+   // fill momentum array
+    new ((*fMomentumArray)[goodTrack]) TLorentzVector(p.X(), p.Y(), p.Z(), e);
+    goodTrack++;
+  }
+  printf("\nNumber of good tracks %d \n", goodTrack);
+}
+
+
diff --git a/JETAN/AliJetMCReader.h b/JETAN/AliJetMCReader.h
new file mode 100644 (file)
index 0000000..670d7fc
--- /dev/null
@@ -0,0 +1,25 @@
+#ifndef ALIJETMCREADER_H
+#define ALIJETMCREADER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+// Jet MC Reader 
+// MC reader for jet analysis
+// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+
+#include "AliJetESDReader.h"
+
+class AliJetMCReader : public AliJetESDReader
+{
+ public: 
+    AliJetMCReader();
+    virtual ~AliJetMCReader();
+    void FillMomentumArray(Int_t event);
+    
+ protected:
+    Float_t fPdgC;   // Pdg code
+    ClassDef(AliJetMCReader,1)
+};
+#endif
diff --git a/JETAN/AliJetMCReaderHeader.cxx b/JETAN/AliJetMCReaderHeader.cxx
new file mode 100644 (file)
index 0000000..0c149d1
--- /dev/null
@@ -0,0 +1,38 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+// Jet MC Reader Header
+// Header for the MC reader in the jet analysis
+// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+
+#include "AliJetMCReaderHeader.h"
+
+ClassImp(AliJetMCReaderHeader)
+
+//____________________________________________________________________________
+
+AliJetMCReaderHeader::AliJetMCReaderHeader():
+ AliJetReaderHeader("AliJetMCReaderHeader") 
+{
+  // Constructor
+  SetPtCut();
+}
+
+//____________________________________________________________________________
+AliJetMCReaderHeader::~AliJetMCReaderHeader()
+{
+  // destructor
+}
diff --git a/JETAN/AliJetMCReaderHeader.h b/JETAN/AliJetMCReaderHeader.h
new file mode 100644 (file)
index 0000000..8b198b8
--- /dev/null
@@ -0,0 +1,37 @@
+#ifndef ALIJETMCREADERHEADER_H
+#define ALIJETMCREADERHEADER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+// Jet MC Reader Header
+// Header for the MC reader in the jet analysis
+// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+
+#include "AliJetReaderHeader.h"
+class AliJetMCReaderHeader : public AliJetReaderHeader
+{
+
+ public:
+  AliJetMCReaderHeader();
+  virtual ~AliJetMCReaderHeader();
+  
+  // Getters
+  Double_t GetPtCut() const {return fPtCut;}
+  Float_t GetDCA() const {return fDCA;} // not working so far..(always 0)
+  
+  // Setters
+  virtual void SetPtCut(const Double_t par=2.0) {fPtCut = par;}
+  virtual void SetDCA(const Float_t dca = 3.0) {fDCA = dca;} // this is for the track!!! (recommended by P.H.)
+  
+ protected:
+  //parameters set by user
+  Double_t fPtCut;
+  Float_t fDCA;
+
+
+  ClassDef(AliJetMCReaderHeader,1);
+};
+#endif
diff --git a/JETAN/AliJetReader.cxx b/JETAN/AliJetReader.cxx
new file mode 100755 (executable)
index 0000000..d60fc23
--- /dev/null
@@ -0,0 +1,67 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * 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
+//---------------------------------------------------------------------
+
+#include <TClonesArray.h>
+
+#include "AliJetReader.h"
+#include "AliJetReaderHeader.h"
+#include "AliESD.h"
+#include "AliHeader.h"
+
+ClassImp(AliJetReader)
+
+////////////////////////////////////////////////////////////////////////
+
+AliJetReader::AliJetReader()
+{
+  // Constructor
+  fChain = 0; 
+  fChainMC = 0;
+  fESD = 0;
+  fMomentumArray = new TClonesArray("TLorentzVector",2000);
+  fArrayMC = 0;
+  fAliHeader = 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+AliJetReader::~AliJetReader()
+{
+  // Destructor
+  delete fChain;
+  delete fChainMC;
+  delete fESD;
+  delete fAliHeader;
+  if (fMomentumArray) {
+      fMomentumArray->Delete();
+      delete fMomentumArray;
+  }
+  delete fArrayMC;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetReader::ClearArray()
+
+{
+  if (fMomentumArray)  fMomentumArray->Clear();
+}
diff --git a/JETAN/AliJetReader.h b/JETAN/AliJetReader.h
new file mode 100755 (executable)
index 0000000..1dc85a0
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef ALIJETREADER_H
+#define ALIJETREADER_H
+/* 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>
+
+class TClonesArray;
+class AliJetReaderHeader;
+class AliESD;
+class AliHeader;
+
+
+class AliJetReader : public TObject 
+{
+ public: 
+  AliJetReader();
+  virtual ~AliJetReader();
+
+  // Getters
+  virtual TClonesArray *GetMomentumArray() {return fMomentumArray;}
+  virtual Int_t GetChainEntries() {return fChain->GetEntries();} 
+  virtual AliJetReaderHeader* GetReaderHeader() { return fReaderHeader;}
+  virtual AliHeader* GetAliHeader() { return fAliHeader;}
+
+  // Setters
+  virtual void FillMomentumArray(Int_t) {}
+  virtual void SetReaderHeader(AliJetReaderHeader* header) {fReaderHeader = header;}
+         
+  // others
+  virtual void OpenInputFiles() {}
+  void ClearArray();
+ protected:
+  TChain                  *fChain;         // chain for reconstructed tracks
+  TChain                  *fChainMC;       // chain for mc information
+  TClonesArray            *fMomentumArray; // array of particle momenta
+  TClonesArray            *fArrayMC;       // array of mc particles
+  AliESD                  *fESD;           // pointer to esd
+  AliJetReaderHeader      *fReaderHeader;  // pointer to header
+  AliHeader               *fAliHeader;     // pointer to event header
+  ClassDef(AliJetReader,1)
+};
+#endif
diff --git a/JETAN/AliJetReaderHeader.cxx b/JETAN/AliJetReaderHeader.cxx
new file mode 100755 (executable)
index 0000000..8521eac
--- /dev/null
@@ -0,0 +1,60 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//---------------------------------------------------------------------
+// base class for Jet Reader Header 
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include "AliJetReaderHeader.h"
+
+ClassImp(AliJetReaderHeader)
+
+////////////////////////////////////////////////////////////////////////
+
+AliJetReaderHeader::AliJetReaderHeader():  
+ TNamed("AliJetReaderHeader", "Jet Reader Header"),
+ fComment("No comment"),fDir(""),fPattern("")
+{
+  //
+  // Constructor
+  //
+  fFirst = 0;
+  fLast = -1;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+AliJetReaderHeader::AliJetReaderHeader(const char * name):  
+ TNamed(name, "Jet Reader Header"),
+ fComment("No comment"),fDir(""),fPattern("")
+{
+  //
+  // Constructor
+  //
+  fFirst = 0;
+  fLast = -1;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+AliJetReaderHeader::~AliJetReaderHeader()
+{
+  // destructor
+  delete fComment;
+  delete fDir;
+  delete fPattern;
+}
diff --git a/JETAN/AliJetReaderHeader.h b/JETAN/AliJetReaderHeader.h
new file mode 100755 (executable)
index 0000000..d458410
--- /dev/null
@@ -0,0 +1,49 @@
+#ifndef ALIJETREADERHEADER_H
+#define ALIJETREADERHEADER_H
+/* 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>
+class AliJetReaderHeader : public TNamed
+{
+
+ public:
+  AliJetReaderHeader(const char* name);
+  AliJetReaderHeader();
+  virtual ~AliJetReaderHeader();
+  
+  // Getters
+  virtual TString GetComment() {return fComment;}
+  virtual const char* GetDirectory() {return fDir.Data();}
+  virtual const char* GetPattern() {return fPattern.Data();}
+  Int_t   GetNEvents() const {return fLast-fFirst;}
+  Int_t   GetLastEvent() const {return fLast;}
+  Int_t   GetFirstEvent() const {return fFirst;}
+  
+  // Setters
+  virtual void SetComment(const char* s) {fComment=TString(s);}
+  virtual void SetPattern(const char* s) {fPattern=TString(s);}
+  virtual void SetDirectory(const char* s) {fDir=TString(s);}
+  virtual void SetFirstEvent(Int_t i=0) {fFirst=i;}
+  virtual void SetLastEvent(Int_t i=-1) {fLast=i;}
+  
+ protected:
+
+  Int_t fFirst;          // First and last events analyzed
+  Int_t fLast;           // in current set of files
+  TString fComment;      // a comment
+  TString fDir;          // directory with input files
+  TString fPattern;      // pattern to look for input files
+  
+  ClassDef(AliJetReaderHeader,1);
+};
+#endif
diff --git a/JETAN/AliLeading.cxx b/JETAN/AliLeading.cxx
new file mode 100644 (file)
index 0000000..0d9afe1
--- /dev/null
@@ -0,0 +1,125 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//---------------------------------------------------------------------
+// Class to find and store the leading particle in event and
+// store its correlation to associated particles
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TMath.h>
+#include <TClonesArray.h>
+#include <TLorentzVector.h>
+
+#include "AliLeading.h"
+#include <AliJetReader.h>
+
+ClassImp(AliLeading)
+
+////////////////////////////////////////////////////////////////////////
+
+AliLeading::AliLeading() 
+{
+  //
+  // Constructor
+  //
+  fNassoc=0;
+  fLeading = new TLorentzVector(0.,0.,0.,0.);
+  fLow = -TMath::Pi()/2.0;
+  fnBin=45;
+  fCorr = TArrayI(fnBin);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+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
+  //
+  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) {
+      ptMax = lv->Pt();
+      idxMax = i;
+    }
+  }
+
+  // fill correlation array
+  fLeading = (TLorentzVector*) lvArray->At(idxMax);
+  for (Int_t i=0; i<nIn;i++){
+    if (i == idxMax) continue;
+    TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+    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()));
+    fCorr.AddAt(fCorr.At(iBin)+1,iBin);
+  }
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliLeading::Reset()
+
+{
+// Reset leding particle information
+  fLeading->SetPxPyPzE(0.,0.,0.,0.);
+  fNassoc=0;
+  fCorr.Reset();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliLeading::PrintLeading()
+
+{
+// Print leading particle information
+  if (fNassoc<0) {
+    cout << " No leading particle in this event" << endl;
+    return;
+  }
+  cout << " Leading particle: " << endl;
+  cout << "    (px,py,pz,e) = (" << fLeading->Px() << ","
+       << fLeading->Py() << "," << fLeading->Pz() << ","
+       << fLeading->E() << ")" << endl;
+  cout << "    (pt,eta,phi) = (" << fLeading->Pt() << ","
+       << fLeading->Eta() << "," << fLeading->Phi() << ")" << endl;
+  cout << "    " << fNassoc << " associated particles." << endl;
+}
diff --git a/JETAN/AliLeading.h b/JETAN/AliLeading.h
new file mode 100644 (file)
index 0000000..784ee86
--- /dev/null
@@ -0,0 +1,51 @@
+#ifndef ALILEADING_H
+#define ALILEADING_H
+/* 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
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include <TObject.h>
+#include <TArrayI.h>
+
+class TLorentzVector;
+class AliJetReader;
+
+
+class AliLeading : public TObject
+{
+ public:
+
+  AliLeading();
+  ~AliLeading();
+
+  // Getters
+  Int_t GetNassoc() const { return fNassoc;}
+  TLorentzVector* GetLeading() const { return fLeading;}
+  TArrayI GetCorr() const { return fCorr;}
+
+  // Setters
+  // Others
+  void FindLeading(AliJetReader* reader);
+  void PrintLeading();
+  void Reset();
+
+ protected:
+
+  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
+
+  ClassDef(AliLeading,1);
+};
+
+#endif
diff --git a/JETAN/AliPxconeJetFinder.cxx b/JETAN/AliPxconeJetFinder.cxx
new file mode 100755 (executable)
index 0000000..becf007
--- /dev/null
@@ -0,0 +1,145 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+   
+  
+//---------------------------------------------------------------------
+// Pxcone Jet finder 
+// manages the search for jets 
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TLorentzVector.h>
+#include "AliPxconeJetFinder.h"
+#include "AliPxconeJetHeader.h"
+#include "AliJetReader.h"
+#include "AliJet.h"
+
+ClassImp(AliPxconeJetFinder)
+
+////////////////////////////////////////////////////////////////////////
+
+AliPxconeJetFinder::AliPxconeJetFinder()
+
+{
+  //
+  // Constructor
+  //
+  fHeader = 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+AliPxconeJetFinder::~AliPxconeJetFinder()
+
+{
+  //
+  // destructor
+  //
+
+  // reset and delete header
+}
+
+////////////////////////////////////////////////////////////////////////
+
+#ifndef WIN32
+# define pxcone pxcone_
+# define type_of_call
+
+#else
+# define pxcone PXCONE
+# define type_of_call _stdcall
+#endif
+
+extern "C" void type_of_call
+pxcone_(int* mode, int* ntrak, int* dim, double* ptrak,
+       double* coner, double* epslon, double* ovlim, 
+       int * maxjet, int* njet, double* pjet, int* ipass, 
+       int* ijmul, int* ierr);
+
+void AliPxconeJetFinder::FindJets()
+
+{
+  // get number of entries
+  // Int_t ne=fReader->GetNEvents();
+  // loops over entries (check user request for number of events)
+  // (this info should be stored in reader-header)
+
+  // test with one event
+  TClonesArray *lvArray = fReader->GetMomentumArray();
+  Int_t nIn = lvArray->GetEntries();
+  
+  // local arrays for input
+  const Int_t kNmaxVec = 30000;
+  if (nIn > kNmaxVec) {
+    cout << " AliPxconeJetFinder::FindJets: Too many input vectors."
+        << endl;
+    cout << " Using only the first " << kNmaxVec << endl;
+    nIn = kNmaxVec;
+  }
+  
+  Double_t pIn[kNmaxVec][4];
+  Double_t pJet[kNmaxVec][5];
+  int ipass[kNmaxVec];
+  int ijmul[kNmaxVec];
+  Int_t ierr;
+
+  // load input vectors
+  for (Int_t i=0; i<nIn;i++){
+    TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+    pIn[i][0]= lv->Px();
+    pIn[i][1]= lv->Py();
+    pIn[i][2]= lv->Pz();
+    pIn[i][3]= lv->E();
+  }
+  fJets->SetNinput(nIn);
+
+  // run the algorithm. Get parameters from header
+  Int_t dim=4;
+  Int_t mode=fHeader->GetMode();
+  Double_t radius=fHeader->GetRadius();
+  Double_t minpt=fHeader->GetMinPt();
+  Double_t ov=fHeader->GetOverlap();
+  Int_t nj;
+  pxcone(&mode,&nIn,&dim,&pIn[0][0],&radius,&minpt,&ov,&nIn,&nj,
+         &pJet[0][0],&ipass[0],&ijmul[0],&ierr);
+
+  // download jets
+  fJets->SetInJet(ipass);
+  for(Int_t i=0; i<nj; i++) 
+    fJets->AddJet(pJet[i][0],pJet[i][1],
+                 pJet[i][2],pJet[i][3]);
+  fJets->SetMultiplicities(ijmul);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliPxconeJetFinder::WriteJHeaderToFile()
+{
+// Write Header to file
+
+  fOut->cd();
+  fHeader->Write();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliPxconeJetFinder::Reset()
+{
+// Reset jet list
+  fJets->ClearJets();
+}
+
diff --git a/JETAN/AliPxconeJetFinder.h b/JETAN/AliPxconeJetFinder.h
new file mode 100755 (executable)
index 0000000..4ab8749
--- /dev/null
@@ -0,0 +1,41 @@
+#ifndef ALIPXCONEJETFINDER_H
+#define ALIPXCONEJETFINDER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Pxcone Jet finder 
+// manages the search for jets 
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include "AliJetFinder.h"
+
+class AliPxconeJetHeader;
+
+class AliPxconeJetFinder : public AliJetFinder 
+{
+ public:
+
+  AliPxconeJetFinder();
+  ~AliPxconeJetFinder();
+
+  // getters
+
+  // setters
+  void SetJetHeader(AliPxconeJetHeader* h) {fHeader= h;}
+  // others
+  void Reset();
+  void FindJets();
+  void WriteJHeaderToFile();
+
+ protected:
+
+  AliPxconeJetHeader* fHeader;         // pointer to jet header
+
+  ClassDef(AliPxconeJetFinder,1)
+};
+
+#endif
diff --git a/JETAN/AliPxconeJetHeader.cxx b/JETAN/AliPxconeJetHeader.cxx
new file mode 100755 (executable)
index 0000000..eb99684
--- /dev/null
@@ -0,0 +1,56 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//---------------------------------------------------------------------
+// Jet header base class 
+// Stores a comment which describes the jet analysis
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include <Riostream.h> 
+#include "AliPxconeJetHeader.h"
+ClassImp(AliPxconeJetHeader)
+////////////////////////////////////////////////////////////////////////
+
+AliPxconeJetHeader::AliPxconeJetHeader():
+  AliJetHeader("AliPxconeJetHeader")
+{
+  //
+  // Constructor
+  //
+  SetMode();
+  SetRadius();
+  SetMinPt();
+  SetOverlap();
+}
+
+////////////////////////////////////////////////////////////////////////
+void AliPxconeJetHeader::PrintParameters()
+
+{
+  //
+  // prints out parameters of jet algorithm
+  //
+
+  cout << " PXCONE jet algorithm " << endl;
+  cout << "  Running mode: " << fMode << endl;
+  cout << "  Cone size: " << fRadius << endl;
+  cout << "  Minimum jet energy: " << fMinPt << endl;
+  cout << "  Overlap fraction: " << fOverlap << endl;
+}
diff --git a/JETAN/AliPxconeJetHeader.h b/JETAN/AliPxconeJetHeader.h
new file mode 100755 (executable)
index 0000000..007b0e6
--- /dev/null
@@ -0,0 +1,48 @@
+#ifndef ALIPXCONEJETHEADER_H
+#define ALIPXCONEJETHEADER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Jet header class for Pxcone algorithm
+// Stores the parameters of the Pxcone jet algorithm
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+#include "AliJetHeader.h"
+class AliPxconeJetHeader : public AliJetHeader
+{
+ public:
+  AliPxconeJetHeader();
+  virtual ~AliPxconeJetHeader() { }
+
+  // Getters
+
+  Int_t    GetMode()    const {return fMode;}
+  Double_t GetRadius()  const {return fRadius;}
+  Double_t GetMinPt()   const {return fMinPt;}
+  Double_t GetOverlap() const {return fOverlap;}
+   
+  // Setters
+  void SetMode(Int_t m=2) {fMode=m;}
+  void SetRadius(Double_t r=0.3) {fRadius=r;}
+  void SetMinPt(Double_t p=10.0) {fMinPt=p;}
+  void SetOverlap(Double_t o=.75) {fOverlap=o;}
+
+  // others
+  void PrintParameters();
+   
+protected:
+  Int_t fMode;           // ee or pp mode
+  Double_t fRadius;      // jet radius
+  Double_t fMinPt;       // min pt of jets  
+  Double_t fOverlap;     // fraction of overlap energy
+
+  ClassDef(AliPxconeJetHeader,1)
+};
+#endif
diff --git a/JETAN/AliUA1JetFinder.cxx b/JETAN/AliUA1JetFinder.cxx
new file mode 100755 (executable)
index 0000000..946bd4d
--- /dev/null
@@ -0,0 +1,230 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+   
+  
+//---------------------------------------------------------------------
+// UA1 Jet finder 
+// manages the search for jets 
+// Author: jgcn@mda.cinvestav.mx
+// (code adapted from EMCAL directory)
+//---------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TLorentzVector.h>
+#include <TH2F.h>
+#include "AliUA1JetFinder.h"
+#include "AliUA1JetHeader.h"
+#include "UA1Common.h"
+#include "AliJetReader.h"
+#include "AliJet.h"
+
+
+ClassImp(AliUA1JetFinder)
+
+////////////////////////////////////////////////////////////////////////
+
+AliUA1JetFinder::AliUA1JetFinder()
+
+{
+  //
+  // Constructor
+  //
+  fHeader = 0x0;
+  fLego   = 0x0;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+AliUA1JetFinder::~AliUA1JetFinder()
+
+{
+  //
+  // destructor
+  //
+
+  delete fLego;            
+  
+  // reset and delete header
+}
+
+////////////////////////////////////////////////////////////////////////
+
+#ifndef WIN32
+# define ua1_jet_finder ua1_jet_finder_
+# define hf1 hf1_
+# define type_of_call
+
+#else
+# define ua1_jet_finder UA1_JET_FINDER
+# define hf1 HF1
+# define type_of_call _stdcall
+#endif
+
+extern "C" void type_of_call
+ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
+               Float_t  etc[60000],  Float_t etac[60000],
+               Float_t  phic[60000],
+               Float_t& min_move, Float_t& max_move, Int_t& mode,
+               Float_t& prec_bg,  Int_t& ierror);
+
+extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
+
+////////////////////////////////////////////////////////////////////////
+
+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
+  
+  // 1) transform input to pt,eta,phi plus lego
+  TClonesArray *lvArray = fReader->GetMomentumArray();
+  Int_t nIn =  lvArray->GetEntries();
+  if (nIn == 0) return;
+    
+  // local arrays for input
+  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());
+      fLego->Fill(etaT[i], phiT[i], ptT[i]);
+  }
+  fJets->SetNinput(nIn);
+  
+  // 2) dump lego
+  // check enough space! *to be done*
+  Float_t etCell[60000];   //! Cell Energy
+  Float_t etaCell[60000];  //! Cell eta
+  Float_t phiCell[60000];  //! Cell phi
+
+  Int_t nCell = 0;
+  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++) {
+         e = fLego->GetBinContent(i,j);
+         if (e > 0.0) e -= fHeader->GetMinCellEt();
+         if (e < 0.0) e = 0.;
+         Float_t eta  = xaxis->GetBinCenter(i);
+         Float_t phi  = yaxis->GetBinCenter(j);            
+         etCell[nCell]  = e;
+         etaCell[nCell] = eta;
+         phiCell[nCell] = phi;
+         nCell++;
+      }
+  }
+
+  // run the algo. Parameters from header
+  Int_t nTot      = (fHeader->GetNbinEta())*(fHeader->GetNbinPhi());
+  Float_t minmove = fHeader->GetMinMove();
+  Float_t maxmove = fHeader->GetMaxMove();
+  Int_t mode      = fHeader->GetMode();
+  Float_t precbg  = fHeader->GetPrecBg();
+  Int_t ierror;
+  ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell, 
+                minmove, maxmove, mode, precbg, ierror);
+
+  // download jet info. 
+  Int_t* mult  = new Int_t[UA1JETS.njet];
+  // sort jets
+  Int_t * idx  = new Int_t[UA1JETS.njet];
+  TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
+  
+  Int_t* injet = new Int_t[nIn];
+  for (Int_t i = 0; i < nIn; i++) injet[i]= -1;
+
+  for(Int_t i = 0; i < UA1JETS.njet; i++) {
+      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
+      mult[i] = 0;
+      for (Int_t j = 0; j < nIn; j++) {
+         Float_t deta = etaT[j] - UA1JETS.etaj[1][idx[i]];
+         Float_t dphi = phiT[j] - UA1JETS.phij[1][idx[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) {
+             injet[j] = i;
+             mult[i]++;
+         }
+      }
+  }
+  fJets->SetMultiplicities(mult);
+  fJets->SetInJet(injet);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliUA1JetFinder::Reset()
+
+{
+  fLego->Reset();
+  fJets->ClearJets();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliUA1JetFinder::WriteJHeaderToFile()
+{
+  fOut->cd();
+  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());
+
+  UA1CELL.etaCellSize = dEta;
+  UA1CELL.phiCellSize = dPhi;
+
+  // jet parameters
+  UA1PARA.coneRad = fHeader->GetRadius();
+  UA1PARA.etSeed  = fHeader->GetEtSeed();
+  UA1PARA.ejMin   = fHeader->GetMinJetEt();
+  UA1PARA.etMin   = fHeader->GetMinCellEt();
+
+  // book lego
+  fLego = new 
+    TH2F("legoH","eta-phi",
+        fHeader->GetNbinEta(), fHeader->GetEtaMin(), fHeader->GetEtaMax(),
+        fHeader->GetNbinPhi(), fHeader->GetPhiMin(), fHeader->GetPhiMax());
+}
diff --git a/JETAN/AliUA1JetFinder.h b/JETAN/AliUA1JetFinder.h
new file mode 100755 (executable)
index 0000000..55fa3bd
--- /dev/null
@@ -0,0 +1,44 @@
+#ifndef ALIUA1JETFINDER_H
+#define ALIUA1JETFINDER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// UA1 Jet finder 
+// manages the search for jets 
+// Author: jgcn@mda.cinvestav.mx
+// (code adapted from EMCAL directory)
+//---------------------------------------------------------------------
+
+#include "AliJetFinder.h"
+class AliUA1JetHeader;
+class TH2F;
+
+class AliUA1JetFinder : public AliJetFinder 
+{
+ public:
+
+  AliUA1JetFinder();
+  ~AliUA1JetFinder();
+
+  // getters
+
+  // setters
+  void SetJetHeader(AliUA1JetHeader* h) {fHeader= h;}
+  // others
+  void FindJets();
+  void Reset();
+  void Init();
+  void WriteJHeaderToFile();
+
+ protected:
+
+  AliUA1JetHeader* fHeader;         // pointer to jet header
+  TH2F           * fLego;           //! Lego Histo
+
+  ClassDef(AliUA1JetFinder,1)
+};
+
+#endif
diff --git a/JETAN/AliUA1JetHeader.cxx b/JETAN/AliUA1JetHeader.cxx
new file mode 100755 (executable)
index 0000000..03566ca
--- /dev/null
@@ -0,0 +1,80 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//---------------------------------------------------------------------
+// UA1 Jet header class 
+// Stores parameters of UA1 jet algo
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+
+#include <Riostream.h> 
+#include <TMath.h>
+#include "AliUA1JetHeader.h"
+ClassImp(AliUA1JetHeader)
+////////////////////////////////////////////////////////////////////////
+
+AliUA1JetHeader::AliUA1JetHeader():
+  AliJetHeader("AliUA1JetHeader")
+{
+  //
+  // Constructor
+  //
+  fConeRadius =  0.3;
+  fEtSeed     =  3.0;
+  fMinJetEt   = 10.0;
+  fMinCellEt  =  0.0;
+  fMode       =  1;
+  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;
+}
+
+////////////////////////////////////////////////////////////////////////
+void AliUA1JetHeader::PrintParameters() const
+
+{
+  //
+  // prints out parameters of jet algorithm
+  //
+
+  cout << " UA1 jet algorithm " << endl;
+  cout << " * Jet parameters: " << endl;
+  cout << "   Cone size: " << fConeRadius<< endl;
+  cout << "   Minimum energy for a seed: " << fEtSeed << endl;
+  cout << "   Minumum energy for a jet: " << fMinJetEt << endl;
+  cout << "   Minumum energy for a cell: " << fMinCellEt << endl;
+  cout << " * Background substraction parameters: " << endl;
+  cout << "   Substraction mode: " << fMode << endl;
+  cout << "   Minimum allowed move: " << fMinMove << endl;
+  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;
+}
diff --git a/JETAN/AliUA1JetHeader.h b/JETAN/AliUA1JetHeader.h
new file mode 100755 (executable)
index 0000000..62eb5aa
--- /dev/null
@@ -0,0 +1,81 @@
+#ifndef ALIUA1JETHEADER_H
+#define ALIUA1JETHEADER_H
+/* 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
+// Author: jgcn@mda.cinvestav.mx
+//---------------------------------------------------------------------
+#include "AliJetHeader.h"
+class AliUA1JetHeader : public AliJetHeader
+{
+ public:
+  AliUA1JetHeader();
+  virtual ~AliUA1JetHeader() { }
+
+  // Getters
+  Float_t GetRadius()    const  {return fConeRadius;}
+  Float_t GetEtSeed()    const  {return fEtSeed;}
+  Float_t GetMinJetEt()  const  {return fMinJetEt;}
+  Float_t GetMinCellEt() const  {return fMinCellEt;}
+  Int_t   GetMode()      const  {return fMode;}
+  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;}
+
+  // Setters
+  void SetRadius(Float_t f) {fConeRadius=f;}
+  void SetEtSeed(Float_t f) {fEtSeed=f;}
+  void SetMinJetEt(Float_t f) {fMinJetEt=f;}
+  void SetMinCellEt(Float_t f) {fMinCellEt=f;}
+  void SetMode(Int_t f) {fMode=f;}
+  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;}
+
+  // others
+  void PrintParameters() const; 
+
+protected:
+
+  // parameters of algorithm
+  Float_t fConeRadius;      //  Cone radius
+  Float_t fEtSeed;          //  Min. Et for seed
+  Float_t fMinJetEt;        //  Min Et of jet
+  Float_t fMinCellEt;       //  Min Et in one cell
+  // parameters of backgound substraction
+  Int_t   fMode;            // key for BG subtraction
+  Float_t fMinMove;         // min cone move 
+  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
+
+  ClassDef(AliUA1JetHeader,1)
+};
+#endif
diff --git a/JETAN/JetAnalysisLinkDef.h b/JETAN/JetAnalysisLinkDef.h
new file mode 100644 (file)
index 0000000..be979f9
--- /dev/null
@@ -0,0 +1,25 @@
+#ifdef __CINT__
+  
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+  
+#pragma link C++ class AliJet+;
+#pragma link C++ class AliJetHeader+;
+#pragma link C++ class AliPxconeJetHeader+;
+#pragma link C++ class AliUA1JetHeader+;
+#pragma link C++ class AliJetFinder+;
+#pragma link C++ class AliPxconeJetFinder+;
+#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+;
+#pragma link C++ class AliLeading+;
+
+#endif
diff --git a/JETAN/UA1Common.h b/JETAN/UA1Common.h
new file mode 100755 (executable)
index 0000000..e7f332d
--- /dev/null
@@ -0,0 +1,47 @@
+#ifndef ROOT_UA1Common
+#define ROOT_UA1Common
+
+#ifndef __CFORTRAN_LOADED
+#include "cfortran.h"
+#endif
+
+extern "C" {
+// COMMON /UA1JETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2), NCELLJ(100)
+
+    typedef struct {
+       Int_t   njet;
+       Float_t etj[100];
+       Float_t etaj[2][100];
+       Float_t phij[2][100];
+       Int_t   ncellj[100];
+    } UA1JetsCommon;
+
+#define UA1JETS COMMON_BLOCK(UA1JETS,ua1jets)
+COMMON_BLOCK_DEF(UA1JetsCommon,UA1JETS);
+
+// COMMON /UA1CELL/ etaCellSize, phiCellSize
+
+    typedef struct {
+       Float_t etaCellSize;
+       Float_t phiCellSize;
+    } UA1CellCommon;
+
+#define UA1CELL COMMON_BLOCK(UA1CELL,ua1cell)
+COMMON_BLOCK_DEF(UA1CellCommon,UA1CELL);
+
+// COMMON /UA1PARA/ cone_rad, et_seed, ej_min, et_min
+    typedef struct {
+       Float_t coneRad;
+       Float_t etSeed;
+       Float_t ejMin;
+       Float_t etMin;
+    } UA1ParaCommon;
+
+#define UA1PARA COMMON_BLOCK(UA1PARA,ua1para)
+COMMON_BLOCK_DEF(UA1ParaCommon,UA1PARA);
+    
+}
+#endif
+
+
+
diff --git a/JETAN/libJETAN.pkg b/JETAN/libJETAN.pkg
new file mode 100644 (file)
index 0000000..5e73a0b
--- /dev/null
@@ -0,0 +1,16 @@
+#-*-Mode: Makefile-*-
+
+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 \
+       AliJetKineReader.cxx AliJetKineReaderHeader.cxx 
+
+FSRCS = pxcone.F ua1_jet_finder.F
+
+HDRS:= $(SRCS:.cxx=.h)
+
+DHDR= JetAnalysisLinkDef.h
+
+EINCLUDE:= STEER PYTHIA6
diff --git a/JETAN/pxcone.F b/JETAN/pxcone.F
new file mode 100755 (executable)
index 0000000..b310121
--- /dev/null
@@ -0,0 +1,867 @@
+      SUBROUTINE PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,
+     +                   MXJET,NJET,PJET,IPASS,IJMUL,IERR)
+*.*********************************************************
+*. ------
+*. PXCONE
+*. ------
+*.
+*.********** Pre Release Version 26.2.93
+*.
+*. Driver for the Cone  Jet finding algorithm of L.A. del Pozo.
+*. Based on algorithm from D.E. Soper.
+*. Finds jets inside cone of half angle CONER with energy > EPSLON.
+*. Jets which receive more than a fraction OVLIM of their energy from
+*. overlaps with other jets are excluded.
+*. Output jets are ordered in energy.
+*. If MODE.EQ.2 momenta are stored as (eta,phi,<empty>,pt)
+*. Usage     :
+*.
+*.      INTEGER  ITKDM,MXTRK
+*.      PARAMETER  (ITKDM=4.or.more,MXTRK=1.or.more)
+*.      INTEGER  MXJET, MXTRAK, MXPROT
+*.      PARAMETER  (MXJET=10,MXTRAK=500,MXPROT=500)
+*.      INTEGER  IPASS (MXTRAK),IJMUL (MXJET)
+*.      INTEGER  NTRAK,NJET,IERR,MODE
+*.      DOUBLE PRECISION  PTRAK (ITKDM,MXTRK),PJET (5,MXJET)
+*.      DOUBLE PRECISION  CONER, EPSLON, OVLIM
+*.      NTRAK = 1.to.MXTRAK
+*.      CONER   = ...
+*.      EPSLON  = ...
+*.      OVLIM   = ...
+*.      CALL PXCONE (MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET,
+*.     +             NJET,PJET,IPASS,IJMUL,IERR)
+*.
+*. INPUT     :  MODE      1=>e+e-, 2=>hadron-hadron
+*. INPUT     :  NTRAK     Number of particles
+*. INPUT     :  ITKDM     First dimension of PTRAK array
+*. INPUT     :  PTRAK     Array of particle 4-momenta (Px,Py,Pz,E)
+*. INPUT     :  CONER     Cone size (half angle) in radians
+*. INPUT     :  EPSLON    Minimum Jet energy (GeV)
+*. INPUT     :  OVLIM     Maximum fraction of overlap energy in a jet
+*. INPUT     :  MXJET     Maximum possible number of jets
+*. OUTPUT    :  NJET      Number of jets found
+*. OUTPUT    :  PJET      5-vectors of jets
+*. OUTPUT    :  IPASS(k)  Particle k belongs to jet number IPASS(k)
+*.                        IPASS = -1 if not assosciated to a jet
+*. OUTPUT    :  IJMUL(i)  Jet i contains IJMUL(i) particles
+*. OUTPUT    :  IERR      = 0 if all is OK ;   = -1 otherwise
+*.
+*. CALLS     : PXSEAR, PXSAME, PXNEW, PXTRY, PXORD, PXUVEC, PXOLAP
+*. CALLED    : User
+*.
+*. AUTHOR    :  L.A. del Pozo
+*. CREATED   :  26-Feb-93
+*. LAST MOD  :   2-Mar-93
+*.
+*. Modification Log.
+*. 2-Jan-97: M Wobisch    - fix bug concerning COS2R in eta phi mode
+*. 4-Apr-93: M H Seymour  - Change 2d arrays to 1d in PXTRY & PXNEW
+*. 2-Apr-93: M H Seymour  - Major changes to add boost-invariant mode
+*. 1-Apr-93: M H Seymour  - Increase all array sizes
+*. 30-Mar-93: M H Seymour - Change all REAL variables to DOUBLE PRECISION
+*. 30-Mar-93: M H Seymour - Change OVLIM into an input parameter
+*. 2-Mar-93: L A del Pozo - Fix Bugs in PXOLAP
+*. 1-Mar-93: L A del Pozo - Remove Cern library routine calls
+*. 1-Mar-93: L A del Pozo - Add Print out of welcome and R and Epsilon
+*.
+*.*********************************************************
+C+SEQ,DECLARE.
+*** External Arrays
+      INTEGER  ITKDM,MXJET,NTRAK,NJET,IERR,MODE
+      INTEGER  IPASS (NTRAK),IJMUL (MXJET)
+      DOUBLE PRECISION  PTRAK (ITKDM,NTRAK),PJET (5,MXJET), 
+     +     CONER, EPSLON, OVLIM
+*** Internal Arrays
+      INTEGER MXPROT, MXTRAK
+      PARAMETER (MXPROT=5000, MXTRAK=5000)
+      DOUBLE PRECISION PP(4,MXTRAK), PU(3,MXTRAK), PJ(4,MXPROT)
+      LOGICAL JETLIS(MXPROT,MXTRAK)
+*** Used in the routine.
+      DOUBLE PRECISION COSR,COS2R, VSEED(3), VEC1(3), VEC2(3),PTSQ,PPSQ,
+     +     COSVAL,PXMDPI
+cMWobisch
+      DOUBLE PRECISION RSEP
+cMWobisch
+      LOGICAL UNSTBL
+      INTEGER I,J,N,MU,N1,N2, ITERR
+      INTEGER NCALL, NPRINT
+      DOUBLE PRECISION ROLD, EPSOLD, OVOLD
+      SAVE NCALL,NPRINT,ROLD, EPSOLD, OVOLD
+      DATA NCALL,NPRINT /0,0/
+      DATA ROLD,EPSOLD,OVOLD/0.,0.,0./
+
+cMWobisch
+c***************************************
+      RSEP  = 2D0
+c***************************************
+cMWobisch
+      IERR=0
+*
+*** INITIALIZE
+      IF(NCALL.LE.0)  THEN
+         ROLD = 0.
+         EPSOLD = 0.
+         OVOLD = 0.
+      ENDIF
+      NCALL = NCALL + 1
+*
+*** Print welcome and Jetfinder parameters
+      IF((CONER.NE.ROLD .OR. EPSLON.NE.EPSOLD .OR. OVLIM.NE.OVOLD)
+     +     .AND. NPRINT.LE.10) THEN
+         WRITE (6,*)
+         WRITE (6,*) ' *********** PXCONE: Cone Jet-finder ***********'
+         WRITE (6,*) '    Written by Luis Del Pozo of OPAL'
+         WRITE (6,*) '    Modified for eta-phi by Mike Seymour'
+         WRITE(6,1000)'   Cone Size R = ',CONER,' Radians'
+         WRITE(6,1001)'   Min Jet energy Epsilon = ',EPSLON,' GeV'
+         WRITE(6,1002)'   Overlap fraction parameter = ',OVLIM
+         WRITE (6,*) ' ***********************************************'
+cMWobisch
+         IF (RSEP .lt. 1.999) THEN
+            WRITE(6,*) ' '
+            WRITE (6,*) ' ******************************************'
+            WRITE (6,*) ' ******************************************'
+            WRITE(6,*) ' M Wobisch: private change !!!!!!!!!!!! '
+            WRITE(6,*) '      Rsep is set to ',RSEP
+            WRITE(6,*) ' this is ONLY meaningful in a NLO calculation'
+            WRITE(6,*) '      ------------------------  '
+            WRITE(6,*) '  please check what you''re doing!!'
+            WRITE(6,*) '   or ask:  Markus.Wobisch@desy.de --'
+            WRITE (6,*) ' ******************************************'
+            WRITE (6,*) ' ******************************************'
+            WRITE (6,*) ' ******************************************'
+            WRITE(6,*) ' '
+            WRITE(6,*) ' '
+         ENDIF
+cMWobisch
+
+         WRITE (6,*)
+1000     FORMAT(A18,F5.2,A10)
+1001     FORMAT(A29,F5.2,A5)
+1002     FORMAT(A33,F5.2)
+         NPRINT = NPRINT + 1
+         ROLD=CONER
+         EPSOLD=EPSLON
+         OVOLD=OVLIM
+      ENDIF
+*
+*** Copy calling array PTRAK  to internal array PP(4,NTRAK)
+*
+      IF (NTRAK .GT. MXTRAK) THEN
+         WRITE (6,*) ' PXCONE: Ntrak too large: ',NTRAK
+         IERR=-1
+         RETURN
+      ENDIF
+      IF (MODE.NE.2) THEN
+         DO  100 I=1, NTRAK
+            DO  101 J=1,4
+               PP(J,I)=PTRAK(J,I)
+101         CONTINUE
+100      CONTINUE
+      ELSE
+*** Converting to eta,phi,pt if necessary
+         DO  104 I=1,NTRAK
+            PTSQ=PTRAK(1,I)**2+PTRAK(2,I)**2
+            PPSQ=(SQRT(PTSQ+PTRAK(3,I)**2)+ABS(PTRAK(3,I)))**2
+            IF (PTSQ.LE.4.25E-18*PPSQ) THEN
+               PP(1,I)=20
+            ELSE
+               PP(1,I)=0.5*LOG(PPSQ/PTSQ)
+            ENDIF
+            PP(1,I)=SIGN(PP(1,I),PTRAK(3,I))
+            IF (PTSQ.EQ.0) THEN
+               PP(2,I)=0
+            ELSE
+               PP(2,I)=ATAN2(PTRAK(2,I),PTRAK(1,I))
+            ENDIF
+            PP(3,I)=0
+            PP(4,I)=SQRT(PTSQ)
+            PU(1,I)=PP(1,I)
+            PU(2,I)=PP(2,I)
+            PU(3,I)=PP(3,I)
+104      CONTINUE
+      ENDIF
+*
+*** Zero output variables
+*
+      NJET=0
+      DO 102 I = 1, NTRAK
+         DO 103 J = 1, MXPROT
+           JETLIS(J,I) = .FALSE.
+103      CONTINUE
+102   CONTINUE
+      CALL PXZERV(4*MXPROT,PJ)
+      CALL PXZERI(MXJET,IJMUL)
+*
+      IF (MODE.NE.2) THEN
+         COSR = COS(CONER)
+         COS2R = COS(CONER)
+      ELSE
+*** Purely for convenience, work in terms of 1-R**2
+         COSR = 1-CONER**2
+cMW -- select Rsep: 1-(Rsep*CONER)**2
+         COS2R =  1-(RSEP*CONER)**2
+cORIGINAL         COS2R =  1-(2*CONER)**2
+      ENDIF
+      UNSTBL = .FALSE.
+      IF (MODE.NE.2) THEN
+         CALL PXUVEC(NTRAK,PP,PU,IERR)
+         IF (IERR .NE. 0) RETURN
+      ENDIF
+*** Look for jets using particle diretions as seed axes
+*
+      DO 110 N = 1,NTRAK
+        DO 120 MU = 1,3
+          VSEED(MU) = PU(MU,N)
+120     CONTINUE
+        CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,
+     &                   NJET,JETLIS,PJ,UNSTBL,IERR)
+         IF (IERR .NE. 0) RETURN
+110   CONTINUE
+
+cMW - for Rsep=1 goto 145
+c      GOTO 145
+
+*** Now look between all pairs of jets as seed axes.
+      DO 140 N1 = 1,NJET-1
+         VEC1(1)=PJ(1,N1)
+         VEC1(2)=PJ(2,N1)
+         VEC1(3)=PJ(3,N1)
+         IF (MODE.NE.2) CALL PXNORV(3,VEC1,VEC1,ITERR)
+         DO 150 N2 = N1+1,NJET
+            VEC2(1)=PJ(1,N2)
+            VEC2(2)=PJ(2,N2)
+            VEC2(3)=PJ(3,N2)
+            IF (MODE.NE.2) CALL PXNORV(3,VEC2,VEC2,ITERR)
+            CALL PXADDV(3,VEC1,VEC2,VSEED,ITERR)
+            IF (MODE.NE.2) THEN
+               CALL PXNORV(3,VSEED,VSEED,ITERR)
+            ELSE
+               VSEED(1)=VSEED(1)/2
+               VSEED(2)=VSEED(2)/2
+            ENDIF
+C---ONLY BOTHER IF THEY ARE BETWEEN 1 AND 2 CONE RADII APART
+            IF (MODE.NE.2) THEN
+              COSVAL=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
+            ELSE
+               IF (ABS(VEC1(1)).GE.20.OR.ABS(VEC2(1)).GE.20) THEN
+                  COSVAL=-1000
+               ELSE
+                  COSVAL=1-
+     +              ((VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2)
+               ENDIF
+            ENDIF
+
+            IF (COSVAL.LE.COSR.AND.COSVAL.GE.COS2R)
+     +           CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
+     +           JETLIS,PJ,UNSTBL,IERR)
+c            CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
+c     +           JETLIS,PJ,UNSTBL,IERR)
+            IF (IERR .NE. 0) RETURN
+150      CONTINUE
+140   CONTINUE
+      IF (UNSTBL) THEN
+        IERR=-1
+        WRITE (6,*) ' PXCONE: Too many iterations to find a proto-jet'
+        RETURN
+      ENDIF
+
+ 145  CONTINUE
+*** Now put the jet list into order by jet energy, eliminating jets
+*** with energy less than EPSLON.
+       CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
+*
+*** Take care of jet overlaps
+       CALL PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
+*
+*** Order jets again as some have been eliminated, or lost energy.
+       CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
+*
+*** All done!, Copy output into output arrays
+      IF (NJET .GT. MXJET) THEN
+         WRITE (6,*) ' PXCONE:  Found more than MXJET jets'
+         IERR=-1
+         GOTO 99
+      ENDIF
+      IF (MODE.NE.2) THEN
+         DO 300 I=1, NJET
+            DO 310 J=1,4
+               PJET(J,I)=PJ(J,I)
+310         CONTINUE
+300      CONTINUE
+      ELSE
+         DO 315 I=1, NJET
+            PJET(1,I)=PJ(4,I)*COS(PJ(2,I))
+            PJET(2,I)=PJ(4,I)*SIN(PJ(2,I))
+            PJET(3,I)=PJ(4,I)*SINH(PJ(1,I))
+            PJET(4,I)=PJ(4,I)*COSH(PJ(1,I))
+ 315     CONTINUE
+      ENDIF
+      DO 320 I=1, NTRAK
+         IPASS(I)=-1
+         DO 330 J=1, NJET
+            IF (JETLIS(J,I)) THEN
+               IJMUL(J)=IJMUL(J)+1
+               IPASS(I)=J
+            ENDIF
+330      CONTINUE
+320   CONTINUE
+99    RETURN
+      END
+*CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+      SUBROUTINE PXNORV(N,A,B,ITERR)
+      INTEGER I,N,ITERR
+      DOUBLE PRECISION A(N),B(N),C
+      C=0
+      DO 10 I=1,N
+        C=C+A(I)**2
+ 10   CONTINUE
+      IF (C.LE.0) RETURN
+      C=1/SQRT(C)
+      DO 20 I=1,N
+        B(I)=A(I)*C
+ 20   CONTINUE
+      END
+*CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
+*CMZ :  1.06/00 15/03/94  12.17.46  by  P. Schleper
+*-- Author :
+*
+C+DECK,PXOLAP.
+      SUBROUTINE PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
+*
+*** Looks for particles assigned to more than 1 jet, and reassigns them
+*** If more than a fraction OVLIM of a jet's energy is contained in
+*** higher energy jets, that jet is neglected.
+*** Particles assigned to the jet closest in angle (a la CDF, Snowmass).
+C+SEQ,DECLARE.
+      INTEGER MXTRAK, MXPROT
+      PARAMETER (MXTRAK=5000,MXPROT=5000)
+      INTEGER NJET, NTRAK, MODE
+      LOGICAL JETLIS(MXPROT,MXTRAK)
+      DOUBLE PRECISION PJ(4,MXPROT),PP(4,MXTRAK),PXMDPI
+      INTEGER I,J,N,MU
+      LOGICAL OVELAP
+      DOUBLE PRECISION EOVER
+      DOUBLE PRECISION OVLIM
+      INTEGER ITERR, IJMIN, IJET(MXPROT), NJ
+      DOUBLE PRECISION VEC1(3), VEC2(3), COST, THET, THMIN
+      DATA IJMIN/0/
+*
+      IF (NJET.LE.1) RETURN
+*** Look for jets with large overlaps with higher energy jets.
+      DO 100 I = 2,NJET
+*** Find overlap energy between jets I and all higher energy jets.
+       EOVER = 0.0
+       DO 110 N = 1,NTRAK
+         OVELAP = .FALSE.
+         DO 120 J= 1,I-1
+           IF (JETLIS(I,N).AND.JETLIS(J,N)) THEN
+            OVELAP = .TRUE.
+           ENDIF
+120      CONTINUE
+         IF (OVELAP) THEN
+           EOVER = EOVER + PP(4,N)
+         ENDIF
+110     CONTINUE
+*** Is the fraction of energy shared larger than OVLIM?
+        IF (EOVER.GT.OVLIM*PJ(4,I)) THEN
+*** De-assign all particles from Jet I
+            DO 130 N = 1,NTRAK
+              JETLIS(I,N) = .FALSE.
+130         CONTINUE
+         ENDIF
+100   CONTINUE
+*** Now there are no big overlaps, assign every particle in
+*** more than 1 jet to the closet jet.
+*** Any particles now in more than 1 jet are assigned to the CLOSET
+*** jet (in angle).
+      DO 140 I=1,NTRAK
+         NJ=0
+         DO 150 J=1, NJET
+         IF(JETLIS(J,I)) THEN
+            NJ=NJ+1
+            IJET(NJ)=J
+         ENDIF
+150      CONTINUE
+         IF (NJ .GT. 1) THEN
+*** Particle in > 1 jet - calc angles...
+            VEC1(1)=PP(1,I)
+            VEC1(2)=PP(2,I)
+            VEC1(3)=PP(3,I)
+            THMIN=0.
+            DO 160 J=1,NJ
+               VEC2(1)=PJ(1,IJET(J))
+               VEC2(2)=PJ(2,IJET(J))
+               VEC2(3)=PJ(3,IJET(J))
+               IF (MODE.NE.2) THEN
+                  CALL PXANG3(VEC1,VEC2,COST,THET,ITERR)
+               ELSE
+                  THET=(VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2
+               ENDIF
+               IF (J .EQ. 1) THEN
+                  THMIN=THET
+                  IJMIN=IJET(J)
+               ELSEIF (THET .LT. THMIN) THEN
+                  THMIN=THET
+                  IJMIN=IJET(J)
+               ENDIF
+160         CONTINUE
+*** Assign track to IJMIN
+            DO 170 J=1,NJET
+               JETLIS(J,I) = .FALSE.
+170         CONTINUE
+            JETLIS(IJMIN,I)=.TRUE.
+         ENDIF
+140   CONTINUE
+*** Recompute PJ
+      DO 200 I = 1,NJET
+        DO 210 MU = 1,4
+          PJ(MU,I) = 0.0
+210     CONTINUE
+        DO 220 N = 1,NTRAK
+          IF( JETLIS(I,N) ) THEN
+             IF (MODE.NE.2) THEN
+                DO 230 MU = 1,4
+                   PJ(MU,I) = PJ(MU,I) + PP(MU,N)
+230             CONTINUE
+             ELSE
+                PJ(1,I)=PJ(1,I)
+     +               + PP(4,N)/(PP(4,N)+PJ(4,I))*(PP(1,N)-PJ(1,I))
+                PJ(2,I)=PJ(2,I)
+     +               + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I))
+                PJ(4,I)=PJ(4,I)+PP(4,N)
+             ENDIF
+          ENDIF
+220     CONTINUE
+200   CONTINUE
+      RETURN
+      END
+*CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
+*CMZ :  1.06/00 14/03/94  15.37.45  by  P. Schleper
+*-- Author :
+*
+C+DECK,PXORD.
+       SUBROUTINE PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
+*
+*** Routine to put jets into order and eliminate tose less than EPSLON
+C+SEQ,DECLARE.
+      INTEGER MXTRAK,MXPROT
+      PARAMETER (MXTRAK=5000,MXPROT=5000)
+       INTEGER I, J, INDEX(MXPROT)
+       DOUBLE PRECISION PTEMP(4,MXPROT), ELIST(MXPROT)
+       INTEGER NJET,NTRAK
+       LOGICAL JETLIS(MXPROT,MXTRAK)
+       LOGICAL LOGTMP(MXPROT,MXTRAK)
+       DOUBLE PRECISION EPSLON,PJ(4,MXPROT)
+*** Puts jets in order of energy: 1 = highest energy etc.
+*** Then Eliminate jets with energy below EPSLON
+*
+*** Copy input arrays.
+      DO 100 I=1,NJET
+         DO 110 J=1,4
+            PTEMP(J,I)=PJ(J,I)
+110      CONTINUE
+         DO 120 J=1,NTRAK
+            LOGTMP(I,J)=JETLIS(I,J)
+120      CONTINUE
+100   CONTINUE
+      DO 150 I=1,NJET
+         ELIST(I)=PJ(4,I)
+150   CONTINUE
+*** Sort the energies...
+      CALL PXSORV(NJET,ELIST,INDEX,'I')
+*** Fill PJ and JETLIS according to sort ( sort is in ascending order!!)
+      DO 200 I=1, NJET
+         DO 210 J=1,4
+            PJ(J,I)=PTEMP(J,INDEX(NJET+1-I))
+210      CONTINUE
+         DO 220 J=1,NTRAK
+            JETLIS(I,J)=LOGTMP(INDEX(NJET+1-I),J)
+220      CONTINUE
+200   CONTINUE
+** Jets are now in order
+*** Now eliminate jets with less than Epsilon energy
+      DO 300, I=1, NJET
+         IF (PJ(4,I) .LT. EPSLON) THEN
+            NJET=NJET-1
+            PJ(4,I)=0.
+         ENDIF
+300   CONTINUE
+      RETURN
+      END
+
+********************************************************************
+*CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
+*CMZ :  1.06/00 14/03/94  15.37.44  by  P. Schleper
+*-- Author :
+C+DECK,PXSEAR.
+      SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
+     +                JETLIS,PJ,UNSTBL,IERR)
+*
+C+SEQ,DECLARE.
+      INTEGER MXTRAK, MXPROT
+      PARAMETER (MXTRAK=5000,MXPROT=5000)
+      INTEGER NTRAK, IERR, MODE
+      DOUBLE PRECISION COSR,PU(3,MXTRAK),PP(4,MXTRAK),VSEED(3)
+      LOGICAL UNSTBL
+      LOGICAL JETLIS(MXPROT,MXTRAK)
+      INTEGER NJET
+      DOUBLE PRECISION  PJ(4,MXPROT)
+*** Using VSEED as a trial axis , look for a stable jet.
+*** Check stable jets against those already found and add to PJ.
+*** Will try up to MXITER iterations to get a stable set of particles
+*** in the cone.
+      INTEGER MU,N,ITER
+      LOGICAL PXSAME,PXNEW,OK
+      LOGICAL NEWLIS(MXTRAK),OLDLIS(MXTRAK)
+      DOUBLE PRECISION OAXIS(3),NAXIS(3),PNEW(4)
+      INTEGER MXITER
+      PARAMETER(MXITER = 30)
+*
+      DO 100 MU=1,3
+        OAXIS(MU) = VSEED(MU)
+100   CONTINUE
+      DO 110 N = 1,NTRAK
+        OLDLIS(N) = .FALSE.
+110   CONTINUE
+      DO 120 ITER = 1,MXITER
+        CALL PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,PNEW,NEWLIS,OK)
+*** Return immediately if there were no particles in the cone.
+       IF (.NOT.OK) THEN
+         RETURN
+       ENDIF
+       IF(PXSAME(NEWLIS,OLDLIS,NTRAK)) THEN
+*** We have a stable jet.
+             IF (PXNEW(NEWLIS,JETLIS,NTRAK,NJET)) THEN
+*** And the jet is a new one. So add it to our arrays.
+*** Check arrays are big anough...
+             IF (NJET .EQ. MXPROT) THEN
+             WRITE (6,*) ' PXCONE:  Found more than MXPROT proto-jets'
+                IERR = -1
+                RETURN
+             ENDIF
+               NJET = NJET + 1
+               DO 130 N = 1,NTRAK
+                 JETLIS(NJET,N) = NEWLIS(N)
+130            CONTINUE
+               DO 140 MU=1,4
+                 PJ(MU,NJET)=PNEW(MU)
+140          CONTINUE
+             ENDIF
+             RETURN
+       ENDIF
+*** The jet was not stable, so we iterate again
+       DO 150 N=1,NTRAK
+         OLDLIS(N)=NEWLIS(N)
+150    CONTINUE
+       DO 160 MU=1,3
+         OAXIS(MU)=NAXIS(MU)
+160    CONTINUE
+120   CONTINUE
+      UNSTBL = .TRUE.
+      RETURN
+      END
+*CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+      SUBROUTINE PXSORV(N,A,K,OPT)
+C     Sort A(N) into ascending order
+C     OPT = 'I' : return index array K only
+C     OTHERWISE : return sorted A and index array K
+C-----------------------------------------------------------------------
+      INTEGER NMAX
+      PARAMETER (NMAX=5000)
+*
+*      INTEGER N,I,J,K(N),IL(NMAX),IR(NMAX)
+*LUND
+      INTEGER N,I,J,K(NMAX),IL(NMAX),IR(NMAX)
+      CHARACTER OPT
+*
+*      DOUBLE PRECISION A(N),B(NMAX)
+      DOUBLE PRECISION A(NMAX),B(NMAX)
+*LUND
+      IF (N.GT.NMAX) STOP 'Sorry, not enough room in Mike''s PXSORV'
+      IL(1)=0
+      IR(1)=0
+      DO 10 I=2,N
+      IL(I)=0
+      IR(I)=0
+      J=1
+   2  IF(A(I).GT.A(J)) GO TO 5
+   3  IF(IL(J).EQ.0) GO TO 4
+      J=IL(J)
+      GO TO 2
+   4  IR(I)=-J
+      IL(J)=I
+      GO TO 10
+   5  IF(IR(J).LE.0) GO TO 6
+      J=IR(J)
+      GO TO 2
+   6  IR(I)=IR(J)
+      IR(J)=I
+  10  CONTINUE
+      I=1
+      J=1
+      GO TO 8
+  20  J=IL(J)
+   8  IF(IL(J).GT.0) GO TO 20
+   9  K(I)=J
+      B(I)=A(J)
+      I=I+1
+      IF(IR(J)) 12,30,13
+  13  J=IR(J)
+      GO TO 8
+  12  J=-IR(J)
+      GO TO 9
+  30  IF(OPT.EQ.'I') RETURN
+      DO 31 I=1,N
+  31  A(I)=B(I)
+ 999  END
+
+*********************************************************************
+*CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
+*CMZ :  1.06/00 14/03/94  15.37.44  by  P. Schleper
+*-- Author :
+*
+C+DECK,PXTRY.
+       SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,
+     +                  PNEW,NEWLIS,OK)
+*
+C+SEQ,DECLARE.
+      INTEGER MXTRAK
+      PARAMETER (MXTRAK=5000)
+       INTEGER NTRAK,MODE
+*** Note that although PU and PP are assumed to be 2d arrays, they
+*** are used as 1d in this routine for efficiency
+       DOUBLE PRECISION COSR,PU(3*MXTRAK),PP(4*MXTRAK),OAXIS(3),PXMDPI
+       LOGICAL OK
+       LOGICAL NEWLIS(MXTRAK)
+       DOUBLE PRECISION NAXIS(3),PNEW(4)
+*** Finds all particles in cone of size COSR about OAXIS direction.
+*** Calculates 4-momentum sum of all particles in cone (PNEW) , and
+*** returns this as new jet axis NAXIS (Both unit Vectors)
+       INTEGER N,MU,NPU,NPP
+       DOUBLE PRECISION COSVAL,NORMSQ,NORM
+*
+       OK = .FALSE.
+       DO 100 MU=1,4
+          PNEW(MU)=0.0
+100    CONTINUE
+       NPU=-3
+       NPP=-4
+       DO 110 N=1,NTRAK
+          NPU=NPU+3
+          NPP=NPP+4
+          IF (MODE.NE.2) THEN
+             COSVAL=0.0
+             DO 120 MU=1,3
+                COSVAL=COSVAL+OAXIS(MU)*PU(MU+NPU)
+120          CONTINUE
+          ELSE
+             IF (ABS(PU(1+NPU)).GE.20.OR.ABS(OAXIS(1)).GE.20) THEN
+                COSVAL=-1000
+             ELSE
+                COSVAL=1-
+     +           ((OAXIS(1)-PU(1+NPU))**2+PXMDPI(OAXIS(2)-PU(2+NPU))**2)
+             ENDIF
+          ENDIF
+          IF (COSVAL.GE.COSR)THEN
+             NEWLIS(N) = .TRUE.
+             OK = .TRUE.
+             IF (MODE.NE.2) THEN
+                DO 130 MU=1,4
+                   PNEW(MU) = PNEW(MU) + PP(MU+NPP)
+130             CONTINUE
+             ELSE
+                PNEW(1)=PNEW(1)
+     +              + PP(4+NPP)/(PP(4+NPP)+PNEW(4))*(PP(1+NPP)-PNEW(1))
+                PNEW(2)=PNEW(2)
+     +              + PP(4+NPP)/(PP(4+NPP)+PNEW(4))
+     +               *PXMDPI(PP(2+NPP)-PNEW(2))
+                PNEW(4)=PNEW(4)+PP(4+NPP)
+             ENDIF
+          ELSE
+             NEWLIS(N)=.FALSE.
+          ENDIF
+110   CONTINUE
+*** If there are particles in the cone, calc new jet axis
+       IF (OK) THEN
+          IF (MODE.NE.2) THEN
+             NORMSQ = 0.0
+             DO 140 MU = 1,3
+                NORMSQ = NORMSQ + PNEW(MU)**2
+140          CONTINUE
+             NORM = SQRT(NORMSQ)
+          ELSE
+             NORM = 1
+          ENDIF
+          DO 150 MU=1,3
+             NAXIS(MU) = PNEW(MU)/NORM
+150       CONTINUE
+       ENDIF
+       RETURN
+       END
+
+*********************************************************************
+*CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
+*CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
+*-- Author :
+C+DECK,PXUVEC.
+*
+       SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR)
+*
+*** Routine to calculate unit vectors PU of all particles PP
+C+SEQ,DECLARE.
+      INTEGER MXTRAK
+      PARAMETER (MXTRAK=5000)
+      INTEGER NTRAK, IERR
+      DOUBLE PRECISION PP(4,MXTRAK)
+      DOUBLE PRECISION PU(3,MXTRAK)
+      INTEGER N,MU
+      DOUBLE PRECISION MAG
+       DO 100 N=1,NTRAK
+          MAG=0.0
+          DO 110 MU=1,3
+             MAG=MAG+PP(MU,N)**2
+110       CONTINUE
+          MAG=SQRT(MAG)
+          IF (MAG.EQ.0.0) THEN
+             WRITE(6,*)' PXCONE: An input particle has zero mod(p)'
+             IERR=-1
+             RETURN
+          ENDIF
+          DO 120 MU=1,3
+           PU(MU,N)=PP(MU,N)/MAG
+120       CONTINUE
+100    CONTINUE
+       RETURN
+       END
+*CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+      SUBROUTINE PXZERI(N,A)
+      INTEGER I,N,A(N)
+      DO 10 I=1,N
+        A(I)=0
+ 10   CONTINUE
+      END
+*CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+C     This is a set of routines written by Mike Seymour to provide the
+C     services presumably normally provided by standard OPAL routines
+C     PXZERV zeroes a vector
+C     PXZERI zeroes a vector of integers
+C     PXNORV normalizes a vector
+C     PXADDV adds two vectors
+C     PXSORV sorts a vector (copied from HERWIG)
+C     PXANG3 finds the angle (and its cosine) between two vectors
+C     PXMDPI moves its argument onto the range [-pi,pi)
+C-----------------------------------------------------------------------
+      SUBROUTINE PXZERV(N,A)
+      INTEGER I,N
+      DOUBLE PRECISION A(N)
+      DO 10 I=1,N
+        A(I)=0
+ 10   CONTINUE
+      END
+*-- Author :
+C-----------------------------------------------------------------------
+      SUBROUTINE PXADDV(N,A,B,C,ITERR)
+      INTEGER I,N,ITERR
+      DOUBLE PRECISION A(N),B(N),C(N)
+      DO 10 I=1,N
+        C(I)=A(I)+B(I)
+ 10   CONTINUE
+      END
+*CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+      SUBROUTINE PXANG3(A,B,COST,THET,ITERR)
+      INTEGER ITERR
+      DOUBLE PRECISION A(3),B(3),C,COST,THET
+      C=(A(1)**2+A(2)**2+A(3)**2)*(B(1)**2+B(2)**2+B(3)**2)
+      IF (C.LE.0) RETURN
+      C=1/SQRT(C)
+      COST=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*C
+      THET=ACOS(COST)
+      END
+*CMZ :  1.06/00 14/03/94  15.41.57  by  P. Schleper
+*-- Author :    P. Schleper   28/02/94
+       LOGICAL FUNCTION PXNEW(TSTLIS,JETLIS,NTRAK,NJET)
+*
+      INTEGER MXTRAK,MXPROT
+      PARAMETER (MXTRAK=5000,MXPROT=5000)
+       INTEGER NTRAK,NJET
+*** Note that although JETLIS is assumed to be a 2d array, it
+*** it is used as 1d in this routine for efficiency
+       LOGICAL TSTLIS(MXTRAK),JETLIS(MXPROT*MXTRAK)
+*** Checks to see if TSTLIS entries correspond to a jet already found
+*** and entered in JETLIS
+       INTEGER N, I, IN
+       LOGICAL MATCH
+*
+       PXNEW = .TRUE.
+       DO 100 I = 1,NJET
+          MATCH = .TRUE.
+          IN=I-MXPROT
+          DO 110 N = 1,NTRAK
+            IN=IN+MXPROT
+            IF(TSTLIS(N).NEQV.JETLIS(IN)) THEN
+             MATCH = .FALSE.
+             GO TO 100
+            ENDIF
+110       CONTINUE
+          IF (MATCH) THEN
+           PXNEW = .FALSE.
+           RETURN
+          ENDIF
+100    CONTINUE
+       RETURN
+       END
+*CMZ :  1.06/00 14/03/94  15.41.57  by  P. Schleper
+*-- Author :    P. Schleper   28/02/94
+       LOGICAL FUNCTION PXSAME(LIST1,LIST2,N)
+*
+       LOGICAL LIST1(*),LIST2(*)
+       INTEGER N
+*** Returns T if the first N elements of LIST1 are the same as the
+*** first N elements of LIST2.
+       INTEGER I
+*
+       PXSAME = .TRUE.
+       DO 100 I = 1,N
+        IF ( LIST1(I).NEQV.LIST2(I) ) THEN
+          PXSAME = .FALSE.
+          RETURN
+        ENDIF
+100    CONTINUE
+       RETURN
+       END
+*CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+      FUNCTION PXMDPI(PHI)
+      IMPLICIT NONE
+C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI)
+      DOUBLE PRECISION PXMDPI,PHI,PI,TWOPI,THRPI,EPS
+      PARAMETER (PI=3.141592654,TWOPI=6.283185307,THRPI=9.424777961)
+      PARAMETER (EPS=1E-15)
+      PXMDPI=PHI
+      IF (PXMDPI.LE.PI) THEN
+        IF (PXMDPI.GT.-PI) THEN
+          GOTO 100
+        ELSEIF (PXMDPI.GT.-THRPI) THEN
+          PXMDPI=PXMDPI+TWOPI
+        ELSE
+          PXMDPI=-MOD(PI-PXMDPI,TWOPI)+PI
+        ENDIF
+      ELSEIF (PXMDPI.LE.THRPI) THEN
+        PXMDPI=PXMDPI-TWOPI
+      ELSE
+        PXMDPI=MOD(PI+PXMDPI,TWOPI)-PI
+      ENDIF
+ 100  IF (ABS(PXMDPI).LT.EPS) PXMDPI=0
+      END
diff --git a/JETAN/read_jets.C b/JETAN/read_jets.C
new file mode 100644 (file)
index 0000000..f7df5e7
--- /dev/null
@@ -0,0 +1,170 @@
+
+void read_jets(const char* fn = "jets.root")
+
+{
+    //
+    // Some histos
+    //
+    TH1F* eH     = new TH1F("eH"  , "Jet Energy", 40.,  0., 200.);
+    TH1F* e1H    = new TH1F("e1H" , "Jet Energy", 40.,  0., 200.);
+    TH1F* e2H    = new TH1F("e2H" , "Jet Energy", 40.,  0., 200.);
+    TH1F* e3H    = new TH1F("e3H" , "Jet Energy", 40.,  0., 200.);
+    TH1F* e4H    = new TH1F("e4H" , "Jet Energy", 40.,  0., 200.);
+    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.);
+    
+
+  // load jet library
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libJETAN");
+
+  // open file
+  TFile *jFile = new TFile(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("AliJetESDReaderHeader"));
+  Int_t first = jReaderH->GetFirstEvent(); 
+  Int_t last  = jReaderH->GetLastEvent();
+  cout << "First event = " << first << "   Last event = " << last << endl;
+
+
+  // loop over events
+  AliJet *jets, *gjets;
+  AliLeading *leading;
+
+  for (Int_t i=first; i< last; i++) {
+      cout << "  Analyzing event " << i << endl;
+      // get 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 
+//
+      Int_t njets = jets->GetNJets();
+      
+      Float_t emax = 0.;
+      Int_t   imax = -1;
+      for (Int_t j = 0; j < njets; j++) {
+         if (jets->GetPt(j) > emax && TMath::Abs(jets->GetEta(j)) < 0.5) {
+             emax = jets->GetPt(j);
+             imax = j;
+         }
+      }
+      
+      if (imax == -1) {
+         e2H->Fill(gjets->GetPt(0));
+      } else {
+         eH->Fill(jets->GetPt(imax));
+         dr1H->Fill(jets->GetEta(imax) - gjets->GetEta(0));
+//
+//    Find the generated jet closest to the reconstructed 
+//
+         
+         Float_t rmin;
+         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++) {
+                 Float_t 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;
+                     igen = j;
+                 }
+             }
+
+             Float_t egen = gjets->GetPt(igen);
+             e1H->Fill(gjets->GetPt(igen));
+             
+             if (egen > 105. && egen < 125.) {
+                 e4H->Fill(emax);
+                 if (rmin > 0.3) etaH->Fill(etaj);
+                 dr2H->Fill(rmin);
+             }
+         }
+      }
+      
+
+//
+//   Leading particle
+//
+    Float_t etal = leading->GetLeading()->Eta();
+    Float_t phil = leading->GetLeading()->Phi();
+    Float_t el   = leading->GetLeading()->E();
+    e3H->Fill(el);
+    
+    Float_t rmin = 1.e6;
+    Int_t igen = -1;
+    Int_t ngen = gjets->GetNJets();
+    for (Int_t j = 0; j < ngen; j++) {
+       Float_t 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;
+       }
+    }
+       dr3H->Fill(rmin);
+
+//    cout << " Generated Jets:" << endl;
+//    gjets->PrintJets();
+//    cout << " Leading particle: " << endl;
+//    leading->PrintLeading();
+  }
+
+  // Get Control Plots
+  TCanvas* c1 = new TCanvas("c1");
+  eH->Draw();
+  e1H->SetLineColor(2);
+  e2H->SetLineColor(4);
+  e3H->SetLineColor(5);
+  e1H->Draw("same");
+  e2H->Draw("same");
+  e3H->Draw("same");
+
+  TCanvas* c2 = new TCanvas("c2");
+//  dr1H->Draw();
+  dr2H->SetLineColor(2);
+  dr2H->Draw("");
+  
+  TCanvas* c3 = new TCanvas("c3");
+  dr3H->Draw();
+  dr2H->Draw("same");
+
+  TCanvas* c4 = new TCanvas("c4");
+  eH->Draw();
+
+  TCanvas* c5 = new TCanvas("c5");
+  etaH->Draw();
+
+  TCanvas* c6 = new TCanvas("c6");
+  e4H->Draw();
+}
diff --git a/JETAN/read_jets_kine.C b/JETAN/read_jets_kine.C
new file mode 100644 (file)
index 0000000..45926a2
--- /dev/null
@@ -0,0 +1,195 @@
+
+void read_jets_kine(const char* fn = "jets.root")
+
+{
+    //
+    // Some histos
+    //
+    TH1F* eH     = new TH1F("eH"  , "Jet Energy", 40.,  0., 200.);
+    TH1F* e1H    = new TH1F("e1H" , "Jet Energy", 40.,  0., 200.);
+    TH1F* e2H    = new TH1F("e2H" , "Jet Energy", 40.,  0., 200.);
+    TH1F* e3H    = new TH1F("e3H" , "Jet Energy", 40.,  0., 200.);
+    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.);
+    
+
+  // load jet library
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libJETAN");
+
+  // open file
+  TFile *jFile = new TFile(fn);
+
+  // get jet header and display parameters
+  AliUA1JetHeader* jHeader = 
+    (AliUA1JetHeader*) (jFile->Get("AliUA1JetHeader"));
+  jHeader->PrintParameters();
+
+  // get reader header and events to be looped over
+  AliJetKineReaderHeader *jReaderH = 
+    (AliJetKineReaderHeader*)(jFile->Get("AliJetKineReaderHeader"));
+  Int_t first = jReaderH->GetFirstEvent(); 
+  Int_t last  = jReaderH->GetLastEvent();
+  cout << "First event = " << first << "   Last event = " << last << endl;
+
+
+  // loop over events
+  AliJet *jets, *gjets;
+  AliLeading *leading;
+
+  for (Int_t i=first; i< last; i++) {
+      cout << "  Analyzing event " << i << endl;
+      // get 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 E_T closest to 100. 
+//
+      Int_t njets = jets->GetNJets();
+      
+      Float_t emax  = 0.;
+      Int_t   imax  = -1;
+      Float_t demin = 1.e5;
+      
+      for (Int_t j = 0; j < njets; j++) {
+         if (TMath::Abs(jets->GetPt(j) - 100.)  < demin 
+             && TMath::Abs(jets->GetEta(j)) < 0.5) {
+             emax  = jets->GetPt(j);
+             imax  = j;
+             demin = TMath::Abs(jets->GetPt(j) - 100.);
+         }
+      }
+      
+      if (emax > 105.) {
+         printf("Strange %d %f %f %f\n", i, emax, jets->GetEta(imax), jets->GetPhi(imax));
+         Int_t ngen = gjets->GetNJets();
+           for (Int_t j = 0; j < ngen; j++) {
+                 Float_t etag = gjets->GetEta(j);
+                 Float_t phig = gjets->GetPhi(j);
+                 Float_t eg   = gjets->GetPt(j);
+                 printf("Generated %d %f %f %f\n", j, eg, etag, phig);
+                 
+           }
+      }
+      
+      if (imax == -1) {
+         e2H->Fill(gjets->GetPt(0));
+      } else {
+         eH->Fill(jets->GetPt(imax));
+         dr1H->Fill(jets->GetEta(imax) - gjets->GetEta(0));
+//
+//    Find the generated jet closest to the reconstructed 
+//
+         
+         Float_t rmin;
+         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++) {
+                 Float_t 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;
+                     igen = j;
+                 }
+             }
+             dr2H->Fill(rmin);
+             e1H->Fill(gjets->GetPt(igen));
+             if (rmin > 0.3) {
+                 etaH->Fill(etaj);
+                 printf("rmin %f %f %f %f %f %f\n", rmin, gjets->GetEta(igen), gjets->GetPhi(igen),
+                        etaj, phij, emax);
+
+                 Int_t ngen = gjets->GetNJets();
+                 for (Int_t j = 0; j < ngen; j++) {
+                     Float_t etag = gjets->GetEta(j);
+                     Float_t phig = gjets->GetPhi(j);
+                     Float_t eg   = gjets->GetPt(j);
+                     printf("   Gen %d %f %f %f\n", j, eg, etag, phig);
+           }
+             }
+             
+         }
+      }
+      
+
+//
+//   Leading particle
+//
+    Float_t etal = leading->GetLeading()->Eta();
+    Float_t phil = leading->GetLeading()->Phi();
+    Float_t el   = leading->GetLeading()->E();
+    e3H->Fill(el);
+    
+    Float_t rmin = 1.e6;
+    Int_t igen = -1;
+    Int_t ngen = gjets->GetNJets();
+    for (Int_t j = 0; j < ngen; j++) {
+       Float_t 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;
+       }
+    }
+       dr3H->Fill(rmin);
+
+//    cout << " Generated Jets:" << endl;
+//    gjets->PrintJets();
+//    cout << " Leading particle: " << endl;
+//    leading->PrintLeading();
+  }
+
+  // Get Control Plots
+  TCanvas* c1 = new TCanvas("c1");
+  eH->SetXTitle("E (GeV)");
+  eH->Draw();
+  
+  e1H->SetLineColor(2);
+  e2H->SetLineColor(4);
+  e3H->SetLineColor(5);
+  e1H->Draw("same");
+  e2H->Draw("same");
+  e3H->Draw("same");
+
+  TCanvas* c2 = new TCanvas("c2");
+//  dr1H->Draw();
+  dr2H->SetXTitle("#DeltaR(rec, gen)");
+  dr2H->SetLineColor(2);
+  dr2H->Draw("");
+  
+  TCanvas* c3 = new TCanvas("c3");
+  dr3H->SetXTitle("#DeltaR(rec, gen)");
+  dr3H->Draw();
+  dr2H->Draw("same");
+
+  TCanvas* c4 = new TCanvas("c4");
+  eH->Draw();
+
+  TCanvas* c5 = new TCanvas("c5");
+  etaH->Draw();
+
+  
+}
diff --git a/JETAN/testpxcone.C b/JETAN/testpxcone.C
new file mode 100755 (executable)
index 0000000..7aecc56
--- /dev/null
@@ -0,0 +1,30 @@
+void testpxcone()
+{
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libJETAN");
+
+  // define reader header
+  AliJetESDReaderHeader *jrh = new AliJetESDReaderHeader();
+  jrh->SetComment("testing");
+  jrh->SetDirectory("/home/work/alice/data/hij_central");
+  jrh->SetPattern("run00");
+  jrh->SetFirstEvent(0);
+  jrh->SetLastEvent(1);
+  jrh->SetPtCut(2);
+
+  // define reader and set its header
+  AliJetESDReader *jr = new AliJetESDReader();
+  jr->SetReaderHeader(jrh);
+
+  // define jet header, take the defaults
+  AliPxconeJetHeader *jh=new AliPxconeJetHeader();
+  jh->SetComment("Testing jet code");
+
+  // define jet finder. Set its header and reader
+  AliPxconeJetFinder *jf = new AliPxconeJetFinder();
+  jf->SetJetHeader(jh);
+  jf->SetJetReader(jr);
+  jf->SetPlotMode(kTRUE);
+  jf->SetOutputFile("pxcone.root");
+  // do the job
+  jf->Run();
+}
diff --git a/JETAN/testua1.C b/JETAN/testua1.C
new file mode 100755 (executable)
index 0000000..9ceb57f
--- /dev/null
@@ -0,0 +1,54 @@
+void testua1()
+{
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libJETAN");
+
+  // define reader header
+  AliJetESDReaderHeader *jrh = new AliJetESDReaderHeader();
+  jrh->SetComment("testing");
+  //jrh->SetDirectory("rfio:///castor/cern.ch/user/p/phristov/cent1_nq/104-125GeV");
+  jrh->SetDirectory("/home/morsch/jets_104_125");
+  jrh->SetPattern("00165");
+  jrh->SetFirstEvent(0);
+  jrh->SetLastEvent(1000);
+  jrh->SetPtCut(2.);
+  jrh->SetReadSignalOnly(kFALSE);
+
+  AliJetKineReaderHeader *krh = new AliJetKineReaderHeader();
+  krh->SetComment("testing");
+  krh->SetDirectory("/home/morsch/AliRoot/newio/FASTSIM/ctest/100/uq");
+  krh->SetFirstEvent(0);
+  krh->SetLastEvent(1000);
+  krh->SetPtCut(2.);
+  krh->SetFastSimTPC(kTRUE);
+  
+
+  // define reader and set its header
+  AliJetESDReader *jr = new AliJetESDReader();
+  jr->SetReaderHeader(jrh);
+
+  // define reader and set its header
+  AliJetKineReader *kr = new AliJetKineReader();
+  kr->SetReaderHeader(krh);
+
+  // define jet header
+  AliUA1JetHeader *jh=new AliUA1JetHeader();
+  jh->SetComment("UA1 jet code with default parameters");
+  jh->SetMode(0);
+  jh->SetRadius(0.4);
+  jh->SetMinCellEt(0.);
+  jh->SetEtSeed(4.);
+  jh->SetNbinPhi(420.);
+  jh->SetNbinEta(120.);
+  jh->SetEtaMin(-0.9);
+  jh->SetEtaMax(+0.9);  
+  
+
+  // define jet finder. Set its header and reader
+  AliUA1JetFinder *jf = new AliUA1JetFinder();
+  jf->SetJetHeader(jh);
+  jf->SetJetReader(jr);
+  jf->SetPlotMode(kTRUE);
+  jf->SetOutputFile("jets.root");
+  // do the job
+  jf->Run();
+}
diff --git a/JETAN/ua1_jet_finder.F b/JETAN/ua1_jet_finder.F
new file mode 100755 (executable)
index 0000000..06eeaa0
--- /dev/null
@@ -0,0 +1,322 @@
+c
+c => Original copy from /star/emc/shester/cl/pams/emc/jet/jet_finer_ua1.F
+c
+c:>------------------------------------------------------------------
+C:ROUTINE:      subroutine jet_finder_ua1
+C:DESCRIPTION:  UA1 jet algorithm from LUND JETSET
+C:RETURN VALUE: ierror=1 on error
+c:>------------------------------------------------------------------
+      subroutine ua1_jet_finder(
+     +        ncell, ncell_tot, etc, etac, phic,
+     +        min_move, max_move, mode, prec_bg, ierror)
+      implicit none ! 5-oct-2001
+
+      real C_PI
+      real C_2PI 
+      real etaCellSize 
+      real phiCellSize 
+      real arg
+      INTEGER NMAX,JMAX
+      parameter(NMAX=60000,JMAX=100) ! 10-oct-2201
+      integer ncell, ierror, mode, ncell_tot
+      real etc(NMAX),etac(NMAX),phic(NMAX)
+      real cone_rad, et_seed, ej_min, et_min
+      real min_move, max_move, prec_bg
+      integer flag(NMAX),index(NMAX)
+      integer n_iter,i,j,k,l,nc
+      real et_sum,et_ave,et_sum_old
+      real et,eta,phi,eta0,phi0,etab,phib,ets,etas,phis
+      real deta,dphi,r_int
+!
+      real occupationAll, occupationInJet  ! 3-oct-2001 by PAI 
+      integer maxTowerInJet                ! 5-oct-2001 by PAI
+      save    maxTowerInJet                ! This is obligatory
+      integer idPerfomance
+      data    idPerfomance /119/
+      save    idPerfomance
+!  print parameter - 27-sep-2001 by PAI
+      integer kpri
+      data    kpri /0/
+      integer njet, ncellj
+      real    etj, etaj, phij
+*    Results
+      COMMON /UA1JETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2), 
+     +     NCELLJ(100)
+*    Cell Geometry
+      COMMON /UA1CELL/ etaCellSize, phiCellSize
+*    Parameters
+      COMMON /UA1PARA/ cone_rad, et_seed, ej_min, et_min
+
+      C_PI  = 3.1415926
+      C_2PI = 2.*C_PI
+
+      if(kpri.eq.0) then
+        kpri = 1
+! for add. correction of jet energy if ocupation in jet not 100% !!
+! may differ from real number because grid 
+        maxTowerInJet=
+     +       int((C_PI*cone_rad*cone_rad)/(etaCellSize*phiCellSize)+0.5)
+        print 1, ncell_tot
+     +       ,cone_rad, et_seed, ej_min, et_min
+     +       ,min_move, max_move, mode, prec_bg
+     +       ,maxTowerInJet
+ 1    format(/
+     +    '    == jet_finder_UA1 ==  '/
+     +    ' ncell_tot                   ', i5/
+     +    ' cone rad                    ', f5.2/
+     +    ' et_seed                     ', f5.2,' GeV/C'/
+     +    ' ej_min                      ', f5.2,' GeV/C'/
+     +    ' et_min(tower after bg sub.) ', f5.2,' GeV/C'/
+     +    ' min_cone_move               ', f5.3/
+     +    ' max_cone_move               ', f5.3/
+     +    ' Mode for BG subtraction     ', i5/
+     +    ' prec_bg                     ', f5.4/
+     +    ' -- PAI"s addition --'/
+     +    ' maxTowerInJet               ', i5/
+     +    ' -- '/
+     +    )
+        if(NMAX .lt. ncell_tot) then
+          print 2, NMAX, ncell_tot
+ 2    format('<E> you must increase memory -> NMAX ',i6
+     +      ,' ncell_tot ', i6)
+        endif 
+      endif
+      call hf1(idPerfomance, 1., 1.)
+      occupationAll = float(ncell) / float(ncell_tot)
+      call hf1(111, occupationAll, 1.)
+! print parameter - 27-sep-2001 by PAI
+      ierror =0
+      n_iter =0
+c*-sort cells by Et decending order, store the order in index
+      call sortzv(etc,index,ncell,-1,1,0)   
+c*-sum up total energy of all cells
+      if(mode .eq. 1) then
+         n_iter=1
+         et_sum=0.0
+         do i=1, ncell
+            et_sum=et_sum+etc(i)
+         enddo
+         et_sum_old=et_sum
+         et_ave=et_sum/float(ncell_tot)
+      else
+         et_ave=0.0
+      endif
+      print *,'Iter ', n_iter, ' et_ave ', et_ave, ' #cells ', ncell 
+c*-Watch out!!! For mode=1, it can jump back to here for next iteration!!!
+ 999  continue
+c*-kill cells (flag=2) with Et below ET_MIN after background subtraction
+cfca      call vzero(flag,ncell)
+      do i=1, ncell
+         flag(i)=0
+         if(etc(i)-et_ave .le. et_min) flag(i)=2
+      enddo
+      njet = 0
+c*-Initiator cell is the one with highest Et of not yet used ones
+      i=1
+      j=index(i)
+      if(i.eq.1. and. etc(j).lt.et_seed) then
+        if(mode.eq.0 .or. (mode.eq.1 .and. n_iter.eq.1)) then
+          call hf1(idPerfomance, 2., 1.)
+          print *,' no cell with Et higher than et_seed ->', etc(j)       
+          return
+        endif
+      endif 
+      do while(etc(j) .ge. et_seed)
+         if(flag(j) .eq. 0) then
+C         write(6,*) j, etc(j), et_seed, etac(j), phic(j), flag(j) 
+         
+            et =etc(j)-et_ave
+            eta=etac(j)
+            phi=phic(j)
+            eta0=eta
+            phi0=phi
+            etab=eta
+            phib=phi
+            ets =0.0
+            etas=0.0
+            phis=0.0
+c*-weighted eta and phi. 
+            do k = 1, ncell
+               l = index(k)
+               if(flag(l).eq.0) then
+                  deta = etac(l)-eta
+c*-Is the cell is in the cone?
+                  if(abs(deta).le.cone_rad)then
+                     dphi=phic(l)-phi
+                     do while(dphi .gt. C_PI)
+                        dphi=dphi-C_2PI
+                     enddo
+                     do while(dphi .le. -C_PI)
+                        dphi=dphi+C_2PI
+                     enddo
+                     if(abs(dphi).le.cone_rad) then
+                        r_int=sqrt(deta**2+dphi**2)
+                        if(r_int.le.cone_rad)then
+c*-calculate offset from initiate cell
+                           deta=etac(l)-eta0
+                           dphi=phic(l)-phi0
+                           do while(dphi .gt. C_PI)
+                              dphi=dphi-C_2PI
+                           enddo
+                           do while(dphi .lt. -C_PI)
+                              dphi=dphi+C_2PI
+                           enddo
+                           et=etc(l)-et_ave
+                           etas=etas+abs(et)*deta
+                           phis=phis+abs(et)*dphi
+                           ets =ets +et
+c*-New weighted eta and phi including this cell
+                           eta=eta0+etas/ets
+                           phi=phi0+phis/ets                          
+c*-If cone does not move much from previous cone, just go next step
+                           r_int=sqrt((eta-etab)**2+(phi-phib)**2)
+                           if(r_int .le. min_move) then
+                              goto 159
+                           endif
+c*-Cone should not move more than MAX_CONE_MOVE from initiator cell
+                           r_int=sqrt((etas/ets)**2+(phis/ets)**2)              
+                           if(r_int .ge. max_move) then
+                              eta=etab
+                              phi=phib
+                              goto 159
+                           endif
+c*-Store this loop information
+                           etab=eta
+                           phib=phi
+                        endif
+                     endif
+                  endif
+               endif
+            enddo
+ 159        continue 
+            
+c*-sum up unused cells within required distance of given eta/phi
+            nc=0
+            ets=0.0
+            etas=0.0
+            phis=0.0
+            do k=1,ncell
+               l=index(k)
+               if(flag(l) .eq. 0) then
+                  deta=etac(l)-eta
+                  if(abs(deta).le.cone_rad)then
+                     dphi=phic(l)-phi
+                     do while(dphi .gt. C_PI)
+                        dphi=dphi-C_2PI
+                     enddo
+                     do while(dphi .le. -C_PI)
+                        dphi=dphi+C_2PI
+                     enddo
+                     if(abs(dphi).le.cone_rad) then
+                        r_int=sqrt(deta**2+dphi**2)
+                        if(r_int.le.cone_rad)then
+                           flag(l)=-1
+                           et  =etc(l)-et_ave
+                           ets =ets +et
+                           etas=etas+et*deta
+                           phis=phis+et*dphi
+                           nc  = nc + 1
+                        endif
+                     endif
+                  endif
+               endif
+            enddo  ! do k=1,ncell
+!  5-oct-2001 by PAI - remove 20-feb-2002 by PAI
+! 20-feb-2002 - it is work if you apply cut on eT before jet finder !!!
+!            if(maxTowerInJet .gt. nc) then
+!              ets = ets - et_ave*(maxTowerInJet - nc) 
+!            endif
+! 5-oct-2001 by PAI
+            
+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
+            
+            if(arg .lt. ej_min) then
+               do k=1,ncell
+                  if(flag(k).le.0) flag(k)=0
+               enddo
+               if(njet.eq.0) call hf1(idPerfomance, 3., 1.)
+            else
+c*-eles, store flags and jet variables
+               do k=1,ncell
+                  if(flag(k).eq.-1) flag(k)=1
+               enddo
+               phi=phi+phis/ets
+               do while(phi .ge. C_2PI)
+                  phi=phi-C_2PI
+               enddo
+               do while(phi .lt. 0.0)
+                  phi=phi+C_2PI
+               enddo
+               njet=njet+1
+               etj(njet) =ets
+               etaj(njet,1)=eta0
+               phij(njet,1)=phi0
+               etaj(njet,2)=etas
+               phij(njet,2)=phi
+               ncellj(njet)=nc
+               call hf1(112, float(nc)/float(maxTowerInJet), 1.) ! 8-oct-2001
+            endif 
+         endif
+         i=i+1
+         j=index(i)        
+      enddo
+      
+c*-recalculate energy sum excluding used cells.
+      if(mode.eq.1)then
+         et_sum=0.0
+         nc=0       ! #cells in jets
+         do i=1,ncell
+            if(flag(i).ne.1) then  ! 1 if cell in jet
+               et_sum=et_sum+etc(i)
+            else
+               nc=nc+1
+            endif
+         enddo
+c*-if background level changes more than prec_bg, go next iteration!!!
+c*-after 10 iteration, stop working and finish
+         if( (et_sum .gt. 0.) 
+     +        .and. (abs(et_sum-et_sum_old)/et_sum.gt.prec_bg 
+     +        .and. n_iter.le.10)
+     +        .or. n_iter.eq.1) then ! minimum 2 iteration - 10-oct-2001 by pai
+            et_ave=et_sum/float(ncell_tot-nc)
+            n_iter=n_iter+1
+            et_sum_old = et_sum
+            print *,'End of iteration : et_ave ', et_ave, ' nc ', nc
+     + , ' Jet energy ', ets     
+            goto 999
+c*-Watch out!!! Here is a big jump!!! 
+         endif
+         occupationInJet = float(ncell) / float(ncell_tot)
+         call hf1(111, occupationInJet, 1.)
+         write(*,*) njet,' jet found in ',n_iter,
+     +     ' iteration(s) : EtSum, EtAve =',
+     +     et_sum,et_sum/float(ncell_tot-nc)
+     + , ' ncell_tot ', ncell_tot, ' #cell in jet ', nc
+     + , ' occupationAll ', occupationAll
+     + , ' occupationInJet ', occupationInJet
+      endif
+      
+      if(njet.gt.100)then
+         write(*,*)'UA1:Problem:More than 100 jets found!'
+         ierror = 1
+      elseif(njet.eq.0)then
+         write(*,*)'UA1:Done:No jet found!'
+      else
+         write(*,*)
+     +'UA1:Done:Found ',njet,' jet(s)'         
+      end if
+      return
+      end
+
+
+