diff --git a/Shower/Core/Base/ShowerParticle.cc b/Shower/Core/Base/ShowerParticle.cc
--- a/Shower/Core/Base/ShowerParticle.cc
+++ b/Shower/Core/Base/ShowerParticle.cc
@@ -1,595 +1,584 @@
// -*- C++ -*-
//
// ShowerParticle.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 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 ShowerParticle class.
//
#include "ShowerParticle.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include <ostream>
using namespace Herwig;
PPtr ShowerParticle::clone() const {
return new_ptr(*this);
}
PPtr ShowerParticle::fullclone() const {
return new_ptr(*this);
}
ClassDescription<ShowerParticle> ShowerParticle::initShowerParticle;
// Definition of the static class description member.
void ShowerParticle::vetoEmission(ShowerPartnerType, Energy scale) {
scales_.QED = min(scale,scales_.QED );
scales_.QED_noAO = min(scale,scales_.QED_noAO );
scales_.QCD_c = min(scale,scales_.QCD_c );
scales_.QCD_c_noAO = min(scale,scales_.QCD_c_noAO );
scales_.QCD_ac = min(scale,scales_.QCD_ac );
scales_.QCD_ac_noAO = min(scale,scales_.QCD_ac_noAO);
}
void ShowerParticle::addPartner(EvolutionPartner in ) {
partners_.push_back(in);
}
namespace {
LorentzRotation boostToShower(Lorentz5Momentum & porig, tShowerBasisPtr basis) {
LorentzRotation output;
assert(basis);
if(basis->frame()==ShowerBasis::BackToBack) {
// we are doing the evolution in the back-to-back frame
// work out the boostvector
Boost boostv(-(basis->pVector()+basis->nVector()).boostVector());
// momentum of the parton
Lorentz5Momentum ptest(basis->pVector());
// construct the Lorentz boost
output = LorentzRotation(boostv);
ptest *= output;
Axis axis(ptest.vect().unit());
// now rotate so along the z axis as needed for the splitting functions
if(axis.perp2()>1e-10) {
double sinth(sqrt(1.-sqr(axis.z())));
output.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
}
else if(axis.z()<0.) {
output.rotate(Constants::pi,Axis(1.,0.,0.));
}
porig = output*basis->pVector();
porig.setX(ZERO);
porig.setY(ZERO);
}
else {
output = LorentzRotation(-basis->pVector().boostVector());
porig = output*basis->pVector();
porig.setX(ZERO);
porig.setY(ZERO);
porig.setZ(ZERO);
}
return output;
}
RhoDMatrix bosonMapping(ShowerParticle & particle,
const Lorentz5Momentum & porig,
VectorSpinPtr vspin,
const LorentzRotation & rot) {
// rotate the original basis
vector<LorentzPolarizationVector> sbasis;
for(unsigned int ix=0;ix<3;++ix) {
sbasis.push_back(vspin->getProductionBasisState(ix));
sbasis.back().transform(rot);
}
// splitting basis
vector<LorentzPolarizationVector> fbasis;
bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
VectorWaveFunction wave(porig,particle.dataPtr(),outgoing);
for(unsigned int ix=0;ix<3;++ix) {
if(massless&&ix==1) {
fbasis.push_back(LorentzPolarizationVector());
}
else {
wave.reset(ix);
fbasis.push_back(wave.wave());
}
}
// work out the mapping
RhoDMatrix mapping=RhoDMatrix(PDT::Spin1,false);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
mapping(ix,iy)= sbasis[iy].dot(fbasis[ix].conjugate());
if(particle.id()<0)
mapping(ix,iy)=conj(mapping(ix,iy));
}
}
- // \todo need to fix this
- mapping = RhoDMatrix(PDT::Spin1,false);
- if(massless) {
- mapping(0,0) = 1.;
- mapping(2,2) = 1.;
- }
- else {
- mapping(0,0) = 1.;
- mapping(1,1) = 1.;
- mapping(2,2) = 1.;
- }
return mapping;
}
RhoDMatrix fermionMapping(ShowerParticle & particle,
const Lorentz5Momentum & porig,
FermionSpinPtr fspin,
const LorentzRotation & rot) {
// extract the original basis states
vector<LorentzSpinor<SqrtEnergy> > sbasis;
for(unsigned int ix=0;ix<2;++ix) {
sbasis.push_back(fspin->getProductionBasisState(ix));
sbasis.back().transform(rot);
}
// calculate the states in the splitting basis
vector<LorentzSpinor<SqrtEnergy> > fbasis;
SpinorWaveFunction wave(porig,particle.dataPtr(),
particle.id()>0 ? incoming : outgoing);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
fbasis.push_back(wave.dimensionedWave());
}
RhoDMatrix mapping=RhoDMatrix(PDT::Spin1Half,false);
for(unsigned int ix=0;ix<2;++ix) {
if(fbasis[0].s2()==complex<SqrtEnergy>()) {
mapping(ix,0) = sbasis[ix].s3()/fbasis[0].s3();
mapping(ix,1) = sbasis[ix].s2()/fbasis[1].s2();
}
else {
mapping(ix,0) = sbasis[ix].s2()/fbasis[0].s2();
mapping(ix,1) = sbasis[ix].s3()/fbasis[1].s3();
}
}
return mapping;
}
FermionSpinPtr createFermionSpinInfo(ShowerParticle & particle,
const Lorentz5Momentum & porig,
const LorentzRotation & rot,
Helicity::Direction dir) {
// calculate the splitting basis for the branching
// and rotate back to construct the basis states
LorentzRotation rinv = rot.inverse();
SpinorWaveFunction wave;
if(particle.id()>0)
wave=SpinorWaveFunction(porig,particle.dataPtr(),incoming);
else
wave=SpinorWaveFunction(porig,particle.dataPtr(),outgoing);
FermionSpinPtr fspin = new_ptr(FermionSpinInfo(particle.momentum(),dir==outgoing));
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
LorentzSpinor<SqrtEnergy> basis = wave.dimensionedWave();
basis.transform(rinv);
fspin->setBasisState(ix,basis);
fspin->setDecayState(ix,basis);
}
particle.spinInfo(fspin);
return fspin;
}
VectorSpinPtr createVectorSpinInfo(ShowerParticle & particle,
const Lorentz5Momentum & porig,
const LorentzRotation & rot,
Helicity::Direction dir) {
// calculate the splitting basis for the branching
// and rotate back to construct the basis states
LorentzRotation rinv = rot.inverse();
bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
VectorWaveFunction wave(porig,particle.dataPtr(),dir);
VectorSpinPtr vspin = new_ptr(VectorSpinInfo(particle.momentum(),dir==outgoing));
for(unsigned int ix=0;ix<3;++ix) {
LorentzPolarizationVector basis;
if(massless&&ix==1) {
basis = LorentzPolarizationVector();
}
else {
wave.reset(ix);
basis = wave.wave();
}
basis *= rinv;
vspin->setBasisState(ix,basis);
vspin->setDecayState(ix,basis);
}
particle.spinInfo(vspin);
vspin-> DMatrix() = RhoDMatrix(PDT::Spin1);
vspin->rhoMatrix() = RhoDMatrix(PDT::Spin1);
if(massless) {
vspin-> DMatrix()(0,0) = 0.5;
vspin->rhoMatrix()(0,0) = 0.5;
vspin-> DMatrix()(2,2) = 0.5;
vspin->rhoMatrix()(2,2) = 0.5;
}
return vspin;
}
}
RhoDMatrix ShowerParticle::extractRhoMatrix(bool forward) {
// get the spin density matrix and the mapping
RhoDMatrix mapping;
SpinPtr inspin;
bool needMapping = getMapping(inspin,mapping);
// set the decayed flag
inspin->decay();
// get the spin density matrix
RhoDMatrix rho = forward ? inspin->rhoMatrix() : inspin->DMatrix();
// map to the shower basis if needed
if(needMapping) {
RhoDMatrix rhop(rho.iSpin(),false);
for(int ixb=0;ixb<rho.iSpin();++ixb) {
for(int iyb=0;iyb<rho.iSpin();++iyb) {
if(mapping(iyb,ixb)==0.)continue;
for(int iya=0;iya<rho.iSpin();++iya) {
if(rho(iya,iyb)==0.)continue;
for(int ixa=0;ixa<rho.iSpin();++ixa) {
rhop(ixa,ixb) += rho(iya,iyb)*mapping(iya,ixa)*conj(mapping(iyb,ixb));
}
}
}
}
rhop.normalize();
rho = rhop;
}
return rho;
}
bool ShowerParticle::getMapping(SpinPtr & output, RhoDMatrix & mapping) {
// if the particle is not from the hard process
if(!this->perturbative()) {
// mapping is the identity
output=this->spinInfo();
mapping=RhoDMatrix(this->dataPtr()->iSpin());
if(output) {
return false;
}
else {
Lorentz5Momentum porig;
LorentzRotation rot = boostToShower(porig,showerBasis());
Helicity::Direction dir = this->isFinalState() ? outgoing : incoming;
if(this->dataPtr()->iSpin()==PDT::Spin0) {
assert(false);
}
else if(this->dataPtr()->iSpin()==PDT::Spin1Half) {
output = createFermionSpinInfo(*this,porig,rot,dir);
}
else if(this->dataPtr()->iSpin()==PDT::Spin1) {
output = createVectorSpinInfo(*this,porig,rot,dir);
}
else {
assert(false);
}
return false;
}
}
// if particle is final-state and is from the hard process
else if(this->isFinalState()) {
assert(this->perturbative()==1 || this->perturbative()==2);
// get transform to shower frame
Lorentz5Momentum porig;
LorentzRotation rot = boostToShower(porig,showerBasis());
// the rest depends on the spin of the particle
PDT::Spin spin(this->dataPtr()->iSpin());
mapping=RhoDMatrix(spin,false);
// do the spin dependent bit
if(spin==PDT::Spin0) {
ScalarSpinPtr sspin=dynamic_ptr_cast<ScalarSpinPtr>(this->spinInfo());
if(!sspin) {
ScalarWaveFunction::constructSpinInfo(this,outgoing,true);
}
output=this->spinInfo();
return false;
}
else if(spin==PDT::Spin1Half) {
FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
// spin info exists get information from it
if(fspin) {
output=fspin;
mapping = fermionMapping(*this,porig,fspin,rot);
return true;
}
// spin info does not exist create it
else {
output = createFermionSpinInfo(*this,porig,rot,outgoing);
return false;
}
}
else if(spin==PDT::Spin1) {
VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
// spin info exists get information from it
if(vspin) {
output=vspin;
mapping = bosonMapping(*this,porig,vspin,rot);
return true;
}
else {
output = createVectorSpinInfo(*this,porig,rot,outgoing);
return false;
}
}
// not scalar/fermion/vector
else
assert(false);
}
// incoming to hard process
else if(this->perturbative()==1 && !this->isFinalState()) {
// get the basis vectors
// get transform to shower frame
Lorentz5Momentum porig;
LorentzRotation rot = boostToShower(porig,showerBasis());
porig *= this->x();
// the rest depends on the spin of the particle
PDT::Spin spin(this->dataPtr()->iSpin());
mapping=RhoDMatrix(spin);
// do the spin dependent bit
if(spin==PDT::Spin0) {
cerr << "testing spin 0 not yet implemented " << endl;
assert(false);
}
// spin-1/2
else if(spin==PDT::Spin1Half) {
FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
// spin info exists get information from it
if(fspin) {
output=fspin;
mapping = fermionMapping(*this,porig,fspin,rot);
return true;
}
// spin info does not exist create it
else {
output = createFermionSpinInfo(*this,porig,rot,incoming);
return false;
}
}
// spin-1
else if(spin==PDT::Spin1) {
VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
// spinInfo exists map it
if(vspin) {
output=vspin;
mapping = bosonMapping(*this,porig,vspin,rot);
return true;
}
// create the spininfo
else {
output = createVectorSpinInfo(*this,porig,rot,incoming);
return false;
}
}
assert(false);
}
// incoming to decay
else if(this->perturbative() == 2 && !this->isFinalState()) {
// get the basis vectors
Lorentz5Momentum porig;
LorentzRotation rot=boostToShower(porig,showerBasis());
// the rest depends on the spin of the particle
PDT::Spin spin(this->dataPtr()->iSpin());
mapping=RhoDMatrix(spin);
// do the spin dependent bit
if(spin==PDT::Spin0) {
cerr << "testing spin 0 not yet implemented " << endl;
assert(false);
}
// spin-1/2
else if(spin==PDT::Spin1Half) {
// FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
// // spin info exists get information from it
// if(fspin) {
// output=fspin;
// mapping = fermionMapping(*this,porig,fspin,rot);
// return true;
// // spin info does not exist create it
// else {
// output = createFermionSpinInfo(*this,porig,rot,incoming);
// return false;
// }
// }
assert(false);
}
// // spin-1
// else if(spin==PDT::Spin1) {
// VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
// // spinInfo exists map it
// if(vspin) {
// output=vspin;
// mapping = bosonMapping(*this,porig,vspin,rot);
// return true;
// }
// // create the spininfo
// else {
// output = createVectorSpinInfo(*this,porig,rot,incoming);
// return false;
// }
// }
// assert(false);
assert(false);
}
else
assert(false);
return true;
}
void ShowerParticle::constructSpinInfo(bool timeLike) {
// now construct the required spininfo and calculate the basis states
PDT::Spin spin(dataPtr()->iSpin());
if(spin==PDT::Spin0) {
ScalarWaveFunction::constructSpinInfo(this,outgoing,timeLike);
}
// calculate the basis states and construct the SpinInfo for a spin-1/2 particle
else if(spin==PDT::Spin1Half) {
// outgoing particle
if(id()>0) {
vector<LorentzSpinorBar<SqrtEnergy> > stemp;
SpinorBarWaveFunction::calculateWaveFunctions(stemp,this,outgoing);
SpinorBarWaveFunction::constructSpinInfo(stemp,this,outgoing,timeLike);
}
// outgoing antiparticle
else {
vector<LorentzSpinor<SqrtEnergy> > stemp;
SpinorWaveFunction::calculateWaveFunctions(stemp,this,outgoing);
SpinorWaveFunction::constructSpinInfo(stemp,this,outgoing,timeLike);
}
}
// calculate the basis states and construct the SpinInfo for a spin-1 particle
else if(spin==PDT::Spin1) {
bool massless(id()==ParticleID::g||id()==ParticleID::gamma);
vector<Helicity::LorentzPolarizationVector> vtemp;
VectorWaveFunction::calculateWaveFunctions(vtemp,this,outgoing,massless);
VectorWaveFunction::constructSpinInfo(vtemp,this,outgoing,timeLike,massless);
}
else {
throw Exception() << "Spins higher than 1 are not yet implemented in "
<< "FS_QtildaShowerKinematics1to2::constructVertex() "
<< Exception::runerror;
}
}
void ShowerParticle::initializeDecay() {
Lorentz5Momentum p, n, ppartner, pcm;
assert(perturbative()!=1);
// this is for the initial decaying particle
if(perturbative()==2) {
ShowerBasisPtr newBasis;
p = momentum();
Lorentz5Momentum ppartner(partner()->momentum());
// removed to make inverse recon work properly
//if(partner->thePEGBase()) ppartner=partner->thePEGBase()->momentum();
pcm=ppartner;
Boost boost(p.findBoostToCM());
pcm.boost(boost);
n = Lorentz5Momentum( ZERO,0.5*p.mass()*pcm.vect().unit());
n.boost( -boost);
newBasis = new_ptr(ShowerBasis());
newBasis->setBasis(p,n,ShowerBasis::Rest);
showerBasis(newBasis,false);
}
else {
showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0])->showerBasis(),true);
}
}
void ShowerParticle::initializeInitialState(PPtr parent) {
// For the time being we are considering only 1->2 branching
Lorentz5Momentum p, n, pthis, pcm;
assert(perturbative()!=2);
if(perturbative()==1) {
ShowerBasisPtr newBasis;
// find the partner and its momentum
if(!partner()) return;
if(partner()->isFinalState()) {
Lorentz5Momentum pa = -momentum()+partner()->momentum();
Lorentz5Momentum pb = momentum();
Energy scale=parent->momentum().t();
Lorentz5Momentum pbasis(ZERO,parent->momentum().vect().unit()*scale);
Axis axis(pa.vect().unit());
LorentzRotation rot;
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
if(axis.perp2()>1e-20) {
rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
rot.rotateX(Constants::pi);
}
if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag());
pb *= rot;
if(pb.perp2()/GeV2>1e-20) {
Boost trans = -1./pb.e()*pb.vect();
trans.setZ(0.);
rot.boost(trans);
}
pbasis *=rot;
rot.invert();
n = rot*Lorentz5Momentum(ZERO,-pbasis.vect());
p = rot*Lorentz5Momentum(ZERO, pbasis.vect());
}
else {
pcm = parent->momentum();
p = Lorentz5Momentum(ZERO, pcm.vect());
n = Lorentz5Momentum(ZERO, -pcm.vect());
}
newBasis = new_ptr(ShowerBasis());
newBasis->setBasis(p,n,ShowerBasis::BackToBack);
showerBasis(newBasis,false);
}
else {
showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(children()[0])->showerBasis(),true);
}
}
void ShowerParticle::initializeFinalState() {
// set the basis vectors
Lorentz5Momentum p,n;
if(perturbative()!=0) {
ShowerBasisPtr newBasis;
// find the partner() and its momentum
if(!partner()) return;
Lorentz5Momentum ppartner(partner()->momentum());
// momentum of the emitting particle
p = momentum();
Lorentz5Momentum pcm;
// if the partner is a final-state particle then the reference
// vector is along the partner in the rest frame of the pair
if(partner()->isFinalState()) {
Boost boost=(p + ppartner).findBoostToCM();
pcm = ppartner;
pcm.boost(boost);
n = Lorentz5Momentum(ZERO,pcm.vect());
n.boost( -boost);
}
else if(!partner()->isFinalState()) {
// if the partner is an initial-state particle then the reference
// vector is along the partner which should be massless
if(perturbative()==1) {
n = Lorentz5Momentum(ZERO,ppartner.vect());
}
// if the partner is an initial-state decaying particle then the reference
// vector is along the backwards direction in rest frame of decaying particle
else {
Boost boost=ppartner.findBoostToCM();
pcm = p;
pcm.boost(boost);
n = Lorentz5Momentum( ZERO, -pcm.vect());
n.boost( -boost);
}
}
newBasis = new_ptr(ShowerBasis());
newBasis->setBasis(p,n,ShowerBasis::BackToBack);
showerBasis(newBasis,false);
}
else if(initiatesTLS()) {
showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0]->children()[0])->showerBasis(),true);
}
else {
showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0])->showerBasis(),true);
}
}
void ShowerParticle::setShowerMomentum(bool timeLike) {
Energy m = this->mass() > ZERO ? this->mass() : this->data().mass();
// calculate the momentum of the assuming on-shell
Energy2 pt2 = sqr(this->showerParameters().pt);
double alpha = timeLike ? this->showerParameters().alpha : this->x();
double beta = 0.5*(sqr(m) + pt2 - sqr(alpha)*showerBasis()->pVector().m2())/(alpha*showerBasis()->p_dot_n());
Lorentz5Momentum porig=showerBasis()->sudakov2Momentum(alpha,beta,
this->showerParameters().ptx,
this->showerParameters().pty);
porig.setMass(m);
this->set5Momentum(porig);
}
diff --git a/Shower/Core/Base/ShowerParticle.h b/Shower/Core/Base/ShowerParticle.h
--- a/Shower/Core/Base/ShowerParticle.h
+++ b/Shower/Core/Base/ShowerParticle.h
@@ -1,530 +1,528 @@
// -*- C++ -*-
//
// ShowerParticle.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 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_ShowerParticle_H
#define HERWIG_ShowerParticle_H
//
// This is the declaration of the ShowerParticle class.
//
#include "ThePEG/EventRecord/Particle.h"
#include "Herwig/Shower/Core/SplittingFunctions/SplittingFunction.fh"
#include "Herwig/Shower/Core/ShowerConfig.h"
#include "ShowerBasis.h"
#include "ShowerKinematics.h"
#include "ShowerParticle.fh"
#include <iosfwd>
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
* This class represents a particle in the showering process.
* It inherits from the Particle class of ThePEG and has some
* specifics information useful only during the showering process.
*
* Notice that:
* - for forward evolution, it is clear what is meant by parent/child;
* for backward evolution, however, it depends whether we want
* to keep a physical picture or a Monte-Carlo effective one.
* In the former case, an incoming particle (emitting particle)
* splits into an emitted particle and the emitting particle after
* the emission: the latter two are then children of the
* emitting particle, the parent. In the Monte-Carlo effective
* picture, we have that the particle close to the hard subprocess,
* with higher (space-like) virtuality, splits into an emitted particle
* and the emitting particle at lower virtuality: the latter two are,
* in this case, the children of the first one, the parent. However we
* choose a more physical picture where the new emitting particle is the
* parented of the emitted final-state particle and the original emitting
* particle.
* - the pointer to a SplitFun object is set only in the case
* that the particle has undergone a shower emission. This is similar to
* the case of the decay of a normal Particle where
* the pointer to a Decayer object is set only in the case
* that the particle has undergone to a decay.
* In the case of particle connected directly to the hard subprocess,
* there is no pointer to the hard subprocess, but there is a method
* isFromHardSubprocess() which returns true only in this case.
*
* @see Particle
* @see ShowerConfig
* @see ShowerKinematics
*/
class ShowerParticle: public Particle {
public:
/**
* Struct for all the info on an evolution partner
*/
struct EvolutionPartner {
/**
* Constructor
*/
EvolutionPartner(tShowerParticlePtr p,double w, ShowerPartnerType t,
Energy s) : partner(p), weight(w), type(t), scale(s)
{}
/**
* The partner
*/
tShowerParticlePtr partner;
/**
* Weight
*/
double weight;
/**
* Type
*/
ShowerPartnerType type;
/**
* The assoicated evolution scale
*/
Energy scale;
};
/**
* Struct to store the evolution scales
*/
struct EvolutionScales {
/**
* Constructor
*/
EvolutionScales() : QED(),QCD_c(),QCD_ac(),
QED_noAO(),QCD_c_noAO(),QCD_ac_noAO(),
Max_Q2(Constants::MaxEnergy2)
{}
/**
* QED scale
*/
Energy QED;
/**
* QCD colour scale
*/
Energy QCD_c;
/**
* QCD anticolour scale
*/
Energy QCD_ac;
/**
* QED scale
*/
Energy QED_noAO;
/**
* QCD colour scale
*/
Energy QCD_c_noAO;
/**
* QCD anticolour scale
*/
Energy QCD_ac_noAO;
/**
* Maximum allowed virtuality of the particle
*/
Energy2 Max_Q2;
};
/** @name Construction and descruction functions. */
//@{
/**
* Standard Constructor. Note that the default constructor is
* private - there is no particle without a pointer to a
* ParticleData object.
* @param x the ParticleData object
* @param fs Whether or not the particle is an inital or final-state particle
* @param tls Whether or not the particle initiates a time-like shower
*/
ShowerParticle(tcEventPDPtr x, bool fs, bool tls=false)
: Particle(x), _isFinalState(fs),
_perturbative(0), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
_vMass(ZERO), _thePEGBase() {}
/**
* Copy constructor from a ThePEG Particle
* @param x ThePEG particle
* @param pert Where the particle came from
* @param fs Whether or not the particle is an inital or final-state particle
* @param tls Whether or not the particle initiates a time-like shower
*/
ShowerParticle(const Particle & x, unsigned int pert, bool fs, bool tls=false)
: Particle(x), _isFinalState(fs),
_perturbative(pert), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
_vMass(ZERO), _thePEGBase(&x) {}
//@}
public:
/**
* Set a preliminary momentum for the particle
*/
void setShowerMomentum(bool timelike);
/**
* Construct the spin info object for a shower particle
*/
void constructSpinInfo(bool timelike);
/**
* Perform any initial calculations needed after the branching has been selected
*/
void initializeDecay();
/**
* Perform any initial calculations needed after the branching has been selected
* @param parent The beam particle
*/
void initializeInitialState(PPtr parent);
/**
* Perform any initial calculations needed after the branching has been selected
*/
void initializeFinalState();
/**
* Access/Set various flags about the state of the particle
*/
//@{
/**
* Access the flag that tells if the particle is final state
* or initial state.
*/
bool isFinalState() const { return _isFinalState; }
/**
* Access the flag that tells if the particle is initiating a
* time like shower when it has been emitted in an initial state shower.
*/
bool initiatesTLS() const { return _initiatesTLS; }
/**
* Access the flag which tells us where the particle came from
* This is 0 for a particle produced in the shower, 1 if the particle came
* from the hard sub-process and 2 is it came from a decay.
*/
unsigned int perturbative() const { return _perturbative; }
//@}
/**
* Set/Get the momentum fraction for initial-state particles
*/
//@{
/**
* For an initial state particle get the fraction of the beam momentum
*/
void x(double x) { _x = x; }
/**
* For an initial state particle set the fraction of the beam momentum
*/
double x() const { return _x; }
//@}
/**
* Set/Get methods for the ShowerKinematics objects
*/
//@{
/**
* Access/ the ShowerKinematics object.
*/
const ShoKinPtr & showerKinematics() const { return _showerKinematics; }
/**
* Set the ShowerKinematics object.
*/
void showerKinematics(const ShoKinPtr in) { _showerKinematics = in; }
//@}
/**
* Set/Get methods for the ShowerBasis objects
*/
//@{
/**
* Access/ the ShowerBasis object.
*/
const ShowerBasisPtr & showerBasis() const { return _showerBasis; }
/**
* Set the ShowerBasis object.
*/
void showerBasis(const ShowerBasisPtr in, bool copy) {
if(!copy)
_showerBasis = in;
else {
_showerBasis = new_ptr(ShowerBasis());
_showerBasis->setBasis(in->pVector(),in->nVector(),in->frame());
}
}
//@}
/**
* Members relating to the initial evolution scale and partner for the particle
*/
//@{
/**
* Veto emission at a given scale
*/
void vetoEmission(ShowerPartnerType type, Energy scale);
/**
* Access to the evolution scales
*/
const EvolutionScales & scales() const {return scales_;}
/**
* Access to the evolution scales
*/
EvolutionScales & scales() {return scales_;}
/**
* Return the virtual mass\f$
*/
Energy virtualMass() const { return _vMass; }
/**
* Set the virtual mass
*/
void virtualMass(Energy mass) { _vMass = mass; }
/**
* Return the partner
*/
tShowerParticlePtr partner() const { return _partner; }
/**
* Set the partner
*/
void partner(const tShowerParticlePtr partner) { _partner = partner; }
/**
* Get the possible partners
*/
vector<EvolutionPartner> & partners() { return partners_; }
/**
* Add a possible partners
*/
void addPartner(EvolutionPartner in );
/**
* Clear the evolution partners
*/
void clearPartners() { partners_.clear(); }
/**
* Return the progenitor of the shower
*/
tShowerParticlePtr progenitor() const { return _progenitor; }
/**
* Set the progenitor of the shower
*/
void progenitor(const tShowerParticlePtr progenitor) { _progenitor = progenitor; }
//@}
/**
* Members to store and provide access to variables for a specific
* shower evolution scheme
*/
//@{
struct Parameters {
Parameters() : alpha(1.), beta(), ptx(), pty(), pt() {}
double alpha;
double beta;
Energy ptx;
Energy pty;
Energy pt;
};
/**
* Set the vector containing dimensionless variables
*/
Parameters & showerParameters() { return _parameters; }
//@}
/**
* If this particle came from the hard process get a pointer to ThePEG particle
* it came from
*/
const tcPPtr thePEGBase() const { return _thePEGBase; }
public:
/**
* Extract the rho matrix including mapping needed in the shower
*/
RhoDMatrix extractRhoMatrix(bool forward);
-protected:
-
/**
* For a particle which came from the hard process get the spin density and
* the mapping required to the basis used in the Shower
* @param rho The \f$\rho\f$ matrix
* @param mapping The mapping
* @param showerkin The ShowerKinematics object
*/
bool getMapping(SpinPtr &, RhoDMatrix & map);
protected:
/**
* Standard clone function.
*/
virtual PPtr clone() const;
/**
* Standard clone function.
*/
virtual PPtr fullclone() const;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<ShowerParticle> initShowerParticle;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerParticle & operator=(const ShowerParticle &);
private:
/**
* Whether the particle is in the final or initial state
*/
bool _isFinalState;
/**
* Whether the particle came from
*/
unsigned int _perturbative;
/**
* Does a particle produced in the backward shower initiate a time-like shower
*/
bool _initiatesTLS;
/**
* Dimensionless parameters
*/
Parameters _parameters;
/**
* The beam energy fraction for particle's in the initial state
*/
double _x;
/**
* The shower kinematics for the particle
*/
ShoKinPtr _showerKinematics;
/**
* The shower basis for the particle
*/
ShowerBasisPtr _showerBasis;
/**
* Storage of the evolution scales
*/
EvolutionScales scales_;
/**
* Virtual mass
*/
Energy _vMass;
/**
* Partners
*/
tShowerParticlePtr _partner;
/**
* Pointer to ThePEG Particle this ShowerParticle was created from
*/
const tcPPtr _thePEGBase;
/**
* Progenitor
*/
tShowerParticlePtr _progenitor;
/**
* Partners
*/
vector<EvolutionPartner> partners_;
};
inline ostream & operator<<(ostream & os, const ShowerParticle::EvolutionScales & es) {
os << "Scales: QED=" << es.QED / GeV
<< " QCD_c=" << es.QCD_c / GeV
<< " QCD_ac=" << es.QCD_ac / GeV
<< " QED_noAO=" << es.QED_noAO / GeV
<< " QCD_c_noAO=" << es.QCD_c_noAO / GeV
<< " QCD_ac_noAO=" << es.QCD_ac_noAO / GeV
<< '\n';
return os;
}
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of ShowerParticle. */
template <>
struct BaseClassTrait<Herwig::ShowerParticle,1> {
/** Typedef of the first base class of ShowerParticle. */
typedef Particle NthBase;
};
/** This template specialization informs ThePEG about the name of
* the ShowerParticle class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::ShowerParticle>
: public ClassTraitsBase<Herwig::ShowerParticle> {
/** Return a platform-independent class name */
static string className() { return "Herwig::ShowerParticle"; }
/** Create a Event object. */
static TPtr create() { return TPtr::Create(Herwig::ShowerParticle(tcEventPDPtr(),true)); }
};
/** @endcond */
}
#endif /* HERWIG_ShowerParticle_H */
diff --git a/Shower/Core/Base/ShowerVertex.cc b/Shower/Core/Base/ShowerVertex.cc
--- a/Shower/Core/Base/ShowerVertex.cc
+++ b/Shower/Core/Base/ShowerVertex.cc
@@ -1,76 +1,95 @@
// -*- C++ -*-
//
// ShowerVertex.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 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 ShowerVertex class.
//
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/EventRecord/SpinInfo.h"
#include "ShowerVertex.h"
using namespace Herwig;
using namespace Herwig::Helicity;
using namespace ThePEG;
DescribeNoPIOClass<ShowerVertex,HelicityVertex>
describeShowerVertex ("Herwig::ShowerVertex","");
void ShowerVertex::Init() {
static ClassDocumentation<ShowerVertex> documentation
("The ShowerVertex class is the implementation of a "
"vertex for a shower for the Herwig spin correlation algorithm");
}
// method to get the rho matrix for a given outgoing particle
RhoDMatrix ShowerVertex::getRhoMatrix(int i, bool) const {
assert(matrixElement_->nOut()==2);
// calculate the incoming spin density matrix
RhoDMatrix input=incoming()[0]->rhoMatrix();
if(convertIn_) input = mapIncoming(input);
// get the rho matrices for the outgoing particles
vector<RhoDMatrix> rhoout;
for(unsigned int ix=0,N=outgoing().size();ix<N;++ix) {
if(int(ix)!=i)
rhoout.push_back(outgoing()[ix]->DMatrix());
}
// calculate the spin density matrix
return matrixElement_->calculateRhoMatrix(i,input,rhoout);
}
// method to get the D matrix for an incoming particle
RhoDMatrix ShowerVertex::getDMatrix(int) const {
assert(matrixElement_->nOut()==2);
// get the decay matrices for the outgoing particles
vector<RhoDMatrix> Dout;
for(unsigned int ix=0,N=outgoing().size();ix<N;++ix) {
Dout.push_back(outgoing()[ix]->DMatrix());
}
- // calculate the spin density matrix and return the answer
- return matrixElement_->calculateDMatrix(Dout);
+ // calculate the spin density matrix
+ RhoDMatrix rho = matrixElement_->calculateDMatrix(Dout);
+ // map if needed
+ if(convertIn_) {
+ RhoDMatrix rhop(rho.iSpin(),false);
+ for(int ixb=0;ixb<rho.iSpin();++ixb) {
+ for(int iyb=0;iyb<rho.iSpin();++iyb) {
+ if(inMatrix_(iyb,ixb)==0.)continue;
+ for(int iya=0;iya<rho.iSpin();++iya) {
+ if(rho(iya,iyb)==0.)continue;
+ for(int ixa=0;ixa<rho.iSpin();++ixa) {
+ rhop(ixa,ixb) += rho(iya,iyb)*std::conj(inMatrix_(iya,ixa))*inMatrix_(iyb,ixb);
+ }
+ }
+ }
+ }
+ rhop.normalize();
+ rho = rhop;
+ }
+ // return the answer
+ return rho;
}
RhoDMatrix ShowerVertex::mapIncoming(RhoDMatrix rho) const {
RhoDMatrix output(rho.iSpin());
for(int ixa=0;ixa<rho.iSpin();++ixa) {
for(int ixb=0;ixb<rho.iSpin();++ixb) {
for(int iya=0;iya<rho.iSpin();++iya) {
for(int iyb=0;iyb<rho.iSpin();++iyb) {
output(ixa,ixb) += rho(iya,iyb)*inMatrix_(iya,ixa)*conj(inMatrix_(iyb,ixb));
}
}
}
}
output.normalize();
return output;
}
diff --git a/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.cc b/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.cc
--- a/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.cc
+++ b/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.cc
@@ -1,200 +1,204 @@
// -*- C++ -*-
//
// FS_QTildeShowerKinematics1to2.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 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 FS_QTildeShowerKinematics1to2 class.
//
#include "FS_QTildeShowerKinematics1to2.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "Herwig/Shower/Core/SplittingFunctions/SplittingFunction.h"
#include "Herwig/Shower/Core/Base/ShowerParticle.h"
#include "ThePEG/Utilities/Debug.h"
#include "Herwig/Shower/QTilde/QTildeShowerHandler.h"
#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
#include "Herwig/Shower/QTilde/Base/ShowerModel.h"
#include "Herwig/Shower/QTilde/Base/KinematicsReconstructor.h"
#include "Herwig/Shower/Core/Base/ShowerVertex.h"
using namespace Herwig;
void FS_QTildeShowerKinematics1to2::
updateParameters(tShowerParticlePtr theParent,
tShowerParticlePtr theChild0,
tShowerParticlePtr theChild1,
bool setAlpha) const {
const ShowerParticle::Parameters & parent = theParent->showerParameters();
ShowerParticle::Parameters & child0 = theChild0->showerParameters();
ShowerParticle::Parameters & child1 = theChild1->showerParameters();
// determine alphas of children according to interpretation of z
if ( setAlpha ) {
child0.alpha = z() * parent.alpha;
child1.alpha = (1.-z()) * parent.alpha;
}
// set the values
double cphi = cos(phi());
double sphi = sin(phi());
child0.ptx = pT() * cphi + z() * parent.ptx;
child0.pty = pT() * sphi + z() * parent.pty;
child0.pt = sqrt( sqr(child0.ptx) + sqr(child0.pty) );
child1.ptx = -pT() * cphi + (1.-z())* parent.ptx;
child1.pty = -pT() * sphi + (1.-z())* parent.pty;
child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) );
}
void FS_QTildeShowerKinematics1to2::
updateChildren(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType partnerType,
bool massVeto) const {
assert(children.size()==2);
// calculate the scales
splittingFn()->evaluateFinalStateScales(partnerType,scale(),z(),parent,
children[0],children[1]);
// set the maximum virtual masses
if(massVeto) {
Energy2 q2 = z()*(1.-z())*sqr(scale());
IdList ids(3);
ids[0] = parent->dataPtr();
ids[1] = children[0]->dataPtr();
ids[2] = children[1]->dataPtr();
const vector<Energy> & virtualMasses = SudakovFormFactor()->virtualMasses(ids);
if(ids[0]->id()!=ParticleID::g && ids[0]->id()!=ParticleID::gamma ) {
q2 += sqr(virtualMasses[0]);
}
// limits on further evolution
children[0]->scales().Max_Q2 = z() *(q2-sqr(virtualMasses[2])/(1.-z()));
children[1]->scales().Max_Q2 = (1.-z())*(q2-sqr(virtualMasses[1])/ z() );
}
// update the parameters
updateParameters(parent, children[0], children[1], true);
// set up the colour connections
splittingFn()->colourConnection(parent,children[0],children[1],partnerType,false);
// make the products children of the parent
parent->addChild(children[0]);
parent->addChild(children[1]);
// set the momenta of the children
for(ShowerParticleVector::const_iterator pit=children.begin();
pit!=children.end();++pit) {
(**pit).showerBasis(parent->showerBasis(),true);
(**pit).setShowerMomentum(true);
}
// sort out the helicity stuff
if(! dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->correlations()) return;
SpinPtr pspin(parent->spinInfo());
if(!pspin || !dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->spinCorrelations() ) return;
Energy2 t = sqr(scale())*z()*(1.-z());
IdList ids;
ids.push_back(parent->dataPtr());
ids.push_back(children[0]->dataPtr());
ids.push_back(children[1]->dataPtr());
// create the vertex
SVertexPtr vertex(new_ptr(ShowerVertex()));
// set the matrix element
vertex->ME(splittingFn()->matrixElement(z(),t,ids,phi(),true));
+ RhoDMatrix mapping;
+ SpinPtr inspin;
+ bool needMapping = parent->getMapping(inspin,mapping);
+ if(needMapping) vertex->incomingBasisTransform(mapping);
// set the incoming particle for the vertex
parent->spinInfo()->decayVertex(vertex);
for(ShowerParticleVector::const_iterator pit=children.begin();
pit!=children.end();++pit) {
// construct the spin info for the children
(**pit).constructSpinInfo(true);
// connect the spinInfo object to the vertex
(*pit)->spinInfo()->productionVertex(vertex);
}
}
void FS_QTildeShowerKinematics1to2::
reconstructParent(const tShowerParticlePtr parent,
const ParticleVector & children ) const {
assert(children.size() == 2);
ShowerParticlePtr c1 = dynamic_ptr_cast<ShowerParticlePtr>(children[0]);
ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(children[1]);
parent->showerParameters().beta=
c1->showerParameters().beta + c2->showerParameters().beta;
Lorentz5Momentum pnew = c1->momentum() + c2->momentum();
Energy2 m2 = sqr(pT())/z()/(1.-z()) + sqr(c1->mass())/z()
+ sqr(c2->mass())/(1.-z());
pnew.setMass(sqrt(m2));
parent->set5Momentum( pnew );
}
void FS_QTildeShowerKinematics1to2::reconstructLast(const tShowerParticlePtr last,
Energy mass) const {
// set beta component and consequently all missing data from that,
// using the nominal (i.e. PDT) mass.
Energy theMass = mass > ZERO ? mass : last->data().constituentMass();
Lorentz5Momentum pVector = last->showerBasis()->pVector();
ShowerParticle::Parameters & lastParam = last->showerParameters();
Energy2 denom = 2. * lastParam.alpha * last->showerBasis()->p_dot_n();
if(abs(denom)/(sqr(pVector.e())+pVector.rho2())<1e-10) {
throw KinematicsReconstructionVeto();
}
lastParam.beta = ( sqr(theMass) + sqr(lastParam.pt)
- sqr(lastParam.alpha) * pVector.m2() )
/ denom;
// set that new momentum
Lorentz5Momentum newMomentum = last->showerBasis()->
sudakov2Momentum( lastParam.alpha, lastParam.beta,
lastParam.ptx , lastParam.pty);
newMomentum.setMass(theMass);
newMomentum.rescaleEnergy();
if(last->data().stable()) {
last->set5Momentum( newMomentum );
}
else {
last->boost(last->momentum().findBoostToCM());
last->boost(newMomentum.boostVector());
}
}
void FS_QTildeShowerKinematics1to2::updateParent(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType) const {
IdList ids(3);
ids[0] = parent->dataPtr();
ids[1] = children[0]->dataPtr();
ids[2] = children[1]->dataPtr();
const vector<Energy> & virtualMasses = SudakovFormFactor()->virtualMasses(ids);
if(children[0]->children().empty()) children[0]->virtualMass(virtualMasses[1]);
if(children[1]->children().empty()) children[1]->virtualMass(virtualMasses[2]);
// compute the new pT of the branching
Energy2 pt2=sqr(z()*(1.-z()))*sqr(scale())
- sqr(children[0]->virtualMass())*(1.-z())
- sqr(children[1]->virtualMass())* z() ;
if(ids[0]->id()!=ParticleID::g) pt2 += z()*(1.-z())*sqr(virtualMasses[0]);
if(pt2>ZERO) {
pT(sqrt(pt2));
}
else {
pt2=ZERO;
pT(ZERO);
}
Energy2 q2 =
sqr(children[0]->virtualMass())/z() +
sqr(children[1]->virtualMass())/(1.-z()) +
pt2/z()/(1.-z());
parent->virtualMass(sqrt(q2));
}
void FS_QTildeShowerKinematics1to2::
resetChildren(const tShowerParticlePtr parent,
const ShowerParticleVector & children) const {
updateParameters(parent, children[0], children[1], false);
for(unsigned int ix=0;ix<children.size();++ix) {
if(children[ix]->children().empty()) continue;
ShowerParticleVector newChildren;
for(unsigned int iy=0;iy<children[ix]->children().size();++iy)
newChildren.push_back(dynamic_ptr_cast<ShowerParticlePtr>
(children[ix]->children()[iy]));
children[ix]->showerKinematics()->resetChildren(children[ix],newChildren);
}
}