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