Differential D90 Diff 388 src/EvtGenModels/EvtFourBodyPhsp.cpp

Changeset View

Standalone View

src/EvtGenModels/EvtFourBodyPhsp.cpp

Show All 12 Lines
* but WITHOUT ANY WARRANTY; without even the implied warranty of *		* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *		* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *		* GNU General Public License for more details. *
* *		* *
* You should have received a copy of the GNU General Public License *		* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *		* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/		***********************************************************************/

#include <cmath>		#include "EvtGenModels/EvtFourBodyPhsp.hh"

#include "EvtGenBase/EvtKine.hh"		#include "EvtGenBase/EvtKine.hh"
#include "EvtGenBase/EvtPDL.hh"		#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"		#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"		#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"		#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtFourBodyPhsp.hh"
		#include <cmath>

std::string EvtFourBodyPhsp::getName()		std::string EvtFourBodyPhsp::getName()
{		{
return "FOURBODYPHSP";		return "FOURBODYPHSP";
}		}

EvtDecayBase* EvtFourBodyPhsp::clone()		EvtDecayBase* EvtFourBodyPhsp::clone()
{		{
▲ Show 20 Lines • Show All 121 Lines • ▼ Show 20 Lines	void EvtFourBodyPhsp::initProbMax()
double startM34 = m_m34Min + ( m_m34Max - m_m34Min ) / 20.;		double startM34 = m_m34Min + ( m_m34Max - m_m34Min ) / 20.;
bool contCond = true;		bool contCond = true;
int iteration = 0;		int iteration = 0;

auto parent = getParentId();		auto parent = getParentId();
double mMother = EvtPDL::getMaxMass( parent );		double mMother = EvtPDL::getMaxMass( parent );

double funcValue = 0;		double funcValue = 0;
while (contCond){		while ( contCond ) {
++iteration;		++iteration;
double currentM12 = startM12;		double currentM12 = startM12;
double currentM34 = startM34;		double currentM34 = startM34;
funcValue = phspFactor( mMother, currentM12, currentM34,		funcValue = phspFactor( mMother, currentM12, currentM34,
m_daughterMasses )[0];		m_daughterMasses )[0];
// Maximum along m12		// Maximum along m12
double step = ( m_m12Max - m_m12Min ) / 100.;		double step = ( m_m12Max - m_m12Min ) / 100.;
while ( step > 1e-4 ) {		while ( step > 1e-4 ) {
▲ Show 20 Lines • Show All 45 Lines • ▼ Show 20 Lines	while ( contCond ) {
} else if ( value2 > funcValue ) {		} else if ( value2 > funcValue ) {
currentM34 = point2;		currentM34 = point2;
funcValue = value2;		funcValue = value2;
}		}
step /= 2.;		step /= 2.;
}		}

// Check termination condition		// Check termination condition
double m12Diff = currentM12 - startM12;		double m12Diff = currentM12 - startM12;
double m34Diff = currentM34 - startM34;		double m34Diff = currentM34 - startM34;
double distSq = m12Diff * m12Diff + m34Diff * m34Diff;		double distSq = m12Diff * m12Diff + m34Diff * m34Diff;
if (distSq < 1e-8 \|\| iteration > 50){		if ( distSq < 1e-8 \|\| iteration > 50 ) {
contCond = false;		contCond = false;
}		}
startM12 = currentM12;		startM12 = currentM12;
startM34 = currentM34;		startM34 = currentM34;
}		}

setProbMax( funcValue * 1.05 );		setProbMax( funcValue * 1.05 );
}		}

void EvtFourBodyPhsp::decay( EvtParticle* parent )		void EvtFourBodyPhsp::decay( EvtParticle* parent )
{		{

parent->makeDaughters( getNDaug(), getDaugs() );		parent->makeDaughters( getNDaug(), getDaugs() );
bool massTreeStatus = parent->generateMassTree();		bool massTreeStatus = parent->generateMassTree();
if ( !massTreeStatus ) {		if ( !massTreeStatus ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )		EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Failed to generate daughters masses in EvtFourBodyPhsp."		<< "Failed to generate daughters masses in EvtFourBodyPhsp."
<< std::endl;		<< std::endl;
::abort();		::abort();
}		}
▲ Show 20 Lines • Show All 42 Lines • ▼ Show 20 Lines	void EvtFourBodyPhsp::decay( EvtParticle* parent )
const double m34 = masses.second;		const double m34 = masses.second;

// calculate probability, it will return array with 4 elements with		// calculate probability, it will return array with 4 elements with
// probability, q, p1 and p3		// probability, q, p1 and p3
auto probEval = phspFactor( mMother, m12, m34, m_daughterMasses );		auto probEval = phspFactor( mMother, m12, m34, m_daughterMasses );
setProb( probEval[0] );		setProb( probEval[0] );

// initialise kinematics		// initialise kinematics
const double cosTheta1 = EvtRandom::Flat(-1.0, 1.0);		const double cosTheta1 = EvtRandom::Flat( -1.0, 1.0 );
const double sinTheta1 = std::sqrt( 1 - cosTheta1 * cosTheta1 );		const double sinTheta1 = std::sqrt( 1 - cosTheta1 * cosTheta1 );
const double cosTheta3 = EvtRandom::Flat(-1.0, 1.0);		const double cosTheta3 = EvtRandom::Flat( -1.0, 1.0 );
const double sinTheta3 = std::sqrt( 1 - cosTheta3 * cosTheta3 );		const double sinTheta3 = std::sqrt( 1 - cosTheta3 * cosTheta3 );
const double phi = EvtRandom::Flat( 0., EvtConst::twoPi );		const double phi = EvtRandom::Flat( 0., EvtConst::twoPi );
// m12 and m34 are put along z-axis, 1 and 2 go to x-z plane and 3-4		// m12 and m34 are put along z-axis, 1 and 2 go to x-z plane and 3-4
// plane is rotated by phi compared to 1-2 plane. All momenta are set		// plane is rotated by phi compared to 1-2 plane. All momenta are set
// in 12 and 34 rest frames and then boosted to parent rest frame		// in 12 and 34 rest frames and then boosted to parent rest frame
const double p1x = probEval[2] * sinTheta1;		const double p1x = probEval[2] * sinTheta1;
const double p1z = probEval[2] * cosTheta1;		const double p1z = probEval[2] * cosTheta1;
const double p1Sq = probEval[2] * probEval[2];		const double p1Sq = probEval[2] * probEval[2];
▲ Show 20 Lines • Show All 134 Lines • Show Last 20 Lines