No OneTemporary
Actions

Size

34 KB

Subscribers

None

View Options

	diff --git a/Shower/Dipole/Utility/ConstituentReshuffler.cc b/Shower/Dipole/Utility/ConstituentReshuffler.cc
	--- a/Shower/Dipole/Utility/ConstituentReshuffler.cc
	+++ b/Shower/Dipole/Utility/ConstituentReshuffler.cc
	@@ -1,663 +1,622 @@
	// -- C++ --
	//
	// ConstituentReshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 ConstituentReshuffler class.
	//

	#include <config.h>
	#include "ConstituentReshuffler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"

	#include <limits>

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "DipolePartonSplitter.h"

	#include "Herwig/Utilities/GSLBisection.h"

	#include "Herwig/Shower/Dipole/DipoleShowerHandler.h"

	#include "Herwig/Shower/ShowerHandler.h"

	using namespace Herwig;

	ConstituentReshuffler::ConstituentReshuffler()
	: HandlerBase() {}

	ConstituentReshuffler::~ConstituentReshuffler() {}

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

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

	-double ConstituentReshuffler::ReshuffleEquation::aUnit() {
	- return 1.;
	-}
	-
	-double ConstituentReshuffler::ReshuffleEquation::vUnit() {
	- return 1.;
	-}
	-
	-double ConstituentReshuffler::DecayReshuffleEquation::aUnit() {
	- return 1.;
	-}
	-
	-double ConstituentReshuffler::DecayReshuffleEquation::vUnit() {
	- return 1.;
	-}
	-
	-
	-double ConstituentReshuffler::ReshuffleEquation::operator() (double xi) const {
	-
	- double r = - w/GeV;
	-
	- for (PList::iterator p = p_begin; p != p_end; ++p) {
	- r += sqrt(sqr((**p).dataPtr()->constituentMass()) +
	- xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))) / GeV;
	- }
	- return r;
	-
	-}
	-
	-double ConstituentReshuffler::DecayReshuffleEquation::operator() (double xi) const {
	- double r = - w/GeV;
	-
	- for (PList::iterator pIt = p_begin; pIt != p_end; ++pIt) {
	- r += sqrt(sqr((**pIt).dataPtr()->constituentMass()) +
	- xixi(sqr((pIt).momentum().t())-sqr((pIt).dataPtr()->mass()))) / GeV;
	- }
	-
	- for (PList::iterator rIt = r_begin; rIt != r_end; ++rIt) {
	- r += sqrt(sqr((**rIt).momentum().m()) +
	- xixi(sqr((rIt).momentum().t())-sqr((rIt).momentum().m()))) / GeV;
	- }
	-
	- return r;
	-
	-}
	-
	void ConstituentReshuffler::reshuffle(PList& out,
	PPair& in,
	PList& intermediates,
	const bool decay,
	PList& decayPartons,
	PList& decayRecoilers) {

	assert(ShowerHandler::currentHandler()->retConstituentMasses());

	if ( !decay ) {

	if (out.size() == 0)
	return;

	if (out.size() == 1) {

	PPtr recoiler;
	PPtr parton = out.front();

	if (DipolePartonSplitter::colourConnected(parton,in.first) &&
	DipolePartonSplitter::colourConnected(parton,in.second)) {
	if (UseRandom::rnd() < .5)
	recoiler = in.first;
	else
	recoiler = in.second;
	} else if (DipolePartonSplitter::colourConnected(parton,in.first)) {
	recoiler = in.first;
	} else if (DipolePartonSplitter::colourConnected(parton,in.second)) {
	recoiler = in.second;
	} else assert(false);

	assert(abs(recoiler->momentum().vect().perp2()/GeV2) < 1e-6);

	double sign = recoiler->momentum().z() < 0.*GeV ? -1. : 1.;

	Energy2 qperp2 = parton->momentum().perp2();

	if (qperp2/GeV2 < Constants::epsilon) {
	// no emission off a 2 -> singlet process which
	// needed a single forced splitting: should never happen (?)
	assert(false);
	throw Veto();
	}

	Energy2 m2 = sqr(parton->dataPtr()->constituentMass());

	Energy abs_q = parton->momentum().vect().mag();
	Energy qz = parton->momentum().z();
	Energy abs_pz = recoiler->momentum().t();
	assert(abs_pz > 0.*GeV);

	Energy xi_pz = sign(2.qperp2abs_pz + m2(abs_q + signqz))/(2.qperp2);
	Energy x_qz = (2.qperp2qz + m2(qz+signabs_q))/(2.*qperp2);

	Lorentz5Momentum recoiler_momentum
	(0.GeV,0.GeV,xi_pz,xi_pz < 0.*GeV ? - xi_pz : xi_pz);

	recoiler_momentum.rescaleMass();

	Lorentz5Momentum parton_momentum
	(parton->momentum().x(),parton->momentum().y(),x_qz,sqrt(m2+qperp2+x_qz*x_qz));

	parton_momentum.rescaleMass();

	PPtr n_parton = new_ptr(Particle(parton->dataPtr()));
	n_parton->set5Momentum(parton_momentum);

	DipolePartonSplitter::change(parton,n_parton,false);

	out.pop_front();
	intermediates.push_back(parton);
	out.push_back(n_parton);

	PPtr n_recoiler = new_ptr(Particle(recoiler->dataPtr()));
	n_recoiler->set5Momentum(recoiler_momentum);

	DipolePartonSplitter::change(recoiler,n_recoiler,true);

	intermediates.push_back(recoiler);

	if (recoiler == in.first) {
	in.first = n_recoiler;
	}

	if (recoiler == in.second) {
	in.second = n_recoiler;
	}

	return;

	}

	}

	Energy zero (0.*GeV);
	Lorentz5Momentum Q (zero,zero,zero,zero);

	for (PList::iterator p = out.begin();
	p != out.end(); ++p) {
	Q += (**p).momentum();
	}

	Boost beta = Q.findBoostToCM();

	list<Lorentz5Momentum> mbackup;

	bool need_boost = (beta.mag2() > Constants::epsilon);

	if (need_boost) {

	for (PList::iterator p = out.begin();
	p != out.end(); ++p) {
	Lorentz5Momentum mom = (**p).momentum();
	mbackup.push_back(mom);
	(**p).set5Momentum(mom.boost(beta));
	}

	}

	double xi;

	// Only partons
	if ( decayRecoilers.size()==0 ) {
	- ReshuffleEquation solve (Q.m(),out.begin(),out.end());
	+ list<Energy> masses;
	+ for ( auto p : out )
	+ masses.push_back(p->dataPtr()->constituentMass());
	+ ReshuffleEquation<PList::iterator,list<Energy>::const_iterator>
	+ solve (Q.m(),out.begin(),out.end(),
	+ masses.begin(),masses.end());

	GSLBisection solver(1e-10,1e-8,10000);

	try {
	xi = solver.value(solve,0.0,1.1);
	} catch (GSLBisection::GSLerror) {
	throw DipoleShowerHandler::RedoShower();
	} catch (GSLBisection::IntervalError) {
	throw DipoleShowerHandler::RedoShower();
	}
	}

	// Partons and decaying recoilers
	else {
	DecayReshuffleEquation solve (Q.m(),decayPartons.begin(),decayPartons.end(),decayRecoilers.begin(),decayRecoilers.end());

	GSLBisection solver(1e-10,1e-8,10000);

	try {
	xi = solver.value(solve,0.0,1.1);
	} catch (GSLBisection::GSLerror) {
	throw DipoleShowerHandler::RedoShower();
	} catch (GSLBisection::IntervalError) {
	throw DipoleShowerHandler::RedoShower();
	}
	}


	PList reshuffled;

	list<Lorentz5Momentum>::const_iterator backup_it;
	if (need_boost)
	backup_it = mbackup.begin();


	// Reshuffling of non-decaying partons only
	if ( decayRecoilers.size()==0 ) {

	for (PList::iterator p = out.begin();
	p != out.end(); ++p) {

	PPtr rp = new_ptr(Particle((**p).dataPtr()));

	DipolePartonSplitter::change(*p,rp,false);

	Lorentz5Momentum rm;

	rm = Lorentz5Momentum (xi(*p).momentum().x(),
	xi(*p).momentum().y(),
	xi(*p).momentum().z(),
	sqrt(sqr((**p).dataPtr()->constituentMass()) +
	xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))));

	rm.rescaleMass();

	if (need_boost) {
	(*p).set5Momentum(backup_it);
	++backup_it;
	rm.boost(-beta);
	}

	rp->set5Momentum(rm);

	intermediates.push_back(*p);
	reshuffled.push_back(rp);


	}
	}

	// For the case of a decay process with non-partonic recoilers
	else {
	assert ( decay );

	for (PList::iterator p = out.begin();
	p != out.end(); ++p) {

	// Flag to update spinInfo
	bool updateSpin = false;

	PPtr rp = new_ptr(Particle((**p).dataPtr()));

	DipolePartonSplitter::change(*p,rp,false);

	Lorentz5Momentum rm;

	// If the particle is a parton and not a recoiler
	if ( find( decayRecoilers.begin(), decayRecoilers.end(), *p ) == decayRecoilers.end() ) {
	rm = Lorentz5Momentum (xi(*p).momentum().x(),
	xi(*p).momentum().y(),
	xi(*p).momentum().z(),
	sqrt(sqr((**p).dataPtr()->constituentMass()) +
	xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))));
	}


	// Otherwise the parton is a recoiler
	// and its invariant mass must be preserved
	else {
	if ( (*p)-> spinInfo() )
	updateSpin = true;
	rm = Lorentz5Momentum (xi(*p).momentum().x(),
	xi(*p).momentum().y(),
	xi(*p).momentum().z(),
	sqrt(sqr((**p).momentum().m()) +
	xixi(sqr((p).momentum().t())-sqr((p).momentum().m()))));
	}

	rm.rescaleMass();

	if (need_boost) {
	(*p).set5Momentum(backup_it);
	++backup_it;
	rm.boost(-beta);
	}

	rp->set5Momentum(rm);

	// Update SpinInfo if required
	if ( updateSpin )
	updateSpinInfo(*p, rp);

	intermediates.push_back(*p);
	reshuffled.push_back(rp);


	}
	}

	out.clear();
	out.splice(out.end(),reshuffled);

	}


	void ConstituentReshuffler::hardProcDecayReshuffle(PList& decaying,
	PList& eventOutgoing,
	PList& eventHard,
	PPair& eventIncoming,
	PList& eventIntermediates) {

	// Note, when this function is called, the particle pointers
	// in theDecays/decaying are those prior to the showering.
	// Here we find the newest pointers in the outgoing.
	// The update of the PPtrs in theDecays is done in DipoleShowerHandler::constituentReshuffle()
	// as this needs to be done if ConstituentReshuffling is switched off.


	//Make sure the shower should return constituent masses:
	assert(ShowerHandler::currentHandler()->retConstituentMasses());

	// Find the outgoing decaying particles
	PList recoilers;
	for ( PList::iterator decIt = decaying.begin(); decIt != decaying.end(); ++decIt) {

	// First find the particles in the intermediates
	PList::iterator pos = find(eventIntermediates.begin(),eventIntermediates.end(), *decIt);

	// Colourless particle or coloured particle that did not radiate.
	if(pos==eventIntermediates.end()) {

	// Check that this is not a particle from a subsequent decay.
	// e.g. the W from a top decay from an LHE file.
	if ( find( eventHard.begin(), eventHard.end(), *decIt ) == eventHard.end() &&
	find( eventOutgoing.begin(), eventOutgoing.end(), *decIt ) == eventOutgoing.end() )
	continue;

	else
	recoilers.push_back( *decIt );

	}

	// Coloured decaying particle that radiated
	else {
	PPtr unstable = *pos;
	while(!unstable->children().empty()) {
	unstable = unstable->children()[0];
	}
	assert( find( eventOutgoing.begin(),eventOutgoing.end(), unstable ) != eventOutgoing.end() );

	recoilers.push_back( unstable );
	}
	}

	// Make a list of partons
	PList partons;
	for ( PList::iterator outPos = eventOutgoing.begin(); outPos != eventOutgoing.end(); ++outPos ) {
	if ( find (recoilers.begin(), recoilers.end(), *outPos ) == recoilers.end() ) {
	partons.push_back( *outPos );
	}
	}


	// If no outgoing partons, do nothing
	if ( partons.size() == 0 ){
	return;
	}

	// Otherwise reshuffling needs to be done.

	// If there is only one parton, attempt to reshuffle with
	// the incoming to be consistent with the reshuffle for a
	// hard process with no decays.
	else if ( partons.size() == 1 &&
	( DipolePartonSplitter::colourConnected(partons.front(),eventIncoming.first) \|\|
	DipolePartonSplitter::colourConnected(partons.front(),eventIncoming.second) ) ) {

	// Erase the parton from the event outgoing
	eventOutgoing.erase( find( eventOutgoing.begin(), eventOutgoing.end(), partons.front() ) );
	// Perform the reshuffle, this update the intermediates and the incoming
	reshuffle(partons, eventIncoming, eventIntermediates);
	// Update the outgoing
	eventOutgoing.push_back(partons.front());
	return;
	}

	// If reshuffling amongst the incoming is not possible
	// or if we have multiple outgoing partons.
	else {

	// Create a complete list of the outgoing from the process
	PList out;
	// Make an empty list for storing the new intermediates
	PList intermediates;
	// Empty incoming particles pair
	PPair in;

	// A single parton which cannot be reshuffled
	// with the incoming.
	if ( partons.size() == 1 ) {

	// Populate the out for the reshuffling
	out.insert(out.end(),partons.begin(),partons.end());
	out.insert(out.end(),recoilers.begin(),recoilers.end());
	assert( out.size() > 1 );

	// Perform the reshuffle with the temporary particle lists
	reshuffle(out, in, intermediates, true, partons, recoilers);
	}

	// If there is more than one parton, reshuffle only
	// amongst the partons
	else {
	assert(partons.size() > 1);

	// Populate the out for the reshuffling
	out.insert(out.end(),partons.begin(),partons.end());
	assert( out.size() > 1 );

	// Perform the reshuffle with the temporary particle lists
	reshuffle(out, in, intermediates, true);
	}

	// Update the dipole event record
	updateEvent(intermediates, eventIntermediates, out, eventOutgoing, eventHard );
	return;
	}
	}


	void ConstituentReshuffler::decayReshuffle(PerturbativeProcessPtr& decayProc,
	PList& eventOutgoing,
	PList& eventHard,
	PList& eventIntermediates ) {

	// Separate particles into those to be assigned constituent masses
	// i.e. non-decaying coloured partons
	// and those which must only absorb recoil
	// i.e. non-coloured and decaying particles
	PList partons;
	PList recoilers;

	//Make sure the shower should return constituent masses:
	assert(ShowerHandler::currentHandler()->retConstituentMasses());


	// Populate the particle lists from the outgoing of the decay process
	for( unsigned int ix = 0; ix<decayProc->outgoing().size(); ++ix) {

	// Identify recoilers
	if ( !decayProc->outgoing()[ix].first->coloured() \|\|
	ShowerHandler::currentHandler()->decaysInShower(decayProc->outgoing()[ix].first->id() ) )
	recoilers.push_back(decayProc->outgoing()[ix].first);

	else
	partons.push_back(decayProc->outgoing()[ix].first);
	}

	// If there are no outgoing partons, then no reshuffling
	// needs to be done
	if ( partons.size() == 0 )
	return;

	// Reshuffling needs to be done:
	else {

	// Create a complete list of the outgoing from the process
	PList out;
	// Make an empty list for storing the new intermediates
	PList intermediates;
	// Empty incoming particles pair
	PPair in;


	// SW - 15/06/2018, 31/01/2019 - Always include 'recoilers' in
	// reshuffling, regardless of the number of partons to be put on their
	// constituent mass shell. This is because reshuffling between 2 partons
	// frequently leads to a redoShower exception. This treatment is
	// consistent with the AO shower

	// Populate the out for the reshuffling
	out.insert(out.end(),partons.begin(),partons.end());
	out.insert(out.end(),recoilers.begin(),recoilers.end());
	assert( out.size() > 1 );

	// Perform the reshuffle with the temporary particle lists
	reshuffle(out, in, intermediates, true, partons, recoilers);

	// Update the dipole event record and the decay process
	updateEvent(intermediates, eventIntermediates, out, eventOutgoing, eventHard, decayProc );
	return;
	}
	}


	void ConstituentReshuffler::updateEvent( PList& intermediates,
	PList& eventIntermediates,
	#ifndef NDEBUG
	PList& out,
	#else
	PList&,
	#endif
	PList& eventOutgoing,
	PList& eventHard,
	PerturbativeProcessPtr decayProc ) {

	// Loop over the new intermediates following the reshuffling
	for (PList::iterator p = intermediates.begin();
	p != intermediates.end(); ++p) {

	// Update the event record intermediates
	eventIntermediates.push_back(*p);

	// Identify the reshuffled particle
	assert( (*p)->children().size()==1 );
	PPtr reshuffled = (*p)->children()[0];
	assert( find(out.begin(), out.end(), reshuffled) != out.end() );

	// Update the event record outgoing
	PList::iterator posOut = find(eventOutgoing.begin(), eventOutgoing.end(), *p);

	if ( posOut != eventOutgoing.end() ) {
	eventOutgoing.erase(posOut);
	eventOutgoing.push_back(reshuffled);
	}

	else {
	PList::iterator posHard = find(eventHard.begin(), eventHard.end(), *p);
	assert( posHard != eventHard.end() );
	eventHard.erase(posHard);
	eventHard.push_back(reshuffled);
	}

	// Replace the particle in the the decay process outgoing
	if ( decayProc ) {
	vector<pair<PPtr,PerturbativeProcessPtr> >::iterator decayOutIt = decayProc->outgoing().end();
	for ( decayOutIt = decayProc->outgoing().begin();
	decayOutIt!= decayProc->outgoing().end(); ++decayOutIt ) {

	if ( decayOutIt->first == *p ){
	break;
	}
	}
	assert( decayOutIt != decayProc->outgoing().end() );
	decayOutIt->first = reshuffled;
	}
	}
	}

	void ConstituentReshuffler::updateSpinInfo( PPtr& oldPart,
	PPtr& newPart ) {

	const Lorentz5Momentum& oldMom = oldPart->momentum();
	const Lorentz5Momentum& newMom = newPart->momentum();

	// Rotation from old momentum to +ve z-axis
	LorentzRotation oldToZAxis;
	Axis axisOld(oldMom.vect().unit());
	if( axisOld.perp2() > 1e-12 ) {
	double sinth(sqrt(1.-sqr(axisOld.z())));
	oldToZAxis.rotate( -acos(axisOld.z()),Axis(-axisOld.y()/sinth,axisOld.x()/sinth,0.));
	}

	// Rotation from new momentum to +ve z-axis
	LorentzRotation newToZAxis;
	Axis axisNew(newMom.vect().unit());
	if( axisNew.perp2() > 1e-12 ) {
	double sinth(sqrt(1.-sqr(axisNew.z())));
	newToZAxis.rotate( -acos(axisNew.z()),Axis(-axisNew.y()/sinth,axisNew.x()/sinth,0.));
	}

	// Boost from old momentum to new momentum along z-axis
	Lorentz5Momentum momOldRotated = oldToZAxis*Lorentz5Momentum(oldMom);
	Lorentz5Momentum momNewRotated = newToZAxis*Lorentz5Momentum(newMom);

	Energy2 a = sqr(momOldRotated.z()) + sqr(momNewRotated.t());
	Energy2 b = 2.momOldRotated.t()momOldRotated.z();
	Energy2 c = sqr(momOldRotated.t()) - sqr(momNewRotated.t());
	double beta;

	// The rotated momentum should always lie along the +ve z-axis
	if ( momOldRotated.z() > ZERO )
	beta = (-b + sqrt(sqr(b)-4.ac)) / 2. / a;
	else
	beta = (-b - sqrt(sqr(b)-4.ac)) / 2. / a;

	LorentzRotation boostOldToNew(0., 0., beta);

	// Total transform
	LorentzRotation transform = (newToZAxis.inverse())boostOldToNewoldToZAxis;

	// Assign the same spin info to the old and new particles
	newPart->spinInfo(oldPart->spinInfo());
	newPart->spinInfo()->transform(oldMom, transform);
	}

	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void ConstituentReshuffler::persistentOutput(PersistentOStream &) const {
	}

	void ConstituentReshuffler::persistentInput(PersistentIStream &, int) {
	}

	ClassDescription<ConstituentReshuffler> ConstituentReshuffler::initConstituentReshuffler;
	// Definition of the static class description member.

	void ConstituentReshuffler::Init() {

	static ClassDocumentation<ConstituentReshuffler> documentation
	("The ConstituentReshuffler class implements reshuffling "
	"of partons on their nominal mass shell to their constituent "
	"mass shells.");

	}

	diff --git a/Shower/Dipole/Utility/ConstituentReshuffler.h b/Shower/Dipole/Utility/ConstituentReshuffler.h
	--- a/Shower/Dipole/Utility/ConstituentReshuffler.h
	+++ b/Shower/Dipole/Utility/ConstituentReshuffler.h
	@@ -1,299 +1,274 @@
	// -- C++ --
	//
	// ConstituentReshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ConstituentReshuffler_H
	#define HERWIG_ConstituentReshuffler_H
	//
	// This is the declaration of the ConstituentReshuffler class.
	//

	#include "ThePEG/Handlers/HandlerBase.h"
	#include "ThePEG/Utilities/Exception.h"
	#include "Herwig/Shower/PerturbativeProcess.h"
	+#include "Herwig/Utilities/Reshuffler.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup DipoleShower
	* \author Simon Platzer, Stephen Webster
	*
	* \brief The ConstituentReshuffler class implements reshuffling
	* of partons on their nominal mass shell to their constituent
	* mass shells.
	*
	*/
	-class ConstituentReshuffler: public HandlerBase {
	+class ConstituentReshuffler: public HandlerBase, public Reshuffler {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	ConstituentReshuffler();

	/**
	* The destructor.
	*/
	virtual ~ConstituentReshuffler();
	//@}

	public:

	/**
	* Reshuffle the outgoing partons to constituent
	* masses. Optionally, incoming partons are given
	* to absorb recoils. Add the non-reshuffled partons
	* to the intermediates list. Throw ConstituentReshufflerProblem
	* if a numerical problem prevents the solution of
	* the reshuffling equation.
	*/
	void reshuffle(PList& out,
	PPair& in,
	PList& intermediates,
	const bool decay,
	PList& decayPartons,
	PList& decayRecoilers);

	/**
	* Reshuffle the outgoing partons to constituent
	* masses. Optionally, incoming partons are given
	* to absorb recoils. Add the non-reshuffled partons
	* to the intermediates list. Throw ConstituentReshufflerProblem
	* if a numerical problem prevents the solution of
	* the reshuffling equation.
	*/
	void reshuffle(PList& out,
	PPair& in,
	PList& intermediates,
	const bool decay=false) {

	PList decayPartons;
	PList decayRecoilers;

	reshuffle(out,
	in,
	intermediates,
	decay,
	decayPartons,
	decayRecoilers);
	}


	/**
	* Reshuffle the outgoing partons following the showering
	* of the initial hard interaction to constituent masses,
	* for the case of outgoing decaying particles.
	* Throw ConstituentReshufflerProblem
	* if a numerical problem prevents the solution of
	* the reshuffling equation.
	*/
	void hardProcDecayReshuffle(PList& decaying,
	PList& eventOutgoing,
	PList& eventHard,
	PPair& eventIncoming,
	PList& eventIntermediates) ;

	/**
	* Reshuffle the outgoing partons following the showering
	* of a particle decay to constituent masses.
	* Throw ConstituentReshufflerProblem
	* if a numerical problem prevents the solution of
	* the reshuffling equation.
	*/
	void decayReshuffle(PerturbativeProcessPtr& decayProc,
	PList& eventOutgoing,
	PList& eventHard,
	PList& eventIntermediates) ;

	/**
	* Update the dipole event record and, if appropriate,
	* the relevant decay process.
	**/
	void updateEvent( PList& intermediates,
	PList& eventIntermediates,
	PList& out,
	PList& eventOutgoing,
	PList& eventHard,
	PerturbativeProcessPtr decayProc = PerturbativeProcessPtr() ) ;


	/**
	* Update the spinInfo of a particle following reshuffling
	* to take account of the change in momentum.
	* Used only for unstable particles that need to be dealt with.
	**/
	void updateSpinInfo( PPtr& oldPart,
	PPtr& newPart ) ;

	protected:

	/**
	* The function object defining the equation
	- * to be solved.
	- */
	- struct ReshuffleEquation {
	-
	- ReshuffleEquation (Energy q,
	- PList::iterator m_begin,
	- PList::iterator m_end)
	- : w(q), p_begin(m_begin), p_end(m_end) {}
	-
	- typedef double ArgType;
	- typedef double ValType;
	-
	- static double aUnit();
	- static double vUnit();
	-
	- double operator() (double xi) const;
	-
	- Energy w;
	-
	- PList::iterator p_begin;
	- PList::iterator p_end;
	-
	- };
	-
	- /**
	- * The function object defining the equation
	* to be solved in the case of separate recoilers
	* TODO - refine the whole implementation of separate partons and recoilers
	*/
	struct DecayReshuffleEquation {

	DecayReshuffleEquation (Energy q,
	PList::iterator m_begin,
	PList::iterator m_end,
	PList::iterator n_begin,
	PList::iterator n_end)
	: w(q), p_begin(m_begin), p_end(m_end), r_begin(n_begin), r_end(n_end) {}

	typedef double ArgType;
	typedef double ValType;

	static double aUnit();
	static double vUnit();

	double operator() (double xi) const;

	Energy w;

	PList::iterator p_begin;
	PList::iterator p_end;

	PList::iterator r_begin;
	PList::iterator r_end;
	};



	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	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;
	//@}


	// If needed, insert declarations of virtual function defined in the
	// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<ConstituentReshuffler> initConstituentReshuffler;

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	ConstituentReshuffler & operator=(const ConstituentReshuffler &) = delete;

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of ConstituentReshuffler. */
	template <>
	struct BaseClassTrait<Herwig::ConstituentReshuffler,1> {
	/** Typedef of the first base class of ConstituentReshuffler. */
	typedef HandlerBase NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the ConstituentReshuffler class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::ConstituentReshuffler>
	: public ClassTraitsBase<Herwig::ConstituentReshuffler> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::ConstituentReshuffler"; }
	/**
	* The name of a file containing the dynamic library where the class
	* ConstituentReshuffler is implemented. It may also include several, space-separated,
	* libraries if the class ConstituentReshuffler depends on other classes (base classes
	* excepted). In this case the listed libraries will be dynamically
	* linked in the order they are specified.
	*/
	static string library() { return "HwDipoleShower.so"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_ConstituentReshuffler_H */
	diff --git a/Utilities/Makefile.am b/Utilities/Makefile.am
	--- a/Utilities/Makefile.am
	+++ b/Utilities/Makefile.am
	@@ -1,46 +1,47 @@
	SUBDIRS = XML Statistics
	noinst_LTLIBRARIES = libHwUtils.la

	libHwUtils_la_SOURCES = \
	EnumParticles.h \
	Interpolator.tcc Interpolator.h \
	Kinematics.cc Kinematics.h \
	Progress.h Progress.cc \
	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 \
	+Reshuffler.h Reshuffler.cc \
	expm-1.h \
	HiggsLoopFunctions.h AlphaS.h

	nodist_libHwUtils_la_SOURCES = hgstamp.inc
	BUILT_SOURCES = hgstamp.inc
	CLEANFILES = hgstamp.inc

	HGVERSION := $(shell hg -R $(top_srcdir) parents --template '"Herwig {node\|short} ({branch})"' 2> /dev/null \|\| echo \"$(PACKAGE_STRING)\" \|\| true )

	.PHONY: update_hgstamp
	hgstamp.inc: update_hgstamp
	@[ -f $@ ] \|\| touch $@
	@echo '$(HGVERSION)' \| cmp -s $@ - \|\| echo '$(HGVERSION)' > $@

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


	check_PROGRAMS = utilities_test
	utilities_test_SOURCES = \
	tests/utilitiesTestsMain.cc \
	tests/utilitiesTestsGlobalFixture.h \
	tests/utilitiesTestsKinematics.h \
	tests/utilitiesTestMaths.h \
	tests/utilitiesTestsStatistic.h
	utilities_test_LDADD = $(BOOST_UNIT_TEST_FRAMEWORK_LIBS) $(THEPEGLIB) -ldl libHwUtils.la
	utilities_test_LDFLAGS = $(AM_LDFLAGS) -export-dynamic $(BOOST_UNIT_TEST_FRAMEWORK_LDFLAGS) $(THEPEGLDFLAGS)
	utilities_test_CPPFLAGS = $(AM_CPPFLAGS) $(BOOST_CPPFLAGS)
	TESTS = utilities_test
	diff --git a/Utilities/Reshuffler.cc b/Utilities/Reshuffler.cc
	new file mode 100644
	--- /dev/null
	+++ b/Utilities/Reshuffler.cc
	@@ -0,0 +1,85 @@
	+// -- C++ --
	+//
	+// Reshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	+// Copyright (C) 2002-2019 The Herwig Collaboration
	+//
	+// Herwig is licenced under version 3 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 Reshuffler class.
	+//
	+
	+#include <config.h>
	+#include "Reshuffler.h"
	+#include "Herwig/Utilities/GSLBisection.h"
	+#include "Herwig/Shower/ShowerHandler.h"
	+
	+using namespace Herwig;
	+
	+void Reshuffler::reshuffle(const PVector& particles,
	+ const vector<Energy>& masses) const {
	+
	+ Energy zero (0.*GeV);
	+ Lorentz5Momentum Q (zero,zero,zero,zero);
	+
	+ for (PVector::const_iterator p = particles.begin();
	+ p != particles.end(); ++p) {
	+ Q += (**p).momentum();
	+ }
	+
	+ Boost beta = Q.findBoostToCM();
	+
	+ vector<Lorentz5Momentum> mbackup;
	+
	+ bool need_boost = (beta.mag2() > Constants::epsilon);
	+
	+ if (need_boost) {
	+
	+ for (PVector::const_iterator p = particles.begin();
	+ p != particles.end(); ++p) {
	+ Lorentz5Momentum mom = (**p).momentum();
	+ mbackup.push_back(mom);
	+ (**p).set5Momentum(mom.boost(beta));
	+ }
	+
	+ }
	+
	+ double xi;
	+
	+ ReshuffleEquation<PVector::const_iterator,vector<Energy>::const_iterator>
	+ solve (Q.m(),particles.begin(),particles.end(),
	+ masses.begin(),masses.end());
	+
	+ GSLBisection solver(1e-10,1e-8,10000);
	+
	+ try {
	+ xi = solver.value(solve,0.0,1.1);
	+ } catch (GSLBisection::GSLerror) {
	+ throw Exception("Failed to reshuffle.",Exception::eventerror);
	+ } catch (GSLBisection::IntervalError) {
	+ throw Exception("Failed to reshuffle.",Exception::eventerror);
	+ }
	+
	+ PVector::const_iterator p = particles.begin();
	+ vector<Energy>::const_iterator m = masses.begin();
	+ for (; p != particles.end(); ++p, ++m) {
	+
	+ Lorentz5Momentum rm = Lorentz5Momentum (xi(*p).momentum().x(),
	+ xi(*p).momentum().y(),
	+ xi(*p).momentum().z(),
	+ sqrt(sqr(*m) +
	+ xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))));
	+
	+ rm.rescaleMass();
	+
	+ if (need_boost) {
	+ rm.boost(-beta);
	+ }
	+
	+ (**p).set5Momentum(rm);
	+
	+ }
	+
	+}
	diff --git a/Utilities/Reshuffler.h b/Utilities/Reshuffler.h
	new file mode 100644
	--- /dev/null
	+++ b/Utilities/Reshuffler.h
	@@ -0,0 +1,86 @@
	+// -- C++ --
	+//
	+// Reshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	+// Copyright (C) 2002-2019 The Herwig Collaboration
	+//
	+// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	+// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	+//
	+#ifndef HERWIG_Reshuffler_H
	+#define HERWIG_Reshuffler_H
	+//
	+// This is the declaration of the Reshuffler class.
	+//
	+
	+#include "ThePEG/Config/ThePEG.h"
	+
	+namespace Herwig {
	+
	+using namespace ThePEG;
	+
	+/**
	+ * \author Simon Platzer, Stephen Webster
	+ *
	+ * \brief The Reshuffler class implements reshuffling
	+ * of partons on their nominal mass shell to their constituent
	+ * mass shells.
	+ *
	+ */
	+class Reshuffler {
	+
	+protected:
	+
	+ /**
	+ * The function object defining the equation
	+ * to be solved.
	+ */
	+ template<class PIterator, class MIterator>
	+ struct ReshuffleEquation {
	+
	+ ReshuffleEquation (Energy q,
	+ PIterator pp_begin,
	+ PIterator pp_end,
	+ MIterator mm_begin,
	+ MIterator mm_end)
	+ : w(q), p_begin(pp_begin), p_end(pp_end),
	+ m_begin(mm_begin), m_end(mm_end) {}
	+
	+ typedef double ArgType;
	+ typedef double ValType;
	+
	+ static double aUnit() { return 1.; }
	+ static double vUnit() { return 1.; }
	+
	+ double operator() (double xi) const {
	+ double r = - w/GeV;
	+ PIterator p = p_begin;
	+ MIterator m = m_begin;
	+ for (; p != p_end; ++p, ++m) {
	+ r += sqrt(sqr(*m) +
	+ xixi(sqr((p).momentum().t())-sqr((p).dataPtr()->mass()))) / GeV;
	+ }
	+ return r;
	+ }
	+
	+ Energy w;
	+
	+ PIterator p_begin;
	+ PIterator p_end;
	+
	+ MIterator m_begin;
	+ MIterator m_end;
	+
	+ };
	+
	+ /**
	+ * Reshuffle to consitutent masses
	+ */
	+ void reshuffle(const PVector& particles,
	+ const vector<Energy>& masses) const;
	+
	+};
	+
	+}
	+
	+#endif /* HERWIG_Reshuffler_H */
	+

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:59 PM (1 d, 7 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3799670
Default Alt Text: (34 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions