Page MenuHomeHEPForge
Files F8309017

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc
--- a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc
+++ b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc
@@ -1,842 +1,843 @@
// -*- C++ -*-
//
// MadGraphAmplitude.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MadGraphAmplitude class.
//
#include "MadGraphAmplitude.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/EnumParticles.h"
#include <boost/lexical_cast.hpp>
#include <boost/filesystem.hpp>
#include <cstdlib>
#include <dlfcn.h>
#include <errno.h>
#include <sstream>
using namespace Herwig;
#ifndef HERWIG_BINDIR
#error Makefile.am needs to define HERWIG_BINDIR
#endif
#ifndef HERWIG_PKGDATADIR
#error Makefile.am needs to define HERWIG_PKGDATADIR
#endif
#ifndef MADGRAPH_PREFIX
#error Makefile.am needs to define MADGRAPH_PREFIX
#endif
extern "C" void mginitproc_(char *i,int);
extern "C" void MG_Calculate_wavefunctions_virt(int* proc,double*,double*);
extern "C" void MG_Calculate_wavefunctions_born(int* proc,double*, int*);
extern "C" void MG_Jamp (int* proc,int*, double*);
extern "C" void MG_LNJamp (int* proc,int*, double*);
extern "C" void MG_Virt (int* proc,double*);
extern "C" void MG_NCol (int* proc,int*);
extern "C" void MG_vxxxxx (double* p,double* n,int* inc,double* );
extern "C" void MG_Colour (int* proc,int* i,int* j ,int* color);
MadGraphAmplitude::MadGraphAmplitude()
: theMGmodel("loop_sm"),keepinputtopmass(false),
bindir_(HERWIG_BINDIR), includedir_(HERWIG_INCLUDEDIR), pkgdatadir_(HERWIG_PKGDATADIR), madgraphPrefix_(MADGRAPH_PREFIX)
{}
MadGraphAmplitude::~MadGraphAmplitude() {
}
IBPtr MadGraphAmplitude::clone() const {
return new_ptr(*this);
}
IBPtr MadGraphAmplitude::fullclone() const {
return new_ptr(*this);
}
bool MadGraphAmplitude::initializedMad=false;
vector<string> MadGraphAmplitude::BornAmplitudes=vector<string>();
vector<string> MadGraphAmplitude::VirtAmplitudes=vector<string>();
void MadGraphAmplitude::initProcess(const cPDVector& ) {
if ( lastMatchboxXComb()->initialized() )
return;
if ( !DynamicLoader::load(mgProcLibPath()+"InterfaceMadGraph.so") )
throw Exception() << "MadGraphAmplitude: Failed to load MadGraph amplitudes\n"
<< DynamicLoader::lastErrorMessage
<< Exception::runerror;
if (!initializedMad){
string mstr=(factory()->runStorage()+"MadGraphAmplitudes"+"/param_card"+((theMGmodel=="loop_sm")?"":("_"+theMGmodel))+".dat");
+ if( theMGmodel[0]=='/')mstr="param_card.dat";
size_t len = mstr.size();
mginitproc_(const_cast<char*>(mstr.c_str()),len);
initializedMad=true;
}
lastMatchboxXComb()->isInitialized();
}
bool MadGraphAmplitude::writeAmplitudesDat(){
bool res=false;
string born= mgProcLibPath()+"BornAmplitudes.dat";
if ( !boost::filesystem::exists(born) ) {
ofstream borns(born.c_str());
for (vector<string>::iterator amps=BornAmplitudes.begin();amps!=BornAmplitudes.end();amps++)
borns<<*amps<<endl;
borns.close();
res=true;
}
string virt= mgProcLibPath()+"VirtAmplitudes.dat";
if ( !boost::filesystem::exists(virt) ) {
ofstream virts(virt.c_str());
for (vector<string>::iterator amps=VirtAmplitudes.begin();amps!=VirtAmplitudes.end();amps++)
virts<<*amps<<endl;
virts.flush();
virts.close();
res=true;
}
return res;
}
bool MadGraphAmplitude::checkAmplitudes(){
string born= mgProcLibPath()+"BornAmplitudes.dat";
string virt= mgProcLibPath()+"VirtAmplitudes.dat";
assert ( boost::filesystem::exists(born)|| boost::filesystem::exists(virt));
bool foundallborns=true;
for (vector<string>::iterator amps=BornAmplitudes.begin();amps!=BornAmplitudes.end();amps++){
ifstream borns(born.c_str());
string line;
bool foundthisborn=false;
while (std::getline(borns, line)) {
if(line==*amps)foundthisborn=true;
}
foundallborns&=foundthisborn;
}
bool foundallvirts=true;
for (vector<string>::iterator amps=VirtAmplitudes.begin();amps!=VirtAmplitudes.end();amps++){
ifstream virts(virt.c_str());
string line;
bool foundthisvirt=false;
while (std::getline(virts, line)) {
if(line==*amps)foundthisvirt=true;
}
foundallvirts&=foundthisvirt;
}
if (!foundallborns||!foundallvirts)
throw Exception() << "MadGraphAmplitude: One amplitude has no externalId. Please remove the MadGraphAmplitude-folder and rebuild.\n" << Exception::runerror;
return foundallborns && foundallvirts;
}
string MadGraphAmplitude::mgProcLibPath(){
string res=theProcessPath == "" ? factory()->buildStorage()+"MadGraphAmplitudes" : theProcessPath;
if (res.at(res.length()-1) != '/') res.append("/");
return res;
}
bool MadGraphAmplitude::initializeExternal() {
if ( boost::filesystem::exists(mgProcLibPath()) ) {
if ( !boost::filesystem::is_directory(mgProcLibPath()) )
throw Exception() << "MadGraphAmplitude: MadGraph amplitude storage '"
<< mgProcLibPath() << "' existing but not a directory."
<< Exception::runerror;
} else {
boost::filesystem::create_directories(mgProcLibPath());
}
string runAmplitudes = factory()->runStorage() + "/MadGraphAmplitudes";
if ( boost::filesystem::exists(runAmplitudes) ) {
if ( !boost::filesystem::is_directory(runAmplitudes) )
throw Exception() << "MadGraphAmplitude: MadGraph amplitude storage '"
<< runAmplitudes << "' existing but not a directory."
<< Exception::runerror;
} else {
boost::filesystem::create_directories(runAmplitudes);
}
//EW-consistency check:
Energy MW=getParticleData(ParticleID::Wplus)->hardProcessMass();
Energy MZ=getParticleData(ParticleID::Z0)->hardProcessMass();
if( MW!= sqrt(MZ*MZ/2.0+sqrt(MZ*MZ*MZ*MZ/4.0-Constants::pi*SM().alphaEMMZ()*MZ*MZ/ sqrt(2.0)/SM().fermiConstant()))){
generator()->log()<<"\n\n-----!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-----";
generator()->log() << "\nYou are using a EW scheme which is inconsistent with the MadGraph parametisation:\n\n"
<<MW/GeV<< " GeV==MW!= sqrt(MZ^2/2+sqrt(MZ^4/4.0-pi*alphaEMMZ*MZ^2/ sqrt(2)/G_f))=="<<
sqrt(MZ*MZ/2.0+sqrt(MZ*MZ*MZ*MZ/4.0-Constants::pi*SM().alphaEMMZ()*MZ*MZ/ sqrt(2.0)/SM().fermiConstant()))/GeV
<<" GeV\n\n-----!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-----\n";
}
string para= factory()->runStorage()+"/MadGraphAmplitudes"+"/MG-Parameter.dat";
ofstream params(para.c_str());
params<<"$WZ$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Z0)->hardProcessWidth() /GeV;
params<<"\n$WW$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Wplus)->hardProcessWidth()/GeV;
params<<"\n$alphas$ " <<std::setiosflags(ios::scientific) <<SM().alphaS();
params<<"\n$GF$ " <<std::setiosflags(ios::scientific) <<SM().fermiConstant()*GeV2 ;
params<<"\n$alphaMZ$ " <<std::setiosflags(ios::scientific) <<1/SM().alphaEMMZ();
params<<"\n$MZ$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Z0)->hardProcessMass() /GeV<<flush;
params<<"\n$MW$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Wplus)->hardProcessMass() /GeV<<flush;
params<<"\n$sw2$ " <<std::setiosflags(ios::scientific) << SM().sin2ThetaW() <<flush;
if(theMGmodel=="heft"&&!keepinputtopmass){
if ( factory()->initVerbose() ) {
generator()->log()<<"\n---------------------------------------------------------------";
generator()->log()<<"\n---------------------------------------------------------------";
generator()->log()<<"\nNote: You are using the Higgs Effective model (heft) in ";
generator()->log()<<"\n Madgraph. We assume you try to calculate NLO with ";
generator()->log()<<"\n the GoSam virtual amplitudes. To match the models we ";
generator()->log()<<"\n therefore set the topmass to 10000000 GeV.";
generator()->log()<<"\n\n For more information see the \\tau parameter in:";
generator()->log()<<"\n https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Models/HiggsEffective";
generator()->log()<<"\n\n The Effective Higgs model in Gosam is using mT=infinity";
generator()->log()<<"\n\n\n If you want to use the LO matrixelements of MadGraph with finite' topmass you need to add: ";
generator()->log()<<"\n\n set Madgraph:KeepInputTopMass True";
generator()->log()<<"\n\n to your input file.";
generator()->log()<<"\n---------------------------------------------------------------";
generator()->log()<<"\n---------------------------------------------------------------\n";
}
params<<"\n$MT$ 10000000." <<flush;
}else{
params<<"\n$MT$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::t)->hardProcessMass() /GeV <<flush;
}
params<<"\n$WT$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::t)->hardProcessWidth() /GeV <<flush;
params<<"\n$MB$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::b)->hardProcessMass() /GeV <<flush;
params<<"\n$MH$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::h0)->hardProcessMass() /GeV <<flush;
params<<"\n$WH$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::h0)->hardProcessWidth() /GeV <<flush;
params<<"\n$MTA$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::tauplus)->hardProcessMass() /GeV <<flush;
string cmd = "python " + bindir_ + "/mg2herwig ";
cmd +=" --buildpath "+mgProcLibPath();
cmd +=" --model "+theMGmodel;
cmd +=" --runpath "+factory()->runStorage()+"/MadGraphAmplitudes ";
cmd +=" --datadir "+pkgdatadir_;
cmd +=" --includedir "+includedir_;
std::stringstream as,aem;
as << factory()->orderInAlphaS();
cmd +=" --orderas "+as.str() ;
aem <<factory()->orderInAlphaEW();
cmd +=" --orderew "+aem.str();
// TODO move to boost::system
writeAmplitudesDat();
if (boost::filesystem::exists(mgProcLibPath()+"InterfaceMadGraph.so") ){
//set the parameters
checkAmplitudes();
std::system(cmd.c_str());
ranMadGraphInitializeExternal = true;
return true;
}
char cwd[1024];
if ( !getcwd(cwd,sizeof(cwd)) )
throw Exception() << "MadGraphAmplitude: failed to determine current working directory\n"
<< Exception::runerror;
cmd +=" --madgraph " + madgraphPrefix_ + "/bin " ;
cmd +="--build > ";
cmd += mgProcLibPath()+"MG.log 2>&1";
generator()->log() << "\n>>> Compiling MadGraph amplitudes. This may take some time -- please be patient.\n"
<< ">>> In case of problems see " << mgProcLibPath() << "MG.log for details.\n\n"
<< flush;
std::system(cmd.c_str());
cmd = "python " + bindir_ + "/mg2herwig ";
cmd +=" --buildpath "+mgProcLibPath();
cmd +=" --model "+theMGmodel;
cmd +=" --runpath "+factory()->runStorage()+"/MadGraphAmplitudes ";
cmd +=" --datadir "+pkgdatadir_;
as.clear();
aem.clear();
as << factory()->orderInAlphaS();
cmd +=" --orderas "+as.str() ;
aem <<factory()->orderInAlphaEW();
cmd +=" --orderew "+aem.str();
std::system(cmd.c_str());
ranMadGraphInitializeExternal = true;
return boost::filesystem::exists(mgProcLibPath()+"InterfaceMadGraph.so");
}
int MadGraphAmplitude::externalId(const cPDVector& proc) {
for (int i=0;i<100;i++){
colourindex.push_back(-2);
}
assert(!BornAmplitudes.empty()||!VirtAmplitudes.empty());
writeAmplitudesDat();
int res=0;
string amp="";
int k=0;
for (cPDVector::const_iterator it=proc.begin();it!=proc.end();it++,k++){
amp+=boost::lexical_cast<string>( (*it)->id())+" ";if (k==1)amp+=" > ";
}
string born= mgProcLibPath()+"BornAmplitudes.dat";
string virt= mgProcLibPath()+"VirtAmplitudes.dat";
assert ( boost::filesystem::exists(born)|| boost::filesystem::exists(virt));
ifstream borns(born.c_str());
string line;
while (std::getline(borns, line)) {
res+=1;
if(line==amp)return res;
}
ifstream virts(virt.c_str());
while (std::getline(virts, line)) {
res+=1;
if(line==amp)return res;
}
throw Exception() << "MadGraphAmplitude: One amplitude has no externalId. Please remove the MadGraphAmplitude-folder and rebuild.\n" << Exception::runerror;
return res;
}
bool MadGraphAmplitude::ranMadGraphInitializeExternal = false;
void MadGraphAmplitude::doinit() {
if ( !ranMadGraphInitializeExternal ) {
initializeExternal();
}
MatchboxAmplitude::doinit();
}
void MadGraphAmplitude::doinitrun() {
if ( !ranMadGraphInitializeExternal ) {
initializeExternal();
}
MatchboxAmplitude::doinitrun();
}
bool MadGraphAmplitude::canHandle(const PDVector& p,
Ptr<MatchboxFactory>::tptr factory,
bool virt) const {
if ( factory->processData()->diagramMap().find(p) !=
factory->processData()->diagramMap().end() )
return true;
vector<Ptr<Tree2toNDiagram>::ptr> diags =
factory->diagramGenerator()->generate(p,orderInGs(),orderInGem());
if ( diags.empty() )
return false;
factory->processData()->diagramMap()[p] = diags;
string amp="";
int k=0;
for (PDVector::const_iterator it=p.begin();it!=p.end();it++,k++){
amp+=boost::lexical_cast<string>( (*it)->id())+" ";if (k==1)amp+=" > ";
}
if (virt && factory->highestVirt()>=p.size()){
VirtAmplitudes.push_back(amp);
}else{
BornAmplitudes.push_back(amp);
}
return true;
}
void MadGraphAmplitude::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr me) {
useMe();
if ( !calculateTreeAmplitudes() ) {
MatchboxAmplitude::prepareAmplitudes(me);
return;
}
if (colourindex.empty()) {
for (int i=0;i<100;i++){
colourindex.push_back(-2);
}
}
lastMatchboxXComb()->clearheljamp();
lastMatchboxXComb()->clearhelLNjamp();
initProcess(mePartonData());
MatchboxAmplitude::prepareAmplitudes(me);
}
Complex MadGraphAmplitude::evaluate(size_t i, const vector<int>& hel, Complex& largeN) {
//find the colourline:
int ii = -1;
int xx=lastMatchboxXComb()->externalId();
if (colourindex.size()<=i) {
colourindex.clear();
for (size_t l=0;l<=i+10;l++){
colourindex.push_back(-2);
}
}
if(colourindex[i]!=-2){
ii = colourindex[i];
if (ii==-1) {
largeN = Complex(0.0);
return Complex(0.0);
}
} else {
set<vector<size_t> > a = colourOrdering(i);
int ncol=-1;
MG_NCol(&xx,&ncol);
assert(ncol!=-1);
for( int it = 0; it < ncol; it++ ){
int n = 0;
for ( cPDVector::const_iterator nx = mePartonData().begin();
nx != mePartonData().end(); nx++ )
if ( (*nx)->coloured() ) n++;
set<vector<size_t> > tmpset;
vector<size_t> tmpvek;
for ( int it2 = 0; it2 < n; it2++ ) {
int ret=-2;
MG_Colour(&xx,&it,&it2,&ret);
assert(ret !=-2);
if (ret== -1)
break;
if ( ret == 0 ) {
n++;
tmpset.insert(tmpvek);
tmpvek.clear();
} else {
tmpvek.push_back(ret-1);
}
if( it2 == n-1 ) tmpset.insert(tmpvek);
}
bool found_all = true;
for ( set<vector<size_t> >::iterator it3 = a.begin(); it3 != a.end(); it3++ ) {
bool found_it3=false;
for ( set<vector<size_t> >::iterator it4 = tmpset.begin(); it4 != tmpset.end(); it4++ ) {
vector<size_t> it3tmp = gluonsFirst(*it3);
vector<size_t> it4tmp = (*it4);
if ( it3tmp.size() != it4tmp.size() ) continue;
if ( it3tmp == it4tmp ) found_it3 = true;
}
found_all = found_all && found_it3;
}
if ( found_all ) {
colourindex[i]=it;
ii=it;
}
}
}
if ( ii == -1 ){
colourindex[i]=ii;
largeN = Complex(0.0);
return Complex(0.0);
}
const map<vector<int>,vector < complex<double> > >& tmp = lastMatchboxXComb()->heljamp();
const map<vector<int>,vector < complex<double> > >& tmpLN = lastMatchboxXComb()->helLNjamp();
if( tmp.find(hel) != tmp.end()) {
largeN = tmpLN.find(hel)->second[ii];
return tmp.find(hel)->second[ii];;
}
double units = pow(sqrt(lastSHat())/GeV,int(hel.size())-4);
int heltmp[10];
for(size_t j=0;j<hel.size();j++){
int cross=crossingMap()[j];
if( (cross>1&&j<=1)||(cross<=1&&j>1)){
heltmp[cross]=-1*hel[j];}
else{heltmp[cross]=hel[j];}
}
vector<Lorentz5Momentum> reshuffled = meMomenta();
if ( !reshuffleMasses().empty() && reshuffled.size() > 3 ) {
const cPDVector& pdata = mePartonData();
const map<long,Energy>& masses = reshuffleMasses();
lastMatchboxXComb()->reshuffle(reshuffled,pdata,masses);
}
double momenta[50];
size_t j=0;
for (size_t i=0;i<mePartonData().size();i++){
momenta[j]=abs(reshuffled[i].e()/GeV)<1.e-13?0.:double(reshuffled[i].e()/GeV);
momenta[j+1]=abs(reshuffled[i].x()/GeV)<1.e-13?0.:double(reshuffled[i].x()/GeV);
momenta[j+2]=abs(reshuffled[i].y()/GeV)<1.e-13?0.:double(reshuffled[i].y()/GeV);
momenta[j+3]=abs(reshuffled[i].z()/GeV)<1.e-13?0.:double(reshuffled[i].z()/GeV);
if(momenta[j ] == 0. && momenta[j+1] == 0. &&
momenta[j+2] == 0. && momenta[j+3] == 0. )
return 0.;
j+=4;
}
MG_Calculate_wavefunctions_born(&xx, &momenta[0], &heltmp[0]);
int ncol=-1;
MG_NCol(&xx,&ncol);
Complex res;
Complex resLN;
for( int it = 0; it < ncol; it++ ){
double dd[2];
MG_Jamp(&xx,&it,&dd[0]);
Complex d(dd[0],dd[1]);
if(it==ii)res=d*units;
lastMatchboxXComb()->pushheljamp(hel,d*units);
double ddLN[2];
MG_LNJamp(&xx,&it,&ddLN[0]);
Complex dLN(ddLN[0],ddLN[1]);
if(it==ii)resLN=dLN*units;
lastMatchboxXComb()->pushhelLNjamp(hel,dLN*units);
}
largeN = resLN;
return res;
}
vector<unsigned int> MadGraphAmplitude::physicalHelicities(const vector<int>& hel) const {
vector<unsigned int> res(hel.size(),0);
for ( size_t j = 0; j < hel.size(); ++j ) {
int cross = crossingMap()[j];
int xhel = 0;
if ( (cross > 1 && j <= 1) || (cross <= 1 && j > 1) )
xhel = -1*hel[j];
else
xhel = hel[j];
if ( mePartonData()[cross]->iSpin() == PDT::Spin1Half )
res[cross] = (xhel == -1 ? 0 : 1);
else if ( mePartonData()[cross]->iSpin() == PDT::Spin1 )
res[cross] = (unsigned int)(xhel + 1);
else if ( mePartonData()[cross]->iSpin() == PDT::Spin0 )
res[cross] = 0;
else assert(false);
}
return res;
}
LorentzVector<Complex> MadGraphAmplitude::plusPolarization(const Lorentz5Momentum& p,
const Lorentz5Momentum& n,
int i) const {
int tmp=i;
double pg[4],ng[4],poltmp[8];
pg[0]=p.e()/GeV;pg[1]=p.x()/GeV;pg[2]=p.y()/GeV;pg[3]=p.z()/GeV;
ng[0]=n.e()/GeV;ng[1]=n.x()/GeV;ng[2]=n.y()/GeV;ng[3]=n.z()/GeV;
MG_vxxxxx(&pg[0],&ng[0],&tmp,&poltmp[0]);
complex<double> pol[6];
pol[0]=Complex(poltmp[0],poltmp[1]);
pol[1]=Complex(poltmp[2],poltmp[3]);
pol[2]=Complex(poltmp[4],poltmp[5]);
pol[3]=Complex(poltmp[6],poltmp[7]);
LorentzVector<Complex> polarization(pol[1],pol[2],pol[3],pol[0]);
return polarization;
}
bool equalsModulo(unsigned int i, const vector<int>& a, const vector<int>& b) {
assert(a.size()==b.size());
if ( a[i] == b[i] )
return false;
for ( unsigned int k = 0; k < a.size(); ++k ) {
if ( k == i )
continue;
if ( a[k] != b[k] )
return false;
}
return true;
}
vector<size_t> MadGraphAmplitude::gluonsFirst(vector<size_t> vec) {
vector<size_t> vecout;
for(vector<size_t>::iterator it= vec.begin();it!= vec.end();++it)
if ( mePartonData()[crossingMap()[*it]]->id()==21)
vecout.push_back(crossingMap()[*it]);
for(vector<size_t>::iterator it= vec.begin();it!= vec.end();++it)
if ( mePartonData()[crossingMap()[*it]]->id()!=21)
vecout.push_back(crossingMap()[*it]);
return vecout;
}
double MadGraphAmplitude::spinColourCorrelatedME2(pair<int,int> ij,
const SpinCorrelationTensor& c) const {
vector<Lorentz5Momentum> reshuffled = meMomenta();
if ( !reshuffleMasses().empty() && reshuffled.size() > 3 ) {
const cPDVector& pdata = mePartonData();
const map<long,Energy>& masses = reshuffleMasses();
lastMatchboxXComb()->reshuffle(reshuffled,pdata,masses);
}
Lorentz5Momentum p = reshuffled[ij.first];
Lorentz5Momentum n = reshuffled[ij.second];
LorentzVector<Complex> polarization = plusPolarization(p,n,ij.first<2?-1:1);
int iCrossed = -1;
for ( unsigned int k = 0; k < crossingMap().size(); ++k )
if ( crossingMap()[k] == ij.first ) {
iCrossed = k;
break;
}
assert(iCrossed!=-1);
if(ij.first>1) polarization =polarization.conjugate();
if(iCrossed<2) polarization =polarization.conjugate();
Complex pFactor = (polarization*c.momentum())/sqrt(abs(c.scale()));
double avg =
colourCorrelatedME2(ij)*(-c.diagonal()+ (c.scale() > ZERO ? 1. : -1.)*norm(pFactor));
Complex csCorr = 0.0;
if ( calculateColourSpinCorrelator(ij) ) {
set<const CVector*> done;
for ( AmplitudeConstIterator a = lastAmplitudes().begin();
a != lastAmplitudes().end(); ++a ) {
if ( done.find(&(a->second)) != done.end() )
continue;
AmplitudeConstIterator b = lastAmplitudes().begin();
while ( !equalsModulo(iCrossed,a->first,b->first) )
if ( ++b == lastAmplitudes().end() )
break;
if ( b == lastAmplitudes().end() || done.find(&(b->second)) != done.end() )
continue;
done.insert(&(a->second)); done.insert(&(b->second));
if ( a->first[iCrossed] == 1 )
swap(a,b);
csCorr -= colourBasis()->colourCorrelatedInterference(ij,mePartonData(),a->second,b->second);
}
lastColourSpinCorrelator(ij,csCorr);
} else {
csCorr = lastColourSpinCorrelator(ij);
}
double corr =
2.*real(csCorr*sqr(pFactor));
double Nc = generator()->standardModel()->Nc();
double cfac = 1.;
if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) {
cfac = Nc;
} else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 ||
mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) {
cfac = (sqr(Nc)-1.)/(2.*Nc);
} else assert(false);
return
( avg +(c.scale() > ZERO ? 1. : -1.)*corr/cfac);
}
void MadGraphAmplitude::prepareOneLoopAmplitudes(Ptr<MatchboxMEBase>::tcptr ){
assert(false);
}
double MadGraphAmplitude::oneLoopInterference() const {
if ( !calculateOneLoopInterference() )
return lastOneLoopInterference();
evaloneLoopInterference();
return lastOneLoopInterference();
}
void MadGraphAmplitude::evaloneLoopInterference() const {
double units = pow(lastSHat()/GeV2,int(mePartonData().size())-4);
vector<Lorentz5Momentum> reshuffled = meMomenta();
if ( !reshuffleMasses().empty() && reshuffled.size() > 3 ) {
const cPDVector& pdata = mePartonData();
const map<long,Energy>& masses = reshuffleMasses();
lastMatchboxXComb()->reshuffle(reshuffled,pdata,masses);
}
double virt[20];
double momenta[50];
size_t j=0;
for (size_t i=0;i<mePartonData().size();i++){
momenta[j]=abs(reshuffled[i].e()/GeV)<1.e-13?0.:double(reshuffled[i].e()/GeV);
momenta[j+1]=abs(reshuffled[i].x()/GeV)<1.e-13?0.:double(reshuffled[i].x()/GeV);
momenta[j+2]=abs(reshuffled[i].y()/GeV)<1.e-13?0.:double(reshuffled[i].y()/GeV);
momenta[j+3]=abs(reshuffled[i].z()/GeV)<1.e-13?0.:double(reshuffled[i].z()/GeV);
j+=4;
}
int xx=lastMatchboxXComb()->externalId();
MG_Calculate_wavefunctions_virt(&xx,&momenta[0],&virt[0]);
double ifact = 1.;
ifact = 1./4.;
if (lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour3 ||
lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour3bar )
ifact /= SM().Nc();
else if ( lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour8 )
ifact /= (SM().Nc()*SM().Nc()-1.);
if ( lastMatchboxXComb()->matchboxME()->mePartonData()[1]->iColour() == PDT::Colour3 ||
lastMatchboxXComb()->matchboxME()->mePartonData()[1]->iColour() == PDT::Colour3bar )
ifact /= SM().Nc();
else if ( mePartonData()[1]->iColour() == PDT::Colour8 )
ifact /= (SM().Nc()*SM().Nc()-1.);
ifact *= lastMatchboxXComb()->matchboxME()->finalStateSymmetry();
lastOneLoopInterference(virt[1]/ifact*units);
lastOneLoopPoles(pair<double, double>(virt[2]/ifact*units,virt[3]/ifact*units));
}
void MadGraphAmplitude::persistentOutput(PersistentOStream & os) const {
os << theOrderInGs << theOrderInGem << BornAmplitudes << VirtAmplitudes
<< colourindex<<crossing << theProcessPath << theMGmodel << bindir_
<< pkgdatadir_ << madgraphPrefix_;
}
void MadGraphAmplitude::persistentInput(PersistentIStream & is, int) {
is >> theOrderInGs >> theOrderInGem >> BornAmplitudes >> VirtAmplitudes
>> colourindex>>crossing >> theProcessPath >> theMGmodel >> bindir_
>> pkgdatadir_ >> madgraphPrefix_;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<MadGraphAmplitude,MatchboxAmplitude>
describeHerwigMadGraphAmplitude("Herwig::MadGraphAmplitude", "HwMatchboxMadGraph.so");
void MadGraphAmplitude::Init() {
static ClassDocumentation<MadGraphAmplitude>
documentation("MadGraphAmplitude",
"Matrix elements have been calculated using MadGraph5 \\cite{Alwall:2011uj}",
"%\\cite{Alwall:2011uj}\n"
"\\bibitem{Alwall:2011uj}\n"
"J. Alwall et al.,\n"
"``MadGraph 5 : Going Beyond,''\n"
"arXiv:1106.0522 [hep-ph].\n"
"%%CITATION = ARXIV:1106.0522;%%");
static Parameter<MadGraphAmplitude,string> interfaceProcessPath
("ProcessPath",
"The Process Path.",
&MadGraphAmplitude::theProcessPath, "",false, false);
static Parameter<MadGraphAmplitude,string> interfaceModel
("Model",
"The MadGraph-Model.",
&MadGraphAmplitude::theMGmodel, "loop_sm",false, false);
static Switch<MadGraphAmplitude,bool> interfacekeepinputtopmass
("KeepInputTopMass",
"Switch On/Off formopt",
&MadGraphAmplitude::keepinputtopmass, false, false, false);
static SwitchOption interfacekeepinputtopmassTrue
(interfacekeepinputtopmass,
"On",
"On",
true);
static SwitchOption interfacekeepinputtopmassFalse
(interfacekeepinputtopmass,
"Off",
"Off",
false);
static Parameter<MadGraphAmplitude,string> interfaceBinDir
("BinDir",
"The location for the installed executable",
&MadGraphAmplitude::bindir_, string(HERWIG_BINDIR),
false, false);
static Parameter<MadGraphAmplitude,string> interfacePKGDATADIR
("DataDir",
"The location for the installed Herwig data files",
&MadGraphAmplitude::pkgdatadir_, string(HERWIG_PKGDATADIR),
false, false);
static Parameter<MadGraphAmplitude,string> interfaceMadgraphPrefix
("MadgraphPrefix",
"The prefix for the location of MadGraph",
&MadGraphAmplitude::madgraphPrefix_, string(MADGRAPH_PREFIX),
false, false);
}
diff --git a/MatrixElement/Matchbox/External/MadGraph/mg2herwig.py.in b/MatrixElement/Matchbox/External/MadGraph/mg2herwig.py.in
--- a/MatrixElement/Matchbox/External/MadGraph/mg2herwig.py.in
+++ b/MatrixElement/Matchbox/External/MadGraph/mg2herwig.py.in
@@ -1,376 +1,382 @@
#! /usr/bin/env python
import os,sys,glob,errno,shutil,time,fnmatch #argparse
from optparse import OptionParser
# helper to replace all sourceText in fileName with replaceText
def replacetext(fileName, sourceText, replaceText):
file = open(fileName, "r")
text = file.read()
file.close()
file = open(fileName, "w")
file.write(text.replace(sourceText, replaceText))
file.close()
# helper to build recursivly path
def mkdir_p(path):
try:
os.makedirs(path)
except OSError as exc: # Python >2.5
if exc.errno == errno.EEXIST and os.path.isdir(path):
pass
else: raise
# helper to find all files of with name in path
def find(name, path):
for root, dirs, files in os.walk(path):
if name in files:
return os.path.join(root, name)
# helper to find all file paths which contain file names matching filepattern
def finddirs(filepattern, path):
founddirs = []
for root, dirs, files in os.walk(path):
if fnmatch.filter(files, filepattern):
founddirs.append(root)
return founddirs
# fill the proc.dat file from BornAmplitudes.dat and VirtAmplitudes.dat.
def fillprocs(model,oras,orew):
bornlist=[]
virtlist=[]
fileproc=open("proc.dat","w")
fileproc.write("set fortran_compiler @FC@ --no_save\n")
fileproc.write("import model "+model+"\n")
borns="BornAmplitudes.dat"
virts="VirtAmplitudes.dat"
first=True
procnr=0
virtlines=""
bornlines=""
minlegs=100
legs=0
for i in [borns, virts]:
file = open(i, "r")
for line in file:
if (len(line.split(" "))<minlegs):
minlegs=len(line.split(" "))
for it in line.split(" "):
if it.replace("-","").isdigit():
legs+=1
file.close()
#conversion for heft model to go from (2QCD+1QED)->1HIG for each FS HIGGS.
HIG=0
if (model=="heft"):
HIG=(int(oras)+int(orew)-legs+2)/2
if (int(oras)+int(orew)-legs+2)%2!=0:
print "Warning: No possible coupling power:(int(oras)+int(orew)-legs+2)%2!=0 "
exit()
return
file = open(borns, "r")
for line in file:
#this assumes extra QCD emmissions
addalphas=len(line.split(" "))-minlegs
linetmp=line.rstrip()
procnr+=1
bornlist+=[str(procnr)]
if first:
if HIG ==0 :
bornlines+="generate "+linetmp+" QCD="+str(int(oras)+addalphas)+" QED="+str(orew)+" @"+str(procnr)+"\n"
else:
bornlines+="generate "+linetmp+" HIG="+str(HIG)+" QCD="+str(int(oras)+addalphas-2*HIG)+" QED="+str(int(orew)-HIG)+" @"+str(procnr)+"\n"
first=False
else:
if HIG ==0 :
bornlines+="add process "+linetmp+" QCD="+str(int(oras)+addalphas)+" QED="+str(orew)+" @"+str(procnr)+"\n"
else:
bornlines+="add process "+linetmp+" HIG="+str(HIG)+" QCD="+str(int(oras)+addalphas-2*HIG)+" QED="+str(int(orew)-HIG)+" @"+str(procnr)+"\n"
file.close()
first=True
file = open(virts, "r")
for line in file:
addalphas=len(line.split(" "))-minlegs
- linetmp=line.rstrip()+" QED="+str(int(orew))+" [ virt=QCD ]"
+ linetmp=line.rstrip()+" QCD="+str(int(oras)+addalphas)+" QED="+str(int(orew))+" [ virt=QCD ]"
procnr+=1
virtlist+=[str(procnr)]
if first:
virtlines+="generate "+linetmp+" @"+str(procnr)+"\n"
first=False
else:
virtlines+="add process "+linetmp+" @"+str(procnr)+"\n"
file.close()
fileproc.write(bornlines)
if virtlines!="" and bornlines!="":
fileproc.write("output matchbox MG5 --postpone_model\n")
fileproc.write(virtlines)
fileproc.write("output matchbox MG5 -f\n")
fileproc.close()
return bornlist,virtlist
def build_matchbox_tmp(pwd,buildpath):
cwd=os.getcwd()
os.chdir(pwd)
mkdir_p(pwd+"/Herwig-scratch/MG_tmp/")
if not buildpath.startswith("/"):
buildpath=pwd+"/"+buildpath.lstrip("./")
resources=glob.glob(buildpath +"/MG5/SubProcesses/MadLoop5_resources/*")
resources+=glob.glob(buildpath +"/MG5/Cards/*")
resources+=glob.glob(buildpath +"/MG5/Cards/SubProcesses/*")
for i in resources:
if not os.path.isfile( pwd+"/Herwig-scratch/MG_tmp/"+os.path.basename(i)) \
and not os.path.islink( pwd+"/Herwig-scratch/MG_tmp/"+os.path.basename(i)):
os.symlink(i, pwd+"/Herwig-scratch/MG_tmp/"+os.path.basename(i))
os.chdir(cwd)
parser = OptionParser()
parser.add_option("-a", "--buildpath", dest="buildpath",help="Do not use this script. Only for Herwig internal use. ")
parser.add_option("-b", "--build", action="store_true", dest="build", default=True,help="Do not use this script. Only for Herwig internal use.")
parser.add_option("-c", "--madgraph", dest="madgraph",help="Do not use this script. Only for Herwig internal use.")
parser.add_option("-d", "--runpath", dest="runpath",help="Do not use this script. Only for Herwig internal use.")
parser.add_option("-e", "--model", dest="model",help="Do not use this script. Only for Herwig internal use.")
parser.add_option("-f", "--orderas", dest="orderas",help="Do not use this script. Only for Herwig internal use.")
parser.add_option("-g", "--orderew", dest="orderew",help="Do not use this script. Only for Herwig internal use.")
parser.add_option("-i", "--datadir",dest="datadir",help="Do not use this script. Only for Herwig internal use.")
parser.add_option("-I", "--includedir",dest="includedir",help="Do not use this script. Only for Herwig internal use.")
(options, args) = parser.parse_args()
#parser = argparse.ArgumentParser()
#parser.add_argument('--buildpath', help='installpath')
#parser.add_argument('--build', help='build', action="store_true")
#parser.add_argument('--madgraph', help='madgraph_installpath')
#parser.add_argument('--runpath', help='runpath')
#parser.add_argument('--model', help='model')
#parser.add_argument('--orderas', help='orderas')
#parser.add_argument('--orderew', help='orderew')
#parser.add_argument('--datadir', help='datadir')
#args = parser.parse_args()
pwd=os.getcwd()
param_card=""
mkdir_p(pwd+"/Herwig-scratch/MG_tmp/")
if options.model=="loop_sm" or options.model=="heft":
if options.model=="loop_sm":
param_card="param_card.dat"
else:
param_card="param_card_"+options.model+".dat"
file = open("%s/MadGraphInterface/%s.in" % (options.datadir,param_card) , "r")
paramcard = file.read()
file.close()
file = open(options.runpath+"/"+param_card, "w")
params=open(options.runpath+"/MG-Parameter.dat", "r")
for line in params:
a=line.split()
paramcard=paramcard.replace(a[0],a[1])
params.close()
file.write(paramcard)
file.close()
+elif options.model.startswith("/"):
+
+ os.system("python %s/write_param_card.py " % options.model)
+
+
+
else:
print "---------------------------------------------------------------"
print "---------------------------------------------------------------"
print "Warning: The model set for the MadGraph Interface "
print " needs a parameter setting by hand."
print " Please fill the param_card_"+options.model+".dat"
print " with your favourite assumptions."
print " And make sure Herwig uses the same parameters."
print "---------------------------------------------------------------"
print "---------------------------------------------------------------"
if os.path.isfile(options.buildpath +"/MG5/Cards/param_card.dat") and not os.path.isfile(options.runpath+"/"+"param_card_"+options.model+".dat"):
shutil.copyfile(options.buildpath +"/MG5/Cards/param_card.dat", options.runpath+"/"+"param_card_"+options.model+".dat")
time.sleep(1)
if not os.path.isdir(options.buildpath):
print "The MadGraph Install path was not existend. It has been created for you."
print "Just start Herwig read again.."
mkdir_p(options.buildpath)
exit()
os.chdir(options.buildpath)
if os.path.isfile("InterfaceMadGraph.so"):
build_matchbox_tmp(pwd,options.buildpath)
exit()
Bornlist,Virtlist=fillprocs(options.model,options.orderas,options.orderew)
if not options.madgraph and not os.path.isfile("InterfaceMadGraph.so"):
print "*** warning *** MadGraph build failed, check logfile for details"
exit()
os.system("python "+options.madgraph+"/mg5_aMC proc.dat")
routines=[["","BORN(momenta,hel)"],
["","SLOOPMATRIX(momenta,virt)"],
["","GET_JAMP(color,Jamp)"],
["","GET_LNJAMP(color,Jamp)"],
["","GET_NCOL(color)"],
["","GET_NCOLOR(i,j,color)"]]
for routine in routines:
for i in Bornlist + list(set(Virtlist) - set(Bornlist)):
if routine[1]=="Virt(amp)" or routine[1]=="SLOOPMATRIX(momenta,virt)" and i not in Virtlist:
continue
if routine[0]=="":
routine[0]+=" SELECT CASE (proc) \n"
routine[0]+=" CASE("+i+") \n CALL "
routine[0]+= "MG5_"+i+"_"+routine[1]+"\n"
else:
routine[0]+=" CASE("+i+") \n"\
" CALL "
routine[0]+= "MG5_"+i+"_"+routine[1]+"\n"
if routine[0]!="":
routine[0]+=" CASE DEFAULT\n"
routine[0]+=" WRITE(*,*) '##W02A WARNING No id found '\n"
routine[0]+=" END SELECT \n"
shutil.copyfile("%s/MadGraphInterface/InterfaceMadGraph.f.in" % options.datadir, "InterfaceMadGraph.f")
replacetext("InterfaceMadGraph.f","MG_CalculateBORNtxt",routines[0][0])
replacetext("InterfaceMadGraph.f","MG_CalculateVIRTtxt",routines[1][0])
replacetext("InterfaceMadGraph.f","MG_Jamptxt", routines[2][0])
replacetext("InterfaceMadGraph.f","MG_LNJamptxt", routines[3][0])
replacetext("InterfaceMadGraph.f","MG_NColtxt", routines[4][0])
replacetext("InterfaceMadGraph.f","MG_ColourMattxt",routines[5][0])
MG_vxxxxxtxt=""
if routines[1][0]!="":
MG_vxxxxxtxt=""" subroutine MG_vxxxxx(p, n,inc,VC)
$ bind(c, name='MG_vxxxxx')
IMPLICIT NONE
double precision p(0:3)
double precision n(0:3)
INTEGER inc
double precision VC(0:7)
double complex VCtmp(8)
call vxxxxx(p, 0d0,1,inc ,VCtmp)
VC(0)= real(VCtmp(5))
VC(1)=aimag(VCtmp(5))
VC(2)= real(VCtmp(6))
VC(3)=aimag(VCtmp(6))
VC(4)= real(VCtmp(7))
VC(5)=aimag(VCtmp(7))
VC(6)= real(VCtmp(8))
VC(7)=aimag(VCtmp(8))
END"""
else:
MG_vxxxxxtxt=""" subroutine MG_vxxxxx(p, n,inc,VC)
$ bind(c, name='MG_vxxxxx')
IMPLICIT NONE
double precision p(0:3)
double precision n(0:3)
INTEGER inc
double precision VC(0:7)
double complex VCtmp(6)
call vxxxxx(p, 0d0,1,inc ,VCtmp)
VC(0)= real(VCtmp(3))
VC(1)=aimag(VCtmp(3))
VC(2)= real(VCtmp(4))
VC(3)=aimag(VCtmp(4))
VC(4)= real(VCtmp(5))
VC(5)=aimag(VCtmp(5))
VC(6)= real(VCtmp(6))
VC(7)=aimag(VCtmp(6))
END"""
replacetext("InterfaceMadGraph.f","MG_vxxxxxtxt",MG_vxxxxxtxt)
make=" "
fortanfiles=glob.glob('*/*/*.f')+glob.glob('*/*/*/*.f')
for i in fortanfiles:
if "check_sa" not in i:
if not os.path.islink(i):
make += " "+i+"\\\n "
incfiles=glob.glob('*/*/*.inc')+glob.glob('*/*/*/*.inc')
coefdir=""
for i in incfiles:
if "nexternal.inc" in i:
coefdir+=" -I"+i.replace("nexternal.inc"," ")
file=open("makefile","w")
file.write("include MG5/Source/make_opts ")
if Virtlist!=[]:
file.write("\nLIBDIR = MG5/lib\nLINKLIBS = -L$(LIBDIR) -lcts -liregi -L$(LIBDIR)/golem95_lib -lgolem")
file.write("\nLIBS = $(LIBDIR)libcts.$(libext) $(LIBDIR)libgolem.$(libext) $(LIBDIR)libiregi.$(libext)")
file.write("\nPROCESS= InterfaceMadGraph.f "+make+"\n\nall: \n\t @FC@ @FFLAGS@ -w -fbounds-check -ffixed-line-length-132 -fPIC -fno-f2c -shared -s -o InterfaceMadGraph.so -IMG5/SubProcesses/" )
if Virtlist!=[]:
file.write(" -IMG5/lib/golem95_include ")
# Find all .mod files also in /usr/include if golem was build there.
# There can be an error message in the MadGraph output to add the golem include path to the makefiles.
# Usually MadGraph finds the path if its Golem was build in an separate dictionary.
# Our bootstrap script installs golem with gosam beside boost. Here MadGraph creates a link (->errormessage).
# If we can find the modfiles easily the user doesn't need to change the makefiles.
moddirs=finddirs('*.mod',options.includedir)
for moddir in moddirs:
file.write(" -I%s " % moddir)
if os.path.isdir("/usr/include"):
moddirs=finddirs('*.mod',"/usr/include")
for moddir in moddirs:
file.write(" -I%s " % moddir)
if coefdir != "":
file.write(coefdir)
file.write(" $(PROCESS) $(LINKLIBS) ")
file.close()
os.chdir(pwd)
os.chdir(options.buildpath)
replacetext("MG5/Source/MODEL/lha_read.f", "ident_card.dat","Herwig-scratch/MG_tmp/ident_card.dat")
replacetext("MG5/Source/MODEL/lha_read.f", "param.log","Herwig-scratch/MG_tmp/param.log")
if Virtlist!=[]:
replacetext("MG5/SubProcesses/MadLoopCommons.f", "PREFIX='./'","PREFIX='./Herwig-scratch/MG_tmp/'")
os.system("make")
build_matchbox_tmp(pwd,options.buildpath)
diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc
--- a/Shower/Base/Evolver.cc
+++ b/Shower/Base/Evolver.cc
@@ -1,3222 +1,3222 @@
// -*- C++ -*-
//
// Evolver.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Evolver class.
//
#include "Evolver.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/Base/ShowerParticle.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ShowerKinematics.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Utilities/Throw.h"
#include "ShowerTree.h"
#include "ShowerProgenitor.h"
#include "KinematicsReconstructor.h"
#include "PartnerFinder.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ShowerVertex.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "ThePEG/Handlers/StandardXComb.h"
using namespace Herwig;
namespace {
/**
* A struct to order the particles in the same way as in the DecayMode's
*/
struct ParticleOrdering {
/**
* Operator for the ordering
* @param p1 The first ParticleData object
* @param p2 The second ParticleData object
*/
bool operator() (tcPDPtr p1, tcPDPtr p2) {
return abs(p1->id()) > abs(p2->id()) ||
( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
}
};
typedef multiset<tcPDPtr,ParticleOrdering> OrderedParticles;
/**
* Cached lookup of decay modes.
* Generator::findDecayMode() is not efficient.
*/
tDMPtr findDecayMode(const string & tag) {
static map<string,DMPtr> cache;
map<string,DMPtr>::const_iterator pos = cache.find(tag);
if ( pos != cache.end() )
return pos->second;
tDMPtr dm = CurrentGenerator::current().findDecayMode(tag);
cache[tag] = dm;
return dm;
}
}
DescribeClass<Evolver,Interfaced>
describeEvolver ("Herwig::Evolver","HwShower.so");
bool Evolver::_hardEmissionModeWarn = true;
bool Evolver::_missingTruncWarn = true;
IBPtr Evolver::clone() const {
return new_ptr(*this);
}
IBPtr Evolver::fullclone() const {
return new_ptr(*this);
}
void Evolver::persistentOutput(PersistentOStream & os) const {
os << _model << _splittingGenerator << _maxtry
<< _meCorrMode << _hardVetoMode << _hardVetoRead << _hardVetoReadOption
<< _limitEmissions << _spinOpt << _softOpt << _hardPOWHEG
<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
<< _vetoes << _trunc_Mode << _hardEmissionMode << _reconOpt
<< _massVetoOption << isMCatNLOSEvent << isMCatNLOHEvent
<< isPowhegSEvent << isPowhegHEvent
<< theFactorizationScaleFactor << theRenormalizationScaleFactor << ounit(muPt,GeV)
<< interaction_ << _maxTryFSR << _maxFailFSR << _fracFSR << interactions_.size();
for(unsigned int ix=0;ix<interactions_.size();++ix)
os << oenum(interactions_[ix]);
}
void Evolver::persistentInput(PersistentIStream & is, int) {
unsigned int isize;
is >> _model >> _splittingGenerator >> _maxtry
>> _meCorrMode >> _hardVetoMode >> _hardVetoRead >> _hardVetoReadOption
>> _limitEmissions >> _spinOpt >> _softOpt >> _hardPOWHEG
>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
>> _vetoes >> _trunc_Mode >> _hardEmissionMode >> _reconOpt
>> _massVetoOption >> isMCatNLOSEvent >> isMCatNLOHEvent
>> isPowhegSEvent >> isPowhegHEvent
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor >> iunit(muPt,GeV)
>> interaction_ >> _maxTryFSR >> _maxFailFSR >> _fracFSR >> isize;
interactions_.resize(isize);
for(unsigned int ix=0;ix<interactions_.size();++ix)
is >> ienum(interactions_[ix]);
}
void Evolver::doinit() {
Interfaced::doinit();
// interactions may have been changed through a setup file so we
// clear it up here
interactions_.clear();
if(interaction_==0) {
interactions_.push_back(ShowerInteraction::QCD);
interactions_.push_back(ShowerInteraction::QED);
}
else if(interaction_==1) {
interactions_.push_back(ShowerInteraction::QCD);
}
else if(interaction_==2) {
interactions_.push_back(ShowerInteraction::QED);
interactions_.push_back(ShowerInteraction::QCD);
}
else if(interaction_==3) {
interactions_.push_back(ShowerInteraction::QED);
}
else if(interaction_==4) {
interactions_.push_back(ShowerInteraction::Both);
}
// calculate max no of FSR vetos
- _maxFailFSR = max(_maxFailFSR,int(_fracFSR*double(generator()->N())));
+ _maxFailFSR = max(int(_maxFailFSR), int(_fracFSR*double(generator()->N())));
}
void Evolver::Init() {
static ClassDocumentation<Evolver> documentation
("This class is responsible for carrying out the showering,",
"including the kinematics reconstruction, in a given scale range,"
"including the option of the POWHEG approach to simulated next-to-leading order"
" radiation\\cite{Nason:2004rx}.",
"%\\cite{Nason:2004rx}\n"
"\\bibitem{Nason:2004rx}\n"
" P.~Nason,\n"
" ``A new method for combining NLO QCD with shower Monte Carlo algorithms,''\n"
" JHEP {\\bf 0411} (2004) 040\n"
" [arXiv:hep-ph/0409146].\n"
" %%CITATION = JHEPA,0411,040;%%\n");
static Reference<Evolver,SplittingGenerator>
interfaceSplitGen("SplittingGenerator",
"A reference to the SplittingGenerator object",
&Herwig::Evolver::_splittingGenerator,
false, false, true, false);
static Reference<Evolver,ShowerModel> interfaceShowerModel
("ShowerModel",
"The pointer to the object which defines the shower evolution model.",
&Evolver::_model, false, false, true, false, false);
static Parameter<Evolver,unsigned int> interfaceMaxTry
("MaxTry",
"The maximum number of attempts to generate the shower from a"
" particular ShowerTree",
&Evolver::_maxtry, 100, 1, 1000,
false, false, Interface::limited);
static Switch<Evolver, unsigned int> ifaceMECorrMode
("MECorrMode",
"Choice of the ME Correction Mode",
&Evolver::_meCorrMode, 1, false, false);
static SwitchOption off
(ifaceMECorrMode,"No","MECorrections off", 0);
static SwitchOption on
(ifaceMECorrMode,"Yes","hard+soft on", 1);
static SwitchOption hard
(ifaceMECorrMode,"Hard","only hard on", 2);
static SwitchOption soft
(ifaceMECorrMode,"Soft","only soft on", 3);
static Switch<Evolver, unsigned int> ifaceHardVetoMode
("HardVetoMode",
"Choice of the Hard Veto Mode",
&Evolver::_hardVetoMode, 1, false, false);
static SwitchOption HVoff
(ifaceHardVetoMode,"No","hard vetos off", 0);
static SwitchOption HVon
(ifaceHardVetoMode,"Yes","hard vetos on", 1);
static SwitchOption HVIS
(ifaceHardVetoMode,"Initial", "only IS emissions vetoed", 2);
static SwitchOption HVFS
(ifaceHardVetoMode,"Final","only FS emissions vetoed", 3);
static Switch<Evolver, unsigned int> ifaceHardVetoRead
("HardVetoScaleSource",
"If hard veto scale is to be read",
&Evolver::_hardVetoRead, 0, false, false);
static SwitchOption HVRcalc
(ifaceHardVetoRead,"Calculate","Calculate from hard process", 0);
static SwitchOption HVRread
(ifaceHardVetoRead,"Read","Read from XComb->lastScale", 1);
static Switch<Evolver, bool> ifaceHardVetoReadOption
("HardVetoReadOption",
"Apply read-in scale veto to all collisions or just the primary one?",
&Evolver::_hardVetoReadOption, false, false, false);
static SwitchOption AllCollisions
(ifaceHardVetoReadOption,
"AllCollisions",
"Read-in pT veto applied to primary and secondary collisions.",
false);
static SwitchOption PrimaryCollision
(ifaceHardVetoReadOption,
"PrimaryCollision",
"Read-in pT veto applied to primary but not secondary collisions.",
true);
static Parameter<Evolver, Energy> ifaceiptrms
("IntrinsicPtGaussian",
"RMS of intrinsic pT of Gaussian distribution:\n"
"2*(1-Beta)*exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
&Evolver::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV,
false, false, Interface::limited);
static Parameter<Evolver, double> ifacebeta
("IntrinsicPtBeta",
"Proportion of inverse quadratic distribution in generating intrinsic pT.\n"
"(1-Beta) is the proportion of Gaussian distribution",
&Evolver::_beta, 0, 0, 1,
false, false, Interface::limited);
static Parameter<Evolver, Energy> ifacegamma
("IntrinsicPtGamma",
"Parameter for inverse quadratic:\n"
"2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT))",
&Evolver::_gamma,GeV, ZERO, ZERO, 100000.0*GeV,
false, false, Interface::limited);
static Parameter<Evolver, Energy> ifaceiptmax
("IntrinsicPtIptmax",
"Upper bound on intrinsic pT for inverse quadratic",
&Evolver::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV,
false, false, Interface::limited);
static RefVector<Evolver,ShowerVeto> ifaceVetoes
("Vetoes",
"The vetoes to be checked during showering",
&Evolver::_vetoes, -1,
false,false,true,true,false);
static Switch<Evolver,unsigned int> interfaceLimitEmissions
("LimitEmissions",
"Limit the number and type of emissions for testing",
&Evolver::_limitEmissions, 0, false, false);
static SwitchOption interfaceLimitEmissionsNoLimit
(interfaceLimitEmissions,
"NoLimit",
"Allow an arbitrary number of emissions",
0);
static SwitchOption interfaceLimitEmissionsOneInitialStateEmission
(interfaceLimitEmissions,
"OneInitialStateEmission",
"Allow one emission in the initial state and none in the final state",
1);
static SwitchOption interfaceLimitEmissionsOneFinalStateEmission
(interfaceLimitEmissions,
"OneFinalStateEmission",
"Allow one emission in the final state and none in the initial state",
2);
static SwitchOption interfaceLimitEmissionsHardOnly
(interfaceLimitEmissions,
"HardOnly",
"Only allow radiation from the hard ME correction",
3);
static SwitchOption interfaceLimitEmissionsOneEmission
(interfaceLimitEmissions,
"OneEmission",
"Allow one emission in either the final state or initial state, but not both",
4);
static Switch<Evolver,bool> interfaceTruncMode
("TruncatedShower", "Include the truncated shower?",
&Evolver::_trunc_Mode, 1, false, false);
static SwitchOption interfaceTruncMode0
(interfaceTruncMode,"No","Truncated Shower is OFF", 0);
static SwitchOption interfaceTruncMode1
(interfaceTruncMode,"Yes","Truncated Shower is ON", 1);
static Switch<Evolver,int> interfaceHardEmissionMode
("HardEmissionMode",
"Whether to use ME corrections or POWHEG for the hardest emission",
&Evolver::_hardEmissionMode, 0, false, false);
static SwitchOption interfaceHardEmissionModeDecayMECorrection
(interfaceHardEmissionMode,
"DecayMECorrection",
"Old fashioned ME correction for decays only",
-1);
static SwitchOption interfaceHardEmissionModeMECorrection
(interfaceHardEmissionMode,
"MECorrection",
"Old fashioned ME correction",
0);
static SwitchOption interfaceHardEmissionModePOWHEG
(interfaceHardEmissionMode,
"POWHEG",
"Powheg style hard emission using internal matrix elements",
1);
static SwitchOption interfaceHardEmissionModeMatchboxPOWHEG
(interfaceHardEmissionMode,
"MatchboxPOWHEG",
"Powheg style emission for the hard process using Matchbox",
2);
static SwitchOption interfaceHardEmissionModeFullPOWHEG
(interfaceHardEmissionMode,
"FullPOWHEG",
"Powheg style emission for the hard process using Matchbox "
"and decays using internal matrix elements",
3);
static Switch<Evolver,unsigned int > interfaceInteractions
("Interactions",
"The interactions to be used in the shower",
&Evolver::interaction_, 1, false, false);
static SwitchOption interfaceInteractionsQCDFirst
(interfaceInteractions,
"QCDFirst",
"QCD first then QED",
0);
static SwitchOption interfaceInteractionsQCDOnly
(interfaceInteractions,
"QCDOnly",
"Only QCD",
1);
static SwitchOption interfaceInteractionsQEDFirst
(interfaceInteractions,
"QEDFirst",
"QED first then QCD",
2);
static SwitchOption interfaceInteractionsQEDOnly
(interfaceInteractions,
"QEDOnly",
"Only QED",
3);
static SwitchOption interfaceInteractionsBothAtOnce
(interfaceInteractions,
"BothAtOnce",
"Generate both at the same time",
4);
static Switch<Evolver,unsigned int> interfaceReconstructionOption
("ReconstructionOption",
"Treatment of the reconstruction of the transverse momentum of "
"a branching from the evolution scale.",
&Evolver::_reconOpt, 0, false, false);
static SwitchOption interfaceReconstructionOptionCutOff
(interfaceReconstructionOption,
"CutOff",
"Use the cut-off masses in the calculation",
0);
static SwitchOption interfaceReconstructionOptionOffShell
(interfaceReconstructionOption,
"OffShell",
"Use the off-shell masses in the calculation",
1);
static SwitchOption interfaceReconstructionOptionOffShell2
(interfaceReconstructionOption,
"OffShell2",
"Use the off-shell masses in the calculation but only locally for each branching",
2);
static Switch<Evolver,unsigned int> interfaceMassVetoOption
("MassVetoOption",
"Option for the handling of the mass vetos",
&Evolver::_massVetoOption, 1, false, false);
static SwitchOption interfaceMassVetoOptionReset
(interfaceMassVetoOption,
"Reset",
"Try another branching without resetting the starting scale",
0);
static SwitchOption interfaceMassVetoOptionInclude
(interfaceMassVetoOption,
"Include",
"Include the veto in the scale generation via the veto algorithm",
1);
static Switch<Evolver,unsigned int> interfaceSpinCorrelations
("SpinCorrelations",
"Treatment of spin correlations in the parton shower",
&Evolver::_spinOpt, 1, false, false);
static SwitchOption interfaceSpinCorrelationsOff
(interfaceSpinCorrelations,
"No",
"No spin correlations",
0);
static SwitchOption interfaceSpinCorrelationsSpin
(interfaceSpinCorrelations,
"Yes",
"Include the azimuthal spin correlations only",
1);
static Switch<Evolver,unsigned int> interfaceSoftCorrelations
("SoftCorrelations",
"Option for the treatment of soft correlations in the parton shower",
&Evolver::_softOpt, 2, false, false);
static SwitchOption interfaceSoftCorrelationsNone
(interfaceSoftCorrelations,
"No",
"No soft correlations",
0);
static SwitchOption interfaceSoftCorrelationsFull
(interfaceSoftCorrelations,
"Full",
"Use the full eikonal",
1);
static SwitchOption interfaceSoftCorrelationsSingular
(interfaceSoftCorrelations,
"Singular",
"Use original Webber-Marchisini form",
2);
static Switch<Evolver,bool> interfaceHardPOWHEG
("HardPOWHEG",
"Treatment of powheg emissions which are too hard to have a shower interpretation",
&Evolver::_hardPOWHEG, false, false, false);
static SwitchOption interfaceHardPOWHEGAsShower
(interfaceHardPOWHEG,
"AsShower",
"Still interpret as shower emissions",
false);
static SwitchOption interfaceHardPOWHEGRealEmission
(interfaceHardPOWHEG,
"RealEmission",
"Generate shower from the real emmission configuration",
true);
static Parameter<Evolver,unsigned int> interfaceMaxTryFSR
("MaxTryFSR",
"The maximum number of attempted FSR emissions in"
" the generation of the FSR",
&Evolver::_maxTryFSR, 100000, 10, 100000000,
false, false, Interface::limited);
static Parameter<Evolver,unsigned int> interfaceMaxFailFSR
("MaxFailFSR",
"Maximum number of failures generating the FSR",
&Evolver::_maxFailFSR, 100, 1, 100000000,
false, false, Interface::limited);
static Parameter<Evolver,double> interfaceFSRFailureFraction
("FSRFailureFraction",
"Maximum fraction of events allowed to fail due to too many FSR emissions",
&Evolver::_fracFSR, 0.001, 1e-10, 1,
false, false, Interface::limited);
}
void Evolver::generateIntrinsicpT(vector<ShowerProgenitorPtr> particlesToShower) {
_intrinsic.clear();
if ( !ipTon() || !isISRadiationON() ) return;
// don't do anything for the moment for secondary scatters
if( !ShowerHandler::currentHandler()->firstInteraction() ) return;
// generate intrinsic pT
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
// only consider initial-state particles
if(particlesToShower[ix]->progenitor()->isFinalState()) continue;
if(!particlesToShower[ix]->progenitor()->dataPtr()->coloured()) continue;
Energy ipt;
if(UseRandom::rnd() > _beta) {
ipt=_iptrms*sqrt(-log(UseRandom::rnd()));
}
else {
ipt=_gamma*sqrt(pow(1.+sqr(_iptmax/_gamma), UseRandom::rnd())-1.);
}
pair<Energy,double> pt = make_pair(ipt,UseRandom::rnd(Constants::twopi));
_intrinsic[particlesToShower[ix]] = pt;
}
}
void Evolver::setupMaximumScales(const vector<ShowerProgenitorPtr> & p,
XCPtr xcomb) {
// let POWHEG events radiate freely
if(_hardEmissionMode>0&&hardTree()) {
vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(Constants::MaxEnergy);
return;
}
// return if no vetos
if (!hardVetoOn()) return;
// find out if hard partonic subprocess.
bool isPartonic(false);
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit = _currenttree->incomingLines().begin();
Lorentz5Momentum pcm;
for(; cit!=currentTree()->incomingLines().end(); ++cit) {
pcm += cit->first->progenitor()->momentum();
isPartonic |= cit->first->progenitor()->coloured();
}
// find minimum pt from hard process, the maximum pt from all outgoing
// coloured lines (this is simpler and more general than
// 2stu/(s^2+t^2+u^2)). Maximum scale for scattering processes will
// be transverse mass.
Energy ptmax = generator()->maximumCMEnergy();
// general case calculate the scale
if (!hardVetoXComb()||
(hardVetoReadOption()&&
!ShowerHandler::currentHandler()->firstInteraction())) {
// scattering process
if(currentTree()->isHard()) {
assert(xcomb);
// coloured incoming particles
if (isPartonic) {
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cjt = currentTree()->outgoingLines().begin();
for(; cjt!=currentTree()->outgoingLines().end(); ++cjt) {
if (cjt->first->progenitor()->coloured())
ptmax = min(ptmax,cjt->first->progenitor()->momentum().mt());
}
}
if (ptmax == generator()->maximumCMEnergy() ) ptmax = pcm.m();
if(hardVetoXComb()&&hardVetoReadOption()&&
!ShowerHandler::currentHandler()->firstInteraction()) {
ptmax=min(ptmax,sqrt(xcomb->lastCentralScale()));
}
}
// decay, incoming() is the decaying particle.
else {
ptmax = currentTree()->incomingLines().begin()->first
->progenitor()->momentum().mass();
}
}
// hepeup.SCALUP is written into the lastXComb by the
// LesHouchesReader itself - use this by user's choice.
// Can be more general than this.
else {
if(currentTree()->isHard()) {
assert(xcomb);
ptmax = sqrt( xcomb->lastCentralScale() );
}
else {
ptmax = currentTree()->incomingLines().begin()->first
->progenitor()->momentum().mass();
}
}
ptmax *= ShowerHandler::currentHandler()->hardScaleFactor();
// set maxHardPt for all progenitors. For partonic processes this
// is now the max pt in the FS, for non-partonic processes or
// processes with no coloured FS the invariant mass of the IS
vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(ptmax);
}
void Evolver::setupHardScales(const vector<ShowerProgenitorPtr> & p,
XCPtr xcomb) {
if ( hardVetoXComb() &&
(!hardVetoReadOption() ||
ShowerHandler::currentHandler()->firstInteraction()) ) {
Energy hardScale = ZERO;
if(currentTree()->isHard()) {
assert(xcomb);
hardScale = sqrt( xcomb->lastCentralScale() );
}
else {
hardScale = currentTree()->incomingLines().begin()->first
->progenitor()->momentum().mass();
}
hardScale *= ShowerHandler::currentHandler()->hardScaleFactor();
vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
for (; ckt != p.end(); ckt++) (*ckt)->hardScale(hardScale);
muPt = hardScale;
}
}
void Evolver::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {
isMCatNLOSEvent = false;
isMCatNLOHEvent = false;
isPowhegSEvent = false;
isPowhegHEvent = false;
Ptr<SubtractedME>::tptr subme;
Ptr<MatchboxMEBase>::tptr me;
Ptr<SubtractionDipole>::tptr dipme;
Ptr<StandardXComb>::ptr sxc = dynamic_ptr_cast<Ptr<StandardXComb>::ptr>(xcomb);
if ( sxc ) {
subme = dynamic_ptr_cast<Ptr<SubtractedME>::tptr>(sxc->matrixElement());
me = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(sxc->matrixElement());
dipme = dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(sxc->matrixElement());
}
if ( subme ) {
if ( subme->showerApproximation() ) {
theShowerApproximation = subme->showerApproximation();
// separate MCatNLO and POWHEG-type corrections
if ( !subme->showerApproximation()->needsSplittingGenerator() ) {
if ( subme->realShowerSubtraction() )
isMCatNLOHEvent = true;
else if ( subme->virtualShowerSubtraction() )
isMCatNLOSEvent = true;
}
else {
if ( subme->realShowerSubtraction() )
isPowhegHEvent = true;
else if ( subme->virtualShowerSubtraction() || subme->loopSimSubtraction() )
isPowhegSEvent = true;
}
}
} else if ( me ) {
if ( me->factory()->showerApproximation() ) {
theShowerApproximation = me->factory()->showerApproximation();
if ( !me->factory()->showerApproximation()->needsSplittingGenerator() )
isMCatNLOSEvent = true;
else
isPowhegSEvent = true;
}
}
string error = "Inconsistent hard emission set-up in Evolver::showerHardProcess(). ";
if ( ( isMCatNLOSEvent || isMCatNLOHEvent ) ){
if (_hardEmissionMode > 1)
throw Exception() << error
<< "Cannot generate POWHEG matching with MC@NLO shower "
<< "approximation. Add 'set Evolver:HardEmissionMode 0' to input file."
<< Exception::runerror;
if ( ShowerHandler::currentHandler()->canHandleMatchboxTrunc())
throw Exception() << error
<< "Cannot use truncated qtilde shower with MC@NLO shower "
<< "approximation. Set LHCGenerator:EventHandler"
<< ":CascadeHandler to '/Herwig/Shower/ShowerHandler' or "
<< "'/Herwig/DipoleShower/DipoleShowerHandler'."
<< Exception::runerror;
}
else if ( ((isPowhegSEvent || isPowhegHEvent) || dipme) &&
_hardEmissionMode < 2){
if ( ShowerHandler::currentHandler()->canHandleMatchboxTrunc())
throw Exception() << error
<< "Unmatched events requested for POWHEG shower "
<< "approximation. Set Evolver:HardEmissionMode to "
<< "'MatchboxPOWHEG' or 'FullPOWHEG'."
<< Exception::runerror;
else if (_hardEmissionModeWarn){
_hardEmissionModeWarn = false;
_hardEmissionMode+=2;
throw Exception() << error
<< "Unmatched events requested for POWHEG shower "
<< "approximation. Changing Evolver:HardEmissionMode from "
<< _hardEmissionMode-2 << " to " << _hardEmissionMode
<< Exception::warning;
}
}
if ( isPowhegSEvent || isPowhegHEvent) {
if (theShowerApproximation->needsTruncatedShower() &&
!ShowerHandler::currentHandler()->canHandleMatchboxTrunc() )
throw Exception() << error
<< "Current shower handler cannot generate truncated shower. "
<< "Set Generator:EventHandler:CascadeHandler to "
<< "'/Herwig/Shower/PowhegShowerHandler'."
<< Exception::runerror;
}
else if ( dipme && _missingTruncWarn){
_missingTruncWarn=false;
throw Exception() << "Warning: POWHEG shower approximation used without "
<< "truncated shower. Set Generator:EventHandler:"
<< "CascadeHandler to '/Herwig/Shower/PowhegShowerHandler' and "
<< "'MEMatching:TruncatedShower Yes'."
<< Exception::warning;
}
else if ( !dipme && _hardEmissionMode > 1 &&
ShowerHandler::currentHandler()->firstInteraction())
throw Exception() << error
<< "POWHEG matching requested for LO events. Include "
<< "'set Factory:ShowerApproximation MEMatching' in input file."
<< Exception::runerror;
_hardme = HwMEBasePtr();
// extract the matrix element
tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(xcomb);
if(lastXC) {
_hardme = dynamic_ptr_cast<HwMEBasePtr>(lastXC->matrixElement());
}
_decayme = HwDecayerBasePtr();
// set the current tree
currentTree(hard);
hardTree(HardTreePtr());
// number of attempts if more than one interaction switched on
unsigned int interactionTry=0;
do {
try {
// generate the showering
doShowering(true,xcomb);
// if no vetos return
return;
}
catch (InteractionVeto) {
currentTree()->clear();
++interactionTry;
}
}
while(interactionTry<=5);
throw Exception() << "Too many tries for shower in "
<< "Evolver::showerHardProcess()"
<< Exception::eventerror;
}
void Evolver::hardMatrixElementCorrection(bool hard) {
// set the initial enhancement factors for the soft correction
_initialenhance = 1.;
_finalenhance = 1.;
// if hard matrix element switched off return
if(!MECOn(hard)) return;
// see if we can get the correction from the matrix element
// or decayer
if(hard) {
if(_hardme&&_hardme->hasMECorrection()) {
_hardme->initializeMECorrection(_currenttree,
_initialenhance,_finalenhance);
if(hardMEC(hard))
_hardme->applyHardMatrixElementCorrection(_currenttree);
}
}
else {
if(_decayme&&_decayme->hasMECorrection()) {
_decayme->initializeMECorrection(_currenttree,
_initialenhance,_finalenhance);
if(hardMEC(hard))
_decayme->applyHardMatrixElementCorrection(_currenttree);
}
}
}
Branching Evolver::selectTimeLikeBranching(tShowerParticlePtr particle,
ShowerInteraction::Type type) {
Branching fb;
while (true) {
fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
// no emission return
if(!fb.kinematics) return fb;
// if emission OK break
if(!timeLikeVetoed(fb,particle)) break;
// otherwise reset scale and continue - SO IS involved in veto algorithm
particle->vetoEmission(fb.type,fb.kinematics->scale());
if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
}
return fb;
}
ShowerParticleVector Evolver::createTimeLikeChildren(tShowerParticlePtr particle, IdList ids) {
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
tcPDPtr pdata[2];
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(ids[ix+1]);
if(particle->id()!=ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
ShowerParticleVector children;
for(unsigned int ix=0;ix<2;++ix) {
children.push_back(new_ptr(ShowerParticle(pdata[ix],true)));
if(children[ix]->id()==_progenitor->id()&&!pdata[ix]->stable())
children[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
else
children[ix]->set5Momentum(Lorentz5Momentum(pdata[ix]->mass()));
}
return children;
}
bool Evolver::timeLikeShower(tShowerParticlePtr particle,
ShowerInteraction::Type type,
Branching fb, bool first) {
// don't do anything if not needed
if(_limitEmissions == 1 || hardOnly() ||
( _limitEmissions == 2 && _nfs != 0) ||
( _limitEmissions == 4 && _nfs + _nis != 0) ) {
if(particle->spinInfo()) particle->spinInfo()->develop();
return false;
}
// too many tries
if(_nFSR>=_maxTryFSR) {
++_nFailedFSR;
// too many failed events
if(_nFailedFSR>=_maxFailFSR)
throw Exception() << "Too many events have failed due to too many shower emissions, in\n"
<< "Evolver::timeLikeShower(). Terminating run\n"
<< Exception::runerror;
throw Exception() << "Too many attempted emissions in Evolver::timeLikeShower()\n"
<< Exception::eventerror;
}
// generate the emission
ShowerParticleVector children;
int ntry=0;
while (ntry<50) {
++ntry;
// generate the emission
if(!fb.kinematics)
fb = selectTimeLikeBranching(particle,type);
// no emission, return
if(!fb.kinematics) {
if(particle->spinInfo()) particle->spinInfo()->develop();
return false;
}
// has emitted
// Assign the shower kinematics to the emitting particle.
++_nFSR;
particle->showerKinematics(fb.kinematics);
// generate phi
fb.kinematics->phi(fb.sudakov->generatePhiForward(*particle,fb.ids,fb.kinematics));
// check highest pT
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// create the children
children = createTimeLikeChildren(particle,fb.ids);
// update the children
particle->showerKinematics()->
updateChildren(particle, children,fb.type);
// update number of emissions
++_nfs;
if(_limitEmissions!=0) {
if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
// select branchings for children
Branching fc[2] = {selectTimeLikeBranching(children[0],type),
selectTimeLikeBranching(children[1],type)};
// old recon option
if(_reconOpt==0) {
// shower the first particle
if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
// shower the second particle
if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
}
else if(_reconOpt==1) {
// shower the first particle
if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
// shower the second particle
if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
// branching has happened
particle->showerKinematics()->
updateParent(particle, children,fb.type);
// clean up the vetoed emission
if(particle->virtualMass()==ZERO) {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<children.size();++ix)
particle->abandonChild(children[ix]);
children.clear();
if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
if(_massVetoOption==1) particle->vetoEmission(fb.type,fb.kinematics->scale());
fb = Branching();
continue;
}
}
else if(_reconOpt==2) {
// cut-off masses for the branching
const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
// compute the masses of the children
Energy masses[3];
for(unsigned int ix=0;ix<2;++ix) {
if(fc[ix].kinematics) {
const vector<Energy> & vm = fc[ix].sudakov->virtualMasses(fc[ix].ids);
Energy2 q2 =
fc[ix].kinematics->z()*(1.-fc[ix].kinematics->z())*sqr(fc[ix].kinematics->scale());
if(fc[ix].ids[0]!=ParticleID::g) q2 += sqr(vm[0]);
masses[ix+1] = sqrt(q2);
}
else {
masses[ix+1] = virtualMasses[ix+1];
}
}
masses[0] = fb.ids[0]!=ParticleID::g ? virtualMasses[0] : ZERO;
double z = fb.kinematics->z();
Energy2 pt2 = z*(1.-z)*(z*(1.-z)*sqr(fb.kinematics->scale())
+sqr(masses[0]))
- sqr(masses[1])*(1.-z) - sqr(masses[2])*z;
if(pt2>=ZERO) {
// branching has happened
particle->showerKinematics()->
updateParent(particle, children,fb.type);
// shower the first particle
if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
// shower the second particle
if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
}
else {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<children.size();++ix)
particle->abandonChild(children[ix]);
children.clear();
if(_massVetoOption==1) particle->vetoEmission(fb.type,fb.kinematics->scale());
if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
fb = Branching();
continue;
}
}
break;
};
if(first&&!children.empty())
particle->showerKinematics()->resetChildren(particle,children);
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
bool
Evolver::spaceLikeShower(tShowerParticlePtr particle, PPtr beam,
ShowerInteraction::Type type) {
//using the pdf's associated with the ShowerHandler assures, that
//modified pdf's are used for the secondary interactions via
//CascadeHandler::resetPDFs(...)
tcPDFPtr pdf;
if(ShowerHandler::currentHandler()->firstPDF().particle() == _beam)
pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
if(ShowerHandler::currentHandler()->secondPDF().particle() == _beam)
pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
// don't do anything if not needed
if(_limitEmissions == 2 || hardOnly() ||
( _limitEmissions == 1 && _nis != 0 ) ||
( _limitEmissions == 4 && _nis + _nfs != 0 ) ) {
if(particle->spinInfo()) particle->spinInfo()->develop();
return false;
}
Branching bb;
// generate branching
while (true) {
bb=_splittingGenerator->chooseBackwardBranching(*particle,beam,
_initialenhance,
_beam,type,
pdf,freeze);
// return if no emission
if(!bb.kinematics) {
if(particle->spinInfo()) particle->spinInfo()->develop();
return false;
}
// if not vetoed break
if(!spaceLikeVetoed(bb,particle)) break;
// otherwise reset scale and continue
particle->vetoEmission(bb.type,bb.kinematics->scale());
if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
}
// assign the splitting function and shower kinematics
particle->showerKinematics(bb.kinematics);
if(bb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(bb.kinematics->pT());
// For the time being we are considering only 1->2 branching
// particles as in Sudakov form factor
tcPDPtr part[2]={getParticleData(bb.ids[0]),
getParticleData(bb.ids[2])};
if(particle->id()!=bb.ids[1]) {
if(part[0]->CC()) part[0]=part[0]->CC();
if(part[1]->CC()) part[1]=part[1]->CC();
}
// Now create the actual particles, make the otherChild a final state
// particle, while the newParent is not
ShowerParticlePtr newParent=new_ptr(ShowerParticle(part[0],false));
ShowerParticlePtr otherChild = new_ptr(ShowerParticle(part[1],true,true));
ShowerParticleVector theChildren;
theChildren.push_back(particle);
theChildren.push_back(otherChild);
//this updates the evolution scale
particle->showerKinematics()->
updateParent(newParent, theChildren,bb.type);
// update the history if needed
_currenttree->updateInitialStateShowerProduct(_progenitor,newParent);
_currenttree->addInitialStateBranching(particle,newParent,otherChild);
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
++_nis;
bool emitted = _limitEmissions==0 ?
spaceLikeShower(newParent,beam,type) : false;
if(newParent->spinInfo()) newParent->spinInfo()->develop();
// now reconstruct the momentum
if(!emitted) {
if(_intrinsic.find(_progenitor)==_intrinsic.end()) {
bb.kinematics->updateLast(newParent,ZERO,ZERO);
}
else {
pair<Energy,double> kt=_intrinsic[_progenitor];
bb.kinematics->updateLast(newParent,
kt.first*cos(kt.second),
kt.first*sin(kt.second));
}
}
particle->showerKinematics()->
updateChildren(newParent, theChildren,bb.type);
if(_limitEmissions!=0) {
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
// perform the shower of the final-state particle
timeLikeShower(otherChild,type,Branching(),true);
updateHistory(otherChild);
if(theChildren[1]->spinInfo()) theChildren[1]->spinInfo()->develop();
// return the emitted
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
void Evolver::showerDecay(ShowerTreePtr decay) {
_decayme = HwDecayerBasePtr();
_hardme = HwMEBasePtr();
// find the decayer
// try the normal way if possible
tDMPtr dm = decay->incomingLines().begin()->first->original() ->decayMode();
if(!dm) dm = decay->incomingLines().begin()->first->copy() ->decayMode();
if(!dm) dm = decay->incomingLines().begin()->first->progenitor()->decayMode();
// otherwise make a string and look it up
if(!dm) {
string tag = decay->incomingLines().begin()->first->original()->dataPtr()->name()
+ "->";
OrderedParticles outgoing;
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
if(abs(decay->incomingLines().begin()->first->original()->id()) == ParticleID::t &&
abs(it->first->original()->id())==ParticleID::Wplus &&
decay->treelinks().size() == 1) {
ShowerTreePtr Wtree = decay->treelinks().begin()->first;
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
it2=Wtree->outgoingLines().begin();it2!=Wtree->outgoingLines().end();++it2) {
outgoing.insert(it2->first->original()->dataPtr());
}
}
else {
outgoing.insert(it->first->original()->dataPtr());
}
}
for(OrderedParticles::const_iterator it=outgoing.begin(); it!=outgoing.end();++it) {
if(it!=outgoing.begin()) tag += ",";
tag +=(**it).name();
}
tag += ";";
dm = findDecayMode(tag);
}
if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
// set the ShowerTree to be showered
currentTree(decay);
decay->applyTransforms();
hardTree(HardTreePtr());
unsigned int interactionTry=0;
do {
try {
// generate the showering
doShowering(false,XCPtr());
// if no vetos
// force calculation of spin correlations
SpinPtr spInfo = decay->incomingLines().begin()->first->progenitor()->spinInfo();
if(spInfo) {
if(!spInfo->developed()) spInfo->needsUpdate();
spInfo->develop();
}
// and then return
return;
}
catch (InteractionVeto) {
currentTree()->clear();
++interactionTry;
}
}
while(interactionTry<=5);
throw Exception() << "Too many tries for QED shower in Evolver::showerDecay()"
<< Exception::eventerror;
}
bool Evolver::spaceLikeDecayShower(tShowerParticlePtr particle,
const ShowerParticle::EvolutionScales & maxScales,
Energy minmass,ShowerInteraction::Type type) {
Branching fb;
while (true) {
fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
_initialenhance,type);
// return if no radiation
if(!fb.kinematics) return false;
// if not vetoed break
if(!spaceLikeDecayVetoed(fb,particle)) break;
// otherwise reset scale and continue
particle->vetoEmission(fb.type,fb.kinematics->scale());
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
tcPDPtr pdata[2];
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
ShowerParticleVector theChildren;
for(unsigned int ix=0;ix<2;++ix) {
theChildren.push_back(new_ptr(ShowerParticle(pdata[ix],true)));
if(theChildren[ix]->id()==_progenitor->id()&&!pdata[ix]->stable())
theChildren[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
else
theChildren[ix]->set5Momentum(Lorentz5Momentum(pdata[ix]->mass()));
}
// some code moved to updateChildren
particle->showerKinematics()->
updateChildren(particle, theChildren, fb.type);
// In the case of splittings which involves coloured particles,
// set properly the colour flow of the branching.
// update the history if needed
_currenttree->updateInitialStateShowerProduct(_progenitor,theChildren[0]);
_currenttree->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the first particle
spaceLikeDecayShower(theChildren[0],maxScales,minmass,type);
// shower the second particle
timeLikeShower(theChildren[1],type,Branching(),true);
updateHistory(theChildren[1]);
// branching has happened
return true;
}
vector<ShowerProgenitorPtr> Evolver::setupShower(bool hard) {
// generate POWHEG hard emission if needed
if(_hardEmissionMode>0) hardestEmission(hard);
ShowerInteraction::Type inter = interactions_[0];
if(_hardtree&&inter!=ShowerInteraction::Both) {
inter = _hardtree->interaction();
}
// set the initial colour partners
setEvolutionPartners(hard,inter,false);
// generate hard me if needed
if(_hardEmissionMode==0 ||
(!hard && _hardEmissionMode==-1)) hardMatrixElementCorrection(hard);
// get the particles to be showered
vector<ShowerProgenitorPtr> particlesToShower =
currentTree()->extractProgenitors();
// remake the colour partners if needed
if(_currenttree->hardMatrixElementCorrection()) {
setEvolutionPartners(hard,interactions_[0],true);
_currenttree->resetShowerProducts();
}
// return the answer
return particlesToShower;
}
void Evolver::setEvolutionPartners(bool hard,ShowerInteraction::Type type,
bool clear) {
// match the particles in the ShowerTree and hardTree
if(hardTree() && !hardTree()->connect(currentTree()))
throw Exception() << "Can't match trees in "
<< "Evolver::setEvolutionPartners()"
<< Exception::eventerror;
// extract the progenitors
vector<ShowerParticlePtr> particles =
currentTree()->extractProgenitorParticles();
// clear the partners if needed
if(clear) {
for(unsigned int ix=0;ix<particles.size();++ix) {
particles[ix]->partner(ShowerParticlePtr());
particles[ix]->clearPartners();
}
}
// sort out the colour partners
if(hardTree()) {
// find the partner
for(unsigned int ix=0;ix<particles.size();++ix) {
tHardBranchingPtr partner =
hardTree()->particles()[particles[ix]]->colourPartner();
if(!partner) continue;
for(map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
it=hardTree()->particles().begin();
it!=hardTree()->particles().end();++it) {
if(it->second==partner) particles[ix]->partner(it->first);
}
if(!particles[ix]->partner())
throw Exception() << "Can't match partners in "
<< "Evolver::setEvolutionPartners()"
<< Exception::eventerror;
}
}
// Set the initial evolution scales
showerModel()->partnerFinder()->
setInitialEvolutionScales(particles,!hard,type,!_hardtree);
if(hardTree() && _hardPOWHEG) {
bool tooHard=false;
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end();
for(unsigned int ix=0;ix<particles.size();++ix) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
mit = hardTree()->particles().find(particles[ix]);
Energy hardScale(ZERO);
ShowerPartnerType::Type type(ShowerPartnerType::Undefined);
// final-state
if(particles[ix]->isFinalState()) {
if(mit!= eit && !mit->second->children().empty()) {
hardScale = mit->second->scale();
type = mit->second->type();
}
}
// initial-state
else {
if(mit!= eit && mit->second->parent()) {
hardScale = mit->second->parent()->scale();
type = mit->second->parent()->type();
}
}
if(type!=ShowerPartnerType::Undefined) {
if(type==ShowerPartnerType::QED) {
tooHard |= particles[ix]->scales().QED_noAO<hardScale;
}
else if(type==ShowerPartnerType::QCDColourLine) {
tooHard |= particles[ix]->scales().QCD_c_noAO<hardScale;
}
else if(type==ShowerPartnerType::QCDAntiColourLine) {
tooHard |= particles[ix]->scales().QCD_ac_noAO<hardScale;
}
}
}
if(tooHard) convertHardTree(hard,type);
}
}
void Evolver::updateHistory(tShowerParticlePtr particle) {
if(!particle->children().empty()) {
ShowerParticleVector theChildren;
for(unsigned int ix=0;ix<particle->children().size();++ix) {
ShowerParticlePtr part = dynamic_ptr_cast<ShowerParticlePtr>
(particle->children()[ix]);
theChildren.push_back(part);
}
// update the history if needed
if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
_currenttree->updateFinalStateShowerProduct(_progenitor,
particle,theChildren);
_currenttree->addFinalStateBranching(particle,theChildren);
for(unsigned int ix=0;ix<theChildren.size();++ix)
updateHistory(theChildren[ix]);
}
}
bool Evolver::startTimeLikeShower(ShowerInteraction::Type type) {
_nFSR = 0;
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && !mit->second->children().empty() ) {
bool output=truncatedTimeLikeShower(progenitor()->progenitor(),
mit->second ,type,true);
if(output) updateHistory(progenitor()->progenitor());
return output;
}
}
bool output = hardOnly() ? false :
timeLikeShower(progenitor()->progenitor() ,type,Branching(),true) ;
if(output) updateHistory(progenitor()->progenitor());
return output;
}
bool Evolver::startSpaceLikeShower(PPtr parent, ShowerInteraction::Type type) {
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit =hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && mit->second->parent() ) {
return truncatedSpaceLikeShower( progenitor()->progenitor(),
parent, mit->second->parent(), type );
}
}
return hardOnly() ? false :
spaceLikeShower(progenitor()->progenitor(),parent,type);
}
bool Evolver::
startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
Energy minimumMass,ShowerInteraction::Type type) {
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit =hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && mit->second->parent() ) {
HardBranchingPtr branch=mit->second;
while(branch->parent()) branch=branch->parent();
return truncatedSpaceLikeDecayShower(progenitor()->progenitor(),maxScales,
minimumMass, branch ,type);
}
}
return hardOnly() ? false :
spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type);
}
bool Evolver::timeLikeVetoed(const Branching & fb,
ShowerParticlePtr particle) {
// work out type of interaction
ShowerInteraction::Type type = fb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
// check whether emission was harder than largest pt of hard subprocess
if ( hardVetoFS() && fb.kinematics->pT() > _progenitor->maxHardPt() )
return true;
// soft matrix element correction veto
if( softMEC()) {
if(_hardme && _hardme->hasMECorrection()) {
if(_hardme->softMatrixElementVeto(_progenitor,particle,fb))
return true;
}
else if(_decayme && _decayme->hasMECorrection()) {
if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
return true;
}
}
// veto on maximum pt
if(fb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
// general vetos
if (fb.kinematics && !_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoTimeLike(_progenitor,particle,fb);
switch((**v).vetoType()) {
case ShowerVeto::Emission:
vetoed |= test;
break;
case ShowerVeto::Shower:
if(test) throw VetoShower();
break;
case ShowerVeto::Event:
if(test) throw Veto();
break;
}
}
if(vetoed) return true;
}
if ( ShowerHandler::currentHandler()->firstInteraction() &&
ShowerHandler::currentHandler()->profileScales() ) {
double weight =
ShowerHandler::currentHandler()->profileScales()->
hardScaleProfile(_progenitor->hardScale(),fb.kinematics->pT());
if ( UseRandom::rnd() > weight )
return true;
}
return false;
}
bool Evolver::spaceLikeVetoed(const Branching & bb,
ShowerParticlePtr particle) {
// work out type of interaction
ShowerInteraction::Type type = bb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
// check whether emission was harder than largest pt of hard subprocess
if (hardVetoIS() && bb.kinematics->pT() > _progenitor->maxHardPt())
return true;
// apply the soft correction
if( softMEC() && _hardme && _hardme->hasMECorrection() ) {
if(_hardme->softMatrixElementVeto(_progenitor,particle,bb))
return true;
}
// the more general vetos
// check vs max pt for the shower
if(bb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
if (!_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoSpaceLike(_progenitor,particle,bb);
switch ((**v).vetoType()) {
case ShowerVeto::Emission:
vetoed |= test;
break;
case ShowerVeto::Shower:
if(test) throw VetoShower();
break;
case ShowerVeto::Event:
if(test) throw Veto();
break;
}
}
if (vetoed) return true;
}
if ( ShowerHandler::currentHandler()->firstInteraction() &&
ShowerHandler::currentHandler()->profileScales() ) {
double weight =
ShowerHandler::currentHandler()->profileScales()->
hardScaleProfile(_progenitor->hardScale(),bb.kinematics->pT());
if ( UseRandom::rnd() > weight )
return true;
}
return false;
}
bool Evolver::spaceLikeDecayVetoed( const Branching & fb,
ShowerParticlePtr particle) {
// work out type of interaction
ShowerInteraction::Type type = fb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
// apply the soft correction
if( softMEC() && _decayme && _decayme->hasMECorrection() ) {
if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
return true;
}
// veto on hardest pt in the shower
if(fb.kinematics->pT()> _progenitor->maximumpT(type)) return true;
// general vetos
if (!_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoSpaceLike(_progenitor,particle,fb);
switch((**v).vetoType()) {
case ShowerVeto::Emission:
vetoed |= test;
break;
case ShowerVeto::Shower:
if(test) throw VetoShower();
break;
case ShowerVeto::Event:
if(test) throw Veto();
break;
}
if (vetoed) return true;
}
}
return false;
}
void Evolver::hardestEmission(bool hard) {
HardTreePtr ISRTree;
if( ( _hardme && _hardme->hasPOWHEGCorrection()!=0 && _hardEmissionMode< 2) ||
( _decayme && _decayme->hasPOWHEGCorrection()!=0 && _hardEmissionMode!=2) ) {
if(_hardme) {
assert(hard);
if(interaction_==4) {
vector<ShowerInteraction::Type> inter(2);
inter[0] = ShowerInteraction::QCD;
inter[1] = ShowerInteraction::QED;
_hardtree = _hardme->generateHardest( currentTree(),inter );
}
else {
_hardtree = _hardme->generateHardest( currentTree(),interactions_ );
}
}
else {
assert(!hard);
_hardtree = _decayme->generateHardest( currentTree() );
}
// store initial state POWHEG radiation
if(_hardtree && _hardme && _hardme->hasPOWHEGCorrection()==1)
ISRTree=_hardtree;
}
else if (_hardEmissionMode>1 && hard) {
// Get minimum pT cutoff used in shower approximation
Energy maxpt = 1.*GeV;
int colouredIn = 0;
int colouredOut = 0;
for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
= currentTree()->outgoingLines().begin();
it != currentTree()->outgoingLines().end(); ++it ) {
if( it->second->coloured() ) colouredOut+=1;
}
for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
= currentTree()->incomingLines().begin();
it != currentTree()->incomingLines().end(); ++it ) {
if( ! it->second->coloured() ) colouredIn+=1;
}
if ( theShowerApproximation ){
if ( theShowerApproximation->ffPtCut() == theShowerApproximation->fiPtCut() &&
theShowerApproximation->ffPtCut() == theShowerApproximation->iiPtCut() )
maxpt = theShowerApproximation->ffPtCut();
else if ( colouredIn == 2 && colouredOut == 0 )
maxpt = theShowerApproximation->iiPtCut();
else if ( colouredIn == 0 && colouredOut > 1 )
maxpt = theShowerApproximation->ffPtCut();
else if ( colouredIn == 2 && colouredOut == 1 )
maxpt = min(theShowerApproximation->iiPtCut(), theShowerApproximation->fiPtCut());
else if ( colouredIn == 1 && colouredOut > 1 )
maxpt = min(theShowerApproximation->ffPtCut(), theShowerApproximation->fiPtCut());
else
maxpt = min(min(theShowerApproximation->iiPtCut(), theShowerApproximation->fiPtCut()),
theShowerApproximation->ffPtCut());
}
// Generate hardtree from born and real emission subprocesses
_hardtree = ShowerHandler::currentHandler()->generateCKKW(currentTree());
// Find transverse momentum of hardest emission
if (_hardtree){
for(set<HardBranchingPtr>::iterator it=_hardtree->branchings().begin();
it!=_hardtree->branchings().end();++it) {
if ((*it)->parent() && (*it)->status()==HardBranching::Incoming)
maxpt=(*it)->branchingParticle()->momentum().perp();
if ((*it)->children().size()==2 && (*it)->status()==HardBranching::Outgoing){
if ((*it)->branchingParticle()->id()!=21 &&
abs((*it)->branchingParticle()->id())>5 ){
if ((*it)->children()[0]->branchingParticle()->id()==21 ||
abs((*it)->children()[0]->branchingParticle()->id())<6)
maxpt=(*it)->children()[0]->branchingParticle()->momentum().perp();
else if ((*it)->children()[1]->branchingParticle()->id()==21 ||
abs((*it)->children()[1]->branchingParticle()->id())<6)
maxpt=(*it)->children()[1]->branchingParticle()->momentum().perp();
}
else {
if ( abs((*it)->branchingParticle()->id())<6){
if (abs((*it)->children()[0]->branchingParticle()->id())<6)
maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
else
maxpt = (*it)->children()[0]->branchingParticle()->momentum().perp();
}
else maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
}
}
}
}
// Hardest (pt) emission should be the first powheg emission.
maxpt=min(sqrt(ShowerHandler::currentHandler()->lastXCombPtr()->lastCentralScale()),maxpt);
// Set maxpt to pT of emission when showering POWHEG real-emission subprocesses
if (!isPowhegSEvent && !isPowhegHEvent){
vector<int> outGluon;
vector<int> outQuark;
map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it;
for( it = currentTree()->outgoingLines().begin();
it != currentTree()->outgoingLines().end(); ++it ) {
if ( abs(it->second->id())< 6) outQuark.push_back(it->second->id());
if ( it->second->id()==21 ) outGluon.push_back(it->second->id());
}
if (outGluon.size() + outQuark.size() == 1){
for( it = currentTree()->outgoingLines().begin();
it != currentTree()->outgoingLines().end(); ++it ) {
if ( abs(it->second->id())< 6 || it->second->id()==21 )
maxpt = it->second->momentum().perp();
}
}
else if (outGluon.size() + outQuark.size() > 1){
// assume qqbar pair from a Z/gamma
if (outGluon.size()==1 && outQuark.size() == 2 && outQuark[0]==-outQuark[1]){
for( it = currentTree()->outgoingLines().begin();
it != currentTree()->outgoingLines().end(); ++it ) {
if ( it->second->id()==21 )
maxpt = it->second->momentum().perp();
}
}
// otherwise take the lowest pT avoiding born DY events
else {
maxpt = generator()->maximumCMEnergy();
for( it = currentTree()->outgoingLines().begin();
it != currentTree()->outgoingLines().end(); ++it ) {
if ( abs(it->second->id())< 6 || it->second->id()==21 )
maxpt = min(maxpt,it->second->momentum().perp());
}
}
}
}
// set maximum pT for subsequent emissions from S events
if ( isPowhegSEvent || (!isPowhegSEvent && !isPowhegHEvent)){
for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
= currentTree()->outgoingLines().begin();
it != currentTree()->outgoingLines().end(); ++it ) {
if( ! it->second->coloured() ) continue;
it->first->maximumpT(maxpt, ShowerInteraction::QCD );
}
for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
= currentTree()->incomingLines().begin();
it != currentTree()->incomingLines().end(); ++it ) {
if( ! it->second->coloured() ) continue;
it->first->maximumpT(maxpt, ShowerInteraction::QCD );
}
}
}
else
_hardtree = ShowerHandler::currentHandler()->generateCKKW(currentTree());
// if hard me doesn't have a FSR powheg
// correction use decay powheg correction
if (_hardme && _hardme->hasPOWHEGCorrection()<2) {
// check for intermediate colour singlet resonance
const ParticleVector inter = _hardme->subProcess()->intermediates();
if (inter.size()!=1 ||
inter[0]->momentum().m2()/GeV2 < 0 ||
inter[0]->dataPtr()->iColour()!=PDT::Colour0){
if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
return;
}
map<ShowerProgenitorPtr, tShowerParticlePtr > out = currentTree()->outgoingLines();
// ignore cases where outgoing particles are not coloured
if (out.size()!=2 ||
out. begin()->second->dataPtr()->iColour()==PDT::Colour0 ||
out.rbegin()->second->dataPtr()->iColour()==PDT::Colour0) {
if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
return;
}
// look up decay mode
tDMPtr dm;
string tag;
string inParticle = inter[0]->dataPtr()->name() + "->";
vector<string> outParticles;
outParticles.push_back(out.begin ()->first->progenitor()->dataPtr()->name());
outParticles.push_back(out.rbegin()->first->progenitor()->dataPtr()->name());
for (int it=0; it<2; ++it){
tag = inParticle + outParticles[it] + "," + outParticles[(it+1)%2] + ";";
dm = generator()->findDecayMode(tag);
if(dm) break;
}
// get the decayer
HwDecayerBasePtr decayer;
if(dm) decayer = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
// check if decayer has a FSR POWHEG correction
if (!decayer || decayer->hasPOWHEGCorrection()<2){
if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
return;
}
// generate the hardest emission
ShowerDecayMap decay;
PPtr in = new_ptr(*inter[0]);
ShowerTreePtr decayTree = new_ptr(ShowerTree(in, decay));
HardTreePtr FSRTree = decayer->generateHardest(decayTree);
if (!FSRTree) {
if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
return;
}
// if there is no ISRTree make _hardtree from FSRTree
if (!ISRTree){
vector<HardBranchingPtr> inBranch,hardBranch;
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit =currentTree()->incomingLines().begin();
cit!=currentTree()->incomingLines().end();++cit ) {
inBranch.push_back(new_ptr(HardBranching(cit->second,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Incoming)));
inBranch.back()->beam(cit->first->original()->parents()[0]);
hardBranch.push_back(inBranch.back());
}
if(inBranch[0]->branchingParticle()->dataPtr()->coloured()) {
inBranch[0]->colourPartner(inBranch[1]);
inBranch[1]->colourPartner(inBranch[0]);
}
for(set<HardBranchingPtr>::iterator it=FSRTree->branchings().begin();
it!=FSRTree->branchings().end();++it) {
if((**it).branchingParticle()->id()!=in->id())
hardBranch.push_back(*it);
}
hardBranch[2]->colourPartner(hardBranch[3]);
hardBranch[3]->colourPartner(hardBranch[2]);
HardTreePtr newTree = new_ptr(HardTree(hardBranch,inBranch,
ShowerInteraction::QCD));
_hardtree = newTree;
}
// Otherwise modify the ISRTree to include the emission in FSRTree
else {
vector<tShowerParticlePtr> FSROut, ISROut;
set<HardBranchingPtr>::iterator itFSR, itISR;
// get outgoing particles
for(itFSR =FSRTree->branchings().begin();
itFSR!=FSRTree->branchings().end();++itFSR){
if ((**itFSR).status()==HardBranching::Outgoing)
FSROut.push_back((*itFSR)->branchingParticle());
}
for(itISR =ISRTree->branchings().begin();
itISR!=ISRTree->branchings().end();++itISR){
if ((**itISR).status()==HardBranching::Outgoing)
ISROut.push_back((*itISR)->branchingParticle());
}
// find COM frame formed by outgoing particles
LorentzRotation eventFrameFSR, eventFrameISR;
eventFrameFSR = ((FSROut[0]->momentum()+FSROut[1]->momentum()).findBoostToCM());
eventFrameISR = ((ISROut[0]->momentum()+ISROut[1]->momentum()).findBoostToCM());
// find rotation between ISR and FSR frames
int j=0;
if (ISROut[0]->id()!=FSROut[0]->id()) j=1;
eventFrameISR.rotateZ( (eventFrameFSR*FSROut[0]->momentum()).phi()-
(eventFrameISR*ISROut[j]->momentum()).phi() );
eventFrameISR.rotateY( (eventFrameFSR*FSROut[0]->momentum()).theta()-
(eventFrameISR*ISROut[j]->momentum()).theta() );
eventFrameISR.invert();
for (itFSR=FSRTree->branchings().begin();
itFSR!=FSRTree->branchings().end();++itFSR){
if ((**itFSR).branchingParticle()->id()==in->id()) continue;
for (itISR =ISRTree->branchings().begin();
itISR!=ISRTree->branchings().end();++itISR){
if ((**itISR).status()==HardBranching::Incoming) continue;
if ((**itFSR).branchingParticle()->id()==
(**itISR).branchingParticle()->id()){
// rotate FSRTree particle to ISRTree event frame
(**itISR).branchingParticle()->setMomentum(eventFrameISR*
eventFrameFSR*
(**itFSR).branchingParticle()->momentum());
(**itISR).branchingParticle()->rescaleMass();
// add the children of the FSRTree particles to the ISRTree
if(!(**itFSR).children().empty()){
(**itISR).addChild((**itFSR).children()[0]);
(**itISR).addChild((**itFSR).children()[1]);
// rotate momenta to ISRTree event frame
(**itISR).children()[0]->branchingParticle()->setMomentum(eventFrameISR*
eventFrameFSR*
(**itFSR).children()[0]->branchingParticle()->momentum());
(**itISR).children()[1]->branchingParticle()->setMomentum(eventFrameISR*
eventFrameFSR*
(**itFSR).children()[1]->branchingParticle()->momentum());
}
}
}
}
_hardtree = ISRTree;
}
}
if(_hardtree){
connectTrees(currentTree(),_hardtree,hard);
}
}
bool Evolver::truncatedTimeLikeShower(tShowerParticlePtr particle,
HardBranchingPtr branch,
ShowerInteraction::Type type,bool first) {
int ntry=0;
do {
++ntry;
Branching fb;
unsigned int iout=0;
tcPDPtr pdata[2];
while (true) {
// no truncated shower break
if(!isTruncatedShowerON()||hardOnly()) break;
// generate emission
fb=splittingGenerator()->chooseForwardBranching(*particle,1.,type);
// no emission break
if(!fb.kinematics) break;
// check haven't evolved too far
if(fb.kinematics->scale() < branch->scale()) {
fb=Branching();
break;
}
// get the particle data objects
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
// find the truncated line
iout=0;
if(pdata[0]->id()!=pdata[1]->id()) {
if(pdata[0]->id()==particle->id()) iout=1;
else if (pdata[1]->id()==particle->id()) iout=2;
}
else if(pdata[0]->id()==particle->id()) {
if(fb.kinematics->z()>0.5) iout=1;
else iout=2;
}
// apply the vetos for the truncated shower
// no flavour changing branchings
if(iout==0) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
// only if same interaction for forced branching
ShowerInteraction::Type type2 = fb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
// and evolution
if(type2==branch->sudakov()->interactionType()) {
if(zsplit < 0.5 || // hardest line veto
fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
// should do base class vetos as well
if(timeLikeVetoed(fb,particle)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
break;
}
// if no branching force truncated emission
if(!fb.kinematics) {
// construct the kinematics for the hard emission
ShoKinPtr showerKin=
branch->sudakov()->createFinalStateBranching(branch->scale(),
branch->children()[0]->z(),
branch->phi(),
branch->children()[0]->pT());
showerKin->initialize( *particle,PPtr() );
IdList idlist(3);
idlist[0] = particle->id();
idlist[1] = branch->children()[0]->branchingParticle()->id();
idlist[2] = branch->children()[1]->branchingParticle()->id();
fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() );
// Assign the shower kinematics to the emitting particle.
++_nFSR;
particle->showerKinematics( fb.kinematics );
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
for(unsigned int ix=0;ix<2;++ix) {
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[ix]->
branchingParticle()->dataPtr(),true)));
if(theChildren[ix]->id()==_progenitor->id()&&!theChildren[ix]->dataPtr()->stable())
theChildren[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
else
theChildren[ix]->set5Momentum(Lorentz5Momentum(theChildren[ix]->dataPtr()->mass()));
}
particle->showerKinematics()->
updateChildren(particle, theChildren,fb.type);
for(unsigned int ix=0;ix<2;++ix) {
theChildren[ix]->scales().QED = min(theChildren[ix]->scales().QED ,particle->scales().QED );
theChildren[ix]->scales().QED_noAO = min(theChildren[ix]->scales().QED_noAO ,particle->scales().QED_noAO );
theChildren[ix]->scales().QCD_c = min(theChildren[ix]->scales().QCD_c ,particle->scales().QCD_c );
theChildren[ix]->scales().QCD_c_noAO = min(theChildren[ix]->scales().QCD_c_noAO ,particle->scales().QCD_c_noAO );
theChildren[ix]->scales().QCD_ac = min(theChildren[ix]->scales().QCD_ac ,particle->scales().QCD_ac );
theChildren[ix]->scales().QCD_ac_noAO = min(theChildren[ix]->scales().QCD_ac_noAO,particle->scales().QCD_ac_noAO);
}
// shower the first particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() )
timeLikeShower(theChildren[0],type,Branching(),false);
}
else {
truncatedTimeLikeShower( theChildren[0],branch->children()[0],type,false);
}
// shower the second particle
if( branch->children()[1]->children().empty() ) {
if( ! hardOnly() )
timeLikeShower( theChildren[1] , type,Branching(),false);
}
else {
truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type,false);
}
// that's if for old approach
if(_reconOpt==0) return true;
// branching has happened
particle->showerKinematics()->updateParent(particle, theChildren,fb.type);
// clean up the vetoed emission
if(particle->virtualMass()==ZERO) {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<theChildren.size();++ix)
particle->abandonChild(theChildren[ix]);
theChildren.clear();
continue;
}
else {
if(first&&!theChildren.empty())
particle->showerKinematics()->resetChildren(particle,theChildren);
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
}
// has emitted
// Assign the shower kinematics to the emitting particle.
++_nFSR;
particle->showerKinematics(fb.kinematics);
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
for(unsigned int ix=0;ix<2;++ix) {
theChildren.push_back( new_ptr( ShowerParticle( pdata[ix], true ) ) );
if(theChildren[ix]->id()==_progenitor->id()&&!pdata[ix]->stable())
theChildren[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
else
theChildren[ix]->set5Momentum(Lorentz5Momentum(pdata[ix]->mass()));
}
particle->showerKinematics()->
updateChildren( particle, theChildren , fb.type);
// shower the first particle
if( iout == 1 ) truncatedTimeLikeShower( theChildren[0], branch , type ,false);
else timeLikeShower( theChildren[0] , type,Branching(),false);
// shower the second particle
if( iout == 2 ) truncatedTimeLikeShower( theChildren[1], branch , type ,false);
else timeLikeShower( theChildren[1] , type,Branching(),false);
// that's if for old approach
if(_reconOpt==0) return true;
// branching has happened
particle->showerKinematics()->updateParent(particle, theChildren,fb.type);
// clean up the vetoed emission
if(particle->virtualMass()==ZERO) {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<theChildren.size();++ix)
particle->abandonChild(theChildren[ix]);
theChildren.clear();
}
else {
if(first&&!theChildren.empty())
particle->showerKinematics()->resetChildren(particle,theChildren);
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
}
while(ntry<50);
return false;
}
bool Evolver::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
HardBranchingPtr branch,
ShowerInteraction::Type type) {
tcPDFPtr pdf;
if(ShowerHandler::currentHandler()->firstPDF().particle() == beamParticle())
pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
if(ShowerHandler::currentHandler()->secondPDF().particle() == beamParticle())
pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
Branching bb;
// parameters of the force branching
double z(0.);
HardBranchingPtr timelike;
for( unsigned int ix = 0; ix < branch->children().size(); ++ix ) {
if( branch->children()[ix]->status() ==HardBranching::Outgoing) {
timelike = branch->children()[ix];
}
if( branch->children()[ix]->status() ==HardBranching::Incoming )
z = branch->children()[ix]->z();
}
// generate truncated branching
tcPDPtr part[2];
if(z>=0.&&z<=1.) {
while (true) {
if( !isTruncatedShowerON() || hardOnly() ) break;
bb = splittingGenerator()->chooseBackwardBranching( *particle,
beam, 1., beamParticle(),
type , pdf,freeze);
if( !bb.kinematics || bb.kinematics->scale() < branch->scale() ) {
bb = Branching();
break;
}
// particles as in Sudakov form factor
part[0] = getParticleData( bb.ids[0] );
part[1] = getParticleData( bb.ids[2] );
//is emitter anti-particle
if( particle->id() != bb.ids[1]) {
if( part[0]->CC() ) part[0] = part[0]->CC();
if( part[1]->CC() ) part[1] = part[1]->CC();
}
double zsplit = bb.kinematics->z();
// apply the vetos for the truncated shower
// if doesn't carry most of momentum
ShowerInteraction::Type type2 = bb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
if(type2==branch->sudakov()->interactionType() &&
zsplit < 0.5) {
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
// others
if( part[0]->id() != particle->id() || // if particle changes type
bb.kinematics->pT() > progenitor()->maximumpT(type2) || // pt veto
bb.kinematics->scale() < branch->scale()) { // angular ordering veto
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
// and those from the base class
if(spaceLikeVetoed(bb,particle)) {
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
break;
}
}
if( !bb.kinematics ) {
//do the hard emission
ShoKinPtr kinematics =
branch->sudakov()->createInitialStateBranching( branch->scale(), z, branch->phi(),
branch->children()[0]->pT() );
kinematics->initialize( *particle, beam );
// assign the splitting function and shower kinematics
particle->showerKinematics( kinematics );
if(kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(kinematics->pT());
// For the time being we are considering only 1->2 branching
// Now create the actual particles, make the otherChild a final state
// particle, while the newParent is not
ShowerParticlePtr newParent =
new_ptr( ShowerParticle( branch->branchingParticle()->dataPtr(), false ) );
ShowerParticlePtr otherChild =
new_ptr( ShowerParticle( timelike->branchingParticle()->dataPtr(),
true, true ) );
ShowerParticleVector theChildren;
theChildren.push_back( particle );
theChildren.push_back( otherChild );
particle->showerKinematics()->
updateParent( newParent, theChildren, branch->type());
// update the history if needed
currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
currentTree()->addInitialStateBranching( particle, newParent, otherChild );
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
bool emitted=false;
if(!hardOnly()) {
if( branch->parent() ) {
emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type);
}
else {
emitted = spaceLikeShower( newParent, beam , type);
}
}
if( !emitted ) {
if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
kinematics->updateLast( newParent, ZERO, ZERO );
}
else {
pair<Energy,double> kt = intrinsicpT()[progenitor()];
kinematics->updateLast( newParent,
kt.first*cos( kt.second ),
kt.first*sin( kt.second ) );
}
}
particle->showerKinematics()->
updateChildren( newParent, theChildren,bb.type);
if(hardOnly()) return true;
// perform the shower of the final-state particle
if( timelike->children().empty() ) {
timeLikeShower( otherChild , type,Branching(),true);
}
else {
truncatedTimeLikeShower( otherChild, timelike , type,true);
}
updateHistory(otherChild);
// return the emitted
return true;
}
// assign the splitting function and shower kinematics
particle->showerKinematics( bb.kinematics );
if(bb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(bb.kinematics->pT());
// For the time being we are considering only 1->2 branching
// Now create the actual particles, make the otherChild a final state
// particle, while the newParent is not
ShowerParticlePtr newParent = new_ptr( ShowerParticle( part[0], false ) );
ShowerParticlePtr otherChild = new_ptr( ShowerParticle( part[1], true, true ) );
ShowerParticleVector theChildren;
theChildren.push_back( particle );
theChildren.push_back( otherChild );
particle->showerKinematics()->
updateParent( newParent, theChildren, bb.type);
// update the history if needed
currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
currentTree()->addInitialStateBranching( particle, newParent, otherChild );
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type);
// now reconstruct the momentum
if( !emitted ) {
if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
bb.kinematics->updateLast( newParent, ZERO, ZERO );
}
else {
pair<Energy,double> kt = intrinsicpT()[ progenitor() ];
bb.kinematics->updateLast( newParent,
kt.first*cos( kt.second ),
kt.first*sin( kt.second ) );
}
}
particle->showerKinematics()->
updateChildren( newParent, theChildren, bb.type);
// perform the shower of the final-state particle
timeLikeShower( otherChild , type,Branching(),true);
updateHistory(otherChild);
// return the emitted
return true;
}
bool Evolver::
truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
const ShowerParticle::EvolutionScales & maxScales,
Energy minmass, HardBranchingPtr branch,
ShowerInteraction::Type type) {
Branching fb;
unsigned int iout=0;
tcPDPtr pdata[2];
while (true) {
// no truncated shower break
if(!isTruncatedShowerON()||hardOnly()) break;
fb=splittingGenerator()->chooseDecayBranching(*particle,maxScales,minmass,1.,type);
// return if no radiation
if(!fb.kinematics) break;
// check haven't evolved too far
if(fb.kinematics->scale() < branch->scale()) {
fb=Branching();
break;
}
// get the particle data objects
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
// find the truncated line
iout=0;
if(pdata[0]->id()!=pdata[1]->id()) {
if(pdata[0]->id()==particle->id()) iout=1;
else if (pdata[1]->id()==particle->id()) iout=2;
}
else if(pdata[0]->id()==particle->id()) {
if(fb.kinematics->z()>0.5) iout=1;
else iout=2;
}
// apply the vetos for the truncated shower
// no flavour changing branchings
if(iout==0) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
ShowerInteraction::Type type2 = fb.type==ShowerPartnerType::QED ?
ShowerInteraction::QED : ShowerInteraction::QCD;
double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
if(type2==branch->sudakov()->interactionType()) {
if(zsplit < 0.5 || // hardest line veto
fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
// should do base class vetos as well
// if not vetoed break
if(!spaceLikeDecayVetoed(fb,particle)) break;
// otherwise reset scale and continue
particle->vetoEmission(fb.type,fb.kinematics->scale());
}
// this may not be currently used but in principle could be
// and should be included
if (!fb.kinematics) {
// construct the kinematics for the hard emission
ShoKinPtr showerKin=
branch->sudakov()->createDecayBranching(branch->scale(),
branch->children()[0]->z(),
branch->phi(),
branch->children()[0]->pT());
showerKin->initialize( *particle,PPtr() );
IdList idlist(3);
idlist[0] = particle->id();
idlist[1] = branch->children()[0]->branchingParticle()->id();
idlist[2] = branch->children()[1]->branchingParticle()->id();
// create the branching
fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine );
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics( fb.kinematics );
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[0]->
branchingParticle()->dataPtr(),true)));
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[1]->
branchingParticle()->dataPtr(),true)));
particle->showerKinematics()->
updateChildren(particle, theChildren,fb.type);
if(theChildren[0]->id()==particle->id()) {
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the space-like particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() ) spaceLikeDecayShower(theChildren[0],maxScales,minmass,type);
}
else {
truncatedSpaceLikeDecayShower( theChildren[0],maxScales,minmass,
branch->children()[0],type);
}
// shower the second particle
if( branch->children()[1]->children().empty() ) {
if( ! hardOnly() ) timeLikeShower( theChildren[1] , type,Branching(), true);
}
else {
truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type,true);
}
updateHistory(theChildren[1]);
}
else {
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[1]);
currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the space-like particle
if( branch->children()[1]->children().empty() ) {
if( ! hardOnly() ) spaceLikeDecayShower(theChildren[1],maxScales,minmass,type);
}
else {
truncatedSpaceLikeDecayShower( theChildren[1],maxScales,minmass,
branch->children()[1],type);
}
// shower the second particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() ) timeLikeShower( theChildren[0] , type, Branching(),true);
}
else {
truncatedTimeLikeShower( theChildren[0],branch->children()[0] ,type,true);
}
updateHistory(theChildren[0]);
}
return true;
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
particle->showerKinematics()->updateChildren(particle, theChildren,fb.type);
// In the case of splittings which involves coloured particles,
// set properly the colour flow of the branching.
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the first particle
truncatedSpaceLikeDecayShower(theChildren[0],maxScales,minmass,branch,type);
// shower the second particle
timeLikeShower(theChildren[1],type,Branching(),true);
updateHistory(theChildren[1]);
// branching has happened
return true;
}
bool Evolver::constructDecayTree(vector<ShowerProgenitorPtr> & particlesToShower,
ShowerInteraction::Type inter) {
Energy ptmax(-GeV);
// get the maximum pt is all ready a hard tree
if(hardTree()) {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->maximumpT(inter)>ptmax&&
particlesToShower[ix]->progenitor()->isFinalState())
ptmax = particlesToShower[ix]->maximumpT(inter);
}
}
vector<HardBranchingPtr> spaceBranchings,allBranchings;
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
HardBranchingPtr newBranch;
if(particlesToShower[ix]->hasEmitted()) {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
particlesToShower[ix]->progenitor()->
showerKinematics()->SudakovFormFactor(),
HardBranchingPtr(),HardBranching::Outgoing));
constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
}
else {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Outgoing));
}
allBranchings.push_back(newBranch);
}
else {
HardBranchingPtr newBranch;
if(particlesToShower[ix]->hasEmitted()) {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
particlesToShower[ix]->progenitor()->
showerKinematics()->SudakovFormFactor(),
HardBranchingPtr(),HardBranching::Decay));
constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
HardBranchingPtr last=newBranch;
do {
for(unsigned int ix=0;ix<last->children().size();++ix) {
if(last->children()[ix]->branchingParticle()->id()==
particlesToShower[ix]->id()) {
last = last->children()[ix];
continue;
}
}
}
while(!last->children().empty());
last->status(HardBranching::Incoming);
spaceBranchings.push_back(newBranch);
allBranchings .push_back(last);
}
else {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Incoming));
spaceBranchings.push_back(newBranch);
allBranchings .push_back(newBranch);
}
}
}
HardTreePtr QCDTree = new_ptr(HardTree(allBranchings,spaceBranchings,inter));
// set the charge partners
ShowerParticleVector particles;
particles.push_back(spaceBranchings.back()->branchingParticle());
for(set<HardBranchingPtr>::iterator cit=QCDTree->branchings().begin();
cit!=QCDTree->branchings().end();++cit) {
if((*cit)->status()==HardBranching::Outgoing)
particles.push_back((*cit)->branchingParticle());
}
// get the partners
showerModel()->partnerFinder()->setInitialEvolutionScales(particles,true,inter,true);
// do the inverse recon
if(!showerModel()->kinematicsReconstructor()->
deconstructDecayJets(QCDTree,this,inter)) {
return false;
}
// clear the old shower
currentTree()->clear();
// set the hard tree
hardTree(QCDTree);
// set the charge partners
setEvolutionPartners(false,inter,false);
// get the particles to be showered
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
particlesToShower.clear();
// incoming particles
for(cit=currentTree()->incomingLines().begin();
cit!=currentTree()->incomingLines().end();++cit)
particlesToShower.push_back(((*cit).first));
assert(particlesToShower.size()==1);
// outgoing particles
for(cjt=currentTree()->outgoingLines().begin();
cjt!=currentTree()->outgoingLines().end();++cjt) {
particlesToShower.push_back(((*cjt).first));
if(ptmax>ZERO) particlesToShower.back()->maximumpT(ptmax,inter);
}
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(particlesToShower[ix]->progenitor());
if( mit != eit) {
if(mit->second->status()==HardBranching::Outgoing)
particlesToShower[ix]->progenitor()->set5Momentum(mit->second->pVector());
}
}
return true;
}
bool Evolver::constructHardTree(vector<ShowerProgenitorPtr> & particlesToShower,
ShowerInteraction::Type inter) {
bool noEmission = true;
vector<HardBranchingPtr> spaceBranchings,allBranchings;
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
HardBranchingPtr newBranch;
if(particlesToShower[ix]->hasEmitted()) {
noEmission = false;
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
particlesToShower[ix]->progenitor()->
showerKinematics()->SudakovFormFactor(),
HardBranchingPtr(),HardBranching::Outgoing));
constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
}
else {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Outgoing));
}
allBranchings.push_back(newBranch);
}
else {
HardBranchingPtr first,last;
if(!particlesToShower[ix]->progenitor()->parents().empty()) {
noEmission = false;
constructSpaceLikeLine(particlesToShower[ix]->progenitor(),
first,last,SudakovPtr(),
particlesToShower[ix]->original()->parents()[0]);
}
else {
first = new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Incoming));
if(particlesToShower[ix]->original()->parents().empty())
first->beam(particlesToShower[ix]->original());
else
first->beam(particlesToShower[ix]->original()->parents()[0]);
last = first;
}
spaceBranchings.push_back(first);
allBranchings.push_back(last);
}
}
if(!noEmission) {
HardTreePtr QCDTree = new_ptr(HardTree(allBranchings,spaceBranchings,
inter));
// set the charge partners
ShowerParticleVector particles;
for(set<HardBranchingPtr>::iterator cit=QCDTree->branchings().begin();
cit!=QCDTree->branchings().end();++cit) {
particles.push_back((*cit)->branchingParticle());
}
// get the partners
showerModel()->partnerFinder()->setInitialEvolutionScales(particles,false,
inter,true);
// do the inverse recon
if(!showerModel()->kinematicsReconstructor()->
deconstructHardJets(QCDTree,this,inter))
throw Exception() << "Can't to shower deconstruction for QED shower in"
<< "QEDEvolver::showerHard" << Exception::eventerror;
// set the hard tree
hardTree(QCDTree);
}
// clear the old shower
currentTree()->clear();
// set the charge partners
setEvolutionPartners(true,inter,false);
// get the particles to be showered
particlesToShower = currentTree()->extractProgenitors();
// reset momenta
if(hardTree()) {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(particlesToShower[ix]->progenitor());
if( mit != eit) {
particlesToShower[ix]->progenitor()->set5Momentum(mit->second->showerMomentum());
}
}
}
return true;
}
void Evolver::constructTimeLikeLine(tHardBranchingPtr branch,
tShowerParticlePtr particle) {
for(unsigned int ix=0;ix<particle->children().size();++ix) {
HardBranching::Status status = branch->status();
tShowerParticlePtr child =
dynamic_ptr_cast<ShowerParticlePtr>(particle->children()[ix]);
if(child->children().empty()) {
HardBranchingPtr newBranch =
new_ptr(HardBranching(child,SudakovPtr(),branch,status));
branch->addChild(newBranch);
}
else {
HardBranchingPtr newBranch =
new_ptr(HardBranching(child,child->showerKinematics()->SudakovFormFactor(),
branch,status));
constructTimeLikeLine(newBranch,child);
branch->addChild(newBranch);
}
}
// sort out the type of interaction
if(!branch->children().empty()) {
if(branch->branchingParticle()->id()==ParticleID::gamma ||
branch->children()[0]->branchingParticle()->id()==ParticleID::gamma ||
branch->children()[1]->branchingParticle()->id()==ParticleID::gamma)
branch->type(ShowerPartnerType::QED);
else {
if(branch->branchingParticle()->id()==
branch->children()[0]->branchingParticle()->id()) {
if(branch->branchingParticle()->dataPtr()->iColour()==PDT::Colour8) {
tShowerParticlePtr emittor =
branch->branchingParticle()->showerKinematics()->z()>0.5 ?
branch->children()[0]->branchingParticle() :
branch->children()[1]->branchingParticle();
if(branch->branchingParticle()->colourLine()==emittor->colourLine())
branch->type(ShowerPartnerType::QCDAntiColourLine);
else if(branch->branchingParticle()->antiColourLine()==emittor->antiColourLine())
branch->type(ShowerPartnerType::QCDColourLine);
else
assert(false);
}
else if(branch->branchingParticle()->colourLine()) {
branch->type(ShowerPartnerType::QCDColourLine);
}
else if(branch->branchingParticle()->antiColourLine()) {
branch->type(ShowerPartnerType::QCDAntiColourLine);
}
else
assert(false);
}
else if(branch->branchingParticle()->id()==ParticleID::g &&
branch->children()[0]->branchingParticle()->id()==
-branch->children()[1]->branchingParticle()->id()) {
if(branch->branchingParticle()->showerKinematics()->z()>0.5)
branch->type(ShowerPartnerType::QCDAntiColourLine);
else
branch->type(ShowerPartnerType::QCDColourLine);
}
else
assert(false);
}
}
}
void Evolver::constructSpaceLikeLine(tShowerParticlePtr particle,
HardBranchingPtr & first,
HardBranchingPtr & last,
SudakovPtr sud,PPtr beam) {
if(!particle) return;
if(!particle->parents().empty()) {
tShowerParticlePtr parent =
dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]);
SudakovPtr newSud=particle->showerKinematics()->SudakovFormFactor();
constructSpaceLikeLine(parent,first,last,newSud,beam);
}
HardBranchingPtr newBranch =
new_ptr(HardBranching(particle,sud,last,HardBranching::Incoming));
newBranch->beam(beam);
if(!first) {
first=newBranch;
last =newBranch;
return;
}
last->addChild(newBranch);
tShowerParticlePtr timeChild =
dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]->children()[1]);
HardBranchingPtr timeBranch;
if(!timeChild->children().empty()) {
timeBranch =
new_ptr(HardBranching(timeChild,
timeChild->showerKinematics()->SudakovFormFactor(),
last,HardBranching::Outgoing));
constructTimeLikeLine(timeBranch,timeChild);
}
else {
timeBranch =
new_ptr(HardBranching(timeChild,SudakovPtr(),last,HardBranching::Outgoing));
}
last->addChild(timeBranch);
// sort out the type
if(last->branchingParticle() ->id() == ParticleID::gamma ||
newBranch->branchingParticle() ->id() == ParticleID::gamma ||
timeBranch->branchingParticle()->id() == ParticleID::gamma) {
last->type(ShowerPartnerType::QED);
}
else if(last->branchingParticle()->id()==newBranch->branchingParticle()->id()) {
if(last->branchingParticle()->id()==ParticleID::g) {
if(last->branchingParticle()->colourLine()==
newBranch->branchingParticle()->colourLine()) {
last->type(ShowerPartnerType::QCDAntiColourLine);
}
else {
last->type(ShowerPartnerType::QCDColourLine);
}
}
else if(last->branchingParticle()->hasColour()) {
last->type(ShowerPartnerType::QCDColourLine);
}
else if(last->branchingParticle()->hasAntiColour()) {
last->type(ShowerPartnerType::QCDAntiColourLine);
}
else
assert(false);
}
else if(newBranch->branchingParticle()->id()==ParticleID::g) {
if(last->branchingParticle()->hasColour()) {
last->type(ShowerPartnerType::QCDAntiColourLine);
}
else if(last->branchingParticle()->hasAntiColour()) {
last->type(ShowerPartnerType::QCDColourLine);
}
else
assert(false);
}
else if(newBranch->branchingParticle()->hasColour()) {
last->type(ShowerPartnerType::QCDColourLine);
}
else if(newBranch->branchingParticle()->hasAntiColour()) {
last->type(ShowerPartnerType::QCDAntiColourLine);
}
else {
assert(false);
}
last=newBranch;
}
void Evolver::connectTrees(ShowerTreePtr showerTree,
HardTreePtr hardTree, bool hard ) {
ShowerParticleVector particles;
// find the Sudakovs
for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
cit!=hardTree->branchings().end();++cit) {
// Sudakovs for ISR
if((**cit).parent()&&(**cit).status()==HardBranching::Incoming) {
++_nis;
IdList br(3);
br[0] = (**cit).parent()->branchingParticle()->id();
br[1] = (**cit). branchingParticle()->id();
br[2] = (**cit).parent()->children()[0]==*cit ?
(**cit).parent()->children()[1]->branchingParticle()->id() :
(**cit).parent()->children()[0]->branchingParticle()->id();
BranchingList branchings = splittingGenerator()->initialStateBranchings();
if(br[1]<0&&br[0]==br[1]) {
br[0] = abs(br[0]);
br[1] = abs(br[1]);
}
else if(br[1]<0) {
br[1] = -br[1];
br[2] = -br[2];
}
long index = abs(br[1]);
SudakovPtr sudakov;
for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
cjt != branchings.upper_bound(index); ++cjt ) {
IdList ids = cjt->second.second;
if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
sudakov=cjt->second.first;
break;
}
}
if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
<< "Evolver::connectTrees() for ISR"
<< Exception::runerror;
(**cit).parent()->sudakov(sudakov);
}
// Sudakovs for FSR
else if(!(**cit).children().empty()) {
++_nfs;
IdList br(3);
br[0] = (**cit) .branchingParticle()->id();
br[1] = (**cit).children()[0]->branchingParticle()->id();
br[2] = (**cit).children()[1]->branchingParticle()->id();
BranchingList branchings = splittingGenerator()->finalStateBranchings();
if(br[0]<0) {
br[0] = abs(br[0]);
br[1] = abs(br[1]);
br[2] = abs(br[2]);
}
long index = br[0];
SudakovPtr sudakov;
for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
cjt != branchings.upper_bound(index); ++cjt ) {
IdList ids = cjt->second.second;
if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
sudakov=cjt->second.first;
break;
}
}
if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
<< "Evolver::connectTrees()"
<< Exception::runerror;
(**cit).sudakov(sudakov);
}
}
// calculate the evolution scale
for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
cit!=hardTree->branchings().end();++cit) {
particles.push_back((*cit)->branchingParticle());
}
showerModel()->partnerFinder()->
setInitialEvolutionScales(particles,!hard,hardTree->interaction(),
!hardTree->partnersSet());
hardTree->partnersSet(true);
// inverse reconstruction
if(hard) {
showerModel()->kinematicsReconstructor()->
deconstructHardJets(hardTree,ShowerHandler::currentHandler()->evolver(),
hardTree->interaction());
}
else
showerModel()->kinematicsReconstructor()->
deconstructDecayJets(hardTree,ShowerHandler::currentHandler()->evolver(),
hardTree->interaction());
// now reset the momenta of the showering particles
vector<ShowerProgenitorPtr> particlesToShower;
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=showerTree->incomingLines().begin();
cit!=showerTree->incomingLines().end();++cit )
particlesToShower.push_back(cit->first);
// extract the showering particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=showerTree->outgoingLines().begin();
cit!=showerTree->outgoingLines().end();++cit )
particlesToShower.push_back(cit->first);
// match them
map<ShowerProgenitorPtr,HardBranchingPtr> partners;
for(set<HardBranchingPtr>::const_iterator bit=hardTree->branchings().begin();
bit!=hardTree->branchings().end();++bit) {
Energy2 dmin( 1e30*GeV2 );
ShowerProgenitorPtr partner;
for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
pit!=particlesToShower.end();++pit) {
if(partners.find(*pit)!=partners.end()) continue;
if( (**bit).branchingParticle()->id() != (**pit).progenitor()->id() ) continue;
if( (**bit).branchingParticle()->isFinalState() !=
(**pit).progenitor()->isFinalState() ) continue;
if( (**pit).progenitor()->isFinalState() ) {
Energy2 dtest =
sqr( (**pit).progenitor()->momentum().x() - (**bit).showerMomentum().x() ) +
sqr( (**pit).progenitor()->momentum().y() - (**bit).showerMomentum().y() ) +
sqr( (**pit).progenitor()->momentum().z() - (**bit).showerMomentum().z() ) +
sqr( (**pit).progenitor()->momentum().t() - (**bit).showerMomentum().t() );
// add mass difference for identical particles (e.g. Z0 Z0 production)
dtest += 1e10*sqr((**pit).progenitor()->momentum().m()-(**bit).showerMomentum().m());
if( dtest < dmin ) {
partner = *pit;
dmin = dtest;
}
}
else {
// ensure directions are right
if((**pit).progenitor()->momentum().z()/(**bit).showerMomentum().z()>ZERO) {
partner = *pit;
break;
}
}
}
if(!partner) throw Exception() << "Failed to match shower and hard trees in Evolver::hardestEmission"
<< Exception::eventerror;
partners[partner] = *bit;
}
for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
pit!=particlesToShower.end();++pit) {
HardBranchingPtr partner = partners[*pit];
if((**pit).progenitor()->dataPtr()->stable()) {
(**pit).progenitor()->set5Momentum(partner->showerMomentum());
(**pit).copy()->set5Momentum(partner->showerMomentum());
}
else {
Lorentz5Momentum oldMomentum = (**pit).progenitor()->momentum();
Lorentz5Momentum newMomentum = partner->showerMomentum();
LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
(**pit).progenitor()->transform(boost);
(**pit).copy() ->transform(boost);
boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
(**pit).progenitor()->transform(boost);
(**pit).copy() ->transform(boost);
}
}
// correction boosts for daughter trees
for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit = showerTree->treelinks().begin();
tit != showerTree->treelinks().end();++tit) {
ShowerTreePtr decayTree = tit->first;
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit = decayTree->incomingLines().begin();
// reset the momentum of the decay particle
Lorentz5Momentum oldMomentum = cit->first->progenitor()->momentum();
Lorentz5Momentum newMomentum = tit->second.second->momentum();
LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
decayTree->transform(boost,true);
boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
decayTree->transform(boost,true);
}
}
void Evolver::doShowering(bool hard,XCPtr xcomb) {
// order of the interactions
bool showerOrder(true);
// zero number of emissions
_nis = _nfs = 0;
// if MC@NLO H event and limited emissions
// indicate both final and initial state emission
if ( isMCatNLOHEvent && _limitEmissions != 0 ) {
_nis = _nfs = 1;
}
// extract particles to shower
vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
// setup the maximum scales for the shower
if (hardVetoOn()) setupMaximumScales(particlesToShower,xcomb);
// set the hard scales for the profiles
setupHardScales(particlesToShower,xcomb);
// specific stuff for hard processes and decays
Energy minmass(ZERO), mIn(ZERO);
// hard process generate the intrinsic p_T once and for all
if(hard) {
generateIntrinsicpT(particlesToShower);
}
// decay compute the minimum mass of the final-state
else {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
if(particlesToShower[ix]->progenitor()->dataPtr()->stable())
minmass += particlesToShower[ix]->progenitor()->dataPtr()->constituentMass();
else
minmass += particlesToShower[ix]->progenitor()->mass();
}
else {
mIn = particlesToShower[ix]->progenitor()->mass();
}
}
// throw exception if decay can't happen
if ( minmass > mIn ) {
throw Exception() << "Evolver.cc: Mass of decaying particle is "
<< "below constituent masses of decay products."
<< Exception::eventerror;
}
}
// check if interactions in right order
if(hardTree() && interaction_!=4 &&
hardTree()->interaction()!=interactions_[0]) {
assert(interactions_.size()==2);
showerOrder = false;
swap(interactions_[0],interactions_[1]);
}
// loop over possible interactions
for(unsigned int inter=0;inter<interactions_.size();++inter) {
// set up for second pass if required
if(inter!=0) {
// zero intrinsic pt so only added first time round
intrinsicpT().clear();
// construct the tree and throw veto if not possible
if(!(hard ?
constructHardTree (particlesToShower,interactions_[inter]) :
constructDecayTree(particlesToShower,interactions_[inter])))
throw InteractionVeto();
}
// create random particle vector (only need to do once)
vector<ShowerProgenitorPtr> tmp;
unsigned int nColouredIncoming = 0;
while(particlesToShower.size()>0){
unsigned int xx=UseRandom::irnd(particlesToShower.size());
tmp.push_back(particlesToShower[xx]);
particlesToShower.erase(particlesToShower.begin()+xx);
}
particlesToShower=tmp;
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(!particlesToShower[ix]->progenitor()->isFinalState() &&
particlesToShower[ix]->progenitor()->coloured()) ++nColouredIncoming;
}
bool switchRecon = hard && nColouredIncoming !=1;
// main shower loop
unsigned int ntry(0);
bool reconstructed = false;
do {
// clear results of last attempt if needed
if(ntry!=0) {
currentTree()->clear();
setEvolutionPartners(hard,interactions_[inter],true);
_nis = _nfs = 0;
// if MC@NLO H event and limited emissions
// indicate both final and initial state emission
if ( isMCatNLOHEvent && _limitEmissions != 0 ) {
_nis = _nfs = 1;
}
for(unsigned int ix=0; ix<particlesToShower.size();++ix) {
SpinPtr spin = particlesToShower[ix]->progenitor()->spinInfo();
if(spin && spin->decayVertex() &&
dynamic_ptr_cast<tcSVertexPtr>(spin->decayVertex())) {
spin->decayVertex(VertexPtr());
}
}
}
// loop over particles
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
// extract the progenitor
progenitor(particlesToShower[ix]);
// final-state radiation
if(progenitor()->progenitor()->isFinalState()) {
if(!isFSRadiationON()) continue;
// perform shower
progenitor()->hasEmitted(startTimeLikeShower(interactions_[inter]));
}
// initial-state radiation
else {
if(!isISRadiationON()) continue;
// hard process
if(hard) {
// get the PDF
setBeamParticle(_progenitor->beam());
assert(beamParticle());
// perform the shower
// set the beam particle
tPPtr beamparticle=progenitor()->original();
if(!beamparticle->parents().empty())
beamparticle=beamparticle->parents()[0];
// generate the shower
progenitor()->hasEmitted(startSpaceLikeShower(beamparticle,
interactions_[inter]));
}
// decay
else {
// skip colour and electrically neutral particles
if(!progenitor()->progenitor()->dataPtr()->coloured() &&
!progenitor()->progenitor()->dataPtr()->charged()) {
progenitor()->hasEmitted(false);
continue;
}
// perform shower
// set the scales correctly. The current scale is the maximum scale for
// emission not the starting scale
ShowerParticle::EvolutionScales maxScales(progenitor()->progenitor()->scales());
progenitor()->progenitor()->scales() = ShowerParticle::EvolutionScales();
if(progenitor()->progenitor()->dataPtr()->charged()) {
progenitor()->progenitor()->scales().QED = progenitor()->progenitor()->mass();
progenitor()->progenitor()->scales().QED_noAO = progenitor()->progenitor()->mass();
}
if(progenitor()->progenitor()->hasColour()) {
progenitor()->progenitor()->scales().QCD_c = progenitor()->progenitor()->mass();
progenitor()->progenitor()->scales().QCD_c_noAO = progenitor()->progenitor()->mass();
}
if(progenitor()->progenitor()->hasAntiColour()) {
progenitor()->progenitor()->scales().QCD_ac = progenitor()->progenitor()->mass();
progenitor()->progenitor()->scales().QCD_ac_noAO = progenitor()->progenitor()->mass();
}
// perform the shower
progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
interactions_[inter]));
}
}
}
// do the kinematic reconstruction, checking if it worked
reconstructed = hard ?
showerModel()->kinematicsReconstructor()->
reconstructHardJets (currentTree(),intrinsicpT(),interactions_[inter],
switchRecon && ntry>maximumTries()/2) :
showerModel()->kinematicsReconstructor()->
reconstructDecayJets(currentTree(),interactions_[inter]);
}
while(!reconstructed&&maximumTries()>++ntry);
// check if failed to generate the shower
if(ntry==maximumTries()) {
if(hard)
throw ShowerHandler::ShowerTriesVeto(ntry);
else
throw Exception() << "Failed to generate the shower after "
<< ntry << " attempts in Evolver::showerDecay()"
<< Exception::eventerror;
}
}
// tree has now showered
_currenttree->hasShowered(true);
if(!showerOrder) swap(interactions_[0],interactions_[1]);
hardTree(HardTreePtr());
}
void Evolver:: convertHardTree(bool hard,ShowerInteraction::Type type) {
map<ColinePtr,ColinePtr> cmap;
// incoming particles
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=currentTree()->incomingLines().begin();cit!=currentTree()->incomingLines().end();++cit) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
mit = hardTree()->particles().find(cit->first->progenitor());
// put the colour lines in the map
ShowerParticlePtr oldParticle = cit->first->progenitor();
ShowerParticlePtr newParticle = mit->second->branchingParticle();
ColinePtr cLine = oldParticle-> colourLine();
ColinePtr aLine = oldParticle->antiColourLine();
if(newParticle->colourLine() &&
cmap.find(newParticle-> colourLine())==cmap.end())
cmap[newParticle-> colourLine()] = cLine;
if(newParticle->antiColourLine() &&
cmap.find(newParticle->antiColourLine())==cmap.end())
cmap[newParticle->antiColourLine()] = aLine;
// check whether or not particle emits
bool emission = mit->second->parent();
if(emission) {
if(newParticle->colourLine()) {
ColinePtr ctemp = newParticle-> colourLine();
ctemp->removeColoured(newParticle);
}
if(newParticle->antiColourLine()) {
ColinePtr ctemp = newParticle->antiColourLine();
ctemp->removeAntiColoured(newParticle);
}
newParticle = mit->second->parent()->branchingParticle();
}
// get the new colour lines
ColinePtr newCLine,newALine;
// sort out colour lines
if(newParticle->colourLine()) {
ColinePtr ctemp = newParticle-> colourLine();
ctemp->removeColoured(newParticle);
if(cmap.find(ctemp)!=cmap.end()) {
newCLine = cmap[ctemp];
}
else {
newCLine = new_ptr(ColourLine());
cmap[ctemp] = newCLine;
}
}
// and anticolour lines
if(newParticle->antiColourLine()) {
ColinePtr ctemp = newParticle->antiColourLine();
ctemp->removeAntiColoured(newParticle);
if(cmap.find(ctemp)!=cmap.end()) {
newALine = cmap[ctemp];
}
else {
newALine = new_ptr(ColourLine());
cmap[ctemp] = newALine;
}
}
// remove colour lines from old particle
if(aLine) {
aLine->removeAntiColoured(cit->first->copy());
aLine->removeAntiColoured(cit->first->progenitor());
}
if(cLine) {
cLine->removeColoured(cit->first->copy());
cLine->removeColoured(cit->first->progenitor());
}
// add particle to colour lines
if(newCLine) newCLine->addColoured (newParticle);
if(newALine) newALine->addAntiColoured(newParticle);
// insert new particles
cit->first->copy(newParticle);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,false)));
cit->first->progenitor(sp);
currentTree()->incomingLines()[cit->first]=sp;
cit->first->perturbative(!emission);
// and the emitted particle if needed
if(emission) {
ShowerParticlePtr newOut = mit->second->parent()->children()[1]->branchingParticle();
if(newOut->colourLine()) {
ColinePtr ctemp = newOut-> colourLine();
ctemp->removeColoured(newOut);
assert(cmap.find(ctemp)!=cmap.end());
cmap[ctemp]->addColoured (newOut);
}
if(newOut->antiColourLine()) {
ColinePtr ctemp = newOut->antiColourLine();
ctemp->removeAntiColoured(newOut);
assert(cmap.find(ctemp)!=cmap.end());
cmap[ctemp]->addAntiColoured(newOut);
}
ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
out->perturbative(false);
currentTree()->outgoingLines().insert(make_pair(out,sout));
}
if(hard) {
// sort out the value of x
if(mit->second->beam()->momentum().z()>ZERO) {
sp->x(newParticle->momentum(). plus()/mit->second->beam()->momentum(). plus());
}
else {
sp->x(newParticle->momentum().minus()/mit->second->beam()->momentum().minus());
}
}
}
// outgoing particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=currentTree()->outgoingLines().begin();cit!=currentTree()->outgoingLines().end();++cit) {
map<tShowerTreePtr,pair<tShowerProgenitorPtr,
tShowerParticlePtr> >::const_iterator tit;
for(tit = currentTree()->treelinks().begin();
tit != currentTree()->treelinks().end();++tit) {
if(tit->second.first && tit->second.second==cit->first->progenitor())
break;
}
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
mit = hardTree()->particles().find(cit->first->progenitor());
if(mit==hardTree()->particles().end()) continue;
// put the colour lines in the map
ShowerParticlePtr oldParticle = cit->first->progenitor();
ShowerParticlePtr newParticle = mit->second->branchingParticle();
ShowerParticlePtr newOut;
ColinePtr cLine = oldParticle-> colourLine();
ColinePtr aLine = oldParticle->antiColourLine();
if(newParticle->colourLine() &&
cmap.find(newParticle-> colourLine())==cmap.end())
cmap[newParticle-> colourLine()] = cLine;
if(newParticle->antiColourLine() &&
cmap.find(newParticle->antiColourLine())==cmap.end())
cmap[newParticle->antiColourLine()] = aLine;
// check whether or not particle emits
bool emission = !mit->second->children().empty();
if(emission) {
if(newParticle->colourLine()) {
ColinePtr ctemp = newParticle-> colourLine();
ctemp->removeColoured(newParticle);
}
if(newParticle->antiColourLine()) {
ColinePtr ctemp = newParticle->antiColourLine();
ctemp->removeAntiColoured(newParticle);
}
newParticle = mit->second->children()[0]->branchingParticle();
newOut = mit->second->children()[1]->branchingParticle();
if(newParticle->id()!=oldParticle->id()&&newParticle->id()==newOut->id())
swap(newParticle,newOut);
}
// get the new colour lines
ColinePtr newCLine,newALine;
// sort out colour lines
if(newParticle->colourLine()) {
ColinePtr ctemp = newParticle-> colourLine();
ctemp->removeColoured(newParticle);
if(cmap.find(ctemp)!=cmap.end()) {
newCLine = cmap[ctemp];
}
else {
newCLine = new_ptr(ColourLine());
cmap[ctemp] = newCLine;
}
}
// and anticolour lines
if(newParticle->antiColourLine()) {
ColinePtr ctemp = newParticle->antiColourLine();
ctemp->removeAntiColoured(newParticle);
if(cmap.find(ctemp)!=cmap.end()) {
newALine = cmap[ctemp];
}
else {
newALine = new_ptr(ColourLine());
cmap[ctemp] = newALine;
}
}
// remove colour lines from old particle
if(aLine) {
aLine->removeAntiColoured(cit->first->copy());
aLine->removeAntiColoured(cit->first->progenitor());
}
if(cLine) {
cLine->removeColoured(cit->first->copy());
cLine->removeColoured(cit->first->progenitor());
}
// special for unstable particles
if(newParticle->id()==oldParticle->id() &&
(tit!=currentTree()->treelinks().end()||!oldParticle->dataPtr()->stable())) {
Lorentz5Momentum oldMomentum = oldParticle->momentum();
Lorentz5Momentum newMomentum = newParticle->momentum();
LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
oldParticle->transform(boost);
boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
oldParticle->transform(boost);
if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
newParticle=oldParticle;
}
// add particle to colour lines
if(newCLine) newCLine->addColoured (newParticle);
if(newALine) newALine->addAntiColoured(newParticle);
// insert new particles
cit->first->copy(newParticle);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,true)));
cit->first->progenitor(sp);
currentTree()->outgoingLines()[cit->first]=sp;
cit->first->perturbative(!emission);
// and the emitted particle if needed
if(emission) {
if(newOut->colourLine()) {
ColinePtr ctemp = newOut-> colourLine();
ctemp->removeColoured(newOut);
assert(cmap.find(ctemp)!=cmap.end());
cmap[ctemp]->addColoured (newOut);
}
if(newOut->antiColourLine()) {
ColinePtr ctemp = newOut->antiColourLine();
ctemp->removeAntiColoured(newOut);
assert(cmap.find(ctemp)!=cmap.end());
cmap[ctemp]->addAntiColoured(newOut);
}
ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
out->perturbative(false);
currentTree()->outgoingLines().insert(make_pair(out,sout));
}
// update any decay products
if(tit!=currentTree()->treelinks().end())
currentTree()->updateLink(tit->first,make_pair(cit->first,sp));
}
// reset the tree
currentTree()->resetShowerProducts();
// reextract the particles and set the colour partners
vector<ShowerParticlePtr> particles =
currentTree()->extractProgenitorParticles();
// clear the partners
for(unsigned int ix=0;ix<particles.size();++ix) {
particles[ix]->partner(ShowerParticlePtr());
particles[ix]->clearPartners();
}
// clear the tree
hardTree(HardTreePtr());
// Set the initial evolution scales
showerModel()->partnerFinder()->
setInitialEvolutionScales(particles,!hard,type,!_hardtree);
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 1:57 PM (15 h, 33 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023045
Default Alt Text
(166 KB)

Event Timeline