diff --git a/.gitignore b/.gitignore
index b3cc6e200c94d..e36b99196bf8d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,4 +4,4 @@ __init__.py
 .#*
 #*#
 *~
-*.pb
+
diff --git a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py
index 3f65923eec478..fae52b18c6609 100644
--- a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py
+++ b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py
@@ -20,7 +20,8 @@
        'keep *_ticlTrackstersHFNoseHAD_*_*',
        'keep *_ticlTrackstersHFNoseMerge_*_*',] +
       ['keep *_pfTICL_*_*'] +
-      ['keep *_ticlGraph_*_*']
+      [ 'keep *_ticlGraph_*_*'] +
+      ['keep *_muons1stStep_*_*']
       )
     )
 TICL_RECO.outputCommands.extend(TICL_AOD.outputCommands)
@@ -48,6 +49,7 @@ def cleanOutputAndSet(outputModule, ticl_outputCommads):
                                             'keep SimVertexs_g4SimHits_*_*',
                                             'keep *_layerClusterSimClusterAssociationProducer_*_*',
                                             'keep *_layerClusterCaloParticleAssociationProducer_*_*',
+                                            'keep *_randomEngineStateProducer_*_*',
                                             ])
 
     if hasattr(process, 'FEVTDEBUGEventContent'):
diff --git a/RecoHGCal/TICL/data/tf_models/energy_id_v0.pb b/RecoHGCal/TICL/data/tf_models/energy_id_v0.pb
new file mode 100644
index 0000000000000..6684f4bb2c10a
Binary files /dev/null and b/RecoHGCal/TICL/data/tf_models/energy_id_v0.pb differ
diff --git a/RecoHGCal/TICL/data/tf_models/energy_regression_test.pb b/RecoHGCal/TICL/data/tf_models/energy_regression_test.pb
new file mode 100644
index 0000000000000..519391a165252
Binary files /dev/null and b/RecoHGCal/TICL/data/tf_models/energy_regression_test.pb differ
diff --git a/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h b/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h
index c47aca80a7d92..93fb1bfa61778 100644
--- a/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h
+++ b/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h
@@ -38,6 +38,7 @@ namespace ticl {
       const TILES& tiles;
       const std::vector<TICLSeedingRegion>& regions;
       const tensorflow::Session* tfSession;
+      const tensorflow::Session* tfSessionER;
 
       Inputs(const edm::Event& eV,
              const edm::EventSetup& eS,
@@ -46,8 +47,9 @@ namespace ticl {
              const edm::ValueMap<std::pair<float, float>>& lT,
              const TILES& tL,
              const std::vector<TICLSeedingRegion>& rG,
-             const tensorflow::Session* tS)
-          : ev(eV), es(eS), layerClusters(lC), mask(mS), layerClustersTime(lT), tiles(tL), regions(rG), tfSession(tS) {}
+             const tensorflow::Session* tS,
+             const tensorflow::Session* tSER)
+          : ev(eV), es(eS), layerClusters(lC), mask(mS), layerClustersTime(lT), tiles(tL), regions(rG), tfSession(tS), tfSessionER(tSER) {}
     };
 
     virtual void makeTracksters(const Inputs& input,
diff --git a/RecoHGCal/TICL/interface/commons.h b/RecoHGCal/TICL/interface/commons.h
index ec52182bd49d6..286ffb952b92b 100644
--- a/RecoHGCal/TICL/interface/commons.h
+++ b/RecoHGCal/TICL/interface/commons.h
@@ -4,6 +4,7 @@
 #include "DataFormats/Common/interface/Ref.h"
 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
 #include "DataFormats/HGCalReco/interface/Trackster.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
 
 namespace ticl {
 
@@ -70,6 +71,64 @@ namespace ticl {
     result.emplace_back(tmpTrackster);
   }
 
+  enum LayerType {
+
+    CE_E_120 = 0,
+    CE_E_200 = 1,
+    CE_E_300 = 2,
+    CE_H_120_F = 3,
+    CE_H_200_F = 4,
+    CE_H_300_F = 5,
+    CE_H_120_C = 6,
+    CE_H_200_C = 7,
+    CE_H_300_C = 8,
+    CE_H_SCINT_C = 9,
+    EnumSize = 10
+
+  };
+
+  inline int returnIndex(DetId& lc_seed, const hgcal::RecHitTools& rhtools_) {
+    auto layer_number = rhtools_.getLayerWithOffset(lc_seed);
+    auto thickness = rhtools_.getSiThickIndex(lc_seed);
+    auto isEELayer = (layer_number <= rhtools_.lastLayerEE(false));
+    auto isScintillator = rhtools_.isScintillator(lc_seed);
+    auto isFine = (layer_number <= rhtools_.lastLayerEE(false) + 7);
+
+    if (isEELayer) {
+      if (thickness == 0) {
+        return CE_E_120;
+      } else if (thickness == 1) {
+        return CE_E_200;
+      } else if (thickness == 2) {
+        return CE_E_300;
+      }
+    } else if (!isEELayer) {
+      if (isScintillator) {
+        return CE_H_SCINT_C;
+      } else {
+        if (isFine) {
+          if (thickness == 0) {
+            return CE_H_120_F;
+          } else if (thickness == 1) {
+            return CE_H_200_F;
+          } else if (thickness == 2) {
+            return CE_H_300_F;
+          }
+        } else {
+          if (thickness == 0) {
+            return CE_H_120_C;
+          } else if (thickness == 1) {
+            return CE_H_200_C;
+          }
+            else if (thickness == 2) {
+              return CE_H_300_C;
+          }
+        }
+      }
+    }
+    return -1;
+  };
+
 }  // namespace ticl
 
 #endif
diff --git a/RecoHGCal/TICL/plugins/BuildFile.xml b/RecoHGCal/TICL/plugins/BuildFile.xml
index 5ebee71a3509a..73fe82ec76d9a 100644
--- a/RecoHGCal/TICL/plugins/BuildFile.xml
+++ b/RecoHGCal/TICL/plugins/BuildFile.xml
@@ -30,6 +30,7 @@
 <use name="TrackingTools/GeomPropagators"/>
 <use name="TrackingTools/Records"/>
 <use name="TrackingTools/TrajectoryState"/>
+<use name="CommonTools/UtilAlgos"/>
 <use name="fastjet"/>
 <library file="*.cc" name="RecoHGCalTICLPlugins">
   <flags EDM_PLUGIN="1"/>
diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc
index 16b0ab927f483..f2bca68d0da73 100644
--- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc
+++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc
@@ -9,6 +9,7 @@
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "FWCore/Utilities/interface/Exception.h"
 #include "PatternRecognitionbyCLUE3D.h"
+#include "RecoHGCal/TICL/interface/commons.h"
 
 #include "TrackstersPCA.h"
 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
@@ -38,7 +39,9 @@ PatternRecognitionbyCLUE3D<TILES>::PatternRecognitionbyCLUE3D(const edm::Paramet
       outlierMultiplier_(conf.getParameter<double>("outlierMultiplier")),
       minNumLayerCluster_(conf.getParameter<int>("minNumLayerCluster")),
       eidInputName_(conf.getParameter<std::string>("eid_input_name")),
+      eidInputNameER_(conf.getParameter<std::string>("eid_input_nameER")),
       eidOutputNameEnergy_(conf.getParameter<std::string>("eid_output_name_energy")),
+      eidOutputNameEnergyER_(conf.getParameter<std::string>("eid_output_name_energyER")),
       eidOutputNameId_(conf.getParameter<std::string>("eid_output_name_id")),
       eidMinClusterEnergy_(conf.getParameter<double>("eid_min_cluster_energy")),
       eidNLayers_(conf.getParameter<int>("eid_n_layers")),
@@ -327,7 +330,7 @@ void PatternRecognitionbyCLUE3D<TILES>::makeTracksters(
                               rhtools_.getPositionLayer(rhtools_.lastLayerEE(false), false).z());
 
   // run energy regression and ID
-  energyRegressionAndID(input.layerClusters, input.tfSession, result);
+  energyRegressionAndID(input.layerClusters, input.tfSession, input.tfSessionER, result);
   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > PatternRecognitionAlgoBaseT<TILES>::Advanced) {
     for (auto const &t : result) {
       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Barycenter: " << t.barycenter();
@@ -352,6 +355,7 @@ void PatternRecognitionbyCLUE3D<TILES>::makeTracksters(
 template <typename TILES>
 void PatternRecognitionbyCLUE3D<TILES>::energyRegressionAndID(const std::vector<reco::CaloCluster> &layerClusters,
                                                               const tensorflow::Session *eidSession,
+                                                              const tensorflow::Session *eidSessionEnergyRegression,
                                                               std::vector<Trackster> &tracksters) {
   // Energy regression and particle identification strategy:
   //
@@ -418,7 +422,6 @@ void PatternRecognitionbyCLUE3D<TILES>::energyRegressionAndID(const std::vector<
   // fill input tensor (5)
   for (int i = 0; i < batchSize; i++) {
     const Trackster &trackster = tracksters[tracksterIndices[i]];
-
     // per layer, we only consider the first eidNClusters_ clusters in terms of energy, so in order
     // to avoid creating large / nested structures to do the sorting for an unknown number of total
     // clusters, create a sorted list of layer cluster indices to keep track of the filled clusters
@@ -466,16 +469,6 @@ void PatternRecognitionbyCLUE3D<TILES>::energyRegressionAndID(const std::vector<
   // run the inference (7)
   tensorflow::run(const_cast<tensorflow::Session *>(eidSession), inputList, outputNames, &outputs);
 
-  // store regressed energy per trackster (8)
-  if (!eidOutputNameEnergy_.empty()) {
-    // get the pointer to the energy tensor, dimension is batch x 1
-    float *energy = outputs[0].flat<float>().data();
-
-    for (const int &i : tracksterIndices) {
-      tracksters[i].setRegressedEnergy(*(energy++));
-    }
-  }
-
   // store id probabilities per trackster (8)
   if (!eidOutputNameId_.empty()) {
     // get the pointer to the id probability tensor, dimension is batch x id_probabilities.size()
@@ -487,6 +480,69 @@ void PatternRecognitionbyCLUE3D<TILES>::energyRegressionAndID(const std::vector<
       probs += tracksters[i].id_probabilities().size();
     }
   }
+  // NEW ENERGY REGRESSION!
+  std::vector<int> tracksterIndicesHadronic;
+  auto index_trackster = 0;
+  for (auto &t : tracksters) {
+    t.setRegressedEnergy(t.raw_energy());
+    if ((t.id_probability(ticl::Trackster::ParticleType::photon) +
+         t.id_probability(ticl::Trackster::ParticleType::electron)) >= 0.) {
+      tracksterIndicesHadronic.push_back(index_trackster);
+    }
+    index_trackster++;
+  }
+  batchSize = static_cast<int>(tracksterIndicesHadronic.size());
+  if (batchSize == 0) {
+    return;
+  }
+  // create input and output tensors for energy regression
+  tensorflow::TensorShape shapeER({batchSize, eidNFeaturesER_});
+  tensorflow::Tensor inputEnergyRegression(tensorflow::DT_FLOAT, shapeER);
+  tensorflow::NamedTensorList inputListEnergyRegression = {{eidInputNameER_, inputEnergyRegression}};
+
+  std::vector<tensorflow::Tensor> outputsEnergyRegression;
+  std::vector<std::string> outputNamesEnergyRegression;
+
+  outputNamesEnergyRegression.push_back(eidOutputNameEnergyER_);
+
+  // fill input tensor (5)
+  for (int i = 0; i < batchSize; i++) {
+    const Trackster &trackster = tracksters[tracksterIndicesHadronic[i]]; 
+    float en_total = 0;
+    float *featuresER = &inputEnergyRegression.tensor<float, 2>()(i, 0);
+    std::array<float, eidNFeaturesER_> inputsArrayER;
+    for (size_t j = 0; j < inputsArrayER.size(); j++) {
+      inputsArrayER[j] = 0.;
+    }
+
+    for (int k = 0; k < (int)trackster.vertices().size(); k++) {
+      const reco::CaloCluster &cluster = layerClusters[trackster.vertices(k)];
+      auto lc_seed = cluster.seed();
+      auto index = returnIndex(lc_seed, rhtools_);
+        if(index >= 0){
+          inputsArrayER[index] += cluster.energy() / trackster.raw_energy(); //energy fraction
+        }
+      en_total += cluster.energy();
+    }
+    
+    inputsArrayER[ticl::LayerType::EnumSize] = trackster.barycenter().eta(); //eta
+    inputsArrayER[ticl::LayerType::EnumSize+1] = trackster.raw_energy(); //total trackster energy
+    for(size_t i_p = 0; i_p < trackster.id_probabilities().size(); i_p++){
+      inputsArrayER[ticl::LayerType::EnumSize+2+i_p] = trackster.id_probabilities(i_p); //probabilities
+    }
+    for (size_t j = 0; j < inputsArrayER.size(); j++) {
+      *(featuresER++) = inputsArrayER[j];
+    }
+  }
+
+  tensorflow::run(const_cast<tensorflow::Session *>(eidSessionEnergyRegression),
+                  inputListEnergyRegression,
+                  outputNamesEnergyRegression,
+                  &outputsEnergyRegression);
+  float *energy = outputsEnergyRegression[0].flat<float>().data();
+  for (const int &i : tracksterIndicesHadronic) {
+    tracksters[i].setRegressedEnergy(*(energy++));
+  }
 }
 
 template <typename TILES>
@@ -869,7 +925,9 @@ void PatternRecognitionbyCLUE3D<TILES>::fillPSetDescription(edm::ParameterSetDes
   iDesc.add<double>("outlierMultiplier", 2);
   iDesc.add<int>("minNumLayerCluster", 2)->setComment("Not Inclusive");
   iDesc.add<std::string>("eid_input_name", "input");
+  iDesc.add<std::string>("eid_input_nameER", "input");
   iDesc.add<std::string>("eid_output_name_energy", "output/regressed_energy");
+  iDesc.add<std::string>("eid_output_name_energyER", "output/er_coefficients");
   iDesc.add<std::string>("eid_output_name_id", "output/id_probabilities");
   iDesc.add<double>("eid_min_cluster_energy", 1.);
   iDesc.add<int>("eid_n_layers", 50);
diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h
index 8e77ff55e605e..98a94ea980b12 100644
--- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h
+++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h
@@ -5,6 +5,7 @@
 #define __RecoHGCal_TICL_PRbyCLUE3D_H__
 #include <memory>  // unique_ptr
 #include "RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h"
+#include "DataFormats/DetId/interface/DetId.h"
 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
 
 namespace ticl {
@@ -19,6 +20,7 @@ namespace ticl {
                         std::unordered_map<int, std::vector<int>>& seedToTracksterAssociation) override;
 
     void energyRegressionAndID(const std::vector<reco::CaloCluster>& layerClusters,
+                               const tensorflow::Session*,
                                const tensorflow::Session*,
                                std::vector<Trackster>& result);
 
@@ -129,7 +131,9 @@ namespace ticl {
     const int minNumLayerCluster_;
     const std::vector<int> filter_on_categories_;
     const std::string eidInputName_;
+    const std::string eidInputNameER_;
     const std::string eidOutputNameEnergy_;
+    const std::string eidOutputNameEnergyER_;
     const std::string eidOutputNameId_;
     const float eidMinClusterEnergy_;
     const int eidNLayers_;
@@ -137,8 +141,10 @@ namespace ticl {
 
     hgcal::RecHitTools rhtools_;
     tensorflow::Session* eidSession_;
+    tensorflow::Session* eidSessionER_;
 
     static const int eidNFeatures_ = 3;
+    static const int eidNFeaturesER_ = 20;
   };
 
 }  // namespace ticl
diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc
index d28c80c957095..c7fc9b36eb4ef 100644
--- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc
+++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc
@@ -23,6 +23,7 @@
 #include "DataFormats/Math/interface/Vector3D.h"
 
 #include "RecoHGCal/TICL/interface/GlobalCache.h"
+#include "RecoHGCal/TICL/interface/commons.h"
 
 #include "TrackingTools/Records/interface/TfGraphRecord.h"
 #include "RecoTracker/FinalTrackSelectors/interface/TfGraphDefWrapper.h"
@@ -73,7 +74,8 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<> {
 
   void energyRegressionAndID(const std::vector<reco::CaloCluster> &layerClusters,
                              const tensorflow::Session *,
-                             std::vector<Trackster> &result) const;
+                             const tensorflow::Session *,
+                             std::vector<Trackster> &result);
   void printTrackstersDebug(const std::vector<Trackster> &, const char *label) const;
   void assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates) const;
   void dumpTrackster(const Trackster &) const;
@@ -89,8 +91,11 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<> {
   const edm::EDGetTokenT<edm::ValueMap<float>> tracks_time_err_token_;
   const edm::EDGetTokenT<std::vector<reco::Muon>> muons_token_;
   const std::string tfDnnLabel_;
+  const std::string tfDnnLabelER_;
   const edm::ESGetToken<TfGraphDefWrapper, TfGraphRecord> tfDnnToken_;
+  const edm::ESGetToken<TfGraphDefWrapper, TfGraphRecord> tfDnnTokenER_;
   const tensorflow::Session *tfSession_;
+  const tensorflow::Session *tfSessionER_;
 
   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometry_token_;
   const std::string detector_;
@@ -117,7 +122,9 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<> {
   const double resol_calo_scale_em_;
   const bool debug_;
   const std::string eidInputName_;
+  const std::string eidInputNameER_;
   const std::string eidOutputNameEnergy_;
+  const std::string eidOutputNameEnergyER_;
   const std::string eidOutputNameId_;
   const float eidMinClusterEnergy_;
   const int eidNLayers_;
@@ -128,12 +135,12 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<> {
 
   std::unique_ptr<GeomDet> firstDisk_[2];
 
-  tensorflow::Session *eidSession_;
   hgcal::RecHitTools rhtools_;
 
   static constexpr int eidNFeatures_ = 3;
 
   edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> hdc_token_;
+  static const int eidNFeaturesER_ = 20;
 };
 
 TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps)
@@ -147,14 +154,17 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps)
       tracks_time_err_token_(consumes<edm::ValueMap<float>>(ps.getParameter<edm::InputTag>("tracksTimeErr"))),
       muons_token_(consumes<std::vector<reco::Muon>>(ps.getParameter<edm::InputTag>("muons"))),
       tfDnnLabel_(ps.getParameter<std::string>("tfDnnLabel")),
+      tfDnnLabelER_(ps.getParameter<std::string>("tfDnnLabelER")),
       tfDnnToken_(esConsumes(edm::ESInputTag("", tfDnnLabel_))),
+      tfDnnTokenER_(esConsumes(edm::ESInputTag("", tfDnnLabelER_))),
       tfSession_(nullptr),
+      tfSessionER_(nullptr),
       geometry_token_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
       detector_(ps.getParameter<std::string>("detector")),
       propName_(ps.getParameter<std::string>("propagator")),
       bfield_token_(esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>()),
       propagator_token_(
-          esConsumes<Propagator, TrackingComponentsRecord, edm::Transition::BeginRun>(edm::ESInputTag("", propName_))),
+          esConsumes<Propagator, TrackingComponentsRecord, edm::Transition::BeginRun>(edm::ESInputTag("", propName_))),      
       optimiseAcrossTracksters_(ps.getParameter<bool>("optimiseAcrossTracksters")),
       eta_bin_window_(ps.getParameter<int>("eta_bin_window")),
       phi_bin_window_(ps.getParameter<int>("phi_bin_window")),
@@ -174,12 +184,13 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps)
       resol_calo_scale_em_(ps.getParameter<double>("resol_calo_scale_em")),
       debug_(ps.getParameter<bool>("debug")),
       eidInputName_(ps.getParameter<std::string>("eid_input_name")),
+      eidInputNameER_(ps.getParameter<std::string>("eid_input_nameER")),
       eidOutputNameEnergy_(ps.getParameter<std::string>("eid_output_name_energy")),
+      eidOutputNameEnergyER_(ps.getParameter<std::string>("eid_output_name_energyER")),
       eidOutputNameId_(ps.getParameter<std::string>("eid_output_name_id")),
       eidMinClusterEnergy_(ps.getParameter<double>("eid_min_cluster_energy")),
       eidNLayers_(ps.getParameter<int>("eid_n_layers")),
-      eidNClusters_(ps.getParameter<int>("eid_n_clusters")),
-      eidSession_(nullptr) {
+      eidNClusters_(ps.getParameter<int>("eid_n_clusters")){
   produces<std::vector<Trackster>>();
   produces<std::vector<TICLCandidate>>();
 
@@ -251,6 +262,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es
   auto resultCandidates = std::make_unique<std::vector<TICLCandidate>>();
 
   tfSession_ = es.getData(tfDnnToken_).getSession();
+  tfSessionER_ = es.getData(tfDnnTokenER_).getSession();
 
   edm::Handle<std::vector<reco::Track>> track_h;
   evt.getByToken(tracks_token_, track_h);
@@ -380,8 +392,9 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es
 }
 
 void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::CaloCluster> &layerClusters,
-                                                    const tensorflow::Session *eidSession,
-                                                    std::vector<Trackster> &tracksters) const {
+                                                              const tensorflow::Session *eidSession,
+                                                              const tensorflow::Session *eidSessionEnergyRegression,
+                                                              std::vector<Trackster> &tracksters) {
   // Energy regression and particle identification strategy:
   //
   // 1. Set default values for regressed energy and particle id for each trackster.
@@ -406,14 +419,13 @@ void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::Calo
   // set default values per trackster, determine if the cluster energy threshold is passed,
   // and store indices of hard tracksters
   std::vector<int> tracksterIndices;
-  for (int i = 0; i < (int)tracksters.size(); i++) {
+  for (int i = 0; i < static_cast<int>(tracksters.size()); i++) {
     // calculate the cluster energy sum (2)
-    // note: after the loop, sumClusterEnergy might be just above the threshold
-    // which is enough to decide whether to run inference for the trackster or
-    // not
+    // note: after the loop, sumClusterEnergy might be just above the threshold which is enough to
+    // decide whether to run inference for the trackster or not
     float sumClusterEnergy = 0.;
     for (const unsigned int &vertex : tracksters[i].vertices()) {
-      sumClusterEnergy += (float)layerClusters[vertex].energy();
+      sumClusterEnergy += static_cast<float>(layerClusters[vertex].energy());
       // there might be many clusters, so try to stop early
       if (sumClusterEnergy >= eidMinClusterEnergy_) {
         // set default values (1)
@@ -426,7 +438,7 @@ void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::Calo
   }
 
   // do nothing when no trackster passes the selection (3)
-  int batchSize = (int)tracksterIndices.size();
+  int batchSize = static_cast<int>(tracksterIndices.size());
   if (batchSize == 0) {
     return;
   }
@@ -435,7 +447,6 @@ void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::Calo
   tensorflow::TensorShape shape({batchSize, eidNLayers_, eidNClusters_, eidNFeatures_});
   tensorflow::Tensor input(tensorflow::DT_FLOAT, shape);
   tensorflow::NamedTensorList inputList = {{eidInputName_, input}};
-  static constexpr int inputDimension = 4;
 
   std::vector<tensorflow::Tensor> outputs;
   std::vector<std::string> outputNames;
@@ -449,11 +460,9 @@ void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::Calo
   // fill input tensor (5)
   for (int i = 0; i < batchSize; i++) {
     const Trackster &trackster = tracksters[tracksterIndices[i]];
-
-    // per layer, we only consider the first eidNClusters_ clusters in terms of
-    // energy, so in order to avoid creating large / nested structures to do
-    // the sorting for an unknown number of total clusters, create a sorted
-    // list of layer cluster indices to keep track of the filled clusters
+    // per layer, we only consider the first eidNClusters_ clusters in terms of energy, so in order
+    // to avoid creating large / nested structures to do the sorting for an unknown number of total
+    // clusters, create a sorted list of layer cluster indices to keep track of the filled clusters
     std::vector<int> clusterIndices(trackster.vertices().size());
     for (int k = 0; k < (int)trackster.vertices().size(); k++) {
       clusterIndices[k] = k;
@@ -472,7 +481,7 @@ void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::Calo
       int j = rhtools_.getLayerWithOffset(cluster.hitsAndFractions()[0].first) - 1;
       if (j < eidNLayers_ && seenClusters[j] < eidNClusters_) {
         // get the pointer to the first feature value for the current batch, layer and cluster
-        float *features = &input.tensor<float, inputDimension>()(i, j, seenClusters[j], 0);
+        float *features = &input.tensor<float, 4>()(i, j, seenClusters[j], 0);
 
         // fill features
         *(features++) = float(cluster.energy() / float(trackster.vertex_multiplicity(k)));
@@ -487,7 +496,7 @@ void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::Calo
     // zero-fill features of empty clusters in each layer (6)
     for (int j = 0; j < eidNLayers_; j++) {
       for (int k = seenClusters[j]; k < eidNClusters_; k++) {
-        float *features = &input.tensor<float, inputDimension>()(i, j, k, 0);
+        float *features = &input.tensor<float, 4>()(i, j, k, 0);
         for (int l = 0; l < eidNFeatures_; l++) {
           *(features++) = 0.f;
         }
@@ -498,20 +507,10 @@ void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::Calo
   // run the inference (7)
   tensorflow::run(const_cast<tensorflow::Session *>(eidSession), inputList, outputNames, &outputs);
 
-  // store regressed energy per trackster (8)
-  if (!eidOutputNameEnergy_.empty()) {
-    // get the pointer to the energy tensor, dimension is batch x 1
-    float *energy = outputs[0].flat<float>().data();
-
-    for (const int &i : tracksterIndices) {
-      tracksters[i].setRegressedEnergy(*(energy++));
-    }
-  }
-
   // store id probabilities per trackster (8)
   if (!eidOutputNameId_.empty()) {
     // get the pointer to the id probability tensor, dimension is batch x id_probabilities.size()
-    int probsIdx = !eidOutputNameEnergy_.empty();
+    int probsIdx = eidOutputNameEnergy_.empty() ? 0 : 1;
     float *probs = outputs[probsIdx].flat<float>().data();
 
     for (const int &i : tracksterIndices) {
@@ -519,6 +518,69 @@ void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::Calo
       probs += tracksters[i].id_probabilities().size();
     }
   }
+  // NEW ENERGY REGRESSION!
+  std::vector<int> tracksterIndicesHadronic;
+  auto index_trackster = 0;
+  for (auto &t : tracksters) {
+    t.setRegressedEnergy(t.raw_energy());
+    if ((t.id_probability(ticl::Trackster::ParticleType::photon) +
+         t.id_probability(ticl::Trackster::ParticleType::electron)) >= 0.) {
+      tracksterIndicesHadronic.push_back(index_trackster);
+    }
+    index_trackster++;
+  }
+  batchSize = static_cast<int>(tracksterIndicesHadronic.size());
+  if (batchSize == 0) {
+    return;
+  }
+  // create input and output tensors for energy regression
+  tensorflow::TensorShape shapeER({batchSize, eidNFeaturesER_});
+  tensorflow::Tensor inputEnergyRegression(tensorflow::DT_FLOAT, shapeER);
+  tensorflow::NamedTensorList inputListEnergyRegression = {{eidInputNameER_, inputEnergyRegression}};
+
+  std::vector<tensorflow::Tensor> outputsEnergyRegression;
+  std::vector<std::string> outputNamesEnergyRegression;
+
+  outputNamesEnergyRegression.push_back(eidOutputNameEnergyER_);
+
+  // fill input tensor (5)
+  for (int i = 0; i < batchSize; i++) {
+    const Trackster &trackster = tracksters[tracksterIndicesHadronic[i]]; 
+    float en_total = 0;
+    float *featuresER = &inputEnergyRegression.tensor<float, 2>()(i, 0);
+    std::array<float, eidNFeaturesER_> inputsArrayER;
+    for (size_t j = 0; j < inputsArrayER.size(); j++) {
+      inputsArrayER[j] = 0.;
+    }
+
+    for (int k = 0; k < (int)trackster.vertices().size(); k++) {
+      const reco::CaloCluster &cluster = layerClusters[trackster.vertices(k)];
+      auto lc_seed = cluster.seed();
+      auto index = returnIndex(lc_seed, rhtools_);
+        if(index >= 0){
+          inputsArrayER[index] += cluster.energy() / trackster.raw_energy(); //energy fraction
+        }
+      en_total += cluster.energy();
+    }
+    
+    inputsArrayER[ticl::LayerType::EnumSize] = trackster.barycenter().eta(); //eta
+    inputsArrayER[ticl::LayerType::EnumSize+1] = trackster.raw_energy(); //total trackster energy
+    for(size_t i_p = 0; i_p < trackster.id_probabilities().size(); i_p++){
+      inputsArrayER[ticl::LayerType::EnumSize+2+i_p] = trackster.id_probabilities(i_p); //probabilities
+    }
+    for (size_t j = 0; j < inputsArrayER.size(); j++) {
+      *(featuresER++) = inputsArrayER[j];
+    }
+  }
+
+  tensorflow::run(const_cast<tensorflow::Session *>(eidSessionEnergyRegression),
+                  inputListEnergyRegression,
+                  outputNamesEnergyRegression,
+                  &outputsEnergyRegression);
+  float *energy = outputsEnergyRegression[0].flat<float>().data();
+  for (const int &i : tracksterIndicesHadronic) {
+    tracksters[i].setRegressedEnergy(*(energy++));
+  }
 }
 
 void TrackstersMergeProducer::assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates) const {
@@ -606,8 +668,11 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d
   desc.add<double>("resol_calo_scale_em", 0.15);
   desc.add<bool>("debug", false);
   desc.add<std::string>("tfDnnLabel", "tracksterSelectionTf");
+  desc.add<std::string>("tfDnnLabelER", "tracksterSelectionTfER");
   desc.add<std::string>("eid_input_name", "input");
+  desc.add<std::string>("eid_input_nameER", "input");
   desc.add<std::string>("eid_output_name_energy", "output/regressed_energy");
+  desc.add<std::string>("eid_output_name_energyER", "output/er_coefficients");
   desc.add<std::string>("eid_output_name_id", "output/id_probabilities");
   desc.add<double>("eid_min_cluster_energy", 1.);
   desc.add<int>("eid_n_layers", 50);
diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc
index 9ab49a19662c3..8b69ca3d53dc1 100644
--- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc
+++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc
@@ -47,8 +47,11 @@ class TrackstersProducer : public edm::stream::EDProducer<> {
   std::string detector_;
   bool doNose_;
   const std::string tfDnnLabel_;
+  const std::string tfDnnLabelER_;
   const edm::ESGetToken<TfGraphDefWrapper, TfGraphRecord> tfDnnToken_;
+  const edm::ESGetToken<TfGraphDefWrapper, TfGraphRecord> tfDnnTokenER_;
   const tensorflow::Session* tfSession_;
+  const tensorflow::Session* tfSessionER_;
   std::unique_ptr<PatternRecognitionAlgoBaseT<TICLLayerTiles>> myAlgo_;
   std::unique_ptr<PatternRecognitionAlgoBaseT<TICLLayerTilesHFNose>> myAlgoHFNose_;
 
@@ -68,8 +71,11 @@ TrackstersProducer::TrackstersProducer(const edm::ParameterSet& ps)
     : detector_(ps.getParameter<std::string>("detector")),
       doNose_(detector_ == "HFNose"),
       tfDnnLabel_(ps.getParameter<std::string>("tfDnnLabel")),
+      tfDnnLabelER_(ps.getParameter<std::string>("tfDnnLabelER")),
       tfDnnToken_(esConsumes(edm::ESInputTag("", tfDnnLabel_))),
+      tfDnnTokenER_(esConsumes(edm::ESInputTag("", tfDnnLabelER_))),
       tfSession_(nullptr),
+      tfSessionER_(nullptr),
       clusters_token_(consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("layer_clusters"))),
       filtered_layerclusters_mask_token_(consumes<std::vector<float>>(ps.getParameter<edm::InputTag>("filtered_mask"))),
       original_layerclusters_mask_token_(consumes<std::vector<float>>(ps.getParameter<edm::InputTag>("original_mask"))),
@@ -120,6 +126,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri
   desc.add<std::string>("patternRecognitionBy", "CA");
   desc.add<std::string>("itername", "unknown");
   desc.add<std::string>("tfDnnLabel", "tracksterSelectionTf");
+  desc.add<std::string>("tfDnnLabelER", "tracksterSelectionTfER");
 
   // CA Plugin
   edm::ParameterSetDescription pluginDesc;
@@ -150,6 +157,7 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
   const auto& seeding_regions = evt.get(seeding_regions_token_);
 
   tfSession_ = es.getData(tfDnnToken_).getSession();
+  tfSessionER_ = es.getData(tfDnnTokenER_).getSession();
 
   std::unordered_map<int, std::vector<int>> seedToTrackstersAssociation;
   // if it's regional iteration and there are seeding regions
@@ -169,14 +177,22 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
                                                                                          layerClustersTimes,
                                                                                          layer_clusters_hfnose_tiles,
                                                                                          seeding_regions,
-                                                                                         tfSession_);
+                                                                                         tfSession_,
+                                                                                         tfSessionER_);
 
     myAlgoHFNose_->makeTracksters(inputHFNose, *result, seedToTrackstersAssociation);
 
   } else {
     const auto& layer_clusters_tiles = evt.get(layer_clusters_tiles_token_);
-    const typename PatternRecognitionAlgoBaseT<TICLLayerTiles>::Inputs input(
-        evt, es, layerClusters, inputClusterMask, layerClustersTimes, layer_clusters_tiles, seeding_regions, tfSession_);
+    const typename PatternRecognitionAlgoBaseT<TICLLayerTiles>::Inputs input(evt,
+                                                                             es,
+                                                                             layerClusters,
+                                                                             inputClusterMask,
+                                                                             layerClustersTimes,
+                                                                             layer_clusters_tiles,
+                                                                             seeding_regions,
+                                                                             tfSession_,
+                                                                             tfSessionER_);
 
     myAlgo_->makeTracksters(input, *result, seedToTrackstersAssociation);
   }
diff --git a/RecoHGCal/TICL/python/customiseTICLFromReco.py b/RecoHGCal/TICL/python/customiseTICLFromReco.py
new file mode 100644
index 0000000000000..9695b53bfc8bf
--- /dev/null
+++ b/RecoHGCal/TICL/python/customiseTICLFromReco.py
@@ -0,0 +1,67 @@
+# Reconstruction
+from RecoHGCal.TICL.iterativeTICL_cff import *
+from RecoLocalCalo.HGCalRecProducers.hgcalLayerClusters_cff import hgcalLayerClusters
+# Validation
+from Validation.HGCalValidation.HGCalValidator_cfi import *
+from RecoLocalCalo.HGCalRecProducers.hgcalRecHitMapProducer_cfi import hgcalRecHitMapProducer
+
+# Load DNN ESSource
+from RecoTracker.IterativeTracking.iterativeTk_cff import trackdnn_source
+
+# Automatic addition of the customisation function from RecoHGCal.Configuration.RecoHGCal_EventContent_cff
+from RecoHGCal.Configuration.RecoHGCal_EventContent_cff import customiseHGCalOnlyEventContent
+
+
+
+def customiseTICLFromReco(process):
+
+# TensorFlow ESSource
+    process.TFileService = cms.Service("TFileService",
+                                       # fileName = cms.string('/afs/cern.ch/work/w/wredjeb/public/EnergyRegression/TFInterface/Rebase/CMSSW_12_4_0_pre2/src/38693.0_CloseByParticleGun+2026D86+CloseByParticle_Photon_ERZRanges_GenSimHLBeamSpotHGCALCloseBy+DigiTrigger+RecoGlobal+HARVESTGlobal/ticlTracksterCLUE3DKaonsAllProb.root')
+					fileName = cms.string('/afs/cern.ch/work/w/wredjeb/public/EnergyRegression/TFInterface/Rebase/CMSSW_12_4_0_pre2/src/38694.203_CloseByPGun_CE_E_Front_300um+2026D86_ticl_v4+CE_E_Front_300um_GenSimHLBeamSpotHGCALCloseBy+DigiTrigger+RecoGlobal+HARVESTGlobal/ticlTrackstersCLUE3DElectronAllProb.root')
+                                    )
+    process.TFESSource = cms.Task(process.trackdnn_source)
+
+    process.TICL = cms.Path(process.hgcalLayerClusters,
+                            process.TFESSource,
+                            process.ticlLayerTileTask,
+                            process.ticlIterationsTask,
+                            process.ticlTracksterMergeTask)
+
+    process.trackstersNtuplerCLUE3D = cms.EDAnalyzer('TracksterNtupler',
+        tracksters = cms.InputTag('ticlTrackstersCLUE3DHigh'),
+        caloParticles = cms.InputTag('mix', 'MergedCaloTruth'),
+        layerClusters=cms.InputTag('hgcalLayerClusters'),
+        outfilePath = cms.untracked.string('/afs/cern.ch/work/w/wredjeb/public/EnergyRegression/TFInterface/Rebase/CMSSW_12_4_0_pre2/src/38693.0_CloseByParticleGun+2026D86+CloseByParticle_Photon_ERZRanges_GenSimHLBeamSpotHGCALCloseBy+DigiTrigger+RecoGlobal+HARVESTGlobal/ticlTracksterCLUE3D_newRegression.root')
+    )
+
+
+    process.TFileService = cms.Service("TFileService",
+            fileName = cms.string("ticlTracksterCLUE3D_newRegressionFS.root")
+    )
+
+# Validation
+    process.TICL_ValidationProducers = cms.Task(process.hgcalRecHitMapProducer,
+                                                process.lcAssocByEnergyScoreProducer,
+                                                process.layerClusterCaloParticleAssociationProducer,
+                                                process.scAssocByEnergyScoreProducer,
+                                                process.layerClusterSimClusterAssociationProducer,
+                                               )
+    process.TICL_Validator = cms.Task(process.hgcalValidator)
+    process.TICL_Validation = cms.Path(process.TICL_ValidationProducers,
+                                       process.TICL_Validator
+                                      )
+# Path and EndPath definitions
+    process.FEVTDEBUGHLToutput_step = cms.EndPath(process.FEVTDEBUGHLToutput + process.trackstersNtuplerCLUE3D)
+    process.DQMoutput_step = cms.EndPath(process.DQMoutput)
+
+# Schedule definition
+    process.schedule = cms.Schedule(process.TICL,
+                                    # process.TICL_Validation,
+                                    process.FEVTDEBUGHLToutput_step,
+                                    # process.DQMoutput_step,
+                                    )
+#call to customisation function customiseHGCalOnlyEventContent imported from RecoHGCal.Configuration.RecoHGCal_EventContent_cff
+    process = customiseHGCalOnlyEventContent(process)
+
+    return process
diff --git a/RecoHGCal/TICL/python/tracksterSelectionTf_cfi.py b/RecoHGCal/TICL/python/tracksterSelectionTf_cfi.py
index 8488872a93710..3ea0c14030c04 100644
--- a/RecoHGCal/TICL/python/tracksterSelectionTf_cfi.py
+++ b/RecoHGCal/TICL/python/tracksterSelectionTf_cfi.py
@@ -3,3 +3,8 @@
     ComponentName = "tracksterSelectionTf",
     FileName = "RecoHGCal/TICL/data/tf_models/energy_id_v0.pb"
 )
+
+tracksterSelectionTfER = _tfGraphDefProducer.clone(
+    ComponentName = "tracksterSelectionTfER",
+    FileName = "RecoHGCal/TICL/data/tf_models/energy_regression_tracksters.pb"
+)