diff --git a/core/DecayingParticle.cpp b/core/DecayingParticle.cpp
index 01f11ff..2a76983 100644
--- a/core/DecayingParticle.cpp
+++ b/core/DecayingParticle.cpp
@@ -55,6 +55,19 @@ namespace elSpectro{
       }
     }
 
+    if(Pdg()!=0&&Pdg()!=-2211){//for real particles
+      Double_t productMasses = 0.0;
+      for(auto* prod: products){
+	productMasses+=prod->MinimumMassPossible();
+      }
+      auto mmp=MaximumMassPossible();
+      if(mmp<productMasses){
+	std::cerr<<"DecayingParticle::PostInit  insufficient mass to decay to its products, max mass = "<<mmp<<" while product masses "<<productMasses<<std::endl;
+	Print();
+	exit(0);
+      }
+    }
+    
     if(_decay)_decay->PostInit(info);
     if(_decayer)_decayer->PostInit(info);
  
diff --git a/core/DecayingParticle.h b/core/DecayingParticle.h
index c29644c..154b11c 100644
--- a/core/DecayingParticle.h
+++ b/core/DecayingParticle.h
@@ -45,6 +45,18 @@ namespace elSpectro{
     
     // virtual const CurrentEventInfo* EventInfo() const {return nullptr;}
 
+    double MaximumMassPossible() const  noexcept override {
+
+      Double_t maxMass=0;
+      if(MassDistribution()!=nullptr){
+	maxMass=MassDistribution()->GetMaxX();
+      }
+      else if(Pdg()!=-2211)
+	maxMass = Particle::MaximumMassPossible();
+
+      return maxMass;
+    }
+    
     double MinimumMassPossible() const  noexcept override {
       if(_minMass) return _minMass;
       
diff --git a/core/DistTF1.cpp b/core/DistTF1.cpp
index c314ff6..5df8fb4 100644
--- a/core/DistTF1.cpp
+++ b/core/DistTF1.cpp
@@ -1,4 +1,5 @@
 #include "DistTF1.h"
+#include "TH1.h"
 namespace elSpectro{
 
   DistTF1::DistTF1(const TF1& ff):
@@ -7,5 +8,27 @@ namespace elSpectro{
     _min_val = _tf1.GetMinimum();
     _tf1.SetNpx(1E4);
   }
- 
+  double DistTF1::GetMinX() const noexcept{
+
+    Int_t ib=0;
+    for(ib=1;ib<_tf1.GetHistogram()->GetNbinsX();ib++){
+     auto xval= _tf1.GetHistogram()->GetBinCenter(ib);
+     auto fval=_tf1.Eval(xval);
+     if(fval>0)break;
+    }
+    
+    return _tf1.GetHistogram()->GetBinLowEdge(ib);
+  }
+  double DistTF1::GetMaxX() const noexcept{
+    Int_t ib=0;
+    for(ib=_tf1.GetHistogram()->GetNbinsX(); ib>0 ;ib--){
+     auto xval= _tf1.GetHistogram()->GetBinCenter(ib);
+     auto fval=_tf1.Eval(xval);
+     if(fval>0)break;
+    }
+    double max=_tf1.GetHistogram()->GetBinLowEdge(ib)+_tf1.GetHistogram()->GetBinWidth(ib);
+    return max;
+
+  }
+
 }
diff --git a/core/DistTF1.h b/core/DistTF1.h
index b86d644..641d616 100644
--- a/core/DistTF1.h
+++ b/core/DistTF1.h
@@ -39,8 +39,8 @@ namespace elSpectro{
 
     double GetX() const noexcept { return _x;}
     
-    double GetMinX() const noexcept final{return _tf1.GetXmin();}
-    double GetMaxX() const noexcept final{return _tf1.GetXmax();}
+    double GetMinX() const noexcept final;
+    double GetMaxX() const noexcept final;
 
     double GetWeightFor(double valX) const {return _tf1.Eval(valX)/_max_val;}
     double GetValueFor(double valX,double valY=0) final {return _tf1.Eval(valX);}
diff --git a/core/DistTH1.cpp b/core/DistTH1.cpp
index f5973d4..cb6f181 100644
--- a/core/DistTH1.cpp
+++ b/core/DistTH1.cpp
@@ -7,5 +7,25 @@ namespace elSpectro{
     _min_val = ff.GetMinimum();
  
   }
- 
+  double DistTH1::GetMinX() const noexcept{
+
+    Int_t ib=0;
+    for(ib=1;ib<_th1.GetNbinsX();ib++){
+     auto fval=_th1.GetBinContent(ib);
+     if(fval>0)break;
+    }
+    
+    return _th1.GetBinLowEdge(ib);
+  }
+  double DistTH1::GetMaxX() const noexcept{
+    Int_t ib=0;
+    for(ib=_th1.GetNbinsX(); ib>0 ;ib--){
+     auto fval=_th1.GetBinContent(ib);
+     if(fval>0)break;
+    }
+    
+    return _th1.GetBinLowEdge(ib)+_th1.GetBinWidth(ib);
+
+  }
+
 }
diff --git a/core/DistTH1.h b/core/DistTH1.h
index c37a950..03bf27f 100644
--- a/core/DistTH1.h
+++ b/core/DistTH1.h
@@ -33,8 +33,8 @@ namespace elSpectro{
 
     double GetX() const noexcept { return _x;}
     
-    double GetMinX() const noexcept final{return _th1.GetXaxis()->GetXmin();}
-    double GetMaxX() const noexcept final{return _th1.GetXaxis()->GetXmax();}
+    double GetMinX() const noexcept final;
+    double GetMaxX() const noexcept final;
 
     //  double GetWeightFor(double valX)  {return  (static_cast<TH1D*>(&_th1))->GetBinContent((static_cast<TH1D*>(&_th1))->FindBin(valX))/_max_val;}
     // double GetWeightFor(double valX)  {return (static_cast<TH1D*>(&_th1))->Interpolate(valX)/_max_val;}
diff --git a/core/Particle.h b/core/Particle.h
index ea4f0d7..3b91b41 100644
--- a/core/Particle.h
+++ b/core/Particle.h
@@ -85,6 +85,9 @@ namespace elSpectro{
     virtual double MinimumMassPossible()const  noexcept{
       return  PdgMass();
     }
+    virtual double MaximumMassPossible()const  noexcept{
+      return  PdgMass();
+    }
     
     double MassWeight() const noexcept {
       return _massWeight;
diff --git a/core/ParticleManager.cpp b/core/ParticleManager.cpp
index d00586b..6b25e33 100644
--- a/core/ParticleManager.cpp
+++ b/core/ParticleManager.cpp
@@ -46,9 +46,14 @@ namespace elSpectro{
     auto pdg=p->Pdg();
     //assign mass distribution if exists
     //note maybe this should be just for unstable
-    if(_massDist.count(pdg)!=0)
+    if(_massDist.count(pdg)!=0){
       p->SetMassDist(_massDist.at(pdg).get());
-     
+
+      //Just take minimum mass as pdg value as default in case needed
+      if(dynamic_cast<DecayingParticle* >(p)!=nullptr)
+	dynamic_cast<DecayingParticle* >(p)->SetPdgMass(_massDist.at(pdg)->GetMinX());
+    }
+    
     auto dp=dynamic_cast<DecayingParticle* >( p );
     auto cp=dynamic_cast<CollidingParticle* >( p );