put checks in particle having sufficient mass possibility to decay to…

… its products
dglazier · May 3, 2022 · b26f798 · b26f798
1 parent 75b861e
commit b26f798
Show file tree

Hide file tree

Showing 8 changed files with 84 additions and 8 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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 );