From 471f69dc410a6ae90abbe2c04daa4330449724a9 Mon Sep 17 00:00:00 2001 From: morsch Date: Wed, 16 Jan 2002 17:24:11 +0000 Subject: [PATCH] First commit. --- EMCAL/AliEMCALJet.cxx | 54 +++++ EMCAL/AliEMCALJet.h | 32 +++ EMCAL/AliEMCALJetFinder.cxx | 453 ++++++++++++++++++++++++++++++++++++ EMCAL/AliEMCALJetFinder.h | 71 ++++++ EMCAL/EMCALLinkDef.h | 3 + EMCAL/Makefile | 18 +- EMCAL/jet_finder_ua1.F | 303 ++++++++++++++++++++++++ 7 files changed, 931 insertions(+), 3 deletions(-) create mode 100644 EMCAL/AliEMCALJet.cxx create mode 100644 EMCAL/AliEMCALJet.h create mode 100644 EMCAL/AliEMCALJetFinder.cxx create mode 100644 EMCAL/AliEMCALJetFinder.h create mode 100644 EMCAL/jet_finder_ua1.F diff --git a/EMCAL/AliEMCALJet.cxx b/EMCAL/AliEMCALJet.cxx new file mode 100644 index 00000000000..dd817e81274 --- /dev/null +++ b/EMCAL/AliEMCALJet.cxx @@ -0,0 +1,54 @@ +/************************************************************************** + * 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 +} + + + + + + + + + + + diff --git a/EMCAL/AliEMCALJet.h b/EMCAL/AliEMCALJet.h new file mode 100644 index 00000000000..f5d48b13105 --- /dev/null +++ b/EMCAL/AliEMCALJet.h @@ -0,0 +1,32 @@ +#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 +//*-- 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 diff --git a/EMCAL/AliEMCALJetFinder.cxx b/EMCAL/AliEMCALJetFinder.cxx new file mode 100644 index 00000000000..45da980d08a --- /dev/null +++ b/EMCAL/AliEMCALJetFinder.cxx @@ -0,0 +1,453 @@ +/************************************************************************** + * 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 +#include +#include +#include +#include +#include +#include +#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; njGetHeader()->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; partParticle(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; trackResetHits(); + 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 + ; +} + + diff --git a/EMCAL/AliEMCALJetFinder.h b/EMCAL/AliEMCALJetFinder.h new file mode 100644 index 00000000000..e8d95267399 --- /dev/null +++ b/EMCAL/AliEMCALJetFinder.h @@ -0,0 +1,71 @@ +#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 +#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 + diff --git a/EMCAL/EMCALLinkDef.h b/EMCAL/EMCALLinkDef.h index 34cd59b1e1b..302c75357e2 100644 --- a/EMCAL/EMCALLinkDef.h +++ b/EMCAL/EMCALLinkDef.h @@ -17,4 +17,7 @@ #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 diff --git a/EMCAL/Makefile b/EMCAL/Makefile index 6b858235c6d..b5adce28b46 100644 --- a/EMCAL/Makefile +++ b/EMCAL/Makefile @@ -11,10 +11,17 @@ PACKAGE = EMCAL 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 @@ -29,10 +36,15 @@ DICTO = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(DICT)) 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) @@ -47,9 +59,9 @@ $(LIBDIR)/libEMCAL.$(SL): $(OBJS) $(FOBJS) $(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)) diff --git a/EMCAL/jet_finder_ua1.F b/EMCAL/jet_finder_ua1.F new file mode 100644 index 00000000000..8f293cf4232 --- /dev/null +++ b/EMCAL/jet_finder_ua1.F @@ -0,0 +1,303 @@ +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(' 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 + + + -- 2.43.0