ALICE PPR.
--- /dev/null
+/**************************************************************************
+ * 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;
+ }
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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);
+}
+
+
--- /dev/null
+#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
+
--- /dev/null
+/**************************************************************************
+ * 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);
+}
+
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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();
+ }
+}
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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
+ //
+}
+
+////////////////////////////////////////////////////////////////////////
+
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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);
+}
+
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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);
+}
+
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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();
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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();
+}
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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());
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
--- /dev/null
+#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
--- /dev/null
+#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
+
+
+
--- /dev/null
+#-*-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
--- /dev/null
+ 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
--- /dev/null
+
+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();
+}
--- /dev/null
+
+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();
+
+
+}
--- /dev/null
+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();
+}
--- /dev/null
+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();
+}
--- /dev/null
+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
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff