diff --git a/DataFormats/L1THGCal/interface/HGCalClusterT.h b/DataFormats/L1THGCal/interface/HGCalClusterT.h
index 2dc644d6fcc62..93d949b527c87 100644
--- a/DataFormats/L1THGCal/interface/HGCalClusterT.h
+++ b/DataFormats/L1THGCal/interface/HGCalClusterT.h
@@ -152,6 +152,8 @@ namespace l1t
       float sigmaEtaEtaTot() const { return sigmaEtaEtaTot_; }
       float sigmaPhiPhiTot() const { return sigmaPhiPhiTot_; }
       float sigmaZZ() const { return sigmaZZ_; }
+      float sigmaRRTot() const { return sigmaRRTot_; }
+      float sigmaRRMax() const { return sigmaRRMax_; }
 
       void set_showerLength(int showerLength) { showerLength_ = showerLength;}
       void set_firstLayer(int firstLayer) { firstLayer_ = firstLayer;}
@@ -160,6 +162,8 @@ namespace l1t
       void set_sigmaEtaEtaTot(float sigmaEtaEtaTot) { sigmaEtaEtaTot_ = sigmaEtaEtaTot;}
       void set_sigmaPhiPhiMax(float sigmaPhiPhiMax) { sigmaPhiPhiMax_ = sigmaPhiPhiMax;}
       void set_sigmaPhiPhiTot(float sigmaPhiPhiTot) { sigmaPhiPhiTot_ = sigmaPhiPhiTot;}
+      void set_sigmaRRMax(float sigmaRRMax) { sigmaRRMax_ = sigmaRRMax;}
+      void set_sigmaRRTot(float sigmaRRTot) { sigmaRRTot_ = sigmaRRTot;}
       void set_sigmaZZ(float sigmaZZ) { sigmaZZ_ = sigmaZZ;}
       
       /* operators */
@@ -186,8 +190,10 @@ namespace l1t
       float eMax_;
       float sigmaEtaEtaMax_;
       float sigmaPhiPhiMax_;
+      float sigmaRRMax_;
       float sigmaEtaEtaTot_;
       float sigmaPhiPhiTot_;
+      float sigmaRRTot_;
       float sigmaZZ_;
 
       ClusterShapes shapes_;
diff --git a/DataFormats/L1THGCal/src/classes_def.xml b/DataFormats/L1THGCal/src/classes_def.xml
index 2b823399651ca..8933b357f2327 100644
--- a/DataFormats/L1THGCal/src/classes_def.xml
+++ b/DataFormats/L1THGCal/src/classes_def.xml
@@ -20,7 +20,8 @@
   <class name="edm::Wrapper<l1t::HGCalTowerBxCollection>"/>
 
   <class name="l1t::HGCalClusterT<l1t::HGCalTriggerCell>" />
-  <class name="l1t::HGCalCluster" ClassVersion="11">
+  <class name="l1t::HGCalCluster" ClassVersion="12">
+  <version ClassVersion="12" checksum="623703096"/>
   <version ClassVersion="11" checksum="354623225"/>
   <version ClassVersion="10" checksum="607544984"/>
   </class>
@@ -29,7 +30,8 @@
   <class name="edm::Wrapper<l1t::HGCalClusterBxCollection>"/>
 
   <class name="l1t::HGCalClusterT<l1t::HGCalCluster>" />
-  <class name="l1t::HGCalMulticluster" ClassVersion="11">
+  <class name="l1t::HGCalMulticluster" ClassVersion="12">
+  <version ClassVersion="12" checksum="4221677522"/>
   <version ClassVersion="11" checksum="1508179045"/>
   <version ClassVersion="10" checksum="1878482802"/>
   </class>
diff --git a/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h
index e017769d53173..5b2129507505b 100644
--- a/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h
+++ b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h
@@ -3,10 +3,12 @@
 #include <vector>
 #include <cmath>
 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
 
 class HGCalShowerShape{
 
     public:
+    typedef math::XYZTLorentzVector LorentzVector;
         
     HGCalShowerShape(){}
 
@@ -28,11 +30,15 @@ class HGCalShowerShape{
     float sigmaPhiPhiTot(const l1t::HGCalCluster& c2d) const;       
     float sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const;   
 
+    float sigmaRRTot(const l1t::HGCalMulticluster& c3d) const;
+    float sigmaRRTot(const l1t::HGCalCluster& c2d) const;       
+    float sigmaRRMax(const l1t::HGCalMulticluster& c3d) const;  
+
     private: 
     
-    float sigmaEtaEta(const std::vector<float>& energy, const std::vector<float>& eta) const;
-    float sigmaPhiPhi(const std::vector<float>& energy, const std::vector<float>& phi) const;   
-    float sigmaZZ(const std::vector<float>& energy, const std::vector<float>& z) const;
+    float meanX(const std::vector<pair<float,float> >& energy_X_tc) const;
+    float sigmaXX(const std::vector<pair<float,float> >& energy_X_tc, const float X_cluster) const;
+    float sigmaPhiPhi(const std::vector<pair<float,float> >& energy_phi_tc, const float phi_cluster) const;
 
     static const int kLayersEE_=28;
     static const int kLayersFH_=12;
diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc
index 7972a4c609a2a..a0da63fbb63a3 100644
--- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc
+++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc
@@ -33,6 +33,8 @@ class HGCalTriggerNtupleHGCMulticlusters : public HGCalTriggerNtupleBase
     std::vector<float> cl3d_spptot_;
     std::vector<float> cl3d_sppmax_;
     std::vector<float> cl3d_szz_;
+    std::vector<float> cl3d_srrtot_;
+    std::vector<float> cl3d_srrmax_;
     std::vector<float> cl3d_emaxe_;
     std::vector<std::vector<unsigned>> cl3d_clusters_;   
 };
@@ -66,6 +68,8 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&&
   tree.Branch("cl3d_spptot", &cl3d_spptot_);
   tree.Branch("cl3d_sppmax", &cl3d_sppmax_);
   tree.Branch("cl3d_szz", &cl3d_szz_);
+  tree.Branch("cl3d_srrtot", &cl3d_srrtot_);
+  tree.Branch("cl3d_srrmax", &cl3d_srrmax_);
   tree.Branch("cl3d_emaxe", &cl3d_emaxe_);  
   tree.Branch("cl3d_clusters", &cl3d_clusters_);
 
@@ -102,6 +106,8 @@ fill(const edm::Event& e, const edm::EventSetup& es)
     cl3d_spptot_.emplace_back(cl3d_itr->sigmaPhiPhiTot());
     cl3d_sppmax_.emplace_back(cl3d_itr->sigmaPhiPhiMax());
     cl3d_szz_.emplace_back(cl3d_itr->sigmaZZ());
+    cl3d_srrtot_.emplace_back(cl3d_itr->sigmaRRTot());
+    cl3d_srrmax_.emplace_back(cl3d_itr->sigmaRRMax());
     cl3d_emaxe_.emplace_back(cl3d_itr->eMax()/cl3d_itr->energy());
 
     // Retrieve indices of trigger cells inside cluster
@@ -130,6 +136,8 @@ clear()
   cl3d_spptot_.clear();
   cl3d_sppmax_.clear();
   cl3d_szz_.clear();
+  cl3d_srrtot_.clear();
+  cl3d_srrmax_.clear();
   cl3d_emaxe_.clear();
   cl3d_clusters_.clear();
   
diff --git a/L1Trigger/L1THGCal/src/be_algorithms/HGCalMulticlusteringImpl.cc b/L1Trigger/L1THGCal/src/be_algorithms/HGCalMulticlusteringImpl.cc
index 1389ff1aa11db..3f7e24ad38421 100644
--- a/L1Trigger/L1THGCal/src/be_algorithms/HGCalMulticlusteringImpl.cc
+++ b/L1Trigger/L1THGCal/src/be_algorithms/HGCalMulticlusteringImpl.cc
@@ -126,6 +126,8 @@ void HGCalMulticlusteringImpl::clusterizeDR( const edm::PtrVector<l1t::HGCalClus
             multiclustersTmp.at(i).set_sigmaPhiPhiTot(shape_.sigmaPhiPhiTot(multiclustersTmp.at(i)));
             multiclustersTmp.at(i).set_sigmaPhiPhiMax(shape_.sigmaPhiPhiMax(multiclustersTmp.at(i)));
             multiclustersTmp.at(i).set_sigmaZZ(shape_.sigmaZZ(multiclustersTmp.at(i)));
+            multiclustersTmp.at(i).set_sigmaRRTot(shape_.sigmaRRTot(multiclustersTmp.at(i)));
+            multiclustersTmp.at(i).set_sigmaRRMax(shape_.sigmaRRMax(multiclustersTmp.at(i)));
             multiclustersTmp.at(i).set_eMax(shape_.eMax(multiclustersTmp.at(i)));
 
             multiclusters.push_back( 0, multiclustersTmp.at(i));  
@@ -214,6 +216,8 @@ void HGCalMulticlusteringImpl::clusterizeDBSCAN( const edm::PtrVector<l1t::HGCal
       multiclustersTmp.at(i).set_sigmaPhiPhiTot(shape_.sigmaPhiPhiTot(multiclustersTmp.at(i)));
       multiclustersTmp.at(i).set_sigmaPhiPhiMax(shape_.sigmaPhiPhiMax(multiclustersTmp.at(i)));
       multiclustersTmp.at(i).set_sigmaZZ(shape_.sigmaZZ(multiclustersTmp.at(i)));
+      multiclustersTmp.at(i).set_sigmaRRTot(shape_.sigmaRRTot(multiclustersTmp.at(i)));
+      multiclustersTmp.at(i).set_sigmaRRMax(shape_.sigmaRRMax(multiclustersTmp.at(i)));
       multiclustersTmp.at(i).set_eMax(shape_.eMax(multiclustersTmp.at(i)));
       
       multiclusters.push_back( 0, multiclustersTmp.at(i));  
diff --git a/L1Trigger/L1THGCal/src/be_algorithms/HGCalShowerShape.cc b/L1Trigger/L1THGCal/src/be_algorithms/HGCalShowerShape.cc
index 4a637f080d7fc..ea8af36729224 100644
--- a/L1Trigger/L1THGCal/src/be_algorithms/HGCalShowerShape.cc
+++ b/L1Trigger/L1THGCal/src/be_algorithms/HGCalShowerShape.cc
@@ -4,6 +4,7 @@
 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
 #include "DataFormats/L1THGCal/interface/HGCalCluster.h"
 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
 
 #include <iostream>
 #include <sstream>
@@ -23,83 +24,68 @@ int HGCalShowerShape::HGC_layer(const uint32_t subdet, const uint32_t layer)
 
 }
 
+//Compute energy-weighted mean of any variable X in the cluster
 
+float HGCalShowerShape::meanX(const std::vector<pair<float,float> >& energy_X_tc) const {
 
+  float Etot = 0;
+  float X_sum = 0;
 
+  for(const auto& energy_X : energy_X_tc){
 
-float HGCalShowerShape::sigmaEtaEta(const std::vector<float>& energy, const std::vector<float>& eta) const {
+    X_sum += energy_X.first*energy_X.second;
+    Etot += energy_X.first;
 
-  int ntc=energy.size();
-  //compute weighted eta mean
-  float eta_sum=0;
-  float w_sum=0; //here weight is energy
-  for(int i=0;i<ntc;i++){
-    eta_sum+=energy.at(i)*eta.at(i);
-    w_sum+=energy.at(i);
   }
-  float eta_mean=0;
-  if (w_sum!=0) eta_mean=eta_sum/w_sum;
-  //compute weighted eta RMS
-  float deltaeta2_sum=0;
-  for(int i=0;i<ntc;i++) deltaeta2_sum+=energy.at(i)*pow((eta.at(i)-eta_mean),2);
-  float eta_RMS=0;
-  if (w_sum!=0) eta_RMS=deltaeta2_sum/w_sum;
-  float See=sqrt(eta_RMS);
-  return See;
-    
+
+  float X_mean = 0;
+  if(Etot>0) X_mean = X_sum/Etot;
+  return X_mean;
+
 }
 
+//Compute energy-weighted RMS of any variable X in the cluster
 
-float HGCalShowerShape::sigmaPhiPhi(const std::vector<float>& energy, const std::vector<float>& phi) const {
+float HGCalShowerShape::sigmaXX(const std::vector<pair<float,float> >& energy_X_tc, const float X_cluster) const {
+
+  float Etot = 0;
+  float deltaX2_sum = 0;
+
+  for(const auto& energy_X : energy_X_tc){
+
+    deltaX2_sum += energy_X.first * pow(energy_X.second - X_cluster,2);
+    Etot += energy_X.first;
 
-  int ntc=energy.size();
-  //compute weighted phi mean
-  float phi_sum=0;
-  float w_sum=0; //here weight is energy
-  for(int i=0;i<ntc;i++){
-    if(phi.at(i)>0) phi_sum+=energy.at(i)*phi.at(i);
-    else phi_sum+=energy.at(i)*(phi.at(i)+2*TMath::Pi());
-    w_sum+=energy.at(i);
   }
-  float phi_mean=0;
-  if (w_sum!=0) phi_mean=phi_sum/w_sum;
-  if(phi_mean>TMath::Pi()) phi_mean-=2*TMath::Pi();
-  //compute weighted eta RMS
-  float deltaphi2_sum=0;
-  for(int i=0;i<ntc;i++){
-    float deltaPhi=std::abs(phi.at(i)-phi_mean);
-    if (deltaPhi>TMath::Pi()) deltaPhi=2*TMath::Pi()-deltaPhi;
-    deltaphi2_sum+=energy.at(i)*pow(deltaPhi,2);
-  }        
-  float phi_RMS=0;
-  if (w_sum!=0) phi_RMS=deltaphi2_sum/w_sum;
-  float Spp=sqrt(phi_RMS);
-  return Spp;
-  
+
+  float X_MSE = 0;
+  if (Etot>0) X_MSE=deltaX2_sum/Etot;
+  float X_RMS=sqrt(X_MSE);
+  return X_RMS;
+
 }
 
 
+//Compute energy-weighted RMS of any variable X in the cluster
+//Extra care needed because of deltaPhi
+
+float HGCalShowerShape::sigmaPhiPhi(const std::vector<pair<float,float> >& energy_phi_tc, const float phi_cluster) const {
 
-float HGCalShowerShape::sigmaZZ(const std::vector<float>& energy, const std::vector<float>& z) const {
+  float Etot = 0;
+  float deltaphi2_sum = 0;
+
+  for(const auto& energy_phi : energy_phi_tc){
+
+    deltaphi2_sum += energy_phi.first * pow(deltaPhi(energy_phi.second,phi_cluster),2);
+    Etot += energy_phi.first;
 
-  int ntc=energy.size();
-  //compute weighted eta mean
-  float z_sum=0;
-  float w_sum=0; //here weight is energy
-  for(int i=0;i<ntc;i++){
-    z_sum+=energy.at(i)*z.at(i);
-    w_sum+=energy.at(i);
   }
-  float z_mean=0;
-  if (w_sum!=0) z_mean=z_sum/w_sum;
-  //compute weighted eta RMS
-  float deltaz2_sum=0;
-  for(int i=0;i<ntc;i++) deltaz2_sum+=energy.at(i)*pow((z.at(i)-z_mean),2);
-  float z_RMS=0;
-  if (w_sum!=0) z_RMS=deltaz2_sum/w_sum;
-  float Szz=sqrt(z_RMS);
-  return Szz;
-    
+
+  float phi_MSE = 0;
+  if (Etot>0) phi_MSE=deltaphi2_sum/Etot;
+  float Spp=sqrt(phi_MSE);
+  return Spp;
+
 }
 
 
@@ -148,8 +134,7 @@ float HGCalShowerShape::sigmaEtaEtaTot(const l1t::HGCalMulticluster& c3d) const
 
   const edm::PtrVector<l1t::HGCalCluster>& clustersPtrs = c3d.constituents();
   
-  std::vector<float> tc_energy ; 
-  std::vector<float> tc_eta ;
+  std::vector<std::pair<float,float> > tc_energy_eta ; 
   
   for(const auto& clu : clustersPtrs){
     
@@ -157,17 +142,13 @@ float HGCalShowerShape::sigmaEtaEtaTot(const l1t::HGCalMulticluster& c3d) const
     
     for(const auto& tc : triggerCells){
       
-      tc_energy.emplace_back(tc->energy());
-      tc_eta.emplace_back(tc->eta());
+      tc_energy_eta.emplace_back( std::make_pair(tc->energy(),tc->eta()) );
       
     }
     
   }
   
-  float SeeTot = sigmaEtaEta(tc_energy,tc_eta);
-  
-  tc_energy.clear();
-  tc_eta.clear();
+  float SeeTot = sigmaXX(tc_energy_eta,c3d.eta());
   
   return SeeTot;
   
@@ -180,8 +161,7 @@ float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const
 
   const edm::PtrVector<l1t::HGCalCluster>& clustersPtrs = c3d.constituents();
   
-  std::vector<float> tc_energy ; 
-  std::vector<float> tc_phi ;
+  std::vector<std::pair<float,float> > tc_energy_phi ; 
   
   for(const auto& clu : clustersPtrs){
     
@@ -189,16 +169,13 @@ float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const
     
     for(const auto& tc : triggerCells){
       
-      tc_energy.emplace_back(tc->energy());
-      tc_phi.emplace_back(tc->phi());
+      tc_energy_phi.emplace_back( std::make_pair(tc->energy(),tc->phi()) );
 
     }
+
   }
   
-  float SppTot = sigmaPhiPhi(tc_energy,tc_phi);
-  
-  tc_energy.clear();
-  tc_phi.clear();
+  float SppTot = sigmaPhiPhi(tc_energy_phi,c3d.phi());
   
   return SppTot;
   
@@ -206,10 +183,42 @@ float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const
 
 
 
+
+
+float HGCalShowerShape::sigmaRRTot(const l1t::HGCalMulticluster& c3d) const {
+
+  const edm::PtrVector<l1t::HGCalCluster>& clustersPtrs = c3d.constituents();
+  
+  std::vector<std::pair<float,float> > tc_energy_r ; 
+  
+  for(const auto& clu : clustersPtrs){
+    
+    const edm::PtrVector<l1t::HGCalTriggerCell>& triggerCells = clu->constituents();
+    
+    for(const auto& tc : triggerCells){
+           
+      float r = std::sqrt( pow(tc->position().x(),2) + pow(tc->position().y(),2) )/std::abs(tc->position().z());
+      tc_energy_r.emplace_back( std::make_pair(tc->energy(),r) );
+
+    }
+
+  }
+
+  float r_mean = meanX(tc_energy_r);
+  float Szz = sigmaXX(tc_energy_r,r_mean);
+  
+  return Szz;
+  
+}
+
+
+
+
+
 float HGCalShowerShape::sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const {
 
-  std::unordered_map<int, vector<float> > tc_layer_energy;
-  std::unordered_map<int, vector<float> > tc_layer_eta;
+  std::unordered_map<int, std::vector<std::pair<float,float> > > tc_layer_energy_eta;
+  std::unordered_map<int, LorentzVector> layer_LV;
 
   const edm::PtrVector<l1t::HGCalCluster>& clustersPtrs = c3d.constituents();
 
@@ -217,28 +226,26 @@ float HGCalShowerShape::sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const
     
     int layer = HGC_layer(clu->subdetId(),clu->layer());    
     
-    const edm::PtrVector<l1t::HGCalTriggerCell>& triggerCells = clu->constituents();
+    layer_LV[layer] += clu->p4();
 
-    auto layer_energy_itr = tc_layer_energy.emplace(layer, 0).first;
-    auto layer_eta_itr = tc_layer_eta.emplace(layer, 0).first;
+    const edm::PtrVector<l1t::HGCalTriggerCell>& triggerCells = clu->constituents();
 
     for(const auto& tc : triggerCells){
 
-      layer_energy_itr->second.emplace_back(tc->energy());
-      layer_eta_itr->second.emplace_back(tc->eta());
+      tc_layer_energy_eta[layer].emplace_back( std::make_pair(tc->energy(),tc->eta()) );
 
     }
 
   }
 
+
   float SigmaEtaEtaMax=0;
 
-  for(auto& tc_iter : tc_layer_energy){
+  for(auto& tc_iter : tc_layer_energy_eta){
     
-    const std::vector<float>& energy_layer = tc_iter.second;
-    const std::vector<float>& eta_layer= tc_layer_eta[tc_iter.first];
-
-    float SigmaEtaEtaLayer = sigmaEtaEta(energy_layer,eta_layer); 
+    const std::vector<std::pair<float, float> >& energy_eta_layer = tc_iter.second;
+    const LorentzVector& LV_layer = layer_LV[tc_iter.first];
+    float SigmaEtaEtaLayer = sigmaXX(energy_eta_layer,LV_layer.eta()); //RMS wrt layer eta, not wrt c3d eta  
     if(SigmaEtaEtaLayer > SigmaEtaEtaMax) SigmaEtaEtaMax = SigmaEtaEtaLayer;
 
   }
@@ -252,11 +259,10 @@ float HGCalShowerShape::sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const
 
 
 
-
 float HGCalShowerShape::sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const {
 
-  std::unordered_map<int, vector<float> > tc_layer_energy;
-  std::unordered_map<int, vector<float> > tc_layer_phi;
+  std::unordered_map<int, std::vector<std::pair<float,float> > > tc_layer_energy_phi;
+  std::unordered_map<int, LorentzVector> layer_LV;
 
   const edm::PtrVector<l1t::HGCalCluster>& clustersPtrs = c3d.constituents();
 
@@ -264,32 +270,31 @@ float HGCalShowerShape::sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const
     
     int layer = HGC_layer(clu->subdetId(),clu->layer());    
     
+    layer_LV[layer] += clu->p4();
+
     const edm::PtrVector<l1t::HGCalTriggerCell>& triggerCells = clu->constituents();
-    
-    auto layer_energy_itr = tc_layer_energy.emplace(layer, 0).first;
-    auto layer_phi_itr = tc_layer_phi.emplace(layer, 0).first;
-    
+
     for(const auto& tc : triggerCells){
-      
-      layer_energy_itr->second.emplace_back(tc->energy());
-      layer_phi_itr->second.emplace_back(tc->phi());
-      
+
+      tc_layer_energy_phi[layer].emplace_back( std::make_pair(tc->energy(),tc->phi()) );
+
     }
-    
+
   }
 
+
   float SigmaPhiPhiMax=0;
 
-  for(auto& tc_iter : tc_layer_energy){
+  for(auto& tc_iter : tc_layer_energy_phi){
     
-    const std::vector<float>& energy_layer = tc_iter.second;
-    const std::vector<float>& phi_layer= tc_layer_phi[tc_iter.first];
-    
-    float SigmaPhiPhiLayer = sigmaPhiPhi(energy_layer,phi_layer); 
+    const std::vector<std::pair<float, float> >& energy_phi_layer = tc_iter.second;
+    const LorentzVector& LV_layer = layer_LV[tc_iter.first];
+    float SigmaPhiPhiLayer = sigmaXX(energy_phi_layer,LV_layer.phi()); //RMS wrt layer phi, not wrt c3d phi
     if(SigmaPhiPhiLayer > SigmaPhiPhiMax) SigmaPhiPhiMax = SigmaPhiPhiLayer;
 
   }
 
+
   return SigmaPhiPhiMax;
 
 
@@ -297,6 +302,47 @@ float HGCalShowerShape::sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const
 
 
 
+
+float HGCalShowerShape::sigmaRRMax(const l1t::HGCalMulticluster& c3d) const {
+
+  std::unordered_map<int, std::vector<std::pair<float,float> > > tc_layer_energy_r;
+
+  const edm::PtrVector<l1t::HGCalCluster>& clustersPtrs = c3d.constituents();
+
+  for(const auto& clu : clustersPtrs){
+    
+    int layer = HGC_layer(clu->subdetId(),clu->layer());        
+
+    const edm::PtrVector<l1t::HGCalTriggerCell>& triggerCells = clu->constituents();
+
+    for(const auto& tc : triggerCells){
+
+      float r = std::sqrt( pow(tc->position().x(),2) + pow(tc->position().y(),2) )/std::abs(tc->position().z());
+      tc_layer_energy_r[layer].emplace_back( std::make_pair(tc->energy(),r) );
+
+    }
+
+  }
+
+
+  float SigmaRRMax=0;
+
+  for(auto& tc_iter : tc_layer_energy_r){
+    
+    const std::vector<std::pair<float, float> >& energy_r_layer = tc_iter.second;
+    float r_mean_layer = meanX(energy_r_layer);
+    float SigmaRRLayer = sigmaXX(energy_r_layer,r_mean_layer);
+    if(SigmaRRLayer > SigmaRRMax) SigmaRRMax = SigmaRRLayer;
+
+  }
+
+
+  return SigmaRRMax;
+
+
+}
+
+
 float HGCalShowerShape::eMax(const l1t::HGCalMulticluster& c3d) const {
   
   std::unordered_map<int, float> layer_energy;
@@ -326,56 +372,48 @@ float HGCalShowerShape::eMax(const l1t::HGCalMulticluster& c3d) const {
 
 
 
-
 float HGCalShowerShape::sigmaZZ(const l1t::HGCalMulticluster& c3d) const {
 
-  std::vector<float> tc_energy ; 
-  std::vector<float> tc_z ;
-  
   const edm::PtrVector<l1t::HGCalCluster>& clustersPtrs = c3d.constituents();
   
+  std::vector<std::pair<float,float> > tc_energy_z ; 
+  
   for(const auto& clu : clustersPtrs){
     
     const edm::PtrVector<l1t::HGCalTriggerCell>& triggerCells = clu->constituents();
     
     for(const auto& tc : triggerCells){
-      
-      tc_energy.emplace_back(tc->energy());
-      tc_z.emplace_back(tc->position().z());
-      
+           
+      tc_energy_z.emplace_back( std::make_pair(tc->energy(),tc->position().z()) );
+
     }
-    
+
   }
-  
-  float Szz = sigmaZZ(tc_energy,tc_z);
-  
-  tc_energy.clear();
-  tc_z.clear();
+
+  float z_mean = meanX(tc_energy_z);
+  float Szz = sigmaXX(tc_energy_z,z_mean);
   
   return Szz;
-
+  
 }
 
 
 
+
+
 float HGCalShowerShape::sigmaEtaEtaTot(const l1t::HGCalCluster& c2d) const {
 
   const edm::PtrVector<l1t::HGCalTriggerCell>& cellsPtrs = c2d.constituents();
   
-  std::vector<float> tc_energy ; 
-  std::vector<float> tc_eta ;
+  std::vector<std::pair<float,float> > tc_energy_eta ; 
   
   for(const auto& cell : cellsPtrs){
     
-    tc_energy.emplace_back(cell->energy());
-    tc_eta.emplace_back(cell->eta());
+    tc_energy_eta.emplace_back( std::make_pair(cell->energy(),cell->eta()) ); 
     
   }
   
-  float See = sigmaEtaEta(tc_energy,tc_eta);
-  
-  tc_energy.clear();
-  tc_eta.clear();
+  float See = sigmaXX(tc_energy_eta,c2d.eta());
   
   return See;
   
@@ -387,23 +425,39 @@ float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalCluster& c2d) const {
 
   const edm::PtrVector<l1t::HGCalTriggerCell>& cellsPtrs = c2d.constituents();
   
-  std::vector<float> tc_energy ; 
-  std::vector<float> tc_phi ;
+  std::vector<std::pair<float,float> > tc_energy_phi ; 
   
   for(const auto& cell : cellsPtrs){
     
-    tc_energy.emplace_back(cell->energy());
-    tc_phi.emplace_back(cell->phi());
+    tc_energy_phi.emplace_back( std::make_pair(cell->energy(),cell->phi()) ); 
     
   }
   
-  float Spp = sigmaPhiPhi(tc_energy,tc_phi);
-  
-  tc_energy.clear();
-  tc_phi.clear();
+  float Spp = sigmaXX(tc_energy_phi,c2d.phi());
   
   return Spp;
   
-}  
+}      
+
+
+
 
+float HGCalShowerShape::sigmaRRTot(const l1t::HGCalCluster& c2d) const {
 
+  const edm::PtrVector<l1t::HGCalTriggerCell>& cellsPtrs = c2d.constituents();
+  
+  std::vector<std::pair<float,float> > tc_energy_r ; 
+  
+  for(const auto& cell : cellsPtrs){
+    
+    float r = std::sqrt( pow(cell->position().x(),2) + pow(cell->position().y(),2) )/std::abs(cell->position().z());    
+    tc_energy_r.emplace_back( std::make_pair(cell->energy(),r) ); 
+    
+  }
+  
+  float r_mean = meanX(tc_energy_r);
+  float Srr = sigmaXX(tc_energy_r,r_mean);
+  
+  return Srr;
+  
+}