diff --git a/firmware/data.h b/firmware/data.h
index 607254a..249affe 100644
--- a/firmware/data.h
+++ b/firmware/data.h
@@ -17,6 +17,9 @@ typedef ap_uint<14> tk2em_dr_t;
 typedef ap_uint<14> tk2calo_dr_t;
 typedef ap_uint<10> em2calo_dr_t;
 typedef ap_uint<13> tk2calo_dq_t;
+// FIXME: randomly chosen
+typedef ap_uint<3>  hwflags_t;
+
 
 enum PID { PID_Charged=0, PID_Neutral=1, PID_Photon=2, PID_Electron=3, PID_Muon=4 };
 
@@ -27,11 +30,16 @@ enum PID { PID_Charged=0, PID_Neutral=1, PID_Photon=2, PID_Electron=3, PID_Muon=
     #define NMU 4
     #define NSELCALO 20
     #define NALLNEUTRALS NSELCALO
-    // dummy
-    #define NEMCALO 1
+    // FIXME: determine how big this needs to be
+    #define NEMCALO 20
     #define NPHOTON NEMCALO
     // not used but must be there because used in header files
     #define NNEUTRALS 1
+    // Configuration of EG algo follows
+    // FIXME: determine how big this needs to be
+    #define NEMCALOSEL_EGIN 10
+    #define DOBREMRECOVERY
+    #define NEM_EGOUT 5
 //--------------------------------
 #elif defined(REG_HGCalNoTK)
     #define NCALO 12
@@ -100,6 +108,11 @@ enum PID { PID_Charged=0, PID_Neutral=1, PID_Photon=2, PID_Electron=3, PID_Muon=
        #define NSELCALO 10
        #define NALLNEUTRALS (NPHOTON+NSELCALO)
        #define NNEUTRALS 25
+       // EG configuration
+       #define NEMCALOSEL_EGIN 10
+       // #define DOBREMRECOVERY
+       #define NEM_EGOUT 10
+
    #endif
 
 #endif // region
@@ -141,9 +154,11 @@ struct EmCaloObj {
 	pt_t hwPt, hwPtErr;
 	eta_t hwEta; // relative to the region center, at calo
 	phi_t hwPhi; // relative to the region center, at calo
+  hwflags_t hwFlags;
+
 };
 inline void clear(EmCaloObj & c) {
-    c.hwPt = 0; c.hwPtErr = 0; c.hwEta = 0; c.hwPhi = 0; 
+    c.hwPt = 0; c.hwPtErr = 0; c.hwEta = 0; c.hwPhi = 0; c.hwFlags = 0;
 }
 
 
@@ -252,9 +267,37 @@ inline void fill(PuppiObj &out, const HadCaloObj &src, pt_t puppiPt, puppiWgt_t
     out.setHwPuppiW(puppiWgt);
 }
 
+struct EGIsoEleParticle {
+  pt_t hwPt;
+  eta_t hwEta;  // at calo face
+  phi_t hwPhi;
+  // uint16_t hwQual;
 
+  // track parameters for electrons
+  eta_t hwVtxEta;
+  phi_t hwVtxPhi;
+  z0_t hwZ0;
+  bool hwCharge;
+  // uint16_t hwIso;
+};
 
+inline void clear(EGIsoEleParticle & c) {
+    c.hwPt = 0; c.hwEta = 0; c.hwPhi = 0;
+    c.hwVtxEta = 0; c.hwVtxPhi = 0; c.hwZ0=0;
+    c.hwCharge = false;
+}
 
+struct EGIsoParticle {
+  pt_t hwPt;
+  eta_t hwEta;  // at calo face
+  phi_t hwPhi;
+  // uint16_t hwQual;
+  // uint16_t hwIso;
+};
+
+inline void clear(EGIsoParticle & c) {
+    c.hwPt = 0; c.hwEta = 0; c.hwPhi = 0;
+}
 
 //TMUX
 #define NETA_TMUX 2
diff --git a/firmware/pfalgo_common.icc b/firmware/pfalgo_common.icc
index 00b1611..c906c84 100644
--- a/firmware/pfalgo_common.icc
+++ b/firmware/pfalgo_common.icc
@@ -60,6 +60,7 @@ ap_uint<NB> dr2_plus_dpt_int_cap(int dr2, pt_t pt1, pt_t pt2, PTS_t ptscale, ap_
 
 template<typename T, int NIn, int NOut>
 void ptsort_hwopt(const T in[NIn], T out[NOut]) {
+  #pragma HLS INLINE
     T tmp[NOut];
     #pragma HLS ARRAY_PARTITION variable=tmp complete
 
@@ -218,6 +219,7 @@ struct picker_of_closest<MATCH_T,1> {
 
 template<int DR2MAX, int NCA, typename DR_T>
 ap_uint<NCA>  pick_closest(const DR_T calo_track_drval[NCA]) {
+  #pragma HLS INLINE
     const DR_T eDR2MAX = DR2MAX;
     typedef match<DR_T,ap_uint<ct_log2_ceil<NCA>::value>> match_t;
     match_t candidates[NCA];
@@ -355,5 +357,3 @@ void pfmualgo(const MuObj mu[NMU], const TkObj track[NTRACK], const ap_uint<NMU>
         }
     }
 }
-
-
diff --git a/firmware/pftkegalgo.cpp b/firmware/pftkegalgo.cpp
new file mode 100644
index 0000000..343b88d
--- /dev/null
+++ b/firmware/pftkegalgo.cpp
@@ -0,0 +1,241 @@
+#include "pftkegalgo.h"
+#include <cassert>
+#include "hls_math.h"
+#include "pfalgo_common.icc"
+
+ap_int<pt_t::width+1> ell_dpt_int_cap(eta_t eta1, phi_t phi1, eta_t eta2, phi_t phi2, pt_t pt1, pt_t pt2, ap_int<pt_t::width+1> max) {
+#pragma HLS INLINE
+
+#if defined(REG_HGCal)
+    // const ap_uint<10> cdeta_int = 16;
+    const ap_ufixed<12,10> cdeta = 16;
+#endif
+
+#if defined(REG_Barrel)
+    // const ap_uint<10> cdeta_int = (hls::abs(eta1) > 3) ? 22 : 8;
+    // //FIXME: global absolute eta > 0.9 
+    // const ap_ufixed<10,12> cdeta = (hls::abs(eta1) > 3) ? 21.75 : 7.75;
+    // const ap_uint<10> cdeta = 22;
+    const ap_ufixed<12,10> cdeta =  21.75;
+#endif
+
+    const ap_uint<10> cm = 256;
+    // const ap_uint<10> cm_int = 256;
+
+
+    ap_int<eta_t::width+1> d_eta = (eta1-eta2);
+    ap_int<phi_t::width+1> d_phi = (phi1-phi2);
+
+    // ap_uint<22> ell = d_phi*d_phi + d_eta*d_eta*cdeta;
+    ap_uint<22> ell = d_phi*d_phi + d_eta*d_eta*cdeta;
+    // int ell_int = d_phi*d_phi + d_eta*d_eta*cdeta_int;
+    // int pippo = d_eta*d_eta*cdeta;
+    // int pippo2 = d_eta*d_eta*cdeta_int;
+
+    // std::cout << "cdeta: " << cdeta <<std::endl;
+    // std::cout << "d_eta: " << d_eta << " d_eta*d_eta: " << d_eta*d_eta << " d_eta*d_eta*cdeta: " << d_eta*d_eta*cdeta << " d_eta*d_eta*cdeta_int: "<< d_eta*d_eta*cdeta_int << std::endl;
+    // std::cout << "pippo: " << pippo << " pippo_int: "  << pippo2 << std::endl;
+    // std::cout << "d_phi: " << d_phi << "d_phi*d_phi: " << d_phi*d_phi << std::endl;
+    // std::cout << "ell: " << ell << " ell_int: " << ell_int << " eval: " << (ell <= cm) << " eval_int: " << (ell_int <= int(cm_int)) << std::endl;
+    // FIXME: should be ap_uint<pt_t::width> ?
+    ap_int<pt_t::width+1> d_pt = hls::abs(pt1 - pt2);
+    return (ell <= int(cm)) ? d_pt : max;
+}
+
+
+
+template<int DPTMAX>
+void calo2tk_ellipticdptvals(const EmCaloObj &em, const TkObj track[NTRACK], ap_int<pt_t::width+1> calo_track_dptval[NTRACK]) {
+// #pragma HLS INLINE
+    // ap_int<pt_t::width+1>  ?
+    // FIXME: uint<pt_t::width>
+    const ap_int<pt_t::width+1> eDPTMAX = DPTMAX;
+    const pt_t trkQualityPtMin_ = 40; // 10 GeV
+
+    track_loop: for (int itk = 0; itk < NTRACK; ++itk) {
+      if (track[itk].hwPt < trkQualityPtMin_ || em.hwPt == 0) {
+        calo_track_dptval[itk] = eDPTMAX;
+      } else {
+        calo_track_dptval[itk] = ell_dpt_int_cap(em.hwEta, em.hwPhi, track[itk].hwEta, track[itk].hwPhi, em.hwPt, track[itk].hwPt, eDPTMAX);
+        // std::cout << "[" << itk << "] dpt: " << calo_track_dptval[itk] << std::endl;
+      }
+    }
+
+}
+
+
+
+void link_emCalo2emCalo(const EmCaloObj emcalo[NEMCALOSEL_EGIN], ap_uint<NEMCALOSEL_EGIN> emCalo2emcalo_bit[NEMCALOSEL_EGIN]) {
+  #pragma HLS ARRAY_PARTITION variable=emcalo complete dim=1
+  #pragma HLS ARRAY_PARTITION variable=emCalo2emcalo_bit complete dim=1
+  #pragma HLS INLINE
+
+  const ap_int<eta_t::width+1>  dEtaMaxBrem_ = 5; // 0.02; -> round(0.02*4*180/3.14)
+  const ap_int<phi_t::width+1>  dPhiMaxBrem_ = 23; // 0.1; -> round(0.1*4*180/3.14)
+
+  // NOTE: we assume the input to be sorted!!!
+  brem_reco_outer_loop: for (int ic = 0; ic < NEMCALOSEL_EGIN; ++ic) {
+    auto &calo = emcalo[ic];
+    brem_reco_inner_loop: for (int jc = ic + 1; jc < NEMCALOSEL_EGIN; ++jc) {
+      auto &otherCalo = emcalo[jc];
+      if (calo.hwPt != 0 && otherCalo.hwPt != 0 &&
+        hls::abs(otherCalo.hwEta - calo.hwEta) < dEtaMaxBrem_ &&
+          hls::abs(otherCalo.hwPhi - calo.hwPhi) < dPhiMaxBrem_) {
+            emCalo2emcalo_bit[ic][jc] = 1;
+            emCalo2emcalo_bit[jc][jc] = 1; // use diagonal bit to mark the cluster as already used
+      }
+    }
+  }
+}
+
+
+
+
+
+void link_emCalo2tk(const EmCaloObj emcalo[NEMCALOSEL_EGIN],
+                    const TkObj track[NTRACK],
+                    ap_uint<NTRACK> emCalo2tk_bit[NEMCALOSEL_EGIN]) {
+  #pragma HLS INLINE
+
+  #pragma HLS ARRAY_PARTITION variable=emcalo complete dim=1
+  #pragma HLS ARRAY_PARTITION variable=track complete dim=1
+  #pragma HLS ARRAY_PARTITION variable=emCalo2tk_bit complete dim=1
+
+  // FIXME: parametrize
+  const int DPTMAX = 65535; // DPT = 16+1 bits -> max = 2^16-1 = ((1<<16)-1)
+  calo_loop: for (int ic = 0; ic < NEMCALOSEL_EGIN; ++ic) {
+    ap_int<pt_t::width+1> dptvals[NTRACK];
+    calo2tk_ellipticdptvals<DPTMAX>(emcalo[ic], track, dptvals);
+    emCalo2tk_bit[ic] = pick_closest<DPTMAX,NTRACK,ap_int<pt_t::width+1>>(dptvals);
+  }
+
+}
+
+#if defined(REG_HGCal)
+void sel_emCalo(const EmCaloObj emcalo[NEMCALO], EmCaloObj emcalo_sel[NEMCALOSEL_EGIN]) {
+  #pragma HLS INLINE
+
+  EmCaloObj emcalo_sel_temp[NEMCALO];
+  #pragma HLS ARRAY_PARTITION variable=emcalo_sel complete dim=1
+  #pragma HLS ARRAY_PARTITION variable=emcalo_sel_temp complete dim=1
+  EmCaloObj emcalo_zero;
+  clear(emcalo_zero);
+  in_select_loop: for(int ic = 0; ic < NEMCALO; ++ic) {
+    emcalo_sel_temp[ic] = (emcalo[ic].hwFlags == 4) ? emcalo[ic] : emcalo_zero;
+  }
+  ptsort_hwopt<EmCaloObj,NEMCALO,NEMCALOSEL_EGIN>(emcalo_sel_temp, emcalo_sel);
+}
+#endif
+
+
+#if defined(REG_Barrel)
+// NOTE: for now this is a placeholder more than anything else
+void sel_emCalo(const EmCaloObj emcalo[NEMCALO], EmCaloObj emcalo_sel[NEMCALOSEL_EGIN]) {
+  #pragma HLS INLINE
+
+  EmCaloObj emcalo_sel_temp[NEMCALO];
+  #pragma HLS ARRAY_PARTITION variable=emcalo_sel complete dim=1
+  #pragma HLS ARRAY_PARTITION variable=emcalo_sel_temp complete dim=1
+  EmCaloObj emcalo_zero;
+  clear(emcalo_zero);
+  in_select_loop: for(int ic = 0; ic < NEMCALO; ++ic) {
+    // we require pt>2GeV
+    emcalo_sel_temp[ic] = (emcalo[ic].hwPt > 8) ? emcalo[ic] : emcalo_zero;
+  }
+  ptsort_hwopt<EmCaloObj,NEMCALO,NEMCALOSEL_EGIN>(emcalo_sel_temp, emcalo_sel);
+}
+#endif
+
+
+
+void pftkegalgo(const EmCaloObj emcalo[NCALO], const TkObj track[NTRACK],
+  EGIsoParticle photons[NEM_EGOUT], EGIsoEleParticle eles[NEM_EGOUT]) {
+    #ifdef HLS_pipeline_II
+     #if HLS_pipeline_II == 1
+        #pragma HLS pipeline II=1
+     #elif HLS_pipeline_II == 2
+        #pragma HLS pipeline II=2
+     #elif HLS_pipeline_II == 3
+        #pragma HLS pipeline II=3
+     #elif HLS_pipeline_II == 4
+        #pragma HLS pipeline II=4
+     #elif HLS_pipeline_II == 6
+        #pragma HLS pipeline II=6
+     #endif
+    #else
+        #pragma HLS pipeline II=2
+    #endif
+  // #pragma HLS PIPELINE II=HLS_pipeline_II
+  #pragma HLS ARRAY_PARTITION variable=emcalo complete dim=1
+  #pragma HLS ARRAY_PARTITION variable=track complete dim=1
+
+  #pragma HLS ARRAY_PARTITION variable=photons complete dim=1
+  #pragma HLS ARRAY_PARTITION variable=eles complete dim=1
+
+  EmCaloObj emcalo_sel[NEMCALOSEL_EGIN];
+  sel_emCalo(emcalo, emcalo_sel);
+
+  // FIXME: shall we forseen the same selection step for tracks?
+  ap_uint<NTRACK> emCalo2tk_bit[NEMCALOSEL_EGIN];
+  ap_uint<NEMCALOSEL_EGIN> emCalo2emcalo_bit[NEMCALOSEL_EGIN];
+  #pragma HLS ARRAY_PARTITION variable=emCalo2tk_bit complete dim=1
+  #pragma HLS ARRAY_PARTITION variable=emCalo2emcalo_bit complete dim=1
+
+  // initialize
+  init_loop: for (int ic = 0; ic < NEMCALOSEL_EGIN; ++ic) {
+    emCalo2tk_bit[ic] = 0;
+    emCalo2emcalo_bit[ic] = 0;
+  }
+
+  #if defined(DOBREMRECOVERY)
+    link_emCalo2emCalo(emcalo_sel, emCalo2emcalo_bit);
+  #endif
+  //
+  link_emCalo2tk(emcalo_sel, track, emCalo2tk_bit);
+
+
+  EGIsoParticle photons_temp[NEMCALOSEL_EGIN];
+  EGIsoEleParticle eles_temp[NEMCALOSEL_EGIN];
+  #pragma HLS ARRAY_PARTITION variable=photons_temp complete dim=1
+  #pragma HLS ARRAY_PARTITION variable=eles_temp complete dim=1
+
+  int track_id = -1;
+  loop_calo: for (int ic = 0; ic < NEMCALOSEL_EGIN; ++ic) {
+    loop_track_matched: for(int it = 0; it < NTRACK; ++it) {
+      if(emCalo2tk_bit[ic][it]) {
+        track_id = it;
+        break;
+      }
+    }
+
+    pt_t ptcorr = emcalo_sel[ic].hwPt;
+    #if defined(DOBREMRECOVERY)
+      if(emCalo2emcalo_bit[ic][ic] != 1) { // FIXME: we still save the multiclusters used as brem???
+        // FIXME: we should set the quality bit as "brem-recovery performed"
+        loop_calo_brem_reco: for (int ioc = 0; ioc < NEMCALOSEL_EGIN; ++ioc) {
+          if(emCalo2emcalo_bit[ic][ioc]) {
+            ptcorr += emcalo_sel[ioc].hwPt;
+          }
+        }
+      }
+    #endif
+
+    photons_temp[ic].hwPt = ptcorr;
+    photons_temp[ic].hwEta = emcalo_sel[ic].hwEta;
+    photons_temp[ic].hwPhi = emcalo_sel[ic].hwPhi;
+
+    if(emCalo2tk_bit[ic]) {
+      eles_temp[ic].hwPt = ptcorr;
+      eles_temp[ic].hwEta = emcalo_sel[ic].hwEta;
+      eles_temp[ic].hwPhi = emcalo_sel[ic].hwPhi;
+      // FIXME: add track properties @ vertex using track[track_id]
+      eles_temp[ic].hwZ0 = track[track_id].hwZ0;
+    } else {
+      clear(eles_temp[ic]);
+    }
+
+  }
+  ptsort_hwopt<EGIsoParticle,NEMCALOSEL_EGIN,NEM_EGOUT>(photons_temp, photons);
+  ptsort_hwopt<EGIsoEleParticle,NEMCALOSEL_EGIN,NEM_EGOUT>(eles_temp, eles);
+
+}
diff --git a/firmware/pftkegalgo.h b/firmware/pftkegalgo.h
new file mode 100644
index 0000000..429de4e
--- /dev/null
+++ b/firmware/pftkegalgo.h
@@ -0,0 +1,26 @@
+#ifndef FIRMWARE_PFTKEGALGO_H
+#define FIRMWARE_PFTKEGALGO_H
+
+#include "pfalgo_common.h"
+#include "data.h"
+
+void pftkegalgo(const EmCaloObj emcalo[NEMCALO], const TkObj track[NTRACK],
+  EGIsoParticle photons[NEM_EGOUT],
+  EGIsoEleParticle eles[NEM_EGOUT]
+) ;
+
+void link_emCalo2tk(const EmCaloObj emcalo[NEMCALOSEL_EGIN], const TkObj track[NTRACK], ap_uint<NTRACK> emCalo2tk_bit[NEMCALOSEL_EGIN]);
+
+void link_emCalo2emCalo(const EmCaloObj emcalo[NEMCALOSEL_EGIN], ap_uint<NEMCALOSEL_EGIN> emCalo2emcalo_bit[NEMCALOSEL_EGIN]);
+
+void sel_emCalo(const EmCaloObj emcalo[NEMCALO], EmCaloObj emcalo_sel[NEMCALOSEL_EGIN]);
+
+// FIXME: I added this to have pfalgo_common.icc compile...
+#ifndef CMSSW_GIT_HASH
+#define PFALGO_DR2MAX_TK_CALO 525
+#define PFALGO_TK_MAXINVPT_LOOSE    40
+#define PFALGO_TK_MAXINVPT_TIGHT    80
+#endif
+
+
+#endif
diff --git a/pftkegalgo_test.cpp b/pftkegalgo_test.cpp
new file mode 100644
index 0000000..a0a7f6d
--- /dev/null
+++ b/pftkegalgo_test.cpp
@@ -0,0 +1,200 @@
+#include <cstdio>
+#include "firmware/pftkegalgo.h"
+
+#include "ref/pftkegalgo_ref.h"
+
+#include "utils/DiscretePFInputsReader.h"
+#include "utils/pattern_serializer.h"
+#include "utils/test_utils.h"
+
+#define NTEST 10000
+
+bool emcalo_equal(const EmCaloObj &emcalo1, const EmCaloObj &emcalo2) {
+  return emcalo1.hwPt == emcalo2.hwPt && emcalo1.hwEta == emcalo2.hwEta && emcalo1.hwPhi == emcalo2.hwPhi;
+}
+
+bool compare(const EGIsoEleParticle&ele1, const EGIsoEleParticle&ele2) {
+  // FIXME: add other parameters
+  return (ele1.hwPt == ele2.hwPt &&
+          ele1.hwEta == ele2.hwEta &&
+          ele1.hwPhi == ele2.hwPhi &&
+          ele1.hwZ0 == ele2.hwZ0);
+}
+
+bool compare(const EGIsoParticle&ele1, const EGIsoParticle&ele2) {
+  // FIXME: add other parameters
+  return (ele1.hwPt == ele2.hwPt &&
+          ele1.hwEta == ele2.hwEta &&
+          ele1.hwPhi == ele2.hwPhi);
+}
+
+
+int main() {
+    HumanReadablePatternSerializer debugHR("-", /*zerosuppress=*/true); // this will print on stdout, we'll use it for errors
+#if defined(REG_HGCal)
+std::cout << "REG_Hgcal" << std::endl;
+
+    DiscretePFInputsReader inputs("DoubleElectron_PU200_HGCal.dump");
+    bool filterHwQuality = true;
+    bool doBremRecovery = true;
+#endif
+
+#if defined(REG_Barrel)
+std::cout << "REG_Barrel" << std::endl;
+    DiscretePFInputsReader inputs("DoubleElectron_PU200_Barrel.dump");
+    bool filterHwQuality = false;
+    bool doBremRecovery = false;
+#endif
+    // input TP objects
+    HadCaloObj calo[NCALO]; EmCaloObj emcalo[NEMCALO]; TkObj track[NTRACK]; z0_t hwZPV;
+    MuObj mu[NMU];
+    float reg_eta_center;
+    // output PF objects
+    PFChargedObj outch[NTRACK], outch_ref[NTRACK];
+    PFNeutralObj outne[NSELCALO], outne_ref[NSELCALO];
+    PFChargedObj outmupf[NMU], outmupf_ref[NMU];
+
+    pFTkEGAlgo::pftkegalgo_config cfg(NTRACK, NEMCALO, NEMCALOSEL_EGIN, NEM_EGOUT, filterHwQuality, doBremRecovery);
+    
+    pFTkEGAlgo::PFTkEGAlgo algo(cfg);
+    // -----------------------------------------
+    // run multiple tests
+    int passed_tests = 0;
+    int nEles_tests = 0;
+    int nPhotons_tests = 0;
+
+    for (int test = 1; test <= NTEST; ++test) {
+      // std::cout << "------------ TEST: " << test << std::endl;
+        // get the inputs from the input object
+        if (!inputs.nextRegion(calo, emcalo, track, mu, hwZPV, reg_eta_center)) break;
+
+        // for(auto region: inputs.event().regions) {
+        //   printf("# of EM objects: %u\n", region.emcalo.size());
+        //   for(auto emc: region.emcalo) {
+        //     printf("   - hwFlags: %u\n", emc.hwFlags);
+        //   }
+        // }
+        algo.setRegionEtaCenter(reg_eta_center);
+        std::cout<< "Region eta center: " <<   reg_eta_center << std::endl;      
+        // 0 - input cluster selection
+        EmCaloObj emcalo_sel[NEMCALOSEL_EGIN];
+        EmCaloObj emcalo_sel_ref[NEMCALOSEL_EGIN];
+        sel_emCalo(emcalo, emcalo_sel);
+        algo.sel_emCalo_ref(emcalo, emcalo_sel_ref);
+
+
+        // FIXME: determine output depth
+        int emCalo2tk_ref[NEMCALOSEL_EGIN], emCalo2tk[NEMCALOSEL_EGIN];
+        for(int id = 0; id < NEMCALOSEL_EGIN; id++) {
+          emCalo2tk_ref[id] = -1;
+          emCalo2tk[id] = -1;
+        }
+
+
+        int emCalo2emcalo_ref[NEMCALOSEL_EGIN], emCalo2emcalo[NEMCALOSEL_EGIN];
+        for(int id = 0; id < NEMCALOSEL_EGIN; id++) {
+          emCalo2emcalo_ref[id] = -1;
+          emCalo2emcalo[id] = -1;
+        }
+
+        // pftkegalgo(emcalo, track);
+        EGIsoParticle egphs_ref[NEM_EGOUT];
+        EGIsoEleParticle egele_ref[NEM_EGOUT];
+        algo.pftkegalgo_ref(cfg, emcalo_sel_ref, track, emCalo2tk_ref, emCalo2emcalo_ref, egphs_ref, egele_ref);
+
+        ap_uint<NTRACK> emCalo2tk_bit[NEMCALOSEL_EGIN];
+        ap_uint<NEMCALOSEL_EGIN> emCalo2emcalo_bit[NEMCALOSEL_EGIN];
+        EGIsoParticle egphs[NEM_EGOUT];
+        EGIsoEleParticle egele[NEM_EGOUT];
+
+        for(int i = 0; i < NEMCALOSEL_EGIN; i++) {
+          emCalo2tk_bit[i] = 0;
+          emCalo2emcalo_bit[i] = 0;
+        }
+
+        #ifndef REG_Barrel
+        link_emCalo2emCalo(emcalo_sel, emCalo2emcalo_bit);
+        #endif
+        
+        link_emCalo2tk(emcalo_sel, track, emCalo2tk_bit);
+
+        pftkegalgo(emcalo, track, egphs, egele) ;
+
+        for(int ic = 0; ic < NEMCALOSEL_EGIN; ic++) {
+          for(int it = 0; it < NTRACK; it++) {
+            if(emCalo2tk_bit[ic][it])
+              emCalo2tk[ic] = it;
+          }
+        }
+
+        for(int ic = 0; ic < NEMCALOSEL_EGIN; ic++) {
+          if(emCalo2emcalo_bit[ic][ic]) continue;
+          for(int it = 0; it < NEMCALOSEL_EGIN; it++) {
+            if(emCalo2emcalo_bit[ic][it])
+              emCalo2emcalo[it] = ic;
+          }
+        }
+
+        // -----------------------------------------
+        // validation against the reference algorithm
+        int errors = 0;
+        // int ntot = 0, nch = 0, nneu = 0, nmu = 0;
+        for(int i = 0; i < NEMCALOSEL_EGIN; i++) {
+          if (!emcalo_equal(emcalo_sel_ref[i], emcalo_sel[i])) {
+            std::cout << "emcalo_sel [" << i << "]" << " REF: " << emcalo_sel_ref[i].hwPt << " FW: " << emcalo_sel[i].hwPt << std::endl;
+          errors++;
+          }
+        }
+
+        for(int i = 0; i < NEMCALOSEL_EGIN; i++) {
+          if (emCalo2tk_ref[i] != emCalo2tk[i]) {
+            std::cout << "emCalo2tk [" << i << "]" << " REF: " << emCalo2tk_ref[i] << " FW: " << emCalo2tk[i] << std::endl;
+            errors++;
+          }
+        }
+
+        for(int i = 0; i < NEMCALOSEL_EGIN; i++) {
+          if (emCalo2emcalo_ref[i] != emCalo2emcalo[i]) {
+            std::cout << "emCalo2emcalo [" << i << "]" << " REF: " << emCalo2emcalo_ref[i] << " FW: " << emCalo2emcalo[i] << std::endl;
+            errors++;
+          }
+        }
+
+        for(int it = 0; it < NEM_EGOUT; it++) {
+          if(egele_ref[it].hwPt > 0) {
+            nEles_tests++;
+            if(!compare(egele_ref[it], egele[it])) {
+              errors++;
+              std::cout << "[" << it << "] REF ele pt: " << egele_ref[it].hwPt << " FW ele pt: " << egele[it].hwPt << std::endl;
+              std::cout << "    REF ele Z0: " << egele_ref[it].hwZ0 << " FW ele z0: " << egele[it].hwZ0 << std::endl;
+            }
+          }
+        }
+
+        for(int it = 0; it < NEM_EGOUT; it++) {
+          if(egphs_ref[it].hwPt > 0) {
+            nPhotons_tests++;
+            if(!compare(egphs_ref[it], egphs[it])) {
+              errors++;
+              std::cout << "[" << it << "] REF photon pt: " << egphs_ref[it].hwPt << " FW photon pt: " << egphs[it].hwPt << std::endl;
+            }
+
+          }
+        }
+
+
+        if (errors != 0) {
+            printf("Error in computing test %d (%d)\n", test, errors);
+        //     printf("Inputs: \n"); debugHR.dump_inputs(calo, track, mu);
+        //     printf("Reference output: \n"); debugHR.dump_outputs(outch_ref, outne_ref, outmupf_ref);
+        //     printf("Current output: \n"); debugHR.dump_outputs(outch, outne, outmupf);
+            return 1;
+        } else {
+          passed_tests++;
+            // printf("Passed test %d\n", test);
+        }
+    }
+    printf("Passed %d tests (%d electrons, %d photons)\n", passed_tests, nEles_tests, nPhotons_tests);
+
+    return 0;
+}
diff --git a/ref/pftkegalgo_ref.cpp b/ref/pftkegalgo_ref.cpp
new file mode 100644
index 0000000..253bee9
--- /dev/null
+++ b/ref/pftkegalgo_ref.cpp
@@ -0,0 +1,320 @@
+#include "pftkegalgo_ref.h"
+
+
+#include "../utils/Firmware2DiscretePF.h"
+#include "../utils/DiscretePF2Firmware.h"
+
+#include <cmath>
+#include <cstdio>
+#include <algorithm>
+#include <memory>
+#include <iostream>
+
+
+
+int g_pftkegalgo_debug_ref_ = 0;
+
+using namespace pFTkEGAlgo;
+
+
+float PFTkEGAlgo::deltaPhi(float phi1, float phi2) {
+  // reduce to [-pi,pi]
+  float x = phi1 - phi2;
+  float o2pi = 1. / (2. * M_PI);
+  if (std::abs(x) <= float(M_PI))
+      return x;
+  float n = std::round(x * o2pi);
+  return x - n * float(2. * M_PI);
+}
+
+void PFTkEGAlgo::link_emCalo2emCalo(const std::vector<l1tpf_impl::CaloCluster>& emcalo,
+                        std::vector<int> &emCalo2emCalo) {
+  // NOTE: we assume the input to be sorted!!!
+  for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
+    auto &calo = emcalo[ic];
+    if (filterHwQuality_ && calo.hwFlags != caloHwQual_)
+      continue;
+
+    if (emCalo2emCalo[ic] != -1)
+      continue;
+
+    for (int jc = ic + 1; jc < nc; ++jc) {
+      if (emCalo2emCalo[jc] != -1)
+        continue;
+
+      auto &otherCalo = emcalo[jc];
+      if (filterHwQuality_ && otherCalo.hwFlags != caloHwQual_)
+        continue;
+
+      if (fabs(otherCalo.floatEta() - calo.floatEta()) < dEtaMaxBrem_ &&
+          fabs(deltaPhi(otherCalo.floatPhi(), calo.floatPhi())) < dPhiMaxBrem_) {
+        emCalo2emCalo[jc] = ic;
+      }
+    }
+  }
+}
+
+
+
+void PFTkEGAlgo::link_emCalo2tk(const std::vector<l1tpf_impl::CaloCluster> &emcalo,
+                                const std::vector<l1tpf_impl::PropagatedTrack> &track,
+                                std::vector<int> &emCalo2tk) {
+
+  if (debug_ > 10)
+    std::cout << "[link_emCalo2tk]" << std::endl;
+
+  for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
+    auto &calo = emcalo[ic];
+
+    if (filterHwQuality_ && calo.hwFlags != caloHwQual_)
+      continue;
+
+    if (debug_ > 10) {
+      std::cout << "[REF] EM calo [" << ic << "] hwFlags: " << calo.hwFlags << std::endl;
+      std::cout << "--- calo: pt: " << calo.floatPt() << " eta: " << calo.floatEta() << " phi: " << calo.floatPhi()
+                << std::endl;
+    }
+
+    float dPtMin = 999;
+    for (int itk = 0, ntk = track.size(); itk < ntk; ++itk) {
+      const auto &tk = track[itk];
+      if (tk.floatPt() < trkQualityPtMin_)
+        continue;
+
+      if (debug_ > 10) {
+        std::cout << "[REF] tk [" << itk << "] tk: pt: " << tk.floatPt() << " eta: " << tk.floatEta() << " phi: " << tk.floatPhi() << " z0: " << tk.hwZ0
+        << std::endl;
+      }
+
+      float d_phi = deltaPhi(tk.floatPhi(), calo.floatPhi());
+      float d_eta = tk.floatEta() - calo.floatEta();  // We only use it squared
+
+      auto eta_index =
+          std::distance(
+              absEtaBoundaries_.begin(),
+              std::lower_bound(absEtaBoundaries_.begin(), absEtaBoundaries_.end(), globalAbsEta(calo.floatEta()))) -
+          1;
+      float dEtaMax = dEtaValues_[eta_index];
+      float dPhiMax = dPhiValues_[eta_index];
+      
+      std::cout << "eta: " << calo.floatEta() << " globa eta: " << globalAbsEta(calo.floatEta()) << " deta_max: " << dEtaMax << " dphimax: " << dPhiMax << std::endl;
+      // if (debug_ > 0)
+      //   std::cout << " deta: " << fabs(d_eta) << " dphi: " << d_phi
+      //             << " ell: " << ((d_phi / dPhiMax) * (d_phi / dPhiMax)) + ((d_eta / dEtaMax) * (d_eta / dEtaMax))
+      //             << std::endl;
+      // std::cout << "Global abs eta: " << r.globalAbsEta(calo.floatEta())
+      //           << " abs eta: " << fabs(calo.floatEta()) << std::endl;
+
+      if ((((d_phi / dPhiMax) * (d_phi / dPhiMax)) + ((d_eta / dEtaMax) * (d_eta / dEtaMax))) < 1.) {
+        if (debug_ > 0) {
+              if (debug_ <= 10){
+          std::cout << "[REF] EM calo [" << ic << "] hwFlags: " << calo.hwFlags << std::endl;
+          std::cout << "--- calo: pt: " << calo.floatPt() << " eta: " << calo.floatEta() << " phi: " << calo.floatPhi()
+                    << std::endl;
+          std::cout << "[REF] tk [" << itk << " tk: pt: " << tk.floatPt() << " eta: " << tk.floatEta() << " phi: " << tk.floatPhi() << " z0: " << tk.hwZ0
+                    << std::endl;
+                  }
+          std::cout << "    pass elliptic " << std::endl;
+        }
+        // NOTE: for now we implement only best pt match. This is NOT what is done in the L1TkElectronTrackProducer
+        if (fabs(tk.floatPt() - calo.floatPt()) < dPtMin) {
+          if (debug_ > 0)
+            std::cout << "     best pt match: " << fabs(tk.floatPt() - calo.floatPt()) << std::endl;
+          emCalo2tk[ic] = itk;
+          dPtMin = fabs(tk.floatPt() - calo.floatPt());
+        }
+      }
+    }
+  }
+}
+
+
+
+
+
+// void PFTkEGAlgo::pftkegalgo_ref_set_debug(int debug) { g_pftkegalgo_debug_ref_ = debug; }
+
+bool pFTkEGAlgo::hwpt_sort(const l1tpf_impl::CaloCluster &i, const l1tpf_impl::CaloCluster &j) {
+  return (i.hwPt > j.hwPt);
+}
+
+
+void PFTkEGAlgo::sel_emCalo_ref(const EmCaloObj emcalo[/*cfg.nCALO*/],
+  EmCaloObj emcalo_sel_ref[/*cfg.nCALO*/]) {
+
+
+  std::vector<l1tpf_impl::CaloCluster> emcalos(NEMCALO);
+  for(int ic = 0; ic < NEMCALO; ic++) {
+    // if (emcalo[i].hwPt == 0) continue;
+    l1tpf_impl::CaloCluster c;
+    fw2dpf::convert(emcalo[ic], c);
+    // if(emcalo[ic].hwPt != 0)
+    // std::cout << "[" << ic << "] EmCaloObj pt: " << emcalo[ic].hwPt << " CaloCluster pt: " << c.hwPt << std::endl;
+    emcalos.push_back(c);
+  }
+  std::sort(emcalos.begin(), emcalos.end(), hwpt_sort);
+
+  for(int ic = 0; ic < NEMCALOSEL_EGIN; ic++) {
+    clear(emcalo_sel_ref[ic]);
+  }
+
+  int jc = 0;
+  for(int ic = 0; ic < NEMCALO; ic++) {
+    if(emcalos[ic].hwFlags == 4) {
+      // std::cout << "[" << ic << "] EmCaloObj pt: " << emcalos[ic].hwPt << " hwflags: " << emcalos[ic].hwFlags << std::endl;
+      dpf2fw::convert(emcalos[ic], emcalo_sel_ref[jc++]);
+      if(jc == NEMCALOSEL_EGIN) break;
+    }
+  }
+}
+
+
+
+void PFTkEGAlgo::pftkegalgo_ref(const pftkegalgo_config &cfg,
+                    const EmCaloObj emcalo[/*cfg.nCALO*/],
+                    const TkObj track[/*cfg.nTRACK*/],
+                    int emCalo2tk_ref[],
+                    int emCalo2emcalo_ref[],
+                    EGIsoParticle egphs_ref[],
+                    EGIsoEleParticle egele_ref[]) {
+
+  // printf("[pftkegalgo_ref]\n");
+  //     if (g_pftkegalgo_debug_ref_) {
+  // #ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE
+  pFTkEGAlgo::Region reg;
+
+  for (unsigned int i = 0; i < cfg.nTRACK; ++i) { if (track[i].hwPt == 0) continue;
+      l1tpf_impl::PropagatedTrack tk; fw2dpf::convert(track[i], tk);
+      // printf("FW  \t track %3d: pt %8d [ %7.2f ]  calo eta %+7d [ %+5.2f ]  calo phi %+7d [ %+5.2f ]  calo ptErr %6d [ %7.2f ]   tight %d\n",
+      //                     i, tk.hwPt, tk.floatPt(), tk.hwEta, tk.floatEta(), tk.hwPhi, tk.floatPhi(), tk.hwCaloPtErr, tk.floatCaloPtErr(), int(track[i].hwTightQuality));
+      reg.track.push_back(tk);
+  }
+  for (unsigned int i = 0; i < cfg.nEMCALO; ++i) { if (emcalo[i].hwPt == 0) continue;
+      l1tpf_impl::CaloCluster c; fw2dpf::convert(emcalo[i], c);
+      // printf("FW  \t EM calo  %3d: pt %8d [ %7.2f ]  calo eta %+7d [ %+5.2f ]  calo phi %+7d [ %+5.2f ]  calo emPt %7d [ %7.2f ]   isEM %d \n",
+      //                     i, c.hwPt, c.floatPt(), c.hwEta, c.floatEta(), c.hwPhi, c.floatPhi(), c.hwEmPt, c.floatEmPt(), c.isEM);
+      reg.emcalo.push_back(c);
+  }
+
+  // printf("# tracks: %3d # emcalos: %3d\n", tracks.size(), emcalos.size());
+
+
+  std::vector<int> emCalo2emCalo(reg.emcalo.size(), -1);
+  if (doBremRecovery_)
+    link_emCalo2emCalo(reg.emcalo, emCalo2emCalo);
+  // convert to array
+  for(int id = 0; id < cfg.nEMCALO && id < emCalo2emCalo.size(); id++) {
+      emCalo2emcalo_ref[id] = emCalo2emCalo[id];
+    }
+
+
+  std::vector<int> emCalo2tk(reg.emcalo.size(), -1);
+  link_emCalo2tk(reg.emcalo, reg.track, emCalo2tk);
+
+  // convert to array
+  for(int id = 0; id < cfg.nEMCALO && id < emCalo2tk.size(); id++) {
+    emCalo2tk_ref[id] = emCalo2tk[id];
+  }
+
+
+  eg_algo(reg, emCalo2emCalo, emCalo2tk);
+
+  // std::cout << "# EG photons: " << reg.egphotons.size() << " # eles: " << reg.egeles.size() << std::endl;
+
+  // initialize and fill the output arrays
+  for(int ic = 0; ic != cfg.nEM_EGOUT; ++ic) {
+    clear(egphs_ref[ic]);
+    clear(egele_ref[ic]);
+  }
+
+  for(int ic = 0; ic != reg.egphotons.size(); ++ic) {
+    egphs_ref[ic] = reg.egphotons[ic];
+  }
+  for(int ic = 0; ic != reg.egeles.size(); ++ic) {
+    egele_ref[ic] = reg.egeles[ic];
+  }
+
+}
+
+
+void PFTkEGAlgo::eg_algo(Region &r,
+                         const std::vector<int> &emCalo2emCalo,
+                         const std::vector<int> &emCalo2tk) {
+
+  for (int ic = 0, nc = r.emcalo.size(); ic < nc; ++ic) {
+    auto &calo = r.emcalo[ic];
+    if (filterHwQuality_ && calo.hwFlags != caloHwQual_)
+      continue;
+
+    int itk = emCalo2tk[ic];
+
+
+    // check if brem recovery is on
+    if (!doBremRecovery_) {
+      // 1. create EG objects before brem recovery
+      addEgObjsToPF(r, ic, calo.hwFlags, calo.hwPt, itk);
+      continue;
+    }
+
+    // check if the cluster has already been used in a brem reclustering
+    if (emCalo2emCalo[ic] != -1)
+      continue;
+
+    int ptBremReco = calo.hwPt;
+
+    for (int jc = ic; jc < nc; ++jc) {
+      if (emCalo2emCalo[jc] == ic) {
+        auto &otherCalo = r.emcalo[jc];
+        ptBremReco += otherCalo.hwPt;
+      }
+    }
+
+    // 2. create EG objects with brem recovery
+    // FIXME: duplicating the object is suboptimal but this is done for keeping things as in TDR code...
+    addEgObjsToPF(r, ic, calo.hwFlags + 1, ptBremReco, itk);
+  }
+}
+
+EGIsoParticle &PFTkEGAlgo::addEGIsoToPF(std::vector<EGIsoParticle> &egobjs,
+                                        const l1tpf_impl::CaloCluster &calo,
+                                        const int hwQual,
+                                        const int ptCorr) {
+  EGIsoParticle egiso;
+  egiso.hwPt = ptCorr;
+  egiso.hwEta = calo.hwEta;
+  egiso.hwPhi = calo.hwPhi;
+  // egiso.hwQual = hwQual;
+  // egiso.hwIso = 0;
+  egobjs.push_back(egiso);
+  return egobjs.back();
+}
+
+EGIsoEleParticle &PFTkEGAlgo::addEGIsoEleToPF(std::vector<EGIsoEleParticle> &egobjs,
+                                              const l1tpf_impl::CaloCluster &calo,
+                                              const l1tpf_impl::PropagatedTrack &track,
+                                              const int hwQual,
+                                              const int ptCorr) {
+  EGIsoEleParticle egiso;
+  egiso.hwPt = ptCorr;
+  egiso.hwEta = calo.hwEta;
+  egiso.hwPhi = calo.hwPhi;
+
+  // egiso.hwVtxEta = track.hwVtxEta;
+  // egiso.hwVtxPhi = track.hwVtxPhi;
+  egiso.hwZ0 = track.hwZ0;
+  // egiso.hwCharge = track.hwCharge;
+
+  // egiso.hwQual = hwQual;
+  // egiso.hwIso = 0;
+  egobjs.push_back(egiso);
+  return egobjs.back();
+}
+
+void PFTkEGAlgo::addEgObjsToPF(
+    Region &r, const int calo_idx, const int hwQual, const int ptCorr, const int tk_idx) {
+  EGIsoParticle &egobj = addEGIsoToPF(r.egphotons, r.emcalo[calo_idx], hwQual, ptCorr);
+  if (tk_idx != -1) {
+    // egobj.ele_idx = r.egeles.size();
+    addEGIsoEleToPF(r.egeles, r.emcalo[calo_idx], r.track[tk_idx], hwQual, ptCorr);
+  }
+}
diff --git a/ref/pftkegalgo_ref.h b/ref/pftkegalgo_ref.h
new file mode 100644
index 0000000..c5e78e8
--- /dev/null
+++ b/ref/pftkegalgo_ref.h
@@ -0,0 +1,135 @@
+#ifndef PFTKEGALGO_REF_H
+#define PFTKEGALGO_REF_H
+#ifndef CMSSW_GIT_HASH
+  #include "../DiscretePFInputs.h"
+#else
+  #include "../../interface/DiscretePFInputs.h"
+#endif
+
+#include "../firmware/pfalgo2hgc.h"
+#include "pfalgo_common_ref.h"
+
+
+
+namespace pFTkEGAlgo {
+
+  struct pftkegalgo_config {
+      unsigned int nTRACK;
+      unsigned int nEMCALO;
+      unsigned int nEMCALOSEL_EGIN;
+      unsigned int nEM_EGOUT;
+
+      bool filterHwQuality;
+      bool doBremRecovery;
+      
+      pftkegalgo_config(unsigned int nTrack, 
+                        unsigned int nEmCalo, 
+                        unsigned int nEmCaloSel_in, 
+                        unsigned int nEmOut,
+                        bool filterHwQuality,
+                        bool doBremRecovery) :
+          nTRACK(nTrack), 
+          nEMCALO(nEmCalo), 
+          nEMCALOSEL_EGIN(nEmCaloSel_in), 
+          nEM_EGOUT(nEmOut),
+          filterHwQuality(filterHwQuality),
+          doBremRecovery(doBremRecovery) {}
+  };
+
+  struct Region {
+    std::vector<l1tpf_impl::CaloCluster> emcalo;
+    std::vector<l1tpf_impl::PropagatedTrack> track;
+    std::vector<EGIsoParticle> egphotons;
+    std::vector<EGIsoEleParticle> egeles;
+  };
+
+  // void pftkegalgo_ref_set_debug(int debug) ;
+
+  bool hwpt_sort(const l1tpf_impl::CaloCluster &i, const l1tpf_impl::CaloCluster &j);
+
+
+
+  class PFTkEGAlgo {
+  public:
+    
+    PFTkEGAlgo(const pftkegalgo_config& config) :
+    cfg(config),
+    filterHwQuality_(config.filterHwQuality),
+    doBremRecovery_(config.doBremRecovery) {}
+
+    void pftkegalgo_ref(const pftkegalgo_config &cfg,
+                        const EmCaloObj emcalo[/*cfg.nCALO*/],
+                        const TkObj track[/*cfg.nTRACK*/],
+                        int emCalo2tk_ref[],
+                        int emCalo2emcalo_ref[],
+                        EGIsoParticle egphs_ref[],
+                        EGIsoEleParticle egele_ref[]) ;
+
+    void link_emCalo2emCalo(const std::vector<l1tpf_impl::CaloCluster>& emcalo,
+                            std::vector<int> &emCalo2emCalo);
+
+    void link_emCalo2tk(const std::vector<l1tpf_impl::CaloCluster> &emcalo,
+                        const std::vector<l1tpf_impl::PropagatedTrack> &track,
+                        std::vector<int> &emCalo2tk);
+
+    float deltaPhi(float phi1, float phi2);
+
+
+    void sel_emCalo_ref(const EmCaloObj emcalo[/*cfg.nCALO*/],
+                        EmCaloObj emcalo_sel_ref[/*cfg.nCALO*/]);
+
+    void eg_algo(Region &r,
+                 const std::vector<int> &emCalo2emCalo,
+                 const std::vector<int> &emCalo2tk);
+
+     EGIsoParticle &addEGIsoToPF(std::vector<EGIsoParticle> &egobjs,
+                                 const l1tpf_impl::CaloCluster &calo,
+                                 const int hwQual,
+                                 const int ptCorr);
+
+     EGIsoEleParticle &addEGIsoEleToPF(std::vector<EGIsoEleParticle> &egobjs,
+                                       const l1tpf_impl::CaloCluster &calo,
+                                       const l1tpf_impl::PropagatedTrack &track,
+                                       const int hwQual,
+                                       const int ptCorr);
+
+     void addEgObjsToPF(Region &r, const int calo_idx, const int hwQual, const int ptCorr, const int tk_idx);
+
+     //FIXME: this is hugly but we will soon change the interface
+     void setRegionEtaCenter(float etacenter) {
+       reg_eta_center_ = etacenter;
+     }
+
+   private:
+     
+     
+     float globalAbsEta(float localEta) {
+       return std::abs(localEta + reg_eta_center_);
+     }
+
+     
+     
+     pftkegalgo_config cfg;
+     bool debug_ = 1;
+     bool filterHwQuality_;
+     bool doBremRecovery_;
+
+     int caloHwQual_ = 4;
+     float dEtaMaxBrem_ = 0.02;
+     float dPhiMaxBrem_ = 0.1;
+     // FIXME: need region center in FW, we don't have it for now hence we keep 1 value for the barrel
+     // std::vector<float> absEtaBoundaries_{0.0, 0.9, 1.5};
+     //  // FIXME: should be {0.025, 0.015, 0.01}  but 0.01 it is too small give eta precision. We use  0.0174533 hwEta: 4
+     // std::vector<float> dEtaValues_{0.025, 0.015, 0.0174533};
+     // std::vector<float> dPhiValues_{0.07, 0.07, 0.07};
+     std::vector<float> absEtaBoundaries_{0.0, 1.5};
+      // FIXME: should be {0.025, 0.015, 0.01}  but 0.01 it is too small give eta precision. We use  0.0174533 hwEta: 4
+     std::vector<float> dEtaValues_{0.015, 0.0174533};
+     std::vector<float> dPhiValues_{0.07, 0.07};
+
+     float trkQualityPtMin_ = 10;
+     float reg_eta_center_;
+   };
+}
+
+#endif
diff --git a/run_hls_ptkegalgo.tcl b/run_hls_ptkegalgo.tcl
new file mode 100644
index 0000000..d8dfd77
--- /dev/null
+++ b/run_hls_ptkegalgo.tcl
@@ -0,0 +1,39 @@
+# get the configuration
+#source config_hls_fullpfalgo_mp7.tcl
+set pfBoard "none"
+#set pfBoard "VCU118"
+set pfReg "HGCal"
+set cflags "-std=c++0x -DREG_${pfReg} -DBOARD_${pfBoard} -DHLS_pipeline_II=4"
+
+# open the project, don't forget to reset
+open_project -reset "proj_tkeg${pfReg}_${pfBoard}"
+
+
+if { $pfBoard == "none" } {
+    set_top pftkegalgo
+} else {
+    set_top packed_pftkegalgo
+}
+add_files firmware/pftkegalgo.cpp -cflags "${cflags}"
+add_files -tb pftkegalgo_test.cpp  -cflags "${cflags}"
+add_files -tb ref/pftkegalgo_ref.cpp  -cflags "${cflags}"
+add_files -tb ref/pfalgo_common_ref.cpp  -cflags "${cflags}"
+add_files -tb utils/pattern_serializer.cpp -cflags "${cflags}"
+add_files -tb utils/test_utils.cpp -cflags "${cflags}"
+#add_files -tb data/TTbar_PU200_HGCal.dump
+add_files -tb data/DoubleElectron_PU200_HGCal.dump
+
+# reset the solution
+open_solution -reset "solution"
+set_part {xcvu9p-flga2104-2L-e}
+create_clock -period 3.0 -name default
+
+config_interface -trim_dangling_port
+# do stuff
+csim_design
+csynth_design
+cosim_design -trace_level all
+#export_design -format ip_catalog -vendor "cern-cms" -version ${l1pfIPVersion} -description "${l1pfTopFunc}"
+
+# exit Vivado HLS
+exit
diff --git a/run_hls_ptkegalgo_barrel.tcl b/run_hls_ptkegalgo_barrel.tcl
new file mode 100644
index 0000000..155f807
--- /dev/null
+++ b/run_hls_ptkegalgo_barrel.tcl
@@ -0,0 +1,39 @@
+# get the configuration
+#source config_hls_fullpfalgo_mp7.tcl
+set pfBoard "none"
+#set pfBoard "VCU118"
+set pfReg "Barrel"
+set cflags "-std=c++0x -DREG_${pfReg} -DBOARD_${pfBoard} -DHLS_pipeline_II=2"
+
+# open the project, don't forget to reset
+open_project -reset "proj_tkeg${pfReg}_${pfBoard}"
+
+
+if { $pfBoard == "none" } {
+    set_top pftkegalgo
+} else {
+    set_top packed_pftkegalgo
+}
+add_files firmware/pftkegalgo.cpp -cflags "${cflags}"
+add_files -tb pftkegalgo_test.cpp  -cflags "${cflags}"
+add_files -tb ref/pftkegalgo_ref.cpp  -cflags "${cflags}"
+add_files -tb ref/pfalgo_common_ref.cpp  -cflags "${cflags}"
+add_files -tb utils/pattern_serializer.cpp -cflags "${cflags}"
+add_files -tb utils/test_utils.cpp -cflags "${cflags}"
+#add_files -tb data/TTbar_PU200_HGCal.dump
+add_files -tb data/DoubleElectron_PU200_Barrel.dump
+
+# reset the solution
+open_solution -reset "solution"
+set_part {xcvu9p-flga2104-2L-e}
+create_clock -period 3.0 -name default
+
+config_interface -trim_dangling_port
+# do stuff
+csim_design
+csynth_design
+cosim_design -trace_level all
+#export_design -format ip_catalog -vendor "cern-cms" -version ${l1pfIPVersion} -description "${l1pfTopFunc}"
+
+# exit Vivado HLS
+exit
diff --git a/utils/DiscretePF2Firmware.h b/utils/DiscretePF2Firmware.h
index a438a03..f075fe9 100644
--- a/utils/DiscretePF2Firmware.h
+++ b/utils/DiscretePF2Firmware.h
@@ -36,6 +36,7 @@ namespace dpf2fw {
         out.hwPtErr = in.hwPtErr;
         out.hwEta = in.hwEta;
         out.hwPhi = in.hwPhi;
+        out.hwFlags = in.hwFlags;
     }
     inline void convert(const l1tpf_impl::Muon & in, MuObj & out) {
         out.hwPt = in.hwPt;
diff --git a/utils/DiscretePFInputsReader.h b/utils/DiscretePFInputsReader.h
index 963d591..23456f5 100644
--- a/utils/DiscretePFInputsReader.h
+++ b/utils/DiscretePFInputsReader.h
@@ -48,7 +48,7 @@ class DiscretePFInputsReader {
         }
         ~DiscretePFInputsReader() { fclose(file_); }
         // for region-by-region approach
-        bool nextRegion(HadCaloObj calo[NCALO], EmCaloObj emcalo[NEMCALO], TkObj track[NTRACK], MuObj mu[NMU], z0_t & hwZPV) {
+        bool nextRegion(HadCaloObj calo[NCALO], EmCaloObj emcalo[NEMCALO], TkObj track[NTRACK], MuObj mu[NMU], z0_t & hwZPV, float &reg_eta_center) {
             if (!nextRegion()) return false;
             const Region &r = event_.regions[iregion_];
 
@@ -56,6 +56,7 @@ class DiscretePFInputsReader {
             dpf2fw::convert<NCALO>(r.calo, calo);
             dpf2fw::convert<NEMCALO>(r.emcalo, emcalo);
             dpf2fw::convert<NMU>(r.muon, mu);
+            reg_eta_center = r.etaCenter;
 
 #ifdef __GXX_EXPERIMENTAL_CXX0X__
             hwZPV = event_.z0 * l1tpf_impl::InputTrack::Z0_SCALE;
diff --git a/utils/Firmware2DiscretePF.h b/utils/Firmware2DiscretePF.h
index caf9815..9fe6071 100644
--- a/utils/Firmware2DiscretePF.h
+++ b/utils/Firmware2DiscretePF.h
@@ -128,6 +128,7 @@ namespace fw2dpf {
         out.hwPtErr = in.hwPtErr;
         out.hwEta = in.hwEta;
         out.hwPhi = in.hwPhi;
+        out.hwFlags = in.hwFlags;
         out.src = nullptr;
     }
     inline void convert(const MuObj & in, l1tpf_impl::Muon & out) {