First commit.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 16 Jan 2002 17:24:11 +0000 (17:24 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 16 Jan 2002 17:24:11 +0000 (17:24 +0000)
EMCAL/AliEMCALJet.cxx [new file with mode: 0644]
EMCAL/AliEMCALJet.h [new file with mode: 0644]
EMCAL/AliEMCALJetFinder.cxx [new file with mode: 0644]
EMCAL/AliEMCALJetFinder.h [new file with mode: 0644]
EMCAL/EMCALLinkDef.h
EMCAL/Makefile
EMCAL/jet_finder_ua1.F [new file with mode: 0644]

diff --git a/EMCAL/AliEMCALJet.cxx b/EMCAL/AliEMCALJet.cxx
new file mode 100644 (file)
index 0000000..dd817e8
--- /dev/null
@@ -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 (file)
index 0000000..f5d48b1
--- /dev/null
@@ -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 <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
diff --git a/EMCAL/AliEMCALJetFinder.cxx b/EMCAL/AliEMCALJetFinder.cxx
new file mode 100644 (file)
index 0000000..45da980
--- /dev/null
@@ -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 <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
+    ;
+}
+
+
diff --git a/EMCAL/AliEMCALJetFinder.h b/EMCAL/AliEMCALJetFinder.h
new file mode 100644 (file)
index 0000000..e8d9526
--- /dev/null
@@ -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 <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
+
index 34cd59b..302c753 100644 (file)
@@ -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
index 6b85823..b5adce2 100644 (file)
@@ -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 (file)
index 0000000..8f293cf
--- /dev/null
@@ -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('<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
+
+
+