From 9c8d9257622140e2215bf07e5af4c88e77f49a25 Mon Sep 17 00:00:00 2001
From: dglazier <Derek.Glazier@glasgow.ac.uk>
Date: Tue, 16 Feb 2021 13:49:51 +0000
Subject: [PATCH] fix bug in TCM limits when more than 1 decaying particle

---
 core/DecayModel.cpp | 39 +++++++++++++++++++++++++++++----------
 core/DecayModel.h   |  3 ++-
 core/Interface.h    |  2 +-
 core/Particle.h     |  8 ++++++--
 4 files changed, 38 insertions(+), 14 deletions(-)

diff --git a/core/DecayModel.cpp b/core/DecayModel.cpp
index a9b0b18..c940c01 100644
--- a/core/DecayModel.cpp
+++ b/core/DecayModel.cpp
@@ -20,8 +20,9 @@ namespace elSpectro{
     //store list of unstable particles which decay
     for(auto prod: _products){ 
       auto dp=dynamic_cast<DecayingParticle*>(prod);
-      if(dp!=nullptr)
+      if(dp!=nullptr){
 	_unstables.push_back(dp);
+      }
       else
 	_stables.push_back(prod);
     }
@@ -46,10 +47,22 @@ namespace elSpectro{
   
   void DecayModel::PostInit(ReactionInfo* info){
     if(_unstables.empty()) return;
+    std::vector<double> prodmasses;
     for(auto& p:_unstables){
         p->PostInit(info);
+	prodmasses.push_back(p->MinimumMassPossible());
+    }
+    
+    for(uint i=0;i<_unstables.size();++i){ //for each product
+      //calculate mass reserved for subsequent unstable products
+      float reserve=0;
+      for(uint j=i+1;j<_unstables.size();++j)
+	reserve+=prodmasses[j];
+      
+      _unstableReservedMass.push_back(reserve);
     }
-  };
+
+  }
 
   void  DecayModel::GetStableMasses( std::vector<double >& masses) const{
 
@@ -72,6 +85,7 @@ namespace elSpectro{
   }
   double DecayModel::PhaseSpaceWeightSq(double W){
     //if(Parent()->Pdg()==-2211)std::cout<<GetName()<<" DecayModel::PhaseSpaceWeightSq start "<<MinimumMassPossible()<<" "<<W<<std::endl;
+    //std::cout<<GetName()<<" DecayModel::PhaseSpaceWeightSq start "<<MinimumMassPossible()<<" "<<W<<std::endl;
     //Note use weight squared to reduce sqrt calls
 
     if(_products.size()!=2){
@@ -79,24 +93,29 @@ namespace elSpectro{
       exit(0);
     }
     double result=1;
-    
     double TCM=W;
-    for(auto* p:_stables){
+ 
+   for(auto* p:_stables){
       TCM-=p->Mass();
     }
-
+  
+    uint iu=0; //synch reserved mass vector
     for(auto* p:_unstables){
       //This should be the only call to DetermineDynamicMass in the code.
       //Unless somewhere else uses LockMass in which case
       //this call to DetermineDynamicMass will not change its value
-      p->DetermineDynamicMass(-1,TCM);
+
+      //Note in case there are additional unstable particle we
+      //must subtract off their minimum masses
+      p->DetermineDynamicMass(-1,TCM-_unstableReservedMass[iu++]);
       TCM-=p->Mass();
-      if(TCM<0){ std::cout<<Parent()->Pdg()<<" "<<GetName()<<" "<< MinimumMassPossible()<<" W "<<W<<" T "<<TCM<<" "<<p->Mass()<<" "<<std::endl;return 0;}//below threshold, start again
+      if(TCM<0){ std::cout<<"DecayModel::PhaseSpaceWeightSq "<<Parent()->Pdg()<<" "<<GetName()<<" "<< MinimumMassPossible()<<" W "<<W<<" T "<<TCM<<" "<<p->Mass()<<" "<<std::endl;return 0;}//below threshold, start again
     }
-    //std::cout<<std::endl;
   
-      if(TCM<0){  std::cout<<Parent()->Pdg()<<" "<<GetName()<<" "<< MinimumMassPossible()<<" W "<<W<<" T "<<TCM<<" after stables "<<std::endl; return 0;}//below threshold, start again
-      result  *= kine::PDK2(W,_products[0]->Mass(),_products[1]->Mass());
+    if(TCM<0){  std::cout<<"DecayModel::PhaseSpaceWeightSq "<<Parent()->Pdg()<<" "<<GetName()<<" "<< MinimumMassPossible()<<" W "<<W<<" T "<<TCM<<" after stables "<<std::endl; return 0;}//below threshold, start again
+      
+  
+   result  *= kine::PDK2(W,_products[0]->Mass(),_products[1]->Mass());
     
     for(auto* p:_unstables)
       result*=p->PhaseSpaceWeightSq();
diff --git a/core/DecayModel.h b/core/DecayModel.h
index 6110cf7..7c40334 100644
--- a/core/DecayModel.h
+++ b/core/DecayModel.h
@@ -117,7 +117,8 @@ namespace elSpectro{
     particle_ptrs _products;
     particle_ptrs _stables; //products which are stable
     decaying_ptrs _unstables; //products which decay
-
+    std::vector<double> _unstableReservedMass; //mass reseved for other unstable products
+    
     mutable LorentzVector _parent;
     mutable double _sumOfMasses=0;
 
diff --git a/core/Interface.h b/core/Interface.h
index e65644f..9d0dc90 100644
--- a/core/Interface.h
+++ b/core/Interface.h
@@ -105,5 +105,5 @@ namespace elSpectro{
 	
     return  generator().Reaction();
   }
-  
+ 
 }//namespace elSpectro
diff --git a/core/Particle.h b/core/Particle.h
index 2ff8a6e..2611a83 100644
--- a/core/Particle.h
+++ b/core/Particle.h
@@ -129,17 +129,21 @@ namespace elSpectro{
       
       if(_massDist==nullptr ) return; //stick at pdgMass
       if(_massLocked==true) return; //someone else in charge...
-      
       _dynamicMass=-1;
       auto minposs = MinimumMassPossible();
       while(_dynamicMass<minposs){
 	auto minRange = xmin==-1?minposs:xmin;
 	auto maxRange = xmax==-1?_massDist->GetMaxX():xmax;
+	if(minRange>maxRange){//unphysical
+	  std::cout<<"Warning  Particle::DetermineDynamicMass min "<<minRange<<" greater than max "<<maxRange<<" for "<<_pdg<<std::endl;
+	  _dynamicMass=minRange;
+	  break;
+	}
 	_dynamicMass= _massDist->SampleSingle(minRange,maxRange);
 	
 	//need a weight for "envelope"
 	_massWeight =_massDist->GetCurrentWeight();
-	//	std::cout<<_pdg<<"  DetermineDynamicMass( "<<MinimumMassPossible()<<" "<<_dynamicMass<<" "<<_massWeight<<std::endl;
+	//	std::cout<<_pdg<<"  DetermineDynamicMass( "<<MinimumMassPossible()<<" "<<_dynamicMass<<" "<<_massWeight<<" "<<minRange<<" "<<maxRange<<std::endl;
 	
       }
       SetP4M(_dynamicMass);