--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+//*-- Author: Andreas Morsch (CERN)
+
+#include "AliEMCALJet.h"
+#include "Ecommon.h"
+
+ClassImp(AliEMCALJet)
+
+//____________________________________________________________________________
+AliEMCALJet::AliEMCALJet()
+{
+// Default constructor
+}
+
+AliEMCALJet::AliEMCALJet(Float_t energy, Float_t phi, Float_t eta)
+{
+// Constructor
+ fEnergy = energy;
+ fPhi = phi;
+ fEta = eta;
+}
+
+//____________________________________________________________________________
+AliEMCALJet::~AliEMCALJet()
+{
+// Destructor
+}
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#ifndef ALIEMCALJET_H
+#define ALIEMCALJET_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+#include <TObject.h>
+//*-- Author: Andreas Morsch (CERN)
+
+
+class AliEMCALJet : public TObject {
+ public:
+ AliEMCALJet();
+ AliEMCALJet(Float_t energy, Float_t phi, Float_t eta);
+ virtual ~AliEMCALJet();
+ void SetEnergy(Float_t val) {fEnergy = val;}
+ void SetPhi(Float_t val) {fPhi = val;}
+ void SetEta(Float_t val) {fEta = val;}
+
+ Float_t Energy() {return fEnergy;}
+ Float_t Phi() {return fPhi;}
+ Float_t Eta() {return fEta;}
+protected:
+ Float_t fEnergy; // Jet Energy
+ Float_t fEta; // Jet Phi
+ Float_t fPhi; // Jet Eta
+ ClassDef(AliEMCALJet,1) // Jet for EMCAL
+
+} ;
+
+#endif // ALIEMCALJet_H
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+//*-- Author: Andreas Morsch (CERN)
+#include <TClonesArray.h>
+#include <TTree.h>
+#include <TFile.h>
+#include <TH2.h>
+#include <TAxis.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
+#include "AliEMCALJetFinder.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALHit.h"
+#include "Ecommon.h"
+#include "AliRun.h"
+#include "AliEMCAL.h"
+#include "AliHeader.h"
+
+ClassImp(AliEMCALJetFinder)
+
+//____________________________________________________________________________
+AliEMCALJetFinder::AliEMCALJetFinder()
+{
+// Default constructor
+ fJets = 0;
+ fNjets = 0;
+ fLego = 0;
+}
+
+AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
+ : TTask(name, title)
+{
+// Constructor
+ fJets = new TClonesArray("AliEMCALJet",10000);
+ fNjets = 0;
+ for (Int_t i=0; i<30000; i++)
+ {
+ fEtCell[i] = 0.;
+ fEtaCell[i] = 0.;
+ fPhiCell[i] = 0.;
+ }
+ fLego = 0;
+}
+
+
+
+
+//____________________________________________________________________________
+AliEMCALJetFinder::~AliEMCALJetFinder()
+{
+// Destructor
+//
+ if (fJets){
+ fJets->Delete();
+ delete fJets;
+ }
+}
+
+
+
+#ifndef WIN32
+# define jet_finder_ua1 jet_finder_ua1_
+# define hf1 hf1_
+# define type_of_call
+
+#else
+# define jet_finder_ua1 J
+# define hf1 HF1
+# define type_of_call _stdcall
+#endif
+
+extern "C" void type_of_call
+jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
+ Float_t etc[30000], Float_t etac[30000],
+ Float_t phic[30000],
+ 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 AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
+ Float_t etac[30000], Float_t phic[30000],
+ Float_t min_move, Float_t max_move, Int_t mode,
+ Float_t prec_bg, Int_t ierror)
+{
+// Wrapper for fortran coded jet finder
+// Full argument list
+ jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
+ min_move, max_move, mode, prec_bg, ierror);
+ // Write result to output
+ WriteJets();
+}
+
+void AliEMCALJetFinder::Find()
+{
+// Wrapper for fortran coded jet finder using member data for
+// argument list
+
+ Float_t min_move = 0;
+ Float_t max_move = 0;
+ Int_t mode = 0;
+ Float_t prec_bg = 0.;
+ Int_t ierror = 0;
+
+ jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
+ min_move, max_move, mode, prec_bg, ierror);
+ // Write result to output
+ WriteJets();
+}
+
+
+Int_t AliEMCALJetFinder::Njets()
+{
+// Get number of reconstructed jets
+ return EMCALJETS.njet;
+}
+
+Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
+{
+// Get reconstructed jet energy
+ return EMCALJETS.etj[i];
+}
+
+Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
+{
+// Get reconstructed jet phi from leading particle
+ return EMCALJETS.phij[0][i];
+}
+
+Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
+{
+// Get reconstructed jet phi from weighting
+ return EMCALJETS.phij[1][i];
+}
+
+Float_t AliEMCALJetFinder::JetEtaL(Int_t i)
+{
+// Get reconstructed jet eta from leading particles
+ return EMCALJETS.etaj[0][i];
+}
+
+
+Float_t AliEMCALJetFinder::JetEtaW(Int_t i)
+{
+// Get reconstructed jet eta from weighting
+ return EMCALJETS.etaj[1][i];
+}
+
+void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
+{
+// Set grid cell size
+ EMCALCELLGEO.etaCellSize = eta;
+ EMCALCELLGEO.phiCellSize = phi;
+}
+
+void AliEMCALJetFinder::SetConeRadius(Float_t par)
+{
+// Set jet cone radius
+ EMCALJETPARAM.coneRad = par;
+}
+
+void AliEMCALJetFinder::SetEtSeed(Float_t par)
+{
+// Set et cut for seeds
+ EMCALJETPARAM.etSeed = par;
+}
+
+void AliEMCALJetFinder::SetMinJetEt(Float_t par)
+{
+// Set minimum jet et
+ EMCALJETPARAM.ejMin = par;
+}
+
+void AliEMCALJetFinder::SetMinCellEt(Float_t par)
+{
+// Set et cut per cell
+ EMCALJETPARAM.etMin = par;
+}
+
+
+void AliEMCALJetFinder::Test()
+{
+//
+// Test the finder call
+//
+ const Int_t nmax = 30000;
+ Int_t ncell = 10;
+ Int_t ncell_tot = 100;
+
+ Float_t etc[nmax];
+ Float_t etac[nmax];
+ Float_t phic[nmax];
+ Float_t min_move = 0;
+ Float_t max_move = 0;
+ Int_t mode = 0;
+ Float_t prec_bg = 0;
+ Int_t ierror = 0;
+
+
+ Find(ncell, ncell_tot, etc, etac, phic,
+ min_move, max_move, mode, prec_bg, ierror);
+
+}
+
+//
+// I/O
+//
+
+void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
+{
+ //
+ // Add a jet
+ //
+ TClonesArray &lrawcl = *fJets;
+ new(lrawcl[fNjets++]) AliEMCALJet(jet);
+}
+
+void AliEMCALJetFinder::ResetJets()
+{
+ //
+ // Reset Jet List
+ //
+ fJets->Clear();
+ fNjets = 0;
+}
+
+void AliEMCALJetFinder::WriteJets()
+{
+//
+// Add all jets to the list
+//
+ const Int_t kBufferSize = 4000;
+ TTree *pK = gAlice->TreeK();
+ const char* file = (pK->GetCurrentFile())->GetName();
+// I/O
+ AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
+ printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
+ if (fJets && gAlice->TreeR()) {
+ pEMCAL->MakeBranchInTree(gAlice->TreeR(),
+ "Jets",
+ &fJets,
+ kBufferSize,
+ file);
+ }
+ Int_t njet = Njets();
+ for (Int_t nj=0; nj<njet; nj++)
+ {
+ AliEMCALJet* jet = new AliEMCALJet(JetEnergy(nj),
+ JetPhiW(nj),
+ JetEtaW(nj));
+ AddJet(*jet);
+ delete jet;
+ }
+ Int_t nev = gAlice->GetHeader()->GetEvent();
+ gAlice->TreeR()->Fill();
+ char hname[30];
+ sprintf(hname,"TreeR%d", nev);
+ gAlice->TreeR()->Write(hname);
+ gAlice->TreeR()->Reset();
+ ResetJets();
+}
+
+void AliEMCALJetFinder::BookLego()
+{
+//
+// Book histo for discretisation
+//
+//
+// Get geometry parameters from
+ AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
+ AliEMCALGeometry* geom =
+ AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
+ fNbinEta = geom->GetNZ();
+ fNbinPhi = geom->GetNPhi();
+ const Float_t phiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
+ const Float_t phiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
+ fDphi = (phiMax-phiMin)/fNbinEta;
+ fDeta = 1.4/fNbinEta;
+ fNtot = fNbinPhi*fNbinEta;
+
+//
+ fLego = new TH2F("legoH","eta-phi",
+ fNbinEta, -0.7, 0.7,
+ fNbinPhi, phiMin, phiMax);
+}
+
+void AliEMCALJetFinder::DumpLego()
+{
+//
+// Dump lego histo into array
+//
+ fNcell = 0;
+ for (Int_t i = 1; i < fNbinEta; i++) {
+ for (Int_t j = 1; j < fNbinPhi; j++) {
+ Float_t e = fLego->GetBinContent(i,j);
+ TAxis* Xaxis = fLego->GetXaxis();
+ TAxis* Yaxis = fLego->GetYaxis();
+ Float_t eta = Xaxis->GetBinCenter(i);
+ Float_t phi = Yaxis->GetBinCenter(j);
+ fEtCell[fNcell] = e;
+ fEtaCell[fNcell] = eta;
+ fPhiCell[fNcell] = phi;
+ fNcell++;
+ } // phi
+ } // eta
+ fNcell--;
+}
+
+void AliEMCALJetFinder::ResetMap()
+{
+//
+// Reset eta-phi array
+
+ for (Int_t i=0; i<30000; i++)
+ {
+ fEtCell[i] = 0.;
+ fEtaCell[i] = 0.;
+ fPhiCell[i] = 0.;
+ }
+}
+
+
+void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
+{
+//
+// Fill Cells with hit information
+//
+//
+ ResetMap();
+
+ if (!fLego) BookLego();
+// Reset
+ if (flag == 0) fLego->Reset();
+//
+// Access particle information
+ Int_t npart = (gAlice->GetHeader())->GetNprimary();
+ for (Int_t part=2; part<npart; part++) {
+ TParticle *MPart = gAlice->Particle(part);
+ Int_t mpart = MPart->GetPdgCode();
+ Int_t child1 = MPart->GetFirstDaughter();
+ Float_t pT = MPart->Pt();
+ Float_t phi = MPart->Phi();
+ Float_t theta = MPart->Theta();
+ Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
+// if (part == 6 || part == 7)
+// {
+// printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f",
+// part-5, pT, eta, phi);
+// }
+
+ if (pT == 0.) continue;
+// charged or neutral
+ if (ich == 0) {
+ TParticlePDG* pdgP = MPart->GetPDG();
+ if (pdgP->Charge() == 0) continue;
+ }
+// skip partons
+ if (TMath::Abs(mpart) <= 6 ||
+ mpart == 21 ||
+ mpart == 92) continue;
+// acceptance cut
+ if (TMath::Abs(eta) > 0.7) continue;
+ if (phi*180./TMath::Pi() > 120.) continue;
+// final state only
+ if (child1 >= 0 && child1 < npart) continue;
+// printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
+// part, mpart, child1, eta, phi, pT);
+ fLego->Fill(eta, phi, pT);
+ } // primary loop
+ DumpLego();
+}
+
+void AliEMCALJetFinder::FillFromHits(Int_t flag)
+{
+//
+// Fill Cells with track information
+//
+//
+ ResetMap();
+
+ if (!fLego) BookLego();
+ if (flag == 0) fLego->Reset();
+
+//
+// Access hit information
+ AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
+// AliEMCALGeometry* geom =
+// AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
+
+ TTree *treeH = gAlice->TreeH();
+ Int_t ntracks = (Int_t) treeH->GetEntries();
+//
+// Loop over tracks
+//
+ Int_t nbytes = 0;
+
+
+ for (Int_t track=0; track<ntracks;track++) {
+ gAlice->ResetHits();
+ nbytes += treeH->GetEvent(track);
+//
+// Loop over hits
+//
+ for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
+ mHit;
+ mHit=(AliEMCALHit*) pEMCAL->NextHit())
+ {
+ Float_t x = mHit->X(); // x-pos of hit
+ Float_t y = mHit->Y(); // y-pos
+ Float_t z = mHit->Z(); // z-pos
+ Float_t eloss = mHit->GetEnergy(); // deposited energy
+// Int_t index = mHit->GetId(); // cell index
+// Float_t eta, phi;
+// geom->EtaPhiFromIndex(index, eta, phi);
+ Float_t r = TMath::Sqrt(x*x+y*y);
+ Float_t theta = TMath::ATan2(r,z);
+ Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
+ Float_t phi = TMath::ATan2(y,x);
+ fLego->Fill(eta, phi, eloss);
+// if (eloss > 1.) printf("\nx,y,z %f %f %f %f %f",
+// r, z, eta, phi, eloss);
+// printf("\n Max %f", fLego->GetMaximum());
+ } // Hit Loop
+ } // Track Loop
+ DumpLego();
+}
+
+
+void hf1(Int_t& id, Float_t& x, Float_t& wgt)
+{
+// dummy for hbook calls
+ ;
+}
+
+
--- /dev/null
+#ifndef ALIEMCALJETFINDER_H
+#define ALIEMCALJETFINDER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+//*-- Author:
+//*-- Andreas Morsch (CERN)
+
+#include <TTask.h>
+#include "AliEMCALJet.h"
+
+class TClonesArray;
+class TH2F;
+
+class AliEMCALJetFinder : public TTask {
+ public:
+ AliEMCALJetFinder();
+ AliEMCALJetFinder(const char* name, const char *title);
+ virtual ~AliEMCALJetFinder();
+ virtual void Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
+ Float_t etac[30000], Float_t phic[30000],
+ Float_t min_move, Float_t max_move, Int_t mode,
+ Float_t prec_bg, Int_t ierror);
+ virtual void Find();
+ virtual void Test();
+ // Geometry
+ virtual void SetCellSize(Float_t eta, Float_t phi);
+ // Parameters
+ virtual void SetConeRadius(Float_t par);
+ virtual void SetEtSeed(Float_t par);
+ virtual void SetMinJetEt(Float_t par);
+ virtual void SetMinCellEt(Float_t par);
+ // Results
+ virtual Int_t Njets();
+ virtual Float_t JetEnergy(Int_t);
+ virtual Float_t JetPhiL(Int_t);
+ virtual Float_t JetPhiW(Int_t);
+ virtual Float_t JetEtaL(Int_t);
+ virtual Float_t JetEtaW(Int_t);
+ virtual TH2F* GetLego() {return fLego;}
+ // I/O
+ virtual void FillFromHits(Int_t flag = 0);
+ virtual void FillFromTracks(Int_t flag = 0, Int_t ich = 0);
+ virtual void AddJet(const AliEMCALJet& jet);
+ virtual void WriteJets();
+ virtual void ResetJets();
+ virtual TClonesArray* Jets() {return fJets;}
+ private:
+ virtual void BookLego();
+ virtual void DumpLego();
+ virtual void ResetMap();
+ protected:
+ TClonesArray* fJets; //! List of Jets
+ TH2F* fLego; //! Lego Histo
+ Int_t fNjets; //! Number of Jets
+ Float_t fDeta; //! eta cell size
+ Float_t fDphi; //! phi cell size
+ Int_t fNcell; //! number of cells
+ Int_t fNtot; //! total number of cells
+ Int_t fNbinEta; //! number of cells in eta
+ Int_t fNbinPhi; //! number of cells in phi
+ Float_t fEtCell[30000]; //! Cell Energy
+ Float_t fEtaCell[30000]; //! Cell eta
+ Float_t fPhiCell[30000]; //! Cell phi
+
+ ClassDef(AliEMCALJetFinder,1) // JetFinder for EMCAL
+} ;
+#endif // ALIEMCALJetFinder_H
+
#pragma link C++ class AliEMCALDigitizer+;
//#pragma link C++ class AliEMCALClusterizer+;
#pragma link C++ class AliEMCALGetter+;
+#pragma link C++ class AliEMCALJetFinder+;
+#pragma link C++ class AliEMCALJet+;
+#pragma link C++ class AliEMCALFast+;
#endif
SRCS = AliEMCAL.cxx AliEMCALGeometry.cxx AliEMCALv0.cxx \
AliEMCALv1.cxx AliEMCALHit.cxx AliEMCALDigit.cxx \
- AliEMCALSDigitizer.cxx AliEMCALDigitizer.cxx \
+ AliEMCALSDigitizer.cxx AliEMCALDigitizer.cxx \
+ AliEMCALGetter.cxx \
+ AliEMCALJetFinder.cxx AliEMCALJet.cxx AliEMCALFast.cxx \
AliEMCALGetter.cxx AliEMCALTick.cxx
# AliEMCALRecPoint.cxx AliEMCALEmcRecPoint.cxx AliEMCALClusterizer.cxx
+
+# FORTRAN sources
+
+FSRCS = jet_finder_ua1.F
+
# C++ Headers
HDRS = $(SRCS:.cxx=.h) $(ROOTSYS)/include/TTree.h EMCALLinkDef.h
OBJS = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(SRCS)) $(DICTO)
+# FORTRAN Objects
+FOBJS = $(patsubst %.F,tgt_$(ALICE_TARGET)/%.o,$(FSRCS))
# C++ compilation flags
CXXFLAGS = $(CXXOPTS) -I$(ROOTSYS)/include -I. -I$(ALICE_ROOT)/include/
+
+FFLAGS = $(FOPT)
+
ALSRCS = $(SRCS) $(SHSRCS) $(RCSRCS) $(DUSRCS) dummies.c
ALOBJS = $(SHOBJS) $(RCOBJS) $(DUOBJS)
$(DICT): $(HDRS)
-depend: $(SRCS)
+depend: $(SRCS)
-TOCLEAN = $(OBJS) *Cint.cxx *Cint.h
+TOCLEAN = $(OBJS) $(FOBJS) *Cint.cxx *Cint.h
CHECKS = $(patsubst %.cxx,check/%.viol,$(SRCS))
--- /dev/null
+c
+c => Original copy from /star/emc/shester/cl/pams/emc/jet/jet_finer_ua1.F
+c
+c:>------------------------------------------------------------------
+C:ROUTINE: subroutine jet_finder_ua1
+C:DESCRIPTION: UA1 jet algorithm from LUND JETSET called from EMC_erj
+C:ARGUMENTS: emc_energy_h, emc_energy,
+C:ARGUMENTS: et_tot, EMC_jetpar, ierror
+C:RETURN VALUE: ierror=1 on error
+c:>------------------------------------------------------------------
+ subroutine jet_finder_ua1(
+ + ncell, ncell_tot, etc, etac, phic,
+ + min_move, max_move, mode, prec_bg, ierror)
+ implicit none ! 5-oct-2001
+c#include "math_constants.inc"
+c#include "calorSize.inc"
+ real C_PI
+ real C_2PI
+ real etaCellSize
+ real phiCellSize
+ INTEGER NMAX,JMAX
+ parameter(NMAX=30000,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 /119/
+ save idPerfomance
+! print parameter - 27-sep-2001 by PAI
+ integer kpri/0/
+ integer njet, ncellj
+ real etj, etaj, phij
+* Results
+ COMMON /EMCALJETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2),
+ + NCELLJ(100)
+* Cell Geometry
+ COMMON /EMCALCELLGEO/ etaCellSize, phiCellSize
+* Parameters
+ COMMON /EMCALJETPARAM/ cone_rad, et_seed, ej_min, et_min
+
+ C_PI = 3.1415926
+ C_2PI = 2.*C_PI
+
+ if(kpri.eq.0) then
+ kpri = 1
+! for add. correction of jet energy if ocupation in jet not 100% !!
+! may differ from real number because grid
+ maxTowerInJet=
+ + int((C_PI*cone_rad*cone_rad)/(etaCellSize*phiCellSize)+0.5)
+ print 1, ncell_tot
+ + ,cone_rad, et_seed, ej_min, et_min
+ + ,min_move, max_move, mode, prec_bg
+ + ,maxTowerInJet
+ 1 format(/
+ + ' == jet_finder_UA1 == '/
+ + ' ncell_tot ', i5/
+ + ' cone rad ', f5.2/
+ + ' et_seed ', f5.2,' GeV/C'/
+ + ' ej_min ', f5.2,' GeV/C'/
+ + ' et_min(tower after bg sub.) ', f5.2,' GeV/C'/
+ + ' min_cone_move ', f5.3/
+ + ' max_cone_move ', f5.3/
+ + ' Mode for BG subtraction ', i5/
+ + ' prec_bg ', f5.4/
+ + ' -- PAI"s addition --'/
+ + ' maxTowerInJet ', i5/
+ + ' -- '/
+ + )
+ if(NMAX .lt. ncell_tot) then
+ print 2, NMAX, ncell_tot
+ 2 format('<E> you must increase memory -> NMAX ',i6
+ + ,' ncell_tot ', i6)
+ endif
+ endif
+ call hf1(idPerfomance, 1., 1.)
+ occupationAll = float(ncell) / float(ncell_tot)
+ call hf1(111, occupationAll, 1.)
+! print parameter - 27-sep-2001 by PAI
+ ierror =0
+ n_iter =0
+c*-sort cells by Et decending order, store the order in index
+ call sortzv(etc,index,ncell,-1,1,0)
+c*-sum up total energy of all cells
+ if(mode .eq. 1) then
+ n_iter=1
+ et_sum=0.0
+ do i=1, ncell
+ et_sum=et_sum+etc(i)
+ enddo
+ et_sum_old=et_sum
+ et_ave=et_sum/float(ncell_tot)
+ else
+ et_ave=0.0
+ endif
+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
+ call vzero(flag,ncell)
+ do i=1, ncell
+ 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)
+ endif
+ endif
+ do while(etc(j) .ge. et_seed)
+ if(flag(j) .eq. 0) then
+ 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
+ if(maxTowerInJet .gt. nc) then
+ ets = ets - et_ave*(maxTowerInJet - nc)
+ endif
+! 5-oct-2001 by PAI
+
+c*-reject cluster below minimum Ej_min
+ etas=eta+etas/ets
+ if(ets*cosh(etas/ets).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( (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
+ goto 999
+c*-Watch out!!! Here is a big jump!!!
+ endif
+ occupationInJet = float(ncell) / float(ncell_tot)
+ call hf1(111, occupationInJet, 1.)
+ write(*,*) njet,' jet found in ',n_iter,
+ + ' iteration(s) : EtSum, EtAve =',
+ + et_sum,et_sum/float(ncell_tot-nc)
+ + , ' ncell_tot ', ncell_tot, ' #cell in jet ', nc
+ + , ' occupationAll ', occupationAll
+ + , ' occupationInJet ', occupationInJet
+ endif
+
+ if(njet.gt.100)then
+ write(*,*)'UA1:Problem:More than 100 jets found!'
+ ierror = 1
+ elseif(njet.eq.0)then
+ write(*,*)'UA1:Done:No jet found!'
+ else
+ write(*,*)
+ +'UA1:Done:Found ',njet,' jet(s)'
+ end if
+ return
+ end
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff