Page MenuHomeHEPForge
Files F8308671

No OneTemporary
Actions

View Options

diff --git a/Shower/QTilde/Base/PartnerFinder.cc b/Shower/QTilde/Base/PartnerFinder.cc
--- a/Shower/QTilde/Base/PartnerFinder.cc
+++ b/Shower/QTilde/Base/PartnerFinder.cc
@@ -1,710 +1,710 @@
// -*- C++ -*-
//
// PartnerFinder.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 PartnerFinder class.
//
#include "PartnerFinder.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
-#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
DescribeClass<PartnerFinder,Interfaced>
describePartnerFinder ("Herwig::PartnerFinder","HwShower.so");
// some useful functions to avoid using #define
namespace {
// return bool if final-state particle
inline bool FS(const tShowerParticlePtr a) {
return a->isFinalState();
}
- // return colour line pointer
- inline Ptr<ThePEG::ColourLine>::transient_pointer
+ // return colour line pointer
+ inline Ptr<ThePEG::ColourLine>::transient_pointer
CL(const tShowerParticlePtr a, unsigned int index=0) {
return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() :
const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->colourLines()[index]);
}
// return colour line size
inline size_t
CLSIZE(const tShowerParticlePtr a) {
return a->colourInfo()->colourLines().size();
}
-
+
inline Ptr<ThePEG::ColourLine>::transient_pointer
ACL(const tShowerParticlePtr a, unsigned int index=0) {
return a->colourInfo()->antiColourLines().empty() ? ThePEG::tColinePtr() :
const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->antiColourLines()[index]);
}
inline size_t
ACLSIZE(const tShowerParticlePtr a) {
return a->colourInfo()->antiColourLines().size();
}
}
void PartnerFinder::persistentOutput(PersistentOStream & os) const {
os << partnerMethod_ << QEDPartner_ << scaleChoice_;
}
void PartnerFinder::persistentInput(PersistentIStream & is, int) {
is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_;
}
void PartnerFinder::Init() {
static ClassDocumentation<PartnerFinder> documentation
("This class is responsible for finding the partners for each interaction types ",
"and within the evolution scale range specified by the ShowerVariables ",
"then to determine the initial evolution scales for each pair of partners.");
static Switch<PartnerFinder,int> interfacePartnerMethod
("PartnerMethod",
"Choice of partner finding method for gluon evolution.",
&PartnerFinder::partnerMethod_, 0, false, false);
static SwitchOption interfacePartnerMethodRandom
(interfacePartnerMethod,
"Random",
"Choose partners of a gluon randomly.",
0);
static SwitchOption interfacePartnerMethodMaximum
(interfacePartnerMethod,
"Maximum",
"Choose partner of gluon with largest angle.",
1);
static Switch<PartnerFinder,int> interfaceQEDPartner
("QEDPartner",
"Control of which particles to use as the partner for QED radiation",
&PartnerFinder::QEDPartner_, 0, false, false);
static SwitchOption interfaceQEDPartnerAll
(interfaceQEDPartner,
"All",
"Consider all possible choices which give a positive contribution"
" in the soft limit.",
0);
static SwitchOption interfaceQEDPartnerIIandFF
(interfaceQEDPartner,
"IIandFF",
"Only allow initial-initial or final-final combinations",
1);
static SwitchOption interfaceQEDPartnerIF
(interfaceQEDPartner,
"IF",
"Only allow initial-final combinations",
2);
static Switch<PartnerFinder,int> interfaceScaleChoice
("ScaleChoice",
"The choice of the evolution scales",
&PartnerFinder::scaleChoice_, 0, false, false);
static SwitchOption interfaceScaleChoicePartner
(interfaceScaleChoice,
"Partner",
"Scale of all interactions is that of the evolution partner",
0);
static SwitchOption interfaceScaleChoiceDifferent
(interfaceScaleChoice,
"Different",
"Allow each interaction to have different scales",
1);
}
void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
ShowerInteraction type,
const bool setPartners) {
// clear the existing partners
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) (*cit)->clearPartners();
// set them
if(type==ShowerInteraction::QCD) {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
}
else if(type==ShowerInteraction::QED) {
setInitialQEDEvolutionScales(particles,isDecayCase,setPartners);
}
else if(type==ShowerInteraction::EW) {
setInitialEWEvolutionScales(particles,isDecayCase,false);
}
else if(type==ShowerInteraction::QEDQCD) {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
setInitialQEDEvolutionScales(particles,isDecayCase,false);
}
else if(type==ShowerInteraction::ALL) {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
setInitialQEDEvolutionScales(particles,isDecayCase,false);
setInitialEWEvolutionScales(particles,isDecayCase,false);
}
else
assert(false);
// \todo EW scales here
// print out for debugging
if(Debug::level>=10) {
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
generator()->log() << "Particle: " << **cit << "\n";
if(!(**cit).partner()) continue;
generator()->log() << "Primary partner: " << *(**cit).partner() << "\n";
for(vector<ShowerParticle::EvolutionPartner>::const_iterator it= (**cit).partners().begin();
it!=(**cit).partners().end();++it) {
generator()->log() << static_cast<long>(it->type) << " "
- << it->weight << " "
- << it->scale/GeV << " "
- << *(it->partner)
+ << it->weight << " "
+ << it->scale/GeV << " "
+ << *(it->partner)
<< "\n";
}
}
generator()->log() << flush;
}
}
void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// Loop over particles and consider only coloured particles which don't
// have already their colour partner fixed and that don't have children
- // (the latter requirement is relaxed in the case isDecayCase is true).
- // Build a map which has as key one of these particles (i.e. a pointer
+ // (the latter requirement is relaxed in the case isDecayCase is true).
+ // Build a map which has as key one of these particles (i.e. a pointer
// to a ShowerParticle object) and as a corresponding value the vector
// of all its possible *normal* candidate colour partners, defined as follows:
- // --- have colour, and no children (this is not required in the case
+ // --- have colour, and no children (this is not required in the case
// isDecayCase is true);
- // --- if both are initial (incoming) state particles, then the (non-null) colourLine()
+ // --- if both are initial (incoming) state particles, then the (non-null) colourLine()
// of one of them must match the (non-null) antiColourLine() of the other.
// --- if one is an initial (incoming) state particle and the other is
// a final (outgoing) state particle, then both must have the
// same (non-null) colourLine() or the same (non-null) antiColourLine();
// Notice that this definition exclude the special case of baryon-violating
// processes (as in R-parity Susy), which will show up as particles
// without candidate colour partners, and that we will be treated a part later
// (this means that no modifications of the following loop is needed!
for ( const auto & sp : particles ) {
// Skip colourless particles
if(!sp->data().coloured()) continue;
// find the partners
auto partners = findQCDPartners(sp,particles);
// must have a partner
if(partners.empty()) {
- throw Exception() << "`Failed to make colour connections in "
+ throw Exception() << "`Failed to make colour connections in "
<< "PartnerFinder::setQCDInitialEvolutionScales"
<< *sp
<< Exception::eventerror;
}
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
int position = -1;
for(size_t ix=0; ix< partners.size(); ++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp, partners[ix].second),isDecayCase));
- if (!setPartners && partners[ix].second) position = ix;
+ if (!setPartners && partners[ix].second) position = ix;
}
assert(setPartners || position >= 0);
-
+
// set partners if required
if (setPartners) {
// In the case of more than one candidate colour partners,
// there are now two approaches to choosing the partner. The
// first method is based on two assumptions:
// 1) the choice of which is the colour partner is done
// *randomly* between the available candidates.
// 2) the choice of which is the colour partner is done
// *independently* from each particle: in other words,
- // if for a particle "i" its selected colour partner is
- // the particle "j", then the colour partner of "j"
+ // if for a particle "i" its selected colour partner is
+ // the particle "j", then the colour partner of "j"
// does not have to be necessarily "i".
// The second method always chooses the furthest partner
// for hard gluons and gluinos.
// random choice
if( partnerMethod_ == 0 ) {
// random choice of partner
position = UseRandom::irnd(partners.size());
}
// take the one with largest angle
else if (partnerMethod_ == 1 ) {
- if (sp->perturbative() == 1 &&
+ if (sp->perturbative() == 1 &&
sp->dataPtr()->iColour()==PDT::Colour8 ) {
assert(partners.size()==2);
// Determine largest angle
double maxAngle(0.);
for(unsigned int ix=0;ix<partners.size();++ix) {
double angle = sp->momentum().vect().
angle(partners[ix].second->momentum().vect());
if(angle>maxAngle) {
maxAngle = angle;
position = ix;
}
}
}
else position = UseRandom::irnd(partners.size());
}
else assert(false);
// set the evolution partner
sp->partner(partners[position].second);
}
-
+
// primary partner set, set the others and do the scale
for(size_t ix=0; ix<partners.size(); ++ix) {
sp->addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,1.,partners[ix].first,
scales[ix].first));
}
-
+
// set scales for all interactions to that of the partner, default
Energy scale = scales[position].first;
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first==ShowerPartnerType::QCDColourLine) {
- sp->scales().QCD_c =
- sp->scales().QCD_c_noAO =
+ sp->scales().QCD_c =
+ sp->scales().QCD_c_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
- sp->scales().QCD_ac =
+ sp->scales().QCD_ac =
sp->scales().QCD_ac_noAO =
(scaleChoice_==0 ? scale : scales[ix].first);
}
else assert(false);
}
}
}
void PartnerFinder::setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// loop over all the particles
for(const auto & sp : particles) {
// not charged or photon continue
if(!sp->dataPtr()->charged()) continue;
// find the potential partners
vector<pair<double,tShowerParticlePtr> > partners = findQEDPartners(sp,particles,isDecayCase);
if(partners.empty()) {
throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setQEDInitialEvolutionScales"
<< *sp << Exception::eventerror;
}
// calculate the probabilities
double prob(0.);
for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
// normalise
for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
// set the partner if required
int position(-1);
// use QCD partner if set
if(!setPartners&&sp->partner()) {
for(unsigned int ix=0;ix<partners.size();++ix) {
if(sp->partner()==partners[ix].second) {
position = ix;
break;
}
}
}
// set the partner
if(setPartners||!sp->partner()||position<0) {
prob = UseRandom::rnd();
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first>prob) {
position = ix;
break;
}
prob -= partners[ix].first;
}
if(position>=0&&(setPartners||!sp->partner())) {
sp->partner(partners[position].second);
}
}
// must have a partner
if(position<0) throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setQEDInitialEvolutionScales"
- << *sp << Exception::eventerror;
+ << *sp << Exception::eventerror;
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
for(unsigned int ix=0;ix< partners.size();++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp,partners[ix].second),
isDecayCase));
}
// store all the possible partners
for(unsigned int ix=0;ix<partners.size();++ix) {
sp->addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
partners[ix].first,
ShowerPartnerType::QED,
scales[ix].first));
}
// set scales
sp->scales().QED = scales[position].first;
sp->scales().QED_noAO = scales[position].first;
}
}
pair<Energy,Energy> PartnerFinder::
calculateInitialEvolutionScales(const ShowerPPair &particlePair,
const bool isDecayCase, int key) {
bool FS1=FS(particlePair.first),FS2= FS(particlePair.second);
if(FS1 && FS2){
return calculateFinalFinalScales(particlePair.first->momentum(),particlePair.second->momentum(), key);
}
else if(FS1 && !FS2) {
pair<Energy,Energy> rval = calculateInitialFinalScales(particlePair.second->momentum(),
particlePair.first->momentum(),
isDecayCase);
return { rval.second, rval.first };
}
else if(!FS1 &&FS2)
return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase);
else
return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum());
}
-vector< pair<ShowerPartnerType, tShowerParticlePtr> >
+vector< pair<ShowerPartnerType, tShowerParticlePtr> >
PartnerFinder::findQCDPartners(tShowerParticlePtr particle,
const ShowerParticleVector &particles) {
vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners;
for(const auto & sp : particles) {
if(!sp->data().coloured() || particle==sp) continue;
// one initial-state and one final-state particle
if(FS(particle) != FS(sp)) {
// loop over all the colours of both particles
for(size_t ix=0; ix<CLSIZE(particle); ++ix) {
for(size_t jx=0; jx<CLSIZE(sp); ++jx) {
if((CL(particle,ix) && CL(particle,ix)==CL(sp,jx))) {
partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
}
}
}
//loop over all the anti-colours of both particles
for(size_t ix=0; ix<ACLSIZE(particle); ++ix) {
for(size_t jx=0; jx<ACLSIZE(sp); ++jx) {
if((ACL(particle,ix) && ACL(particle,ix)==ACL(sp,jx))) {
partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
}
}
}
}
// two initial-state or two final-state particles
else {
//loop over the colours of the first particle and the anti-colours of the other
for(size_t ix=0; ix<CLSIZE(particle); ++ix){
for(size_t jx=0; jx<ACLSIZE(sp); ++jx){
if(CL(particle,ix) && CL(particle,ix)==ACL(sp,jx)) {
partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
}
}
}
//loop over the anti-colours of the first particle and the colours of the other
for(size_t ix=0; ix<ACLSIZE(particle); ++ix){
for(size_t jx=0; jx<CLSIZE(sp); jx++){
if(ACL(particle,ix) && ACL(particle,ix)==CL(sp,jx)) {
partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
}
}
}
}
}
// if we haven't found any partners look for RPV
if (partners.empty()) {
// special for RPV
- tColinePtr col = CL(particle);
+ tColinePtr col = CL(particle);
if(FS(particle)&&col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))||
(!FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second ))) {
partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
}
}
}
else if(col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second))||
(!FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))) {
partners.push_back({ ShowerPartnerType:: QCDColourLine, sp });
}
}
}
col = ACL(particle);
if(FS(particle)&&col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second))||
(!FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second ))) {
partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
}
}
}
else if(col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(const auto & sp : particles) {
if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))||
(!FS(sp) && (ACL(sp) == cpair.first ||ACL(sp) == cpair.second))) {
partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
}
}
}
}
// return the partners
return partners;
}
vector< pair<double, tShowerParticlePtr> >
PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
const ShowerParticleVector & particles,
const bool isDecayCase) {
vector< pair<double, tShowerParticlePtr> > partners;
- const double pcharge =
+ const double pcharge =
particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge());
vector< pair<double, tShowerParticlePtr> > photons;
for (const auto & sp : particles) {
if (particle == sp) continue;
if (sp->id()==ParticleID::gamma) photons.push_back(make_pair(1.,sp));
if (!sp->data().charged() ) continue;
double charge = pcharge*double((sp)->data().iCharge());
if ( FS(particle) != FS(sp) ) charge *=-1.;
if ( QEDPartner_ != 0 && !isDecayCase ) {
// only include II and FF as requested
if ( QEDPartner_ == 1 && FS(particle) != FS(sp) )
continue;
// only include IF is requested
else if (QEDPartner_ == 2 && FS(particle) == FS(sp) )
continue;
}
if (particle->id()==ParticleID::gamma) charge = -abs(charge);
// only keep positive dipoles
if (charge<0.) partners.push_back({ -charge, sp });
}
if (particle->id()==ParticleID::gamma && partners.empty())
return photons;
return partners;
}
void PartnerFinder::setInitialEWEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// loop over all the particles
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
// if not weakly interacting continue
if( !weaklyInteracting( (**cit).dataPtr()))
continue;
// find the potential partners
vector<pair<double,tShowerParticlePtr> > partners = findEWPartners(*cit,particles,isDecayCase);
if(partners.empty()) {
throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setEWInitialEvolutionScales"
<< (**cit) << Exception::eventerror;
}
// calculate the probabilities
double prob(0.);
for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
// normalise
for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
// set the partner if required
int position(-1);
// use QCD partner if set
if(!setPartners&&(*cit)->partner()) {
for(unsigned int ix=0;ix<partners.size();++ix) {
if((*cit)->partner()==partners[ix].second) {
position = ix;
break;
}
}
}
// set the partner
if(setPartners||!(*cit)->partner()||position<0) {
prob = UseRandom::rnd();
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first>prob) {
position = ix;
break;
}
prob -= partners[ix].first;
}
if(position>=0&&(setPartners||!(*cit)->partner())) {
(*cit)->partner(partners[position].second);
}
}
// must have a partner
if(position<0) throw Exception() << "Failed to find partner in "
<< "PartnerFinder::setEWInitialEvolutionScales"
<< (**cit) << Exception::eventerror;
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
for(unsigned int ix=0;ix< partners.size();++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
isDecayCase, 0));
}
// store all the possible partners
for(unsigned int ix=0;ix<partners.size();++ix) {
(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
partners[ix].first,
ShowerPartnerType::EW,
scales[ix].first));
}
// set scales
(**cit).scales().EW = scales[position].first;
}
}
vector< pair<double, tShowerParticlePtr> >
PartnerFinder::findEWPartners(tShowerParticlePtr particle,
const ShowerParticleVector &particles,
const bool isDecayCase ) {
vector< pair<double, tShowerParticlePtr> > partners;
ShowerParticleVector::const_iterator cjt;
for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
if(!weaklyInteracting((*cjt)->dataPtr()) ||
particle == *cjt) continue;
if( QEDPartner_ != 0 && !isDecayCase ) {
// only include II and FF as requested
if( QEDPartner_ == 1 && FS(particle) != FS(*cjt) )
continue;
// only include IF is requested
else if(QEDPartner_ == 2 && FS(particle) == FS(*cjt) )
continue;
}
double charge = 1.;
// only keep positive dipoles
if(charge>0.) partners.push_back(make_pair(charge,*cjt));
}
return partners;
}
-pair<Energy,Energy>
+pair<Energy,Energy>
PartnerFinder::calculateFinalFinalScales(
const Lorentz5Momentum & p1,
- const Lorentz5Momentum & p2, int key)
+ const Lorentz5Momentum & p2, int key)
{
static const double eps=1e-7;
// Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda)
// ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2)
// find momenta in rest frame of system
// calculate quantities for the scales
Energy2 Q2 = (p1+p2).m2();
double b = p1.mass2()/Q2;
double c = p2.mass2()/Q2;
if(b<0.) {
if(b<-eps) {
throw Exception() << "Negative Mass squared b = " << b
<< "in PartnerFinder::calculateFinalFinalScales()"
<< Exception::eventerror;
}
b = 0.;
}
if(c<0.) {
if(c<-eps) {
throw Exception() << "Negative Mass squared c = " << c
<< "in PartnerFinder::calculateFinalFinalScales()"
<< Exception::eventerror;
}
c = 0.;
}
- // KMH & PR - 16 May 2008 - swapped lambda calculation from
- // double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)),
+ // KMH & PR - 16 May 2008 - swapped lambda calculation from
+ // double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)),
// which should be identical for p1 & p2 onshell in their COM
// but in the inverse construction for the Nason method, this
- // was not the case, leading to misuse.
+ // was not the case, leading to misuse.
const double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
*(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b)));
-
+
Energy firstQ, secondQ;
double kappab(0.), kappac(0.);
- //key = 0; // symmetric case pre-selection
+ //key = 0; // symmetric case pre-selection
switch(key) {
case 0: // symmetric case
firstQ = sqrt(0.5*Q2*(1.+b-c+lam));
secondQ = sqrt(0.5*Q2*(1.-b+c+lam));
break;
case 1: // maximum emission from both legs
kappab=4.*(1.-2.*sqrt(c)-b+c);
kappac=4.*(1.-2.*sqrt(b)-c+b);
firstQ = sqrt(Q2*kappab);
secondQ = sqrt(Q2*kappac);
break;
default:
assert(false);
- }
+ }
// calculate the scales
return pair<Energy,Energy>(firstQ, secondQ);
}
pair<Energy,Energy>
PartnerFinder::calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
const bool isDecayCase) {
- if(!isDecayCase) {
+ if(!isDecayCase) {
// In this case from JHEP 12(2003)045 we find the conditions
// ktilde_b = (1+c) and ktilde_c = (1+2c)
// We also find that c = m_c^2/Q^2. The process is a+b->c where
// particle a is not colour connected (considered as a colour singlet).
- // Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and
+ // Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and
// q_c = sqrt(Q^2+2 m_c^2)
// We also assume that the first particle in the pair is the initial
- // state particle and the second is the final state one c
+ // state particle and the second is the final state one c
const Energy2 mc2 = sqr(pc.mass());
const Energy2 Q2 = -(pb-pc).m2();
return { sqrt(Q2+mc2), sqrt(Q2+2*mc2) };
}
- else {
+ else {
// In this case from JHEP 12(2003)045 we find, for the decay
// process b->c+a(neutral), the condition
- // (ktilde_b-1)*(ktilde_c-c)=(1/4)*sqr(1-a+c+lambda).
+ // (ktilde_b-1)*(ktilde_c-c)=(1/4)*sqr(1-a+c+lambda).
// We also assume that the first particle in the pair is the initial
// state particle (b) and the second is the final state one (c).
- // - We find maximal phase space coverage through emissions from
+ // - We find maximal phase space coverage through emissions from
// c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c)
// - We find the most 'symmetric' way to populate the phase space
- // occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda)
+ // occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda)
// - We find the most 'smooth' way to populate the phase space
// occurs for...
Energy2 mb2(sqr(pb.mass()));
double a=(pb-pc).m2()/mb2;
double c=sqr(pc.mass())/mb2;
double lambda = 1. + a*a + c*c - 2.*a - 2.*c - 2.*a*c;
lambda = sqrt(max(lambda,0.));
const double PROD = 0.25*sqr(1. - a + c + lambda);
int key = 0;
double ktilde_b, ktilde_c, cosi = 0.;
switch (key) {
case 0: // the 'symmetric' choice
ktilde_c = 0.5*(1-a+c+lambda) + c ;
ktilde_b = 1.+PROD/(ktilde_c-c) ;
break;
case 1: // the 'maximal' choice
ktilde_c = 4.0*(sqr(1.-sqrt(a))-c);
ktilde_b = 1.+PROD/(ktilde_c-c) ;
break;
case 2: // the 'smooth' choice
c = max(c,1.*GeV2/mb2);
cosi = (sqr(1-sqrt(c))-a)/lambda;
ktilde_b = 2.0/(1.0-cosi);
ktilde_c = (1.0-a+c+lambda)*(1.0+c-a-lambda*cosi)/(2.0*(1.0+cosi));
break;
}
return { sqrt(mb2*ktilde_b), sqrt(mb2*ktilde_c) };
}
}
pair<Energy,Energy>
PartnerFinder::calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) {
// This case is quite simple. From JHEP 12(2003)045 we find the condition
// that ktilde_b = ktilde_c = 1. In this case we have the process
// b+c->a so we need merely boost to the CM frame of the two incoming
// particles and then qtilde is equal to the energy in that frame
const Energy Q = sqrt((p1+p2).m2());
return {Q,Q};
}
diff --git a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc
@@ -1,243 +1,242 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HalfHalfZeroEWSplitFn class.
//
#include "HalfHalfZeroEWSplitFn.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/ParticleData.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
#include "Herwig/Models/StandardModel/SMFFHVertex.h"
using namespace Herwig;
IBPtr HalfHalfZeroEWSplitFn::clone() const {
return new_ptr(*this);
}
IBPtr HalfHalfZeroEWSplitFn::fullclone() const {
return new_ptr(*this);
}
void HalfHalfZeroEWSplitFn::persistentOutput(PersistentOStream & os) const {
os << ghqq_;
}
void HalfHalfZeroEWSplitFn::persistentInput(PersistentIStream & is, int) {
is >> ghqq_;
}
-// The following static variable is needed for the type description system in ThePEG.
+// The following static variable is needed for the type description system in ThePEG.
DescribeClass<HalfHalfZeroEWSplitFn,SplittingFunction>
describeHerwigHalfHalfZeroEWSplitFn("Herwig::HalfHalfZeroEWSplitFn", "HwShower.so");
void HalfHalfZeroEWSplitFn::Init() {
static ClassDocumentation<HalfHalfZeroEWSplitFn> documentation
("The HalfHalfZeroEWSplitFn class implements the splitting q->qH");
}
void HalfHalfZeroEWSplitFn::doinit() {
SplittingFunction::doinit();
tcSMPtr sm = generator()->standardModel();
double sw2 = sm->sin2ThetaW();
ghqq_ = 1./sqrt(4.*sw2);
_theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
}
void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids) const {
if(abs(ids[2]->id())==ParticleID::h0) {
//get quark masses
Energy mq;
- if(abs(ids[0]->id())==ParticleID::c)
+ if(abs(ids[0]->id())==ParticleID::c)
mq = getParticleData(ParticleID::c)->mass();
- else if(abs(ids[0]->id())==ParticleID::b)
+ else if(abs(ids[0]->id())==ParticleID::b)
mq = getParticleData(ParticleID::b)->mass();
- else if(abs(ids[0]->id())==ParticleID::t)
+ else if(abs(ids[0]->id())==ParticleID::t)
mq = getParticleData(ParticleID::t)->mass();
Energy mW = getParticleData(ParticleID::Wplus)->mass();
gH = ghqq_*(mq/mW);
}
else
assert(false);
}
void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids, const Energy2 t) const {
if(abs(ids[2]->id())==ParticleID::h0) {
//get quark masses
Energy mq;
- if(abs(ids[0]->id())==ParticleID::c)
+ if(abs(ids[0]->id())==ParticleID::c)
mq = _theSM->mass(t,getParticleData(ParticleID::c));
- else if(abs(ids[0]->id())==ParticleID::b)
+ else if(abs(ids[0]->id())==ParticleID::b)
mq = _theSM->mass(t,getParticleData(ParticleID::b));
- else if(abs(ids[0]->id())==ParticleID::t)
+ else if(abs(ids[0]->id())==ParticleID::t)
mq = _theSM->mass(t,getParticleData(ParticleID::t));
Energy mW = getParticleData(ParticleID::Wplus)->mass();
//Energy mW = _theSM->mass(t,getParticleData(ParticleID::Wplus));
gH = ghqq_*(mq/mW);
}
else
assert(false);
}
double HalfHalfZeroEWSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
double gH(0.);
getCouplings(gH,ids,t);
double val = (1.-z);
Energy mq, mH;
//get masses
if(mass) {
- mq = ids[0]->mass();
- mH = ids[2]->mass();
+ mq = ids[0]->mass();
+ mH = ids[2]->mass();
}
else { // to assure the particle mass in non-zero
- if(abs(ids[0]->id())==ParticleID::c)
+ if(abs(ids[0]->id())==ParticleID::c)
mq = getParticleData(ParticleID::c)->mass();
- else if(abs(ids[0]->id())==ParticleID::b)
+ else if(abs(ids[0]->id())==ParticleID::b)
mq = getParticleData(ParticleID::b)->mass();
- else if(abs(ids[0]->id())==ParticleID::t)
- mq = getParticleData(ParticleID::t)->mass();
- mH = getParticleData(ParticleID::h0)->mass();
+ else if(abs(ids[0]->id())==ParticleID::t)
+ mq = getParticleData(ParticleID::t)->mass();
+ mH = getParticleData(ParticleID::h0)->mass();
}
val += (4.*sqr(mq) - sqr(mH))/(t*(1. - z)*z);
val *= sqr(gH);
return colourFactor(ids)*val;
}
double HalfHalfZeroEWSplitFn::overestimateP(const double z,
const IdList & ids) const {
double gH(0.);
getCouplings(gH,ids);
- return sqr(gH)*colourFactor(ids)*(1.-z);
+ return sqr(gH)*colourFactor(ids)*(1.-z);
}
double HalfHalfZeroEWSplitFn::ratioP(const double z, const Energy2 t,
const IdList & ids, const bool mass,
const RhoDMatrix & rho) const {
double gH(0.);
getCouplings(gH,ids,t);
double val = 1.;
Energy mq, mH;
if(mass) {
mq = ids[0]->mass();
- mH = ids[2]->mass();
+ mH = ids[2]->mass();
}
else { // to assure the particle mass in non-zero
- if(abs(ids[0]->id())==ParticleID::c)
+ if(abs(ids[0]->id())==ParticleID::c)
mq = getParticleData(ParticleID::c)->mass();
- else if(abs(ids[0]->id())==ParticleID::b)
+ else if(abs(ids[0]->id())==ParticleID::b)
mq = getParticleData(ParticleID::b)->mass();
- else if(abs(ids[0]->id())==ParticleID::t)
- mq = getParticleData(ParticleID::t)->mass();
+ else if(abs(ids[0]->id())==ParticleID::t)
+ mq = getParticleData(ParticleID::t)->mass();
mH = getParticleData(ParticleID::h0)->mass();
- }
+ }
val += (4.*sqr(mq) - sqr(mH))/(t*(1. - z)*z);
return val;
-}
+}
double HalfHalfZeroEWSplitFn::integOverP(const double z,
const IdList & ids,
unsigned int PDFfactor) const {
double gH(0.);
getCouplings(gH,ids);
double pre = colourFactor(ids)*sqr(gH);
switch (PDFfactor) {
case 0: //OverP
- return pre*(z-sqr(z)/2.);
+ return pre*(z-sqr(z)/2.);
case 1: //OverP/z
- return pre*(log(z)-z);
+ return pre*(log(z)-z);
case 2: //OverP/(1-z)
- return pre*z;
+ return pre*z;
case 3: //OverP/[z(1-z)]
return pre*log(z);
default:
throw Exception() << "HalfHalfZeroEWSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
- }
+ }
}
double HalfHalfZeroEWSplitFn::invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor) const {
double gH(0.);
getCouplings(gH,ids);
double pre = colourFactor(ids)*sqr(gH);
switch (PDFfactor) {
- case 0:
- return max((-pre+sqrt(sqr(pre)-2.*pre*r))/pre,
- (-pre-sqrt(sqr(pre)-2.*pre*r))/pre);
+ case 0:
+ return 1. - sqrt(1. - 2.*r/pre);
case 1: //OverP/z
case 2: //OverP/(1-z)
- return r/pre;
+ return r/pre;
case 3: //OverP/[z(1-z)]
- return exp(r/pre);
+ return exp(r/pre);
default:
throw Exception() << "HalfHalfZeroEWSplitFn::invIntegOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
- }
+ }
}
bool HalfHalfZeroEWSplitFn::accept(const IdList &ids) const {
if(ids.size()!=3) return false;
if(ids[2]->id()==ParticleID::h0) {
- if(ids[0]->id()==ids[1]->id() && (ids[0]->id()==4 || ids[0]->id()==5 || ids[0]->id()==6))
- return true;
+ if(ids[0]->id()==ids[1]->id() && (ids[0]->id()==4 || ids[0]->id()==5 || ids[0]->id()==6))
+ return true;
}
return false;
}
-vector<pair<int, Complex> >
+vector<pair<int, Complex> >
HalfHalfZeroEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
-vector<pair<int, Complex> >
+vector<pair<int, Complex> >
HalfHalfZeroEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
-DecayMEPtr HalfHalfZeroEWSplitFn::matrixElement(const double z, const Energy2 t,
+DecayMEPtr HalfHalfZeroEWSplitFn::matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool) {
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin0)));
//get masses
Energy mq, mH;
- if(abs(ids[0]->id())==ParticleID::c)
+ if(abs(ids[0]->id())==ParticleID::c)
mq = getParticleData(ParticleID::c)->mass();
- else if(abs(ids[0]->id())==ParticleID::b)
+ else if(abs(ids[0]->id())==ParticleID::b)
mq = getParticleData(ParticleID::b)->mass();
- else if(abs(ids[0]->id())==ParticleID::t)
- mq = getParticleData(ParticleID::t)->mass();
+ else if(abs(ids[0]->id())==ParticleID::t)
+ mq = getParticleData(ParticleID::t)->mass();
mH = getParticleData(ParticleID::h0)->mass();
double gH(0.);
getCouplings(gH,ids,t);
double mqt = mq/sqrt(t);
double mHt = mH/sqrt(t);
double num1 = gH*(1.+z)*mqt;
double num2 = gH*sqrt(-sqr(mqt)*(1.-z) - sqr(mHt)*z + z*(1.-z)*(sqr(mqt)+z*(1.-z))); //watch this
double dnum = sqrt(2.)*sqrt((1.-z)*sqr(z));
Complex phase = exp(Complex(0.,1.)*phi);
Complex cphase = conj(phase);
(*kernal)(0,0,0) = num1/dnum;
(*kernal)(0,1,0) = cphase*num2/dnum;
(*kernal)(1,0,0) = -phase*num2/dnum;
(*kernal)(1,1,0) = num1/dnum;
// return the answer
return kernal;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 12:50 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4020756
Default Alt Text
(38 KB)

Event Timeline