No OneTemporary
Actions

Size

31 KB

Subscribers

None

View Options

	diff --git a/Hadronization/ClusterDecayer.cc b/Hadronization/ClusterDecayer.cc
	--- a/Hadronization/ClusterDecayer.cc
	+++ b/Hadronization/ClusterDecayer.cc
	@@ -1,469 +1,472 @@
	// -- C++ --
	//
	// ClusterDecayer.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ClusterDecayer class.
	//

	#include "ClusterDecayer.h"
	#include <ThePEG/Interface/ClassDocumentation.h>
	#include <ThePEG/Interface/Reference.h>
	#include <ThePEG/Interface/Parameter.h>
	#include <ThePEG/Interface/Switch.h>
	#include <ThePEG/Persistency/PersistentOStream.h>
	#include <ThePEG/Persistency/PersistentIStream.h>
	#include <ThePEG/PDT/EnumParticles.h>
	#include <ThePEG/Repository/EventGenerator.h>
	#include "Herwig++/Utilities/Kinematics.h"
	#include "CheckId.h"
	#include "Herwig++/Utilities/Smearing.h"
	#include "Cluster.h"
	#include <ThePEG/Utilities/DescribeClass.h>
	+#include <ThePEG/Repository/RandomGenerator.h>
	+#include <ThePEG/Repository/StandardRandom.h>

	using namespace Herwig;

	DescribeClass<ClusterDecayer,Interfaced>
	describeClusterDecayer("Herwig::ClusterDecayer","");

	ClusterDecayer::ClusterDecayer() :
	_clDirLight(1),
	_clDirBottom(1),
	_clDirCharm(1),
	_clDirExotic(1),
	_clSmrLight(0.0),
	_clSmrBottom(0.0),
	_clSmrCharm(0.0),
	_clSmrExotic(0.0),
	_onshell(false),
	_masstry(20)
	{}

	IBPtr ClusterDecayer::clone() const {
	return new_ptr(*this);
	}

	IBPtr ClusterDecayer::fullclone() const {
	return new_ptr(*this);
	}

	void ClusterDecayer::persistentOutput(PersistentOStream & os) const
	{
	os << _hadronsSelector << _clDirLight << _clDirBottom
	<< _clDirCharm << _clDirExotic << _clSmrLight << _clSmrBottom
	<< _clSmrCharm << _clSmrExotic << _onshell << _masstry;
	}

	void ClusterDecayer::persistentInput(PersistentIStream & is, int) {
	is >> _hadronsSelector >> _clDirLight >> _clDirBottom
	>> _clDirCharm >> _clDirExotic >> _clSmrLight >> _clSmrBottom
	>> _clSmrCharm >> _clSmrExotic >> _onshell >> _masstry;
	}


	void ClusterDecayer::Init() {

	static ClassDocumentation<ClusterDecayer> documentation
	("This class is responsible for the two-body decays of normal clusters");

	static Reference<ClusterDecayer,HadronSelector>
	interfaceHadronSelector("HadronSelector",
	"A reference to the HadronSelector object",
	&Herwig::ClusterDecayer::_hadronsSelector,
	false, false, true, false);

	//ClDir for light, Bottom, Charm and exotic (e.g Susy) quarks

	static Switch<ClusterDecayer,bool> interfaceClDirLight
	("ClDirLight",
	"Cluster direction for light quarks",
	&ClusterDecayer::_clDirLight, true, false, false);
	static SwitchOption interfaceClDirLightPreserve
	(interfaceClDirLight,
	"Preserve",
	"Preserve the direction of the quark from the perturbative stage"
	" as the direction of the meson produced from it",
	true);
	static SwitchOption interfaceClDirLightIsotropic
	(interfaceClDirLight,
	"Isotropic",
	"Assign the direction of the meson randomly",
	false);

	static Switch<ClusterDecayer,bool> interfaceClDirBottom
	("ClDirBottom",
	"Cluster direction for bottom quarks",
	&ClusterDecayer::_clDirBottom, true, false, false);
	static SwitchOption interfaceClDirBottomPreserve
	(interfaceClDirBottom,
	"Preserve",
	"Preserve the direction of the quark from the perturbative stage"
	" as the direction of the meson produced from it",
	true);
	static SwitchOption interfaceClDirBottomIsotropic
	(interfaceClDirBottom,
	"Isotropic",
	"Assign the direction of the meson randomly",
	false);

	static Switch<ClusterDecayer,bool> interfaceClDirCharm
	("ClDirCharm",
	"Cluster direction for charm quarks",
	&ClusterDecayer::_clDirCharm, true, false, false);
	static SwitchOption interfaceClDirCharmPreserve
	(interfaceClDirCharm,
	"Preserve",
	"Preserve the direction of the quark from the perturbative stage"
	" as the direction of the meson produced from it",
	true);
	static SwitchOption interfaceClDirCharmIsotropic
	(interfaceClDirCharm,
	"Isotropic",
	"Assign the direction of the meson randomly",
	false);

	static Switch<ClusterDecayer,bool> interfaceClDirExotic
	("ClDirExotic",
	"Cluster direction for exotic quarks",
	&ClusterDecayer::_clDirExotic, true, false, false);
	static SwitchOption interfaceClDirExoticPreserve
	(interfaceClDirExotic,
	"Preserve",
	"Preserve the direction of the quark from the perturbative stage"
	" as the direction of the meson produced from it",
	true);
	static SwitchOption interfaceClDirExoticIsotropic
	(interfaceClDirExotic,
	"Isotropic",
	"Assign the direction of the meson randomly",
	false);

	// ClSmr for ligth, Bottom, Charm and exotic (e.g. Susy) quarks
	static Parameter<ClusterDecayer,double>
	interfaceClSmrLight ("ClSmrLight", "cluster direction Gaussian smearing for light quark",
	&ClusterDecayer::_clSmrLight, 0, 0.0 , 0.0 , 2.0,false,false,false);
	static Parameter<ClusterDecayer,double>
	interfaceClSmrBottom ("ClSmrBottom", "cluster direction Gaussian smearing for b quark",
	&ClusterDecayer::_clSmrBottom, 0, 0.0 , 0.0 , 2.0,false,false,false);
	static Parameter<ClusterDecayer,double>
	interfaceClSmrCharm ("ClSmrCharm", "cluster direction Gaussian smearing for c quark",
	&ClusterDecayer::_clSmrCharm, 0, 0.0 , 0.0 , 2.0,false,false,false);
	static Parameter<ClusterDecayer,double>
	interfaceClSmrExotic ("ClSmrExotic", "cluster direction Gaussian smearing for exotic quark",
	&ClusterDecayer::_clSmrExotic, 0, 0.0 , 0.0 , 2.0,false,false,false);

	static Switch<ClusterDecayer,bool> interfaceOnShell
	("OnShell",
	"Whether or not the hadrons produced should by on shell or generated using the"
	" mass generator.",
	&ClusterDecayer::_onshell, false, false, false);
	static SwitchOption interfaceOnShellOnShell
	(interfaceOnShell,
	"Yes",
	"Produce the hadrons on shell",
	true);
	static SwitchOption interfaceOnShellOffShell
	(interfaceOnShell,
	"No",
	"Generate the masses using the mass generator.",
	false);

	static Parameter<ClusterDecayer,unsigned int> interfaceMassTry
	("MassTry",
	"The number attempts to generate the masses of the hadrons produced"
	" in the cluster decay.",
	&ClusterDecayer::_masstry, 20, 1, 50,
	false, false, Interface::limited);

	}


	void ClusterDecayer::decay(const ClusterVector & clusters, tPVector & finalhadrons)
	{
	// Loop over all clusters, and if they are not too heavy (that is
	// intermediate clusters that have undergone to fission) or not
	// too light (that is final clusters that have been already decayed
	// into single hadron) then decay them into two hadrons.
	for (ClusterVector::const_iterator it = clusters.begin();
	it != clusters.end(); ++it) {
	if ((it)->isAvailable() && !(it)->isStatusFinal()
	&& (*it)->isReadyToDecay()) {
	pair<PPtr,PPtr> prod = decayIntoTwoHadrons(*it);
	(*it)->addChild(prod.first);
	(*it)->addChild(prod.second);
	finalhadrons.push_back(prod.first);
	finalhadrons.push_back(prod.second);
	}
	}
	}


	pair<PPtr,PPtr> ClusterDecayer::decayIntoTwoHadrons(tClusterPtr ptr) {
	using Constants::pi;
	using Constants::twopi;
	// To decay the cluster into two hadrons one must distinguish between
	// constituent quarks (or diquarks) that originate from perturbative
	// processes (hard process or parton shower) from those that are generated
	// by the non-perturbative gluon splitting or from fission of heavy clusters.
	// In the latter case the two body decay is assumed to be isotropic.
	// In the former case instead, if proper flag are activated, the two body
	// decay is assumed to "remember" the direction of the constituents inside
	// the cluster, in the cluster frame. The actual smearing of the hadron
	// directions around the direction of the constituents, in the cluster
	// frame, can be different between non-b hadrons and b-hadrons, but is given
	// by the same functional form:
	// cosThetaSmearing = 1 + smearFactor * log( rnd() )
	// (repeated until cosThetaSmearing > -1)
	// where the factor smearFactor is different between b- and non-b hadrons.
	//
	// We need to require (at least at the moment, maybe in the future we
	// could change it) that the cluster has exactly two components.
	// If this is not the case, then send a warning because it is not suppose
	// to happen, and then return.
	if ( ptr->numComponents() != 2 ) {
	generator()->logWarning( Exception("ClusterDecayer::decayIntoTwoHadrons "
	"*Still cluster with not exactly 2 components* ",
	Exception::warning) );
	return pair<PPtr,PPtr>();
	}

	// Extract the id and particle pointer of the two components of the cluster.
	tPPtr ptr1 = ptr->particle(0);
	tPPtr ptr2 = ptr->particle(1);
	tcPDPtr ptr1data = ptr1->dataPtr();
	tcPDPtr ptr2data = ptr2->dataPtr();

	bool isHad1FlavSpecial = false;
	bool cluDirHad1 = _clDirLight;
	double cluSmearHad1 = _clSmrLight;
	bool isHad2FlavSpecial = false;
	bool cluDirHad2 = _clDirLight;
	double cluSmearHad2 = _clSmrLight;

	if (CheckId::isExotic(ptr1data)) {
	isHad1FlavSpecial = true;
	cluDirHad1 = _clDirExotic;
	cluSmearHad1 = _clSmrExotic;
	}
	else if (CheckId::hasBottom(ptr1data)) {
	isHad1FlavSpecial = true;
	cluDirHad1 = _clDirBottom;
	cluSmearHad1 = _clSmrBottom;
	}
	else if (CheckId::hasCharm(ptr1data)) {
	isHad1FlavSpecial = true;
	cluDirHad1 = _clDirCharm;
	cluSmearHad1 = _clSmrCharm;
	}

	if (CheckId::isExotic(ptr2data)) {
	isHad2FlavSpecial = true;
	cluDirHad2 = _clDirExotic;
	cluSmearHad2 = _clSmrExotic;
	}
	else if (CheckId::hasBottom(ptr2data)) {
	isHad2FlavSpecial = true;
	cluDirHad2 = _clDirBottom;
	cluSmearHad2 = _clSmrBottom;
	}
	else if (CheckId::hasCharm(ptr2data)) {
	isHad2FlavSpecial = true;
	cluDirHad2 = _clDirCharm;
	cluSmearHad2 = _clSmrCharm;
	}


	bool isOrigin1Perturbative = ptr->isPerturbative(0);
	bool isOrigin2Perturbative = ptr->isPerturbative(1);

	// We have to decide which, if any, of the two hadrons will have
	// the momentum, in the cluster parent frame, smeared around the
	// direction of its constituent (for Had1 is the one pointed by
	// ptr1, and for Had2 is the one pointed by ptr2).
	// This happens only if the flag _ClDirX is 1 and the constituent is
	// perturbative (that is not coming from nonperturbative gluon splitting
	// or cluster fission). In the case that both the hadrons satisfy this
	// two requirements (of course only one must be treated, because the other
	// one will have the momentum automatically fixed by the momentum
	// conservation) then more priority is given in the case of a b-hadron.
	// Finally, in the case that the two hadrons have same priority, then
	// we choose randomly, with equal probability, one of the two.

	int priorityHad1 = 0;
	if ( cluDirHad1 && isOrigin1Perturbative ) {
	priorityHad1 = isHad1FlavSpecial ? 2 : 1;
	}
	int priorityHad2 = 0;
	if ( cluDirHad2 && isOrigin2Perturbative ) {
	priorityHad2 = isHad2FlavSpecial ? 2 : 1;
	}
	if ( priorityHad2 && priorityHad1 == priorityHad2 && UseRandom::rndbool() ) {
	priorityHad2 = 0;
	}

	Lorentz5Momentum pClu = ptr->momentum();
	bool secondHad = false;
	Axis uSmear_v3;
	if ( priorityHad1 \|\| priorityHad2 ) {

	double cluSmear;
	Lorentz5Momentum pQ;
	if ( priorityHad1 > priorityHad2 ) {
	pQ = ptr1->momentum();
	cluSmear = cluSmearHad1;
	} else {
	pQ = ptr2->momentum();
	cluSmear = cluSmearHad2;
	secondHad = true;
	}

	// To find the momenta of the two hadrons in the parent cluster frame
	// we proceed as follows. First, we determine the unit vector parallel
	// to the direction of the constituent in the cluster frame. Then we
	// have to smear such direction using the following prescription:
	// --- in theta angle w.r.t. such direction (not the z-axis),
	// the drawing of the cosine of such angle is done via:
	// 1.0 + cluSmear*log( rnd() )
	// (repeated until it gives a value above -1.0)
	// --- in phi angle w.r.t. such direction, the drawing is simply flat.
	// Then, given the direction in the parent cluster frame of one of the
	// two hadrons, it is straighforward to get the momenta of both hadrons
	// (always in the same parent cluster frame).

	pQ.boost( -pClu.boostVector() ); // boost from Lab to Cluster frame
	uSmear_v3 = pQ.vect().unit();
	// skip if cluSmear is too small
	if ( cluSmear > 0.001 ) {
	// generate the smearing angle
	double cosSmear;
	do cosSmear = 1.0 + cluSmear*log( UseRandom::rnd() );
	while ( cosSmear < -1.0 );
	double sinSmear = sqrt( 1.0 - sqr(cosSmear) );
	// calculate rotation to z axis
	LorentzRotation rot;
	double sinth(sqrt(1.-sqr(uSmear_v3.z())));
	if(abs(uSmear_v3.perp2()/uSmear_v3.z())>1e-10)
	rot.setRotate(-acos(uSmear_v3.z()),
	Axis(-uSmear_v3.y()/sinth,uSmear_v3.x()/sinth,0.));
	// + random azimuthal rotation
	rot.rotateZ(UseRandom::rnd()*twopi);
	// set direction in rotated frame
	Lorentz5Vector<double> ltemp(0.,sinSmear,cosSmear,0.);
	// rotate back
	rot.invert();
	ltemp *= rot;
	uSmear_v3 = ltemp.vect();
	}
	}
	else {

	// Isotropic decay: flat in cosTheta and phi.
	uSmear_v3 = Axis(1.0, 0.0, 0.0); // just to set the rho to 1
	uSmear_v3.setTheta( acos( UseRandom::rnd( -1.0 , 1.0 ) ) );
	uSmear_v3.setPhi( UseRandom::rnd( -pi , pi ) );

	}

	pair<tcPDPtr,tcPDPtr> dataPair
	= _hadronsSelector->chooseHadronPair(ptr->mass(),
	ptr1data,
	ptr2data);
	if(dataPair.first == tcPDPtr() \|\|
	dataPair.second == tcPDPtr()) return pair<PPtr,PPtr>();

	// Create the two hadron particle objects with the specified id.
	PPtr ptrHad1,ptrHad2;
	// produce the hadrons on mass shell
	if(_onshell) {
	ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass());
	ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass());
	}
	// produce the hadrons with mass given by the mass generator
	else {
	unsigned int ntry(0);
	do {
	ptrHad1 = dataPair.first ->produceParticle();
	ptrHad2 = dataPair.second->produceParticle();
	++ntry;
	}
	while(ntry<_masstry&&ptrHad1->mass()+ptrHad2->mass()>ptr->mass());
	// if fails produce on shell and issue warning (should never happen??)
	if( ptrHad1->mass() + ptrHad2->mass() > ptr->mass() ) {
	generator()->log() << "Failed to produce off-shell hadrons in "
	<< "ClusterDecayer::decayIntoTwoHadrons producing hadrons "
	<< "on-shell" << endl;
	ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass());
	ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass());
	}
	}

	if (!ptrHad1 \|\| !ptrHad2) {
	ostringstream s;
	s << "ClusterDecayer::decayIntoTwoHadrons *Cannot create the two hadrons*"
	<< dataPair.first ->PDGName() << " and "
	<< dataPair.second->PDGName();
	cerr << s.str() << endl;
	generator()->logWarning( Exception(s.str(), Exception::warning) );
	} else {

	Lorentz5Momentum pHad1, pHad2; // 5-momentum vectors that we have to determine
	if ( secondHad ) uSmear_v3 *= -1.0;

	if (pClu.m() < ptrHad1->mass()+ptrHad2->mass() ) {
	throw Exception() << "Impossible Kinematics in ClusterDecayer::decayIntoTwoHadrons()"
	<< Exception::eventerror;
	}
	Kinematics::twoBodyDecay(pClu,ptrHad1->mass(),ptrHad2->mass(),uSmear_v3,
	pHad1,pHad2);
	ptrHad1->set5Momentum(pHad1);
	ptrHad2->set5Momentum(pHad2);

	// Determine the positions of the two children clusters.
	LorentzPoint positionHad1 = LorentzPoint();
	LorentzPoint positionHad2 = LorentzPoint();
	calculatePositions(pClu, ptr->vertex(), pHad1, pHad2, positionHad1, positionHad2);
	ptrHad1->setVertex(positionHad1);
	ptrHad2->setVertex(positionHad2);
	}
	return pair<PPtr,PPtr>(ptrHad1,ptrHad2);
	}


	void ClusterDecayer::
	calculatePositions(const Lorentz5Momentum &pClu,
	const LorentzPoint &positionClu,
	const Lorentz5Momentum &,
	const Lorentz5Momentum &,
	LorentzPoint &positionHad1,
	LorentzPoint &positionHad2 ) const {
	// First, determine the relative positions of the children hadrons
	// with respect to their parent cluster, in the cluster reference frame,
	// assuming gaussian smearing with width inversely proportional to the
	// parent cluster mass.
	+ ThePEG::StandardRandom* randomNumberStandardGenerator = new ThePEG::StandardRandom();
	Length smearingWidth = hbarc / pClu.m();
	LorentzDistance distanceHad[2];
	for ( int i = 0; i < 2; i++ ) { // i==0 is the first hadron; i==1 is the second one
	double delta[4]={0.,0.,0.,0.};
	- // smearing of the four components of the LorentzDistance
	- for ( int j = 0; j < 4; j++ ) {
	- while ( ! Smearing::gaussianSmearing( 0.0, smearingWidth/femtometer, delta[j] ) ) { }
	+ // smearing of the four components of the LorentzDistance, two at the same time to improve speed
	+ for ( int j = 0; j < 3; j += 2 ) {
	+ randomNumberStandardGenerator->rndGaussTwoNumbers(delta[j], delta[j+1], smearingWidth/femtometer, 0.0);
	}
	// set the distance
	delta[0] = abs(delta[0]) +sqrt(sqr(delta[1])+sqr(delta[2])+sqr(delta[3]));
	distanceHad[i] = LorentzVector<double>(delta[1],delta[2],delta[3],delta[0])*femtometer;
	// Boost such relative positions of the children hadrons,
	// with respect to their parent cluster,
	// from the cluster reference frame to the Lab frame.
	distanceHad[i].boost(pClu.boostVector());
	}
	// Finally, determine the absolute positions of the children hadrons
	// in the Lab frame.
	positionHad1 = distanceHad[0] + positionClu;
	positionHad2 = distanceHad[1] + positionClu;
	}

	diff --git a/Hadronization/ClusterDecayer.h b/Hadronization/ClusterDecayer.h
	--- a/Hadronization/ClusterDecayer.h
	+++ b/Hadronization/ClusterDecayer.h
	@@ -1,166 +1,166 @@
	// -- C++ --
	//
	// ClusterDecayer.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ClusterDecayer_H
	#define HERWIG_ClusterDecayer_H

	#include <ThePEG/Interface/Interfaced.h>
	#include <ThePEG/EventRecord/Step.h>
	#include "CluHadConfig.h"
	#include "HadronSelector.h"
	#include "ClusterDecayer.fh"

	namespace Herwig {
	using namespace ThePEG;

	/** \ingroup Hadronization
	* \class ClusterDecayer
	* \brief This class decays the "normal" clusters
	* \author Philip Stephens
	* \author Alberto Ribon
	*
	* This class decays the "normal" clusters, e.g. ones that are not heavy
	* enough for fission, and not too light to decay into one hadron.
	*
	* This class is directs the production of hadrons via 2-body cluster decays.
	* The selection of the hadron flavours is given by Herwig::HadronSelector.
	*
	* @see HadronSelector
	* @see \ref ClusterDecayerInterfaces "The interfaces"
	* defined for ClusterDecayer.
	*/
	class ClusterDecayer: public Interfaced {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Default constructor.
	*/
	ClusterDecayer();
	//@}

	/** Decays all remaining clusters into hadrons.
	* This routine decays the clusters that are left after
	* Herwig::ClusterFissioner::fission and
	* Herwig::LightClusterDecayer::decay have been called. These are all
	* the "normal" clusters which are not forced into hadrons by
	* the other functions.
	*/
	void decay(const ClusterVector & clusters, tPVector & finalhadrons)
	;

	public:

	/**
	* Standard ThePEG function for writing a persistent stream.
	*/
	void persistentOutput(PersistentOStream &) const;

	/**
	* Standard ThePEG function for reading from a persistent stream.
	*/
	void persistentInput(PersistentIStream &, int);

	/**
	* Standard Init function used to initialize the interfaces.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	/**
	* Private and non-existent assignment operator.
	*/
	ClusterDecayer & operator=(const ClusterDecayer &);

	public:

	/** Decays the cluster into two hadrons.
	*
	* This routine is used to take a given cluster and decay it into
	* two hadrons which are returned. If one of the constituents is from
	* the perturbative regime then the direction of the perturbative parton
	* is remembered and the decay is preferentially in that direction. The
	* direction of the decay is given by
	* \f[ \cos \theta = 1 + S \log r_1 \f]
	* where \f$ S \f$ is a parameter of the model and \f$ r_1 \f$ is a random
	* number [0,1].
	*/
	pair<PPtr,PPtr> decayIntoTwoHadrons(tClusterPtr ptr);

	private:

	/** Compute the positions of the new hadrons based on the clusters position.
	*
	* This method calculates the positions of the children hadrons by a
	- * call to Herwig::Smearing::gaussianSmearing with width inversely
	+ * call to ThePEG::RandomGenerator::rndGaussTwoNumbers with width inversely
	* proportional to the cluster mass, around the parent cluster position.
	*/
	void calculatePositions( const Lorentz5Momentum &, const LorentzPoint &,
	const Lorentz5Momentum &, const Lorentz5Momentum &,
	LorentzPoint &, LorentzPoint &) const;

	/**
	* Pointer to a Herwig::HadronSelector for choosing decay types
	*/
	Ptr<HadronSelector>::pointer _hadronsSelector;

	//@{
	/**
	* Whether a cluster decays along the perturbative parton direction.
	*/
	bool _clDirLight;
	bool _clDirBottom;
	bool _clDirCharm;
	bool _clDirExotic;

	/**
	* The S parameter from decayIntoTwoHadrons
	*/
	double _clSmrLight;
	double _clSmrBottom;
	double _clSmrCharm;
	double _clSmrExotic;
	//@}

	/**
	* Whether or not the hadrons produced should be on-shell
	* or generated used the MassGenerator
	*/
	bool _onshell;

	/**
	* Number of tries to generate the masses of the decay products
	*/
	unsigned int _masstry;

	};


	}

	#endif /* HERWIG_ClusterDecayer_H */
	diff --git a/Utilities/Makefile.am b/Utilities/Makefile.am
	--- a/Utilities/Makefile.am
	+++ b/Utilities/Makefile.am
	@@ -1,39 +1,37 @@
	SUBDIRS = XML Statistics
	noinst_LTLIBRARIES = libHwUtils.la
	pkglib_LTLIBRARIES = libHwRunDirectories.la
	libHwUtils_la_SOURCES = \
	EnumParticles.h \
	Interpolator.tcc Interpolator.h \
	Kinematics.cc Kinematics.h \
	-Smearing.cc Smearing.h \
	Maths.h Maths.cc \
	StandardSelectors.cc StandardSelectors.h\
	Histogram.cc Histogram.fh Histogram.h \
	GaussianIntegrator.cc GaussianIntegrator.h \
	GaussianIntegrator.tcc \
	Statistic.h HerwigStrategy.cc HerwigStrategy.h \
	GSLIntegrator.h GSLIntegrator.tcc \
	GSLBisection.h GSLBisection.tcc GSLHelper.h

	libHwRunDirectories_la_SOURCES = \
	RunDirectories.h RunDirectories.cc

	libHwUtils_la_LIBADD = \
	XML/libHwXML.la \
	Statistics/libHwStatistics.la

	CLEANFILES=
	include $(srcdir)/Makefile.am.versionstring

	check_PROGRAMS = utilities_test
	utilities_test_SOURCES = \
	tests/utilitiesTestsMain.cc \
	tests/utilitiesTestsGlobalFixture.h \
	tests/utilitiesTestsKinematics.h \
	tests/utilitiesTestsMaths.h \
	-tests/utilitiesTestsSmearing.h \
	tests/utilitiesTestsStatistic.h
	utilities_test_LDADD = $(BOOST_UNIT_TEST_FRAMEWORK_LIBS) $(BOOST_FILESYSTEM_LIBS) $(BOOST_SYSTEM_LIBS) $(THEPEGLIB) -ldl libHwUtils.la
	utilities_test_LDFLAGS = $(AM_LDFLAGS) -export-dynamic $(BOOST_UNIT_TEST_FRAMEWORK_LDFLAGS) $(BOOST_SYSTEM_LDFLAGS) $(BOOST_FILESYSTEM_LDFLAGS) $(THEPEGLDFLAGS)
	utilities_test_CPPFLAGS = $(AM_CPPFLAGS) $(BOOST_CPPFLAGS) -DHERWIG_PKGDATADIR="\"$(pkgdatadir)\"" -DHERWIG_PKGLIBDIR="\"$(pkglibdir)\"" -DTHEPEG_PKGLIBDIR="\"$(THEPEGLIBPATH)\""
	TESTS = utilities_test
	diff --git a/Utilities/tests/utilitiesTestSmearing.h b/Utilities/tests/utilitiesTestSmearing.h
	new file mode 100644
	--- /dev/null
	+++ b/Utilities/tests/utilitiesTestSmearing.h
	@@ -0,0 +1,102 @@
	+// -- C++ --
	+//
	+// utilitiesTestSmearing.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	+// Copyright (C) 2002-2011 The Herwig Collaboration, 2015 Marco A. Harrendorf
	+//
	+// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	+//
	+#ifndef HERWIG_Utilities_Test_Smearing_H
	+#define HERWIG_Utilities_Test_Smearing_H
	+
	+#include <boost/test/unit_test.hpp>
	+
	+#include "Herwig++/Utilities/Smearing.h"
	+#include <iostream>
	+#include <fstream>
	+
	+/*
	+ * Start of boost unit tests for Smearing.h
	+ *
	+ */
	+BOOST_AUTO_TEST_SUITE(utilitiesSmearingTest)
	+
	+/*
	+ * Boost unit tests
	+ *
	+ */
	+BOOST_AUTO_TEST_CASE(azimuthalSmearing)
	+{
	+ double vectorXvalue, vectorYvalue;
	+ double vectorLengthCalculated;
	+ unsigned int numberOfTrials = 10000;
	+ for(unsigned int vectorLength = 1; vectorLength < 4; ++vectorLength) {
	+ // Variables for statistical analysis
	+ double meanVectorLength = 0, meanVectorX = 0, meanVectorY = 0;
	+ unsigned int inverseSamplingEfficiency = 0;
	+ for(unsigned int i = 0; i < numberOfTrials; ++i) {
	+ if(Herwig::Smearing::azimuthalSmearing(vectorLength, vectorXvalue, vectorYvalue)) {
	+ vectorLengthCalculated = sqrt(vectorXvaluevectorXvalue + vectorYvaluevectorYvalue);
	+ BOOST_CHECK_CLOSE(vectorLengthCalculated, vectorLength, 1e-06);
	+ meanVectorLength += vectorLengthCalculated;
	+ meanVectorX += vectorXvalue;
	+ meanVectorY += vectorYvalue;
	+ }
	+ else {
	+ inverseSamplingEfficiency++;
	+ }
	+ }
	+ // Normalize mean values
	+ meanVectorLength = meanVectorLength/static_cast<double>(numberOfTrials - inverseSamplingEfficiency);
	+ meanVectorX = meanVectorX/static_cast<double>(numberOfTrials - inverseSamplingEfficiency);
	+ meanVectorY = meanVectorY/static_cast<double>(numberOfTrials - inverseSamplingEfficiency);
	+
	+ BOOST_CHECK_CLOSE(meanVectorLength , vectorLength, 0.001);
	+ BOOST_CHECK(meanVectorX < 1e-1 && meanVectorX > -1e-1);
	+ BOOST_CHECK(meanVectorY < 1e-1 && meanVectorY > -1e-1);
	+ BOOST_CHECK(inverseSamplingEfficiency < 0.25 * numberOfTrials);
	+ }
	+}
	+
	+BOOST_AUTO_TEST_CASE(gaussianSmearing)
	+{
	+ std::ofstream myfile;
	+ myfile.open ("out2.txt");
	+
	+ double gaussianValue;
	+ unsigned int numberOfTrials = 1000000;
	+ for(unsigned int gaussianMean = 1; gaussianMean < 2; ++gaussianMean) {
	+ for(double gaussianSigma = 0.25; gaussianSigma < 0.35; gaussianSigma += 0.1) {
	+ // Variables for statistical analysis
	+ double meanValue = 0, meanSigma = 0;
	+ unsigned int inverseSamplingEfficiency = 0;
	+ for(unsigned int i = 0; i < numberOfTrials; ++i) {
	+ if(Herwig::Smearing::gaussianSmearing(gaussianMean, gaussianSigma, gaussianValue)) {
	+ // Check for 8 sigma deviation
	+ BOOST_CHECK(gaussianValue < (gaussianMean + 8 * gaussianSigma));
	+ //std::cout << gaussianValue << " " << gaussianMean << " " << gaussianSigma << " " << (gaussianMean + 8 * gaussianSigma) << std::endl;
	+ myfile << gaussianValue << std::endl;
	+ BOOST_CHECK(gaussianValue > (gaussianMean - 8 * gaussianSigma));
	+ meanValue += gaussianValue;
	+ meanSigma += std::abs(gaussianValue - gaussianMean);
	+ }
	+ else {
	+ inverseSamplingEfficiency++;
	+ }
	+ }
	+ // Normalize mean values
	+ meanValue = meanValue/static_cast<double>(numberOfTrials - inverseSamplingEfficiency);
	+ meanSigma = meanSigma/static_cast<double>(numberOfTrials - inverseSamplingEfficiency);
	+ std::cout << meanValue << " " << meanSigma << std::endl;
	+ BOOST_CHECK_CLOSE(meanValue , gaussianMean, 0.5);
	+ BOOST_CHECK_CLOSE(meanSigma , gaussianSigma, 0.5);
	+ BOOST_CHECK(inverseSamplingEfficiency < 0.25 * numberOfTrials);
	+ }
	+ }
	+}
	+
	+
	+
	+BOOST_AUTO_TEST_SUITE_END()
	+
	+#endif /* HERWIG_Utilities_Test_Smearing_H */
	\ No newline at end of file
	diff --git a/Utilities/tests/utilitiesTestsGlobalFixture.h b/Utilities/tests/utilitiesTestsGlobalFixture.h
	--- a/Utilities/tests/utilitiesTestsGlobalFixture.h
	+++ b/Utilities/tests/utilitiesTestsGlobalFixture.h
	@@ -1,29 +1,27 @@
	// -- C++ --
	//
	// utilitiesTestGlobalFixture.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration, 2015 Marco A. Harrendorf
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#include <boost/test/unit_test.hpp>

	#include "ThePEG/Repository/StandardRandom.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Config/Unitsystem.h"

	-#include <iostream>
	-
	struct FixGlobal1 {
	FixGlobal1() {
	BOOST_TEST_MESSAGE( "setup global fixture for utilitiesTest" );

	// Initialize randomNumberGenerator
	ThePEG::StandardRandom* randomNumberStandardGenerator = new ThePEG::StandardRandom();
	new ThePEG::UseRandom(randomNumberStandardGenerator);
	}

	~FixGlobal1() { BOOST_TEST_MESSAGE( "teardown global fixture for utilitiesTest" ); }
	};

	BOOST_GLOBAL_FIXTURE(FixGlobal1)
	\ No newline at end of file
	diff --git a/Utilities/tests/utilitiesTestsMain.cc b/Utilities/tests/utilitiesTestsMain.cc
	--- a/Utilities/tests/utilitiesTestsMain.cc
	+++ b/Utilities/tests/utilitiesTestsMain.cc
	@@ -1,40 +1,39 @@
	// -- C++ --
	//
	// utilitiesTestMain.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration, 2015 Marco A. Harrendorf
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//

	/**
	* The following part should be included only once.
	*/
	#define BOOST_TEST_DYN_LINK
	#define BOOST_TEST_MAIN
	#include <boost/test/unit_test.hpp>
	#define BOOST_TEST_MODULE utilitiesTest

	/**
	* Include global fixture
	*
	* Global fixture initializes the randomNumber generator
	*/
	#include "Herwig++/Utilities/tests/utilitiesTestsGlobalFixture.h"

	/**
	* Include here the sub tests
	*/
	#include "Herwig++/Utilities/tests/utilitiesTestsKinematics.h"
	#include "Herwig++/Utilities/tests/utilitiesTestMaths.h"
	-#include "Herwig++/Utilities/tests/utilitiesTestSmearing.h"
	#include "Herwig++/Utilities/tests/utilitiesTestsStatistic.h"


	/**
	* Debug and development part
	*/
	BOOST_AUTO_TEST_CASE(fail)
	{
	//BOOST_FAIL("Ende");
	}
	\ No newline at end of file

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:36 PM (1 d, 12 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3802234
Default Alt Text: (31 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions