Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinder.cxx
old mode 100755 (executable)
new mode 100644 (file)
index 3dd4cef..723417d
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
+
+/* $Id */
  
 //---------------------------------------------------------------------
-// UA1 Jet finder 
-// manages the search for jets 
-// Author: jgcn@mda.cinvestav.mx
-// (code adapted from EMCAL directory)
+// UA1 Cone Algorithm Jet finder for jet studies
+// manages the search for jets using charged particle momentum and 
+// neutral cell energy information
+// Authors: Rafael.Diaz.Valdes@cern.ch 
+//          magali.estienne@subatech.in2p3.fr
+//          alexandre.shabetai@cern.ch
+// ** 2011
+// Modified accordingly to reader/finder splitting and new handling of neutral information
+// Versions V1 and V2 merged
 //---------------------------------------------------------------------
 
-#include <Riostream.h>
-#include <TLorentzVector.h>
-#include <TFile.h>
 #include <TH2F.h>
-#include "AliUA1JetFinder.h"
-#include "AliUA1JetHeader.h"
-#include "UA1Common.h"
-#include "AliJetReaderHeader.h"
-#include "AliJetReader.h"
-#include "AliJet.h"
+#include <TMath.h>
 
+#include "AliUA1JetFinder.h"
+#include "AliUA1JetHeaderV1.h"
+#include "AliJetCalTrk.h"
+#include "AliJetBkg.h"
+#include "AliAODJetEventBackground.h"
+#include "AliAODJet.h"
 
 ClassImp(AliUA1JetFinder)
 
 ////////////////////////////////////////////////////////////////////////
 
 AliUA1JetFinder::AliUA1JetFinder():
-  fHeader(0x0),
-  fLego(0x0)
+  AliJetFinder(),
+  fLego(0),  
+  fJetBkg(new AliJetBkg())
 {
   // Default constructor
 }
 
-////////////////////////////////////////////////////////////////////////
-
+//-----------------------------------------------------------------------
 AliUA1JetFinder::~AliUA1JetFinder()
-
 {
-  // destructor
+  // Destructor
+  delete fLego;
+  delete fJetBkg;
+
 }
 
-////////////////////////////////////////////////////////////////////////
+//-----------------------------------------------------------------------
+void AliUA1JetFinder::FindJets()
+{
+  //  Used to find jets using charged particle momentum information 
+  //  & neutral energy from calo cells
+  //
+  //  1) Fill cell map array
+  //  2) calculate total energy and fluctuation level
+  //  3) Run algorithm
+  //     3.1) look centroides in cell map
+  //     3.2) calculate total energy in cones
+  //     3.3) flag as a possible jet
+  //     3.4) reorder cones by energy
+  //  4) subtract backg in accepted jets
+  //  5) fill AliJet list
+  
+  // transform input to pt,eta,phi plus lego
+  
+  AliUA1JetHeader* header  = (AliUA1JetHeader*) fHeader;
+  Int_t              nIn = fCalTrkEvent->GetNCalTrkTracks();
+  fDebug   = fHeader->GetDebug();
+  
+  if (nIn <= 0) return;
+  fJetBkg->SetHeader(fHeader);
+  fJetBkg->SetCalTrkEvent(GetCalTrkEvent());
+  fJetBkg->SetDebug(fDebug);
+  // local arrays for input
+  // ToDo: check memory fragmentation, maybe better to
+  // define them globally and resize as needed
+  // Fragmentation should be worse for low mult...
+  Float_t* ptT      = new Float_t[nIn];
+  Float_t* etaT     = new Float_t[nIn];
+  Float_t* phiT     = new Float_t[nIn];
+  Int_t*   injet    = new Int_t[nIn];
+  Int_t*   injetOk  = new Int_t[nIn];
+
+  memset(ptT,0,sizeof(Float_t)*nIn);
+  memset(etaT,0,sizeof(Float_t)*nIn);
+  memset(phiT,0,sizeof(Float_t)*nIn);
+  memset(injet,0,sizeof(Int_t)*nIn);
+  memset(injetOk,-1,sizeof(Int_t)*nIn);
+
+  // load input vectors and calculate total energy in array
+
+  // total energy in array
+  Float_t etbgTotal = 0.;
+  Float_t npart = 0.;
+  Float_t etbg2 = 0.;
+
+  for (Int_t i = 0; i < fCalTrkEvent->GetNCalTrkTracks(); i++){
+    ptT[i]  = fCalTrkEvent->GetCalTrkTrack(i)->GetPt();
+    etaT[i] = fCalTrkEvent->GetCalTrkTrack(i)->GetEta();
+    phiT[i] = ((fCalTrkEvent->GetCalTrkTrack(i)->GetPhi() < 0) ? (fCalTrkEvent->GetCalTrkTrack(i)->GetPhi()) + 2 * TMath::Pi() : fCalTrkEvent->GetCalTrkTrack(i)->GetPhi());
+    //fCalTrkEvent->GetCalTrkTrack(i)->Print(Form("%d",i));
+    if (fCalTrkEvent->GetCalTrkTrack(i)->GetCutFlag() != 1) continue;
+    fLego->Fill(etaT[i], phiT[i], ptT[i]);
+    npart += 1;
+    etbgTotal+= ptT[i];
+    etbg2 += ptT[i]*ptT[i];
+  }
+  
+  // calculate total energy and fluctuation in map
+  Double_t meanpt = 0.;
+  Double_t ptRMS = 0.;
+  if(npart>0){
+    meanpt = etbgTotal/npart;
+    etbg2 = etbg2/npart;
+    if(etbg2>(meanpt*meanpt)){// prenent NAN, should only happen due to numerical instabilities
+      ptRMS = TMath::Sqrt(etbg2-meanpt*meanpt);
+    }
+  }
+  Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
+  
+  // arrays to hold jets
+  Float_t etaJet[kMaxJets];
+  Float_t phiJet[kMaxJets];
+  Float_t etJet[kMaxJets];
+  Float_t etsigJet[kMaxJets]; //signal et in jet
+  Float_t etallJet[kMaxJets];  // total et in jet (tmp variable)
+  Int_t   ncellsJet[kMaxJets];
+  Int_t   multJetT[kMaxJets];
+  Int_t   multJetC[kMaxJets];
+  Int_t   multJet[kMaxJets];
+  Float_t *areaJet = new Float_t[kMaxJets];
+  // Used for jet reordering at the end of the jet finding procedure
+  Float_t etaJetOk[kMaxJets];
+  Float_t phiJetOk[kMaxJets];
+  Float_t etJetOk[kMaxJets];
+  Float_t etsigJetOk[kMaxJets]; //signal et in jet
+  Float_t etallJetOk[kMaxJets];  // total et in jet (tmp variable)
+  Int_t   ncellsJetOk[kMaxJets];
+  Int_t   multJetOk[kMaxJets];
+  Float_t *areaJetOk = new Float_t[kMaxJets];
+  Int_t   nJets; // to hold number of jets found by algorithm
+  Int_t   nj;    // number of jets accepted
+  Float_t prec  = header->GetPrecBg();
+  Float_t bgprec = 1;
+
+  while(bgprec > prec){
+    //reset jet arrays in memory
+    memset(etaJet,0,sizeof(Float_t)*kMaxJets);
+    memset(phiJet,0,sizeof(Float_t)*kMaxJets);
+    memset(etJet,0,sizeof(Float_t)*kMaxJets);
+    memset(etallJet,0,sizeof(Float_t)*kMaxJets);
+    memset(etsigJet,0,sizeof(Float_t)*kMaxJets);
+    memset(ncellsJet,0,sizeof(Int_t)*kMaxJets);
+    memset(multJetT,0,sizeof(Int_t)*kMaxJets);
+    memset(multJetC,0,sizeof(Int_t)*kMaxJets);
+    memset(multJet,0,sizeof(Int_t)*kMaxJets);
+    memset(areaJet,0,sizeof(Float_t)*kMaxJets);
+    memset(etaJetOk,0,sizeof(Float_t)*kMaxJets);
+    memset(phiJetOk,0,sizeof(Float_t)*kMaxJets);
+    memset(etJetOk,0,sizeof(Float_t)*kMaxJets);
+    memset(etallJetOk,0,sizeof(Float_t)*kMaxJets);
+    memset(etsigJetOk,0,sizeof(Float_t)*kMaxJets);
+    memset(ncellsJetOk,0,sizeof(Int_t)*kMaxJets);
+    memset(multJetOk,0,sizeof(Int_t)*kMaxJets);
+    memset(areaJetOk,0,sizeof(Float_t)*kMaxJets);
+    nJets = 0;
+    nj = 0;
+    // reset particles-jet array in memory
+    memset(injet,-1,sizeof(Int_t)*nIn);
+    //run cone algorithm finder
+    RunAlgoritm(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
+    //run background subtraction
+    if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
+      nj = header->GetNAcceptJets();
+    else
+      nj = nJets;
+
+    //subtract background
+    Float_t etbgTotalN = 0.0; //new background
+    Float_t sigmaN = 0.0; //new background
+    if(header->GetBackgMode() == 1) {// standard
+      fJetBkg->SubtractBackg(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);}
+    if(header->GetBackgMode() == 2) //cone
+      fJetBkg->SubtractBackgCone(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
+    if(header->GetBackgMode() == 3) //ratio
+      fJetBkg->SubtractBackgRatio(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
+    if(header->GetBackgMode() == 4) //statistic
+      fJetBkg->SubtractBackgStat(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
+
+    //calc precision
+    if(etbgTotalN != 0.0)
+      bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
+    else
+      bgprec = 0;
+    etbgTotal = etbgTotalN; // update with new background estimation
+
+  } //end while
+
+  // add tracks to the jet if it wasn't yet done
+  if (header->GetBackgMode() == 0){
+    Float_t rc= header->GetRadius();
+    for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+      for(Int_t ijet=0; ijet<nJets; ijet++){
+        Float_t deta = etaT[jpart] - etaJet[ijet];
+        Float_t dphi = phiT[jpart] - phiJet[ijet];
+        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 <= rc){ // particles inside this cone
+         injet[jpart] = ijet;
+          break;
+        }
+      }// end jets loop
+    } //end particle loop
+  }
+
+  // add jets to list
+  if (fDebug>1) printf("Found %d jets \n", nj);
+  
+  // Reorder jets by et in cone
+  // Sort jets by energy
+  Int_t idx[kMaxJets];
+  TMath::Sort(nJets, etJet, idx);
+  for(Int_t p = 0; p < nJets; p++)
+    {
+      etaJetOk[p]    = etaJet[idx[p]];
+      phiJetOk[p]    = phiJet[idx[p]];
+      etJetOk[p]     = etJet[idx[p]];
+      etallJetOk[p]  = etJet[idx[p]];
+      etsigJetOk[p]  = etsigJet[idx[p]];
+      ncellsJetOk[p] = ncellsJet[idx[p]];
+      multJetOk[p]   = multJet[idx[p]];
+      areaJetOk[p]   = areaJet[idx[p]];
+    }
 
-#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
+  Int_t nTracks = fCalTrkEvent->GetNCalTrkTracks();
+  
+  for(Int_t kj=0; kj<nj; kj++)
+    {
+      if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
+         (etaJetOk[kj] < (header->GetJetEtaMin())) ||
+         (phiJetOk[kj] > (header->GetJetPhiMax())) ||
+         (phiJetOk[kj] < (header->GetJetPhiMin())) ||
+         (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
+      Float_t px=-999, py=-999 ,pz=-999 ,en=-999; // convert to 4-vector
+      px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
+      py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
+      pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
+      en = TMath::Sqrt(px * px + py * py + pz * pz);
+      AliAODJet jet(px, py, pz, en);
 
-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);
+      // Calc jet area if it wasn't yet done
+      if (header->GetBackgMode() == 0){
+        // calculate the area of the jet
+        Float_t rc= header->GetRadius();
+        areaJetOk[kj] = fJetBkg->CalcJetAreaEtaCut(rc,etaJetOk[kj]);
+      }
 
-extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
+      jet.SetEffArea(areaJetOk[kj],0.,0.,0.);
+      jet.SetBgEnergy(etbgTotal,0.);
+      if (fDebug>1) jet.Print(Form("%d",kj));
+      
+      for(Int_t jpart = 0; jpart < nTracks; jpart++) { // loop for all particles in array
+        // Track to jet reordering
+        if(injet[jpart] == idx[kj]){
+          injetOk[jpart] = kj;
+        }
+        // Check if the particle belongs to the jet and add the ref
+        if(injetOk[jpart] == kj && fCalTrkEvent->GetCalTrkTrack(jpart)->GetCutFlag() == 1) {
+          jet.AddTrack(fCalTrkEvent->GetCalTrkTrack(jpart)->GetTrackObject());
+       }
+      }
 
-////////////////////////////////////////////////////////////////////////
+      AddJet(jet);
+      
+    }
 
-void AliUA1JetFinder::FindJetsTPC()
+  //delete
+  delete[] ptT;
+  delete[] etaT;
+  delete[] phiT;
+  delete[] injet;
+  delete[] injetOk;
+  delete[] areaJet;
+  delete[] areaJetOk;
 
-{
-  // 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* enT  = new Float_t[nIn];
-  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);
-    enT[i]  = lv->Energy();
-    ptT[i]  = lv->Pt();
-    etaT[i] = lv->Eta();
-    phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
-    if (fReader->GetCutFlag(i) == 1) 
-      fLego->Fill(etaT[i], phiT[i], ptT[i]);
+//-----------------------------------------------------------------------
+void AliUA1JetFinder::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
+                                 Float_t* const etJet,Float_t* const etaJet, Float_t* const phiJet,
+                                 Float_t* const etallJet, Int_t* const ncellsJet)
+{
+  // Dump lego
+  AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
+  const Int_t nBinsMax = 120000; // we use a fixed array not to fragment memory 
+  const Int_t nBinEta = header->GetLegoNbinEta();
+  const Int_t nBinPhi = header->GetLegoNbinPhi();
+  if((nBinPhi*nBinEta)>nBinsMax){
+    AliError("Too many bins of the ETA-PHI histogram");
   }
-  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
 
+  Float_t etCell[nBinsMax] = {0.};   //! Cell Energy
+  Float_t etaCell[nBinsMax] = {0.};  //! Cell eta
+  Float_t phiCell[nBinsMax] = {0.};  //! Cell phi
+  Short_t flagCell[nBinsMax] = {0}; //! Cell flag
+  
   Int_t nCell = 0;
   TAxis* xaxis = fLego->GetXaxis();
   TAxis* yaxis = fLego->GetYaxis();
   Float_t e = 0.0;
-  for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
-      for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); 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->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
-  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);
-
-  // sort jets
-  Int_t * idx  = new Int_t[UA1JETS.njet];
-  TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
-
-  // download jet info.   
-  for(Int_t i = 0; i < UA1JETS.njet; i++) {
-    // reject events outside acceptable eta range
-    if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
-       || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
-         continue;
-      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);
+  for (Int_t i = 1; i <= nBinEta; i++) {
+    for (Int_t j = 1; j <= nBinPhi; j++) {
+      e = fLego->GetBinContent(i,j);
+      if (e < 0.0) continue; // don't include this cells
+      Float_t eta  = xaxis->GetBinCenter(i);
+      Float_t phi  = yaxis->GetBinCenter(j);
+      etCell[nCell]  = e;
+      etaCell[nCell] = eta;
+      phiCell[nCell] = phi;
+      flagCell[nCell] = 0; //default
+      nCell++;
+    }
   }
+  // Parameters from header
+  Float_t minmove = header->GetMinMove();
+  Float_t maxmove = header->GetMaxMove();
+  Float_t rc      = header->GetRadius();
+  Float_t etseed  = header->GetEtSeed();
+  // Tmp array of jets form algoritm
+  Float_t etaAlgoJet[kMaxJets] = {0.0};
+  Float_t phiAlgoJet[kMaxJets] = {0.0};
+  Float_t etAlgoJet[kMaxJets] = {0.0};
+  Int_t   ncellsAlgoJet[kMaxJets] = {0};
+
+  // Run algorithm//
   
-  // find multiplicities and relationship jet-particle
-  // find also percentage of pt from pythia
-  Int_t* injet = new Int_t[nIn];
-  Int_t* sflag = new Int_t[nIn];
-  for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
-  Int_t* mult  = new Int_t[fJets->GetNJets()];
-  Int_t* ncell  = new Int_t[fJets->GetNJets()];
-  Float_t* percentage  = new Float_t[fJets->GetNJets()];
-
-  for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
-      Float_t pt_sig = 0.0;
-      mult[i] = 0;
-      ncell[i] = UA1JETS.ncellj[i];
-      for (Int_t j = 0; j < nIn; j++) {
-         Float_t deta = etaT[j] - fJets->GetEta(i);
-         Float_t dphi = phiT[j] - fJets->GetPhi(i);
-         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+  // Sort cells by et
+  Int_t index[nBinsMax];
+  TMath::Sort(nCell, etCell, index);
+  // variable used in centroide loop
+  Float_t eta   = 0.0;
+  Float_t phi   = 0.0;
+  Float_t eta0  = 0.0;
+  Float_t phi0  = 0.0;
+  Float_t etab  = 0.0;
+  Float_t phib  = 0.0;
+  Float_t etas  = 0.0;
+  Float_t phis  = 0.0;
+  Float_t ets   = 0.0;
+  Float_t deta  = 0.0;
+  Float_t dphi  = 0.0;
+  Float_t dr    = 0.0;
+  Float_t etsb  = 0.0;
+  Float_t etasb = 0.0;
+  Float_t phisb = 0.0;
+  Float_t dphib = 0.0;
+
+  for(Int_t icell = 0; icell < nCell; icell++)
+    {
+      Int_t jcell = index[icell];
+      if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
+      if(flagCell[jcell] != 0) continue; // if cell was used before
+      
+      eta  = etaCell[jcell];
+      phi  = phiCell[jcell];
+      eta0 = eta;
+      phi0 = phi;
+      etab = eta;
+      phib = phi;
+      ets  = etCell[jcell];
+      etas = 0.0;
+      phis = 0.0;
+      etsb = ets;
+      etasb = 0.0;
+      phisb = 0.0;
+      for(Int_t kcell =0; kcell < nCell; kcell++)
+       {
+         Int_t lcell = index[kcell];
+         if(lcell == jcell) continue; // cell itself
+         if(flagCell[lcell] != 0) continue; // cell used before
+         if(etCell[lcell] > etCell[jcell]) continue; // can this happen
+         //calculate dr
+         deta = etaCell[lcell] - eta;
+         dphi = TMath::Abs(phiCell[lcell] - phi);
          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] == 0) {
-             injet[j] = -(i+1);
-             if ((fReader->GetCutFlag(j)) == 1 &&
-                 (etaT[j] < fHeader->GetLegoEtaMax()) &&
-                 (etaT[j] > fHeader->GetLegoEtaMin())) {
-                 injet[j] = i+1;
-                 mult[i]++;
-                 if (fReader->GetSignalFlag(j) == 1) {
-                   pt_sig+=ptT[j];
-                   sflag[j]=1;
-                 }
+         dr = TMath::Sqrt(deta * deta + dphi * dphi);
+         if(dr <= rc)
+           {
+             // calculate offset from initiate cell
+             deta = etaCell[lcell] - eta0;
+             dphi = phiCell[lcell] - phi0;
+             if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
+             if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
+             etas = etas + etCell[lcell]*deta;
+             phis = phis + etCell[lcell]*dphi;
+             ets = ets + etCell[lcell];
+             //new weighted eta and phi including this cell
+             eta = eta0 + etas/ets;
+             phi = phi0 + phis/ets;
+             // if cone does not move much, just go to next step
+             dphib = TMath::Abs(phi - phib);
+             if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
+             dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
+             if(dr <= minmove) break;
+             // cone should not move more than max_mov
+             dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
+             if(dr > maxmove){
+               eta = etab;
+               phi = phib;
+               ets = etsb;
+               etas = etasb;
+               phis = phisb;
+             } else { // store this loop information
+               etab=eta;
+               phib=phi;
+               etsb = ets;
+               etasb = etas;
+               phisb = phis;
              }
+           } // inside cone
+       }//end of cells loop looking centroide
+      
+      // Avoid cones overloap (to be implemented in the future)
+      
+      // Flag cells in Rc, estimate total energy in cone
+      Float_t etCone   = 0.0;
+      Int_t   nCellIn  = 0;
+      rc = header->GetRadius();
+
+      for(Int_t ncell =0; ncell < nCell; ncell++)
+       {
+         if(flagCell[ncell] != 0) continue; // cell used before
+         //calculate dr
+         deta = etaCell[ncell] - eta;
+         dphi = phiCell[ncell] - phi;
+         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+         if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+         dr = TMath::Sqrt(deta * deta + dphi * dphi);
+         if(dr <= rc){  // cell in cone
+           flagCell[ncell] = -1;
+           etCone+=etCell[ncell];
+           nCellIn++;
          }
-      }
-      percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
-       ((Double_t) fJets->GetPt(i));    
-  }
-  fJets->SetNCells(ncell);
-  fJets->SetPtFromSignal(percentage);
-  fJets->SetMultiplicities(mult);
-  fJets->SetInJet(injet);
-  fJets->SetEtaIn(etaT);
-  fJets->SetPhiIn(phiT);
-  fJets->SetPtIn(ptT);
-  fJets->SetEtAvg(UA1JETS.etavg);
-  delete ncell;
-  delete enT;
-  delete ptT;
-  delete etaT;
-  delete phiT;
-  delete injet;
-  delete idx;
-  delete mult;
-  delete percentage;
-}
-
-void AliUA1JetFinder::FindJets()
-{ // Find jets with TPC or EMCAL or TPC + EMCAL information
-  // initialize event, look for jets, download jet info
-
-  // 1) transform input to pt,eta,phi plus lego
-    AliJetUnitArray *fUnit = fReader->GetUnitArray();
-    Int_t nIn = fUnit->GetUnitEntries();
-    Int_t fOption = fReader->GetReaderHeader()->GetDetector();
-    Int_t fDebug = fReader->GetReaderHeader()->GetDebug();
-    
-    if(fDebug>1){ 
-       printf("Inside FindJets function ! \n");
-       printf("fOption : %d \n",fOption);
-    }
-    
-    if(fDebug>10){ 
-       cout << "fUnit : " << fUnit << endl;
-       printf("nIn = fUnit->GetUnitEntries() : %d \n",nIn);
-    }
-    
-    if(fDebug>10){
-       printf("     -----------------------------------------------------------------\n");
-       printf("     --- All inputs in fUnitArray ---\n");
-       for(Int_t i=0; i<nIn ; i++){
-           if(fUnit[i].GetUnitEnergy()!=0){
-               printf("     -----------------------------------------------------------------\n");
-               cout << "     |  ID : " << fUnit[i].GetUnitID() << "  |  Eta : " << fUnit[i].GetUnitEta() << "  |  Phi : " << fUnit[i].GetUnitPhi() << "  |  Energy : " << fUnit[i].GetUnitEnergy() << " |" <<endl;
-           }
-       }
-    }
-    
-    if (nIn == 0) return;
-    
-    Int_t nCandidateCut = 0;
-    Int_t nCandidateCut2 = 0;
-    Int_t nCandidate = 0;
-    for (Int_t i = 0; i < nIn; i++){
-       if(fUnit[i].GetUnitEnergy()>0. && fUnit[i].GetUnitCutFlag()==1) {
-           // Number of good tracks/digits which have passed the pt cut
-           nCandidateCut += 1;
        }
-    if(fUnit[i].GetUnitEnergy()>0.) {
-       // Number of good tracks/digits without pt cut
-       nCandidate += 1;
-    }
-    }
-    
-    if(fDebug>=3){
-    cout << "nCandidate : " << nCandidate << endl;
-    cout << "nCandidateCut : " << nCandidateCut << endl;
-    cout << "nCandidateCut2 : " << nCandidateCut2 << endl;
-    }
-    
-  // local arrays for input No Cuts
-  // Both pt < ptMin and pt > ptMin
-    Float_t* enT       = new Float_t[nCandidate];
-    Float_t* ptT       = new Float_t[nCandidate];
-    Float_t* etaT      = new Float_t[nCandidate];
-    Float_t* phiT      = new Float_t[nCandidate];
-    Float_t* detaT     = new Float_t[nCandidate];
-    Float_t* dphiT     = new Float_t[nCandidate];
-    Float_t* cFlagT    = new Float_t[nCandidate];
-    Float_t* cClusterT = new Float_t[nCandidate];
-    Float_t* idT       = new Float_t[nCandidate];
-    Int_t loop1 = 0;
-    
-    fJets->SetNinput(nCandidate);
-    
-    if(fDebug>3){
-       cout << "nCandidate : " << nCandidate << endl;
-       //    cout << "nMultCandidate : " << nMultCandidate << endl;
-    }
-    
-    Float_t *etCell = new Float_t[nIn];   //! Cell Energy - Extracted from UnitArray
-    Float_t *etaCell = new Float_t[nIn];  //! Cell eta - Extracted from UnitArray
-    Float_t *phiCell = new Float_t[nIn];  //! Cell phi - Extracted from UnitArray
-    
-    // Information extracted from fUnitArray
-    for(Int_t i=0; i<nIn; i++) 
-    {
-       if(fUnit[i].GetUnitCutFlag()==1){ 
-           etCell[i] = fUnit[i].GetUnitEnergy();
-           if (etCell[i] > 0.0) etCell[i] -= fHeader->GetMinCellEt();
-           if (etCell[i] < 0.0) etCell[i] = 0.;
-           etaCell[i] = fUnit[i].GetUnitEta();
-           phiCell[i] = fUnit[i].GetUnitPhi();
-       }
-       else {
-           etCell[i]  = 0.;
-           etaCell[i] = fUnit[i].GetUnitEta();
-           phiCell[i] = fUnit[i].GetUnitPhi();
+      
+      // Select jets with et > background
+      // estimate max fluctuation of background in cone
+      Double_t ncellin = (Double_t)nCellIn;
+      Double_t ntcell  = (Double_t)nCell;
+      Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
+      // min cone et
+      Double_t etcmin = etCone ;  // could be used etCone - etmin !!
+      //decisions !! etbmax < etcmin
+      
+      for(Int_t mcell =0; mcell < nCell; mcell++){
+       if(flagCell[mcell] == -1){
+         if(etbmax < etcmin)
+           flagCell[mcell] = 1; //flag cell as used
+         else
+           flagCell[mcell] = 0; // leave it free
        }
-       
-       if(fUnit[i].GetUnitEnergy()>0.){
-           ptT[loop1]   = fUnit[i].GetUnitEnergy();
-           enT[loop1]   = fUnit[i].GetUnitEnergy();
-           etaT[loop1]  = fUnit[i].GetUnitEta();
-           phiT[loop1]  = fUnit[i].GetUnitPhi();
-           detaT[loop1] = fUnit[i].GetUnitDeta();
-           dphiT[loop1] = fUnit[i].GetUnitDphi();
-           cFlagT[loop1]= fUnit[i].GetUnitCutFlag();
-           idT[loop1]   = fUnit[i].GetUnitID();
-           loop1++;
-       }
-    }
-    
-    
-    if(fDebug > 40) // For comparison
-    {
-       for(Int_t j=0; j<nIn; j++) {
-           if(etCell[j]>0){
-               cout << "etCell[" << j << "] : " << etCell[j] << endl;
-               cout << "etaCell[" << j << "] : " << etaCell[j] << endl;
-               cout << "phiCell[" << j << "] : " << phiCell[j] << endl;
-           }
+      }
+      //store tmp jet info !!!
+      if(etbmax < etcmin) {
+        if(nJets<kMaxJets){
+         etaAlgoJet[nJets]    = eta;
+         phiAlgoJet[nJets]    = phi;
+         etAlgoJet[nJets]     = etCone;
+         ncellsAlgoJet[nJets] = nCellIn;
+         nJets++;
+        }
+       else{
+         AliError(Form("Too many jets (> %d) found by UA1JetFinder, adapt your cuts",kMaxJets));
+         break;
        }
-    }
-    
-    // Run the algo. Parameters from header
-    //  Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
-    Int_t nTot      = nIn;
-    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(nIn, nTot, etCell, etaCell, phiCell, 
-                  minmove, maxmove, mode, precbg, ierror);
-    
-    // sort jets
-    Int_t * idx  = new Int_t[UA1JETS.njet];
-    TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
-    
-    if(fDebug > 20)
+      }      
+    } // end of cells loop
+
+  //reorder jets by et in cone
+  //sort jets by energy
+  for(Int_t p = 0; p < nJets; p++)
     {
-       for(Int_t i = 0; i < UA1JETS.njet; i++) 
-       {
-           cout << "Number of jets found, UA1JETS.njet : " << UA1JETS.njet << endl;
-           cout << "UA1JETS.etj : " << UA1JETS.etj << endl;
-           cout << "idx[" << i << "] : " << idx[i] << endl;
-           cout << "UA1JETS.etaj[1][" << idx[i] << "] : " << UA1JETS.etaj[1][idx[i]] << endl;
-           cout << "UA1JETS.phij[1][" << idx[i] << "] : " << UA1JETS.phij[1][idx[i]] << endl;
-           cout << "UA1JETS.etj[" << idx[i] << "] : " << UA1JETS.etj[idx[i]] << endl;
-       }
-    }
-    
-    // download jet info.   
-    for(Int_t i = 0; i < UA1JETS.njet; i++) {
-       // reject events outside acceptable eta range
-       if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
-           || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
-       {
-           continue;
-       }
-       
-       Float_t px, py,pz,en,pT; // 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);
-       pT = TMath::Sqrt(px * px + py * py);
-       
-       fJets->AddJet(px, py, pz, en);
-       
+      etaJet[p]    = etaAlgoJet[p];
+      phiJet[p]    = phiAlgoJet[p];
+      etJet[p]     = etAlgoJet[p];
+      etallJet[p]  = etAlgoJet[p];
+      ncellsJet[p] = ncellsAlgoJet[p];
     }
 
-    // find multiplicities and relationship jet-particle
-    // find also percentage of pt from pythia
-    
-    Int_t* injet = new Int_t[nCandidate];
-    Int_t* sflag = new Int_t[nCandidate];
-    for (Int_t i = 0; i < nCandidate; i++) {injet[i]= 0;sflag[i]=0;}
-    Int_t* mult  = new Int_t[fJets->GetNJets()];
-    Int_t* ncell  = new Int_t[fJets->GetNJets()];
-    Float_t* percentage  = new Float_t[fJets->GetNJets()];
-    
-    // Instead of using etaT below, it would be interesting to use the previous fUnitArray object
-    // With the particle ID, it is possible to to have access to its physical properties and one can,
-    // for example, set if the corresponding particle is inside or outside the jet with the flag 
-    // kOutJet/kInJet, other possibilities...
-    
-    for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
-       Float_t pt_sig = 0.0;
-       mult[i] = 0;
-       ncell[i] = UA1JETS.ncellj[i];
-       for (Int_t j = 0; j < nCandidate; j++) {
-           Float_t deta = etaT[j] - fJets->GetEta(i);
-           Float_t dphi = phiT[j] - fJets->GetPhi(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] == 0) {
-               injet[j] = -(i+1);
-               if(cFlagT[j] == 1 &&
-                  (etaT[j] < fHeader->GetLegoEtaMax()) &&
-                  (etaT[j] > fHeader->GetLegoEtaMin())) {
-                   injet[j] = i+1;
-                   mult[i]++;
-                   pt_sig+=enT[j];
-                   sflag[j]=1;
-               }
-           }
-           if(fDebug>10){
-               cout << "mult[" << i << "] : " << mult[i] << endl;
-               cout << "ncell[" << i << "] : " << ncell[i] << endl;
-           }
-       }
-       percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
-           ((Double_t) fJets->GetPt(i));    
-    }
-    
-    fJets->SetNCells(ncell);
-    fJets->SetPtFromSignal(percentage);
-    fJets->SetMultiplicities(mult);
-    fJets->SetInJet(injet);
-    fJets->SetEtaIn(etaT);
-    if(fDebug>10){
-       for(Int_t i=0; i<nCandidate ; i++){
-           cout << "phiT[" << i << "] : " << phiT[i] << endl;
-           cout << "etaT[" << i << "] : " << etaT[i] << endl;
-       }
-    }
-    fJets->SetPhiIn(phiT);
-    fJets->SetPtIn(enT);
-    fJets->SetEtAvg(UA1JETS.etavg);
-    delete etCell;
-    delete etaCell;
-    delete phiCell;
-    delete ncell;
-    delete cFlagT;
-    delete cClusterT;
-    delete enT;
-    delete ptT;
-    delete etaT;
-    delete phiT;
-    delete detaT;
-    delete dphiT;
-    delete injet;
-    delete idx;
-    delete mult;
-    delete percentage;
-    
 }
 
-
-////////////////////////////////////////////////////////////////////////
-
+//-----------------------------------------------------------------------
 void AliUA1JetFinder::Reset()
 {
   fLego->Reset();
-  fJets->ClearJets();
-}
+  AliJetFinder::Reset();
 
-////////////////////////////////////////////////////////////////////////
+}
 
-void AliUA1JetFinder::WriteJHeaderToFile()
+//-----------------------------------------------------------------------
+void AliUA1JetFinder::WriteJHeaderToFile() const
 {
-  fOut->cd();
-  fHeader->Write();
-}
+  AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
+  header->Write();
 
-////////////////////////////////////////////////////////////////////////
+}
 
+//-----------------------------------------------------------------------
 void AliUA1JetFinder::Init()
 {
-  // initializes some variables
-  Float_t dEta, dPhi;
-  dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
-    /((Float_t) fHeader->GetLegoNbinEta());
-  dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
-    /((Float_t) fHeader->GetLegoNbinPhi());
-
-  UA1CELL.etaCellSize = dEta;
-  UA1CELL.phiCellSize = dPhi;
-
-  // jet parameters
-  UA1PARA.coneRad = fHeader->GetRadius();
-  UA1PARA.etSeed  = fHeader->GetEtSeed();
-  UA1PARA.ejMin   = fHeader->GetMinJetEt();
-  UA1PARA.etMin   = fHeader->GetMinCellEt();
 
+  // initializes some variables
+  AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
   // book lego
-  fLego = new 
-    TH2F("legoH","eta-phi",
-        fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(), 
-        fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(), 
-        fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
+  fLego = new TH2F("legoH","eta-phi",
+                  header->GetLegoNbinEta(), header->GetLegoEtaMin(),
+                  header->GetLegoEtaMax(),  header->GetLegoNbinPhi(),
+                  header->GetLegoPhiMin(),  header->GetLegoPhiMax());
+  
 }
+