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,844 +1,849 @@ // -*- 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 #include #include #include #include #include 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 MadGraphAmplitude::BornAmplitudes=vector(); vector MadGraphAmplitude::VirtAmplitudes=vector(); 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(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::iterator amps=BornAmplitudes.begin();amps!=BornAmplitudes.end();amps++) borns<<*amps<::iterator amps=VirtAmplitudes.begin();amps!=VirtAmplitudes.end();amps++) virts<<*amps<::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::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; + + throw Exception() << "MadGraphAmplitude: The MadGraph amplitudes did not match the process.\n" + << " Please remove:"<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" <runStorage()+"/MadGraphAmplitudes"+"/MG-Parameter.dat"; ofstream params(para.c_str()); params<<"$WZ$ " <hardProcessWidth() /GeV; params<<"\n$WW$ " <hardProcessWidth()/GeV; params<<"\n$alphas$ " <hardProcessMass() /GeV<hardProcessMass() /GeV<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." <hardProcessMass() /GeV <hardProcessWidth() /GeV <hardProcessMass() /GeV <hardProcessMass() /GeV <hardProcessWidth() /GeV <hardProcessMass() /GeV <runStorage()+"/MadGraphAmplitudes "; cmd +=" --datadir "+pkgdatadir_; cmd +=" --includedir "+includedir_; std::stringstream as,aem; as << factory()->orderInAlphaS(); cmd +=" --orderas "+as.str() ; aem <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 <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( (*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::tptr factory, bool virt) const { if ( factory->processData()->diagramMap().find(p) != factory->processData()->diagramMap().end() ) return true; vector::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( (*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::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& 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 > 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 > tmpset; vector 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 >::iterator it3 = a.begin(); it3 != a.end(); it3++ ) { bool found_it3=false; for ( set >::iterator it4 = tmpset.begin(); it4 != tmpset.end(); it4++ ) { vector it3tmp = gluonsFirst(*it3); vector 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 < complex > >& tmp = lastMatchboxXComb()->heljamp(); const map,vector < complex > >& 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;j1&&j<=1)||(cross<=1&&j>1)){ heltmp[cross]=-1*hel[j];} else{heltmp[cross]=hel[j];} } vector reshuffled = meMomenta(); if ( !reshuffleMasses().empty() && reshuffled.size() > 3 ) { const cPDVector& pdata = mePartonData(); const map& masses = reshuffleMasses(); lastMatchboxXComb()->reshuffle(reshuffled,pdata,masses); } double momenta[50]; size_t j=0; for (size_t i=0;ipushheljamp(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 MadGraphAmplitude::physicalHelicities(const vector& hel) const { vector 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 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 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 polarization(pol[1],pol[2],pol[3],pol[0]); return polarization; } bool equalsModulo(unsigned int i, const vector& a, const vector& 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 MadGraphAmplitude::gluonsFirst(vector vec) { vector vecout; for(vector::iterator it= vec.begin();it!= vec.end();++it) if ( mePartonData()[crossingMap()[*it]]->id()==21) vecout.push_back(crossingMap()[*it]); for(vector::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 ij, const SpinCorrelationTensor& c) const { vector reshuffled = meMomenta(); if ( !reshuffleMasses().empty() && reshuffled.size() > 3 ) { const cPDVector& pdata = mePartonData(); const map& masses = reshuffleMasses(); lastMatchboxXComb()->reshuffle(reshuffled,pdata,masses); } Lorentz5Momentum p = reshuffled[ij.first]; Lorentz5Momentum n = reshuffled[ij.second]; LorentzVector 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 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::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 reshuffled = meMomenta(); if ( !reshuffleMasses().empty() && reshuffled.size() > 3 ) { const cPDVector& pdata = mePartonData(); const map& masses = reshuffleMasses(); lastMatchboxXComb()->reshuffle(reshuffled,pdata,masses); } double virt[20]; double momenta[50]; size_t j=0; for (size_t i=0;iexternalId(); 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(virt[2]/ifact*units,virt[3]/ifact*units)); } void MadGraphAmplitude::persistentOutput(PersistentOStream & os) const { os << theOrderInGs << theOrderInGem << BornAmplitudes << VirtAmplitudes << colourindex<> 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 describeHerwigMadGraphAmplitude("Herwig::MadGraphAmplitude", "HwMatchboxMadGraph.so"); void MadGraphAmplitude::Init() { static ClassDocumentation 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 interfaceProcessPath ("ProcessPath", "The Process Path.", &MadGraphAmplitude::theProcessPath, "",false, false); static Parameter interfaceModel ("Model", "The MadGraph-Model.", &MadGraphAmplitude::theMGmodel, "loop_sm",false, false); static Switch 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 interfaceBinDir ("BinDir", "The location for the installed executable", &MadGraphAmplitude::bindir_, string(HERWIG_BINDIR), false, false); static Parameter interfacePKGDATADIR ("DataDir", "The location for the installed Herwig data files", &MadGraphAmplitude::pkgdatadir_, string(HERWIG_PKGDATADIR), false, false); static Parameter 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,387 +1,392 @@ #! /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(" "))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()+" 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,absolute_links): cwd=os.getcwd() os.chdir(pwd) mkdir_p(pwd+"/Herwig-scratch/MG_tmp/") if not buildpath.startswith("/"): buildpath=pwd+"/"+buildpath.lstrip("./") if not buildpath.endswith("/"): buildpath=buildpath + "/" 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)): if not absolute_links: source=os.path.dirname(i) dest=pwd+"/Herwig-scratch/MG_tmp/" os.chdir(dest) os.symlink(os.path.relpath(source,dest)+"/"+os.path.basename(i),"./" + os.path.basename(i)) else: 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.") parser.add_option("-l", "--absolute-links",action="store_true", dest="absolute_links", default=False,\ 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,options.absolute_links) 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" + print "Known issue: If this is your first NLO calculation with pure Madgraph Amplitudes" + print " the CutTools compilation can result in a non usable configuration." + print " Please open $HERWIG_ENV/opt/madgraph/vendor/CutTools/makefile" + print " and add a \\ after FC=gfortran... in the ARGS variable." + print " Then run make clean && make in the CutTools folder. " 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,options.absolute_links) diff --git a/MatrixElement/Matchbox/MatchboxFactory.cc b/MatrixElement/Matchbox/MatchboxFactory.cc --- a/MatrixElement/Matchbox/MatchboxFactory.cc +++ b/MatrixElement/Matchbox/MatchboxFactory.cc @@ -1,2109 +1,2226 @@ // -*- C++ -*- // // MatchboxFactory.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 MatchboxFactory class. // #include "MatchboxFactory.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/Utilities/StringUtils.h" #include "ThePEG/Repository/Repository.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Handlers/SamplerBase.h" #include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h" #include "Herwig/MatrixElement/Matchbox/Utility/SU2Helper.h" #include "Herwig/Utilities/RunDirectories.h" #include #include #include using std::ostream_iterator; using namespace Herwig; using std::ostream_iterator; MatchboxFactory::MatchboxFactory() : SubProcessHandler(), theNLight(0), theOrderInAlphaS(0), theOrderInAlphaEW(0), theBornContributions(true), theVirtualContributions(true), theRealContributions(true), theIndependentVirtuals(false), theIndependentPKs(false), theSubProcessGroups(false), theFactorizationScaleFactor(1.0), theRenormalizationScaleFactor(1.0), theFixedCouplings(false), theFixedQEDCouplings(false), theVetoScales(false), theDipoleSet(0), theVerbose(false), theDiagramWeightVerbose(false), theDiagramWeightVerboseNBins(200), theInitVerbose(false), theSubtractionData(""), theSubtractionPlotType(1), theSubtractionScatterPlot(false), thePoleData(""), theRealEmissionScales(false), theAllProcesses(false), theMECorrectionsOnly(false), theLoopSimCorrections(false), ranSetup(false), theFirstPerturbativePDF(true), theSecondPerturbativePDF(true), inProductionMode(false), theSpinCorrelations(false),theAlphaParameter(1.), theEnforceChargeConservation(true), theEnforceColourConservation(false), theEnforceLeptonNumberConservation(false), theEnforceQuarkNumberConservation(false), theLeptonFlavourDiagonal(false), theQuarkFlavourDiagonal(false) {} MatchboxFactory::~MatchboxFactory() {} bool& MatchboxFactory::theIsMatchboxRun() { static bool flag = false; return flag; } IBPtr MatchboxFactory::clone() const { return new_ptr(*this); } IBPtr MatchboxFactory::fullclone() const { return new_ptr(*this); } void MatchboxFactory::prepareME(Ptr::ptr me) { Ptr::ptr amp = dynamic_ptr_cast::ptr>((*me).amplitude()); me->matchboxAmplitude(amp); me->factory(this); if ( phasespace() && !me->phasespace() ) me->phasespace(phasespace()); if ( scaleChoice() && !me->scaleChoice() ) me->scaleChoice(scaleChoice()); if ( !reweighters().empty() ) { for ( vector::const_iterator rw = reweighters().begin(); rw != reweighters().end(); ++rw ) me->addReweighter(*rw); } if ( !preweighters().empty() ) { for ( vector::const_iterator rw = preweighters().begin(); rw != preweighters().end(); ++rw ) me->addPreweighter(*rw); } } string pid(const PDVector& key) { ostringstream res; res << "[" << key[0]->PDGName() << "," << key[1]->PDGName() << "->"; for ( PDVector::const_iterator k = key.begin() + 2; k != key.end(); ++k ) res << (**k).PDGName() << (k != --key.end() ? "," : ""); res << "]"; return res.str(); } vector::ptr> MatchboxFactory:: makeMEs(const vector& proc, unsigned int orderas, bool virt) { generator()->log() << "determining subprocesses for "; copy(proc.begin(),proc.end(),ostream_iterator(generator()->log()," ")); generator()->log() << "\n" << flush; map::ptr,set > ampProcs; map::ptr> > procAmps; set processes = makeSubProcesses(proc); // TODO Fix me for 3.0.x // At the moment we got troubles with processes with no coloured // legs so they will not be supported set colouredProcesses; for ( set::const_iterator pr = processes.begin(); pr != processes.end(); ++pr ) { for ( PDVector::const_iterator pp = pr->begin(); pp != pr->end(); ++pp ) { if ( (**pp).coloured() ) { colouredProcesses.insert(*pr); break; } } } if ( colouredProcesses.size() != processes.size() ) { generator()->log() << "Some or all of the generated subprocesses do not contain coloured legs.\n" << "Processes of this kind are currently not supported.\n" << flush; } if ( colouredProcesses.empty() ) { throw Exception() << "MatchboxFactory::makeMEs(): No processes with coloured legs have been found. " << "This run will be aborted." << Exception::runerror; } processes = colouredProcesses; // end unsupported processes // detect external particles with non-zero width for the hard process bool trouble = false; string troubleMaker; for ( set::const_iterator pr = processes.begin(); pr != processes.end(); ++pr ) { for ( PDVector::const_iterator pp = pr->begin(); pp != pr->end(); ++pp ) { if ( (**pp).hardProcessWidth() != ZERO ) { trouble = true; troubleMaker = (**pp).PDGName(); break; } } } if ( trouble ) { throw Exception() << "MatchboxFactory::makeMEs(): Particle '" << troubleMaker << "' appears as external\nprocess leg with non-zero " << "width to be used in the hard process calculation.\n" << "Please check your setup and consider setting HardProcessWidth to zero." << Exception::runerror; } vector::ptr> matchAmplitudes; unsigned int lowestAsOrder = allProcesses() ? 0 : orderas; unsigned int highestAsOrder = orderas; unsigned int lowestAeOrder = allProcesses() ? 0 : orderInAlphaEW(); unsigned int highestAeOrder = orderInAlphaEW(); for ( unsigned int oas = lowestAsOrder; oas <= highestAsOrder; ++oas ) { for ( unsigned int oae = lowestAeOrder; oae <= highestAeOrder; ++oae ) { for ( vector::ptr>::const_iterator amp = amplitudes().begin(); amp != amplitudes().end(); ++amp ) { if ( !theSelectedAmplitudes.empty() ) { if ( find(theSelectedAmplitudes.begin(),theSelectedAmplitudes.end(),*amp) == theSelectedAmplitudes.end() ) continue; } if ( !theDeselectedAmplitudes.empty() ) { if ( find(theDeselectedAmplitudes.begin(),theDeselectedAmplitudes.end(),*amp) != theDeselectedAmplitudes.end() ) continue; } (**amp).orderInGs(oas); (**amp).orderInGem(oae); if ( (**amp).orderInGs() != oas || (**amp).orderInGem() != oae ) { continue; } matchAmplitudes.push_back(*amp); } } } size_t combinations = processes.size()*matchAmplitudes.size(); size_t procCount = 0; generator()->log() << "building matrix elements." << flush; boost::progress_display * progressBar = new boost::progress_display(combinations,generator()->log()); for ( unsigned int oas = lowestAsOrder; oas <= highestAsOrder; ++oas ) { for ( unsigned int oae = lowestAeOrder; oae <= highestAeOrder; ++oae ) { for ( vector::ptr>::const_iterator amp = matchAmplitudes.begin(); amp != matchAmplitudes.end(); ++amp ) { (**amp).orderInGs(oas); (**amp).orderInGem(oae); for ( set::const_iterator p = processes.begin(); p != processes.end(); ++p ) { ++(*progressBar); if ( !(**amp).canHandle(*p,this,virt) ) continue; if ( (**amp).isExternal() ) externalAmplitudes().insert(*amp); ++procCount; Process proc(*p,oas,oae); ampProcs[*amp].insert(proc); procAmps[proc].insert(*amp); } } } } delete progressBar; generator()->log() << flush; bool clash = false; for ( map::ptr> >::const_iterator check = procAmps.begin(); check != procAmps.end(); ++check ) { if ( check->second.size() > 1 ) { clash = true; generator()->log() << "Several different amplitudes have been found for: " << check->first.legs[0]->PDGName() << " " << check->first.legs[1]->PDGName() << " -> "; for ( PDVector::const_iterator p = check->first.legs.begin() + 2; p != check->first.legs.end(); ++p ) generator()->log() << (**p).PDGName() << " "; generator()->log() << "at alpha_s^" << check->first.orderInAlphaS << " and alpha_ew^" << check->first.orderInAlphaEW << "\n"; generator()->log() << "The following amplitudes claim responsibility:\n"; for ( set::ptr>::const_iterator a = check->second.begin(); a != check->second.end(); ++a ) { generator()->log() << (**a).name() << " "; } generator()->log() << "\n"; } } if ( clash ) { throw Exception() << "MatchboxFactory: Ambiguous amplitude setup - please check your input files.\n" << "To avoid this problem use the SelectAmplitudes or DeselectAmplitudes interfaces.\n" << Exception::runerror; } bool canDoSpinCorrelations = true; vector::ptr> res; for ( map::ptr,set >::const_iterator ap = ampProcs.begin(); ap != ampProcs.end(); ++ap ) { canDoSpinCorrelations &= ap->first->canFillRhoMatrix(); for ( set::const_iterator m = ap->second.begin(); m != ap->second.end(); ++m ) { Ptr::ptr me = ap->first->makeME(m->legs); me->subProcess() = *m; me->amplitude(ap->first); me->matchboxAmplitude(ap->first); prepareME(me); string pname = "ME" + ap->first->name() + pid(m->legs); if ( ! (generator()->preinitRegister(me,pname) ) ) throw Exception() << "MatchboxFactory: Matrix element " << pname << " already existing." << Exception::runerror; if ( me->diagrams().empty() )continue; res.push_back(me); if ( theFirstPerturbativePDF ) theIncoming.insert(m->legs[0]->id()); if ( theSecondPerturbativePDF ) theIncoming.insert(m->legs[1]->id()); } } if ( spinCorrelations() && !canDoSpinCorrelations ) { generator()->log() << "Warning: Spin correlations have been requested, but no amplitude is " << "capable of performing these.\n"; theSpinCorrelations = false; } generator()->log() << "created " << procCount << " subprocesses.\n"; generator()->log() << "--------------------------------------------------------------------------------\n" << flush; return res; } int MatchboxFactory::orderOLPProcess(const Process& proc, Ptr::tptr amp, int type) { map,int>& procs = olpProcesses()[amp]; map,int>::const_iterator it = procs.find(make_pair(proc,type)); if ( it != procs.end() ) return it->second; int id = procs.size(); procs[make_pair(proc,type)] = id + 1; return id + 1; } void MatchboxFactory::productionMode() { if ( inProductionMode ) return; if ( !bornContributions() && !virtualContributions() && !realContributions() ) throw Exception() << "MatchboxFactory: At least one cross section contribution needs to be enabled.\n" << "Please check your setup.\n" << Exception::runerror; bool needTrueVirtuals = virtualContributions() && !meCorrectionsOnly() && !loopSimCorrections(); for ( vector::ptr>::iterator amp = amplitudes().begin(); amp != amplitudes().end(); ++amp ) { if ( !needTrueVirtuals && (**amp).oneLoopAmplitude() ) { Repository::clog() << "One-loop contributions from '" << (**amp).name() << "' are not required and will be disabled.\n" << flush; (**amp).disableOneLoop(); } } if ( subtractionData() != "" && !subProcessGroups() ) { throw Exception() << "MatchboxFactory: Plain NLO settings are required for subtraction checks.\n" << "Please check your setup.\n" << Exception::runerror; } if ( showerApproximation() && !virtualContributions() && !realContributions() ) { Repository::clog() << "Warning: Matching requested for LO run. Matching disabled.\n" << flush; showerApproximation(Ptr::tptr()); } if ( showerApproximation() && (subtractionData() != "" || subProcessGroups()) ) { Repository::clog() << "Warning: Matching requested for plain NLO run. Matching disabled.\n" << flush; showerApproximation(Ptr::tptr()); } if ( showerApproximation() ) { if ( spinCorrelations() && !showerApproximation()->hasSpinCorrelations() ) { Repository::clog() << "Warning: Spin correlations have been requested but the matching " << "object is not capable of these. Spin correlations will be turned of.\n" << flush; theSpinCorrelations = false; } } inProductionMode = true; } void MatchboxFactory::setup() { useMe(); if ( !ranSetup ) { if ( !inProductionMode ) throw Exception() << "MatchboxFactory: The MatchboxFactory object '" << name() << "' has not been switched to production mode.\n" << "Did you use 'do " << name() << ":ProductionMode' before isolating the event generator?\n" << Exception::runerror; olpProcesses().clear(); externalAmplitudes().clear(); theHighestVirtualsize = 0; theIncoming.clear(); bool needTrueVirtuals = virtualContributions() && !meCorrectionsOnly() && !loopSimCorrections(); for ( vector::ptr>::iterator amp = amplitudes().begin(); amp != amplitudes().end(); ++amp ) (**amp).factory(this); if ( bornMEs().empty() ) { if ( particleGroups().find("j") == particleGroups().end() ) throw Exception() << "MatchboxFactory: Could not find a jet particle group named 'j'" << Exception::runerror; // rebind the particle data objects for ( map::iterator g = particleGroups().begin(); g != particleGroups().end(); ++g ) for ( PDVector::iterator p = g->second.begin(); p != g->second.end(); ++p ) { #ifndef NDEBUG long checkid = (**p).id(); #endif *p = getParticleData((**p).id()); assert((**p).id() == checkid); } const PDVector& partons = particleGroups()["j"]; unsigned int nl = 0; for ( PDVector::const_iterator p = partons.begin(); p != partons.end(); ++p ) { if ( abs((**p).id()) < 7 && (**p).hardProcessMass() == ZERO ) ++nl; if ( (**p).id() > 0 && (**p).id() < 7 && (**p).hardProcessMass() == ZERO ) nLightJetVec( (**p).id() ); if ( (**p).id() > 0 && (**p).id() < 7 && (**p).hardProcessMass() != ZERO ) nHeavyJetVec( (**p).id() ); } nLight(nl/2); if ( particleGroups().find("p") == particleGroups().end() ) throw Exception() << "MatchboxFactory: Could not find a hadron particle group named 'p'" << Exception::runerror; const PDVector& partonsInP = particleGroups()["p"]; for ( PDVector::const_iterator pip = partonsInP.begin(); pip != partonsInP.end(); ++pip ) { if ( (**pip).id() > 0 && (**pip).id() < 7 && (**pip).hardProcessMass() == ZERO ) nLightProtonVec( (**pip).id() ); } vector::ptr> mes; for ( vector >::const_iterator p = processes.begin(); p != processes.end(); ++p ) { if( needTrueVirtuals ) { theHighestVirtualsize = max(theHighestVirtualsize,(int((*p).size()))); } mes = makeMEs(*p,orderInAlphaS(),needTrueVirtuals); copy(mes.begin(),mes.end(),back_inserter(bornMEs())); if ( realContributions() || meCorrectionsOnly() || (showerApproximation() && virtualContributions()) || (showerApproximation() && loopSimCorrections()) ) { if ( realEmissionProcesses.empty() ) { vector rproc = *p; rproc.push_back("j"); mes = makeMEs(rproc,orderInAlphaS()+1,false); copy(mes.begin(),mes.end(),back_inserter(realEmissionMEs())); } } } if ( realContributions() || meCorrectionsOnly() || (showerApproximation() && virtualContributions()) || (showerApproximation() && loopSimCorrections()) ) { if ( !realEmissionProcesses.empty() ) { for ( vector >::const_iterator q = realEmissionProcesses.begin(); q != realEmissionProcesses.end(); ++q ) { mes = makeMEs(*q,orderInAlphaS()+1,false); copy(mes.begin(),mes.end(),back_inserter(realEmissionMEs())); } } } } if ( loopInducedMEs().empty() ) { for ( vector >::const_iterator p = loopInducedProcesses.begin(); p != loopInducedProcesses.end(); ++p ) { vector::ptr> mes = makeMEs(*p,orderInAlphaS(),false); copy(mes.begin(),mes.end(),back_inserter(loopInducedMEs())); } } if( bornMEs().empty() && realEmissionMEs().empty() && loopInducedMEs().empty() ) throw Exception() << "MatchboxFactory: No matrix elements have been found.\n\ Please check if your order of Alpha_s and Alpha_ew have the right value.\n" << Exception::runerror; // check if we have virtual contributions bool haveVirtuals = true; // check DR conventions of virtual contributions bool virtualsAreDR = false; bool virtualsAreCDR = false; // check finite term conventions of virtual contributions bool virtualsAreCS = false; bool virtualsAreBDK = false; bool virtualsAreExpanded = false; // renormalization scheme bool virtualsAreDRbar = false; // check and prepare the Born and virtual matrix elements for ( vector::ptr>::iterator born = bornMEs().begin(); born != bornMEs().end(); ++born ) { prepareME(*born); haveVirtuals &= (**born).haveOneLoop(); if ( needTrueVirtuals ) { if ( (**born).haveOneLoop() ) { virtualsAreDRbar |= (**born).isDRbar(); virtualsAreDR |= (**born).isDR(); virtualsAreCDR |= !(**born).isDR(); virtualsAreCS |= (**born).isCS(); virtualsAreBDK |= (**born).isBDK(); virtualsAreExpanded |= (**born).isExpanded(); } } } // prepare the loop induced matrix elements for ( vector::ptr>::iterator looped = loopInducedMEs().begin(); looped != loopInducedMEs().end(); ++looped ) { prepareME(*looped); } if ( needTrueVirtuals ) { // check the additional insertion operators if ( !virtuals().empty() ) haveVirtuals = true; for ( vector::ptr>::const_iterator virt = virtuals().begin(); virt != virtuals().end(); ++virt ) { virtualsAreDRbar |= (**virt).isDRbar(); virtualsAreDR |= (**virt).isDR(); virtualsAreCDR |= !(**virt).isDR(); virtualsAreCS |= (**virt).isCS(); virtualsAreBDK |= (**virt).isBDK(); virtualsAreExpanded |= (**virt).isExpanded(); } // check for consistent conventions on virtuals, if we are to include them if ( virtualContributions() ) { if ( !haveVirtuals ) { throw Exception() << "MatchboxFactory: Could not find amplitudes for all virtual contributions needed.\n" << Exception::runerror; } if ( virtualsAreDR && virtualsAreCDR ) { throw Exception() << "MatchboxFactory: Virtual corrections use inconsistent regularization schemes.\n" << Exception::runerror; } if ( (virtualsAreCS && virtualsAreBDK) || (virtualsAreCS && virtualsAreExpanded) || (virtualsAreBDK && virtualsAreExpanded) || (!virtualsAreCS && !virtualsAreBDK && !virtualsAreExpanded) ) { throw Exception() << "MatchboxFactory: Virtual corrections use inconsistent conventions on finite terms.\n" << Exception::runerror; } } // prepare dipole insertion operators if ( virtualContributions() ) { for ( vector::ptr>::const_iterator virt = DipoleRepository::insertionIOperators(dipoleSet()).begin(); virt != DipoleRepository::insertionIOperators(dipoleSet()).end(); ++virt ) { (**virt).factory(this); if ( virtualsAreDRbar ) (**virt).useDRbar(); if ( virtualsAreDR ) (**virt).useDR(); else (**virt).useCDR(); if ( virtualsAreCS ) (**virt).useCS(); if ( virtualsAreBDK ) (**virt).useBDK(); if ( virtualsAreExpanded ) (**virt).useExpanded(); } for ( vector::ptr>::const_iterator virt = DipoleRepository::insertionPKOperators(dipoleSet()).begin(); virt != DipoleRepository::insertionPKOperators(dipoleSet()).end(); ++virt ) { (**virt).factory(this); if ( virtualsAreDRbar ) (**virt).useDRbar(); if ( virtualsAreDR ) (**virt).useDR(); else (**virt).useCDR(); if ( virtualsAreCS ) (**virt).useCS(); if ( virtualsAreBDK ) (**virt).useBDK(); if ( virtualsAreExpanded ) (**virt).useExpanded(); } } } // prepare the real emission matrix elements if ( realContributions() || meCorrectionsOnly() || (showerApproximation() && virtualContributions()) || (showerApproximation() && loopSimCorrections()) ) { for ( vector::ptr>::iterator real = realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) { prepareME(*real); } } // start creating matrix elements MEs().clear(); // setup born and virtual contributions if ( bornContributions() || virtualContributions() ) { generator()->log() << "preparing Born" << (virtualContributions() ? " and virtual" : "") << " matrix elements.\n" << flush; } if ( (bornContributions() && !virtualContributions()) || (bornContributions() && meCorrectionsOnly()) || (bornContributions() && virtualContributions() && independentVirtuals()) ) { for ( vector::ptr>::iterator born = bornMEs().begin(); born != bornMEs().end(); ++born ) { if ( (**born).onlyOneLoop() ) continue; Ptr::ptr bornme = (**born).cloneMe(); string pname = fullName() + "/" + (**born).name(); if ( virtualContributions() && independentVirtuals() ) pname += ".Born"; if ( ! (generator()->preinitRegister(bornme,pname) ) ) throw Exception() << "MatchboxFactory: Matrix element " << pname << " already existing." << Exception::runerror; if ( bornme->isOLPTree() ) { int id = orderOLPProcess(bornme->subProcess(), (**born).matchboxAmplitude(), ProcessType::treeME2); bornme->olpProcess(ProcessType::treeME2,id); } bornme->needsNoCorrelations(); bornme->cloneDependencies(); MEs().push_back(bornme); } } if ( bornContributions() && !loopInducedMEs().empty() ) { for ( vector::ptr>::iterator looped = loopInducedMEs().begin(); looped != loopInducedMEs().end(); ++looped ) { Ptr::ptr loopme = (**looped).cloneMe(); string pname = fullName() + "/" + (**looped).name() + ".LoopInduced"; if ( ! (generator()->preinitRegister(loopme,pname) ) ) throw Exception() << "MatchboxFactory: Matrix element " << pname << " already existing." << Exception::runerror; if ( loopme->isOLPTree() ) { int id = orderOLPProcess(loopme->subProcess(), (**looped).matchboxAmplitude(), ProcessType::loopInducedME2); loopme->olpProcess(ProcessType::loopInducedME2,id); } loopme->needsNoCorrelations(); loopme->cloneDependencies(); MEs().push_back(loopme); } } if ( needTrueVirtuals ) { bornVirtualMEs().clear(); boost::progress_display * progressBar = new boost::progress_display(bornMEs().size(),generator()->log()); if ( thePoleData != "" ) if ( thePoleData[thePoleData.size()-1] != '/' ) thePoleData += "/"; for ( vector::ptr>::iterator born = bornMEs().begin(); born != bornMEs().end(); ++born ) { Ptr::ptr nlo = (**born).cloneMe(); string pname = fullName() + "/" + (**born).name(); if ( !independentVirtuals() && !(!bornContributions() && virtualContributions()) ) pname += ".BornVirtual"; else if ( independentPKs() && !nlo->onlyOneLoop() ) pname += ".VirtualVI"; else pname += ".Virtual"; if ( ! (generator()->preinitRegister(nlo,pname) ) ) throw Exception() << "MatchboxFactory: NLO ME " << pname << " already existing." << Exception::runerror; nlo->virtuals().clear(); if ( !nlo->onlyOneLoop() ) { for ( vector::ptr>::const_iterator virt = virtuals().begin(); virt != virtuals().end(); ++virt ) { if ( (**virt).apply((**born).diagrams().front()->partons()) ) nlo->virtuals().push_back(*virt); } for ( vector::ptr>::const_iterator virt = DipoleRepository::insertionIOperators(dipoleSet()).begin(); virt != DipoleRepository::insertionIOperators(dipoleSet()).end(); ++virt ) { if ( (**virt).apply((**born).diagrams().front()->partons()) ) nlo->virtuals().push_back(*virt); } if ( !independentVirtuals() || ( independentVirtuals() && !independentPKs() ) ) { for ( vector::ptr>::const_iterator virt = DipoleRepository::insertionPKOperators(dipoleSet()).begin(); virt != DipoleRepository::insertionPKOperators(dipoleSet()).end(); ++virt ) { if ( (**virt).apply((**born).diagrams().front()->partons()) ) nlo->virtuals().push_back(*virt); } } if ( nlo->virtuals().empty() ) throw Exception() << "MatchboxFactory: No insertion operators have been found for " << (**born).name() << "\n" << Exception::runerror; if ( checkPoles() ) { if ( !virtualsAreExpanded ) { throw Exception() << "MatchboxFactory: Cannot check epsilon poles if virtuals are not in `expanded' convention.\n" << Exception::runerror; } } } if ( !bornContributions() || independentVirtuals() ) { nlo->doOneLoopNoBorn(); } else { nlo->doOneLoop(); } if ( nlo->isOLPLoop() ) { int id = orderOLPProcess(nlo->subProcess(), (**born).matchboxAmplitude(), ProcessType::oneLoopInterference); nlo->olpProcess(ProcessType::oneLoopInterference,id); if ( !nlo->onlyOneLoop() && nlo->needsOLPCorrelators() ) { id = orderOLPProcess(nlo->subProcess(), (**born).matchboxAmplitude(), ProcessType::colourCorrelatedME2); nlo->olpProcess(ProcessType::colourCorrelatedME2,id); } } nlo->needsCorrelations(); nlo->cloneDependencies(); bornVirtualMEs().push_back(nlo); MEs().push_back(nlo); if ( independentVirtuals() && independentPKs() && !nlo->onlyOneLoop() ) { Ptr::ptr nlopk = (**born).cloneMe(); string pnamepk = fullName() + "/" + (**born).name(); pnamepk += ".VirtualPK"; if ( ! (generator()->preinitRegister(nlopk,pnamepk) ) ) throw Exception() << "MatchboxFactory: NLO ME " << pnamepk << " already existing." << Exception::runerror; nlopk->virtuals().clear(); for ( vector::ptr>::const_iterator virt = DipoleRepository::insertionPKOperators(dipoleSet()).begin(); virt != DipoleRepository::insertionPKOperators(dipoleSet()).end(); ++virt ) { if ( (**virt).apply((**born).diagrams().front()->partons()) ) nlopk->virtuals().push_back(*virt); } if ( !nlopk->virtuals().empty() ) { nlopk->doOneLoopNoBorn(); nlopk->doOneLoopNoLoops(); if ( nlopk->isOLPLoop() ) { int id = orderOLPProcess(nlopk->subProcess(), (**born).matchboxAmplitude(), ProcessType::treeME2); nlopk->olpProcess(ProcessType::treeME2,id); if ( nlopk->needsOLPCorrelators() ) { id = orderOLPProcess(nlopk->subProcess(), (**born).matchboxAmplitude(), ProcessType::colourCorrelatedME2); nlopk->olpProcess(ProcessType::colourCorrelatedME2,id); } } nlopk->needsCorrelations(); nlopk->cloneDependencies(); bornVirtualMEs().push_back(nlopk); MEs().push_back(nlopk); } } ++(*progressBar); } delete progressBar; generator()->log() << "--------------------------------------------------------------------------------\n" << flush; } theSplittingDipoles.clear(); set bornProcs; if ( showerApproximation() ) { if ( showerApproximation()->needsSplittingGenerator() ) { for ( vector::ptr>::iterator born = bornMEs().begin(); born != bornMEs().end(); ++born ) for ( MEBase::DiagramVector::const_iterator d = (**born).diagrams().begin(); d != (**born).diagrams().end(); ++d ) bornProcs.insert((**d).partons()); } } if ( realContributions() || meCorrectionsOnly() || (showerApproximation() && virtualContributions()) || (showerApproximation() && loopSimCorrections()) ) { generator()->log() << "preparing subtracted matrix elements.\n" << flush; if ( theSubtractionData != "" ) if ( theSubtractionData[theSubtractionData.size()-1] != '/' ) theSubtractionData += "/"; subtractedMEs().clear(); for ( vector::ptr>::iterator born = bornMEs().begin(); born != bornMEs().end(); ++born ) { if ( (**born).onlyOneLoop() ) continue; (**born).needsCorrelations(); if ( (**born).isOLPTree() ) { int id = orderOLPProcess((**born).subProcess(), (**born).matchboxAmplitude(), ProcessType::colourCorrelatedME2); (**born).olpProcess(ProcessType::colourCorrelatedME2,id); bool haveGluon = false; for ( PDVector::const_iterator p = (**born).subProcess().legs.begin(); p != (**born).subProcess().legs.end(); ++p ) if ( (**p).id() == 21 ) { haveGluon = true; break; } if ( haveGluon ) { id = orderOLPProcess((**born).subProcess(), (**born).matchboxAmplitude(), ProcessType::spinColourCorrelatedME2); (**born).olpProcess(ProcessType::spinColourCorrelatedME2,id); } if ( showerApproximation() ) { id = orderOLPProcess((**born).subProcess(), (**born).matchboxAmplitude(), ProcessType::treeME2); (**born).olpProcess(ProcessType::treeME2,id); } } } boost::progress_display * progressBar = new boost::progress_display(realEmissionMEs().size(),generator()->log()); for ( vector::ptr>::iterator real = realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) { Ptr::ptr sub = new_ptr(SubtractedME()); string pname = fullName() + "/" + (**real).name() + ".SubtractedReal"; if ( ! (generator()->preinitRegister(sub,pname) ) ) throw Exception() << "MatchboxFactory: Subtracted ME " << pname << " already existing." << Exception::runerror; sub->factory(this); (**real).needsNoCorrelations(); if ( (**real).isOLPTree() ) { int id = orderOLPProcess((**real).subProcess(), (**real).matchboxAmplitude(), ProcessType::treeME2); (**real).olpProcess(ProcessType::treeME2,id); } sub->head(*real); sub->dependent().clear(); sub->getDipoles(); if ( sub->dependent().empty() ) { // finite real contribution if ( realContributions() ) { Ptr::ptr fme = dynamic_ptr_cast::ptr>(sub->head())->cloneMe(); string qname = fullName() + "/" + (**real).name() + ".FiniteReal"; if ( ! (generator()->preinitRegister(fme,qname) ) ) throw Exception() << "MatchboxFactory: ME " << qname << " already existing." << Exception::runerror; MEs().push_back(fme); finiteRealMEs().push_back(fme); } sub->head(tMEPtr()); ++(*progressBar); continue; } if ( realEmissionScales() ) sub->doRealEmissionScales(); subtractedMEs().push_back(sub); if ( realContributions() ) if ( !showerApproximation() || (showerApproximation() && showerApproximation()->hasHEvents()) ) MEs().push_back(sub); if ( showerApproximation() ) { if ( virtualContributions() && !meCorrectionsOnly() && !loopSimCorrections() ) { Ptr::ptr subv = new_ptr(*sub); string vname = sub->fullName() + ".SubtractionIntegral"; if ( ! (generator()->preinitRegister(subv,vname) ) ) throw Exception() << "MatchboxFactory: Subtracted ME " << vname << " already existing." << Exception::runerror; subv->cloneDependencies(vname); subv->doVirtualShowerSubtraction(); subtractedMEs().push_back(subv); MEs().push_back(subv); } if ( loopSimCorrections() ) { Ptr::ptr subv = new_ptr(*sub); string vname = sub->fullName() + ".SubtractionIntegral"; if ( ! (generator()->preinitRegister(subv,vname) ) ) throw Exception() << "MatchboxFactory: Subtracted ME " << vname << " already existing." << Exception::runerror; subv->cloneDependencies(vname); subv->doLoopSimSubtraction(); subtractedMEs().push_back(subv); MEs().push_back(subv); } sub->doRealShowerSubtraction(); if ( showerApproximation()->needsSplittingGenerator() ) for ( set::const_iterator p = bornProcs.begin(); p != bornProcs.end(); ++p ) { vector::ptr> sdip = sub->splitDipoles(*p); set::ptr>& dips = theSplittingDipoles[*p]; copy(sdip.begin(),sdip.end(),inserter(dips,dips.begin())); } } ++(*progressBar); } delete progressBar; generator()->log() << "--------------------------------------------------------------------------------\n" << flush; } if ( !theSplittingDipoles.empty() ) { map::ptr,Ptr::ptr> cloneMap; for ( map::ptr> >::const_iterator sd = theSplittingDipoles.begin(); sd != theSplittingDipoles.end(); ++sd ) { for ( set::ptr>::const_iterator d = sd->second.begin(); d != sd->second.end(); ++d ) { cloneMap[*d] = Ptr::ptr(); } } for ( map::ptr,Ptr::ptr>::iterator cd = cloneMap.begin(); cd != cloneMap.end(); ++cd ) { Ptr::ptr cloned = cd->first->cloneMe(); string dname = cd->first->fullName() + ".splitting"; if ( ! (generator()->preinitRegister(cloned,dname)) ) throw Exception() << "MatchboxFactory: Dipole '" << dname << "' already existing." << Exception::runerror; cloned->cloneDependencies(); cloned->showerApproximation(Ptr::tptr()); cloned->doSplitting(); cd->second = cloned; } for ( map::ptr> >::iterator sd = theSplittingDipoles.begin(); sd != theSplittingDipoles.end(); ++sd ) { set::ptr> cloned; for ( set::ptr>::iterator d = sd->second.begin(); d != sd->second.end(); ++d ) { cloned.insert(cloneMap[*d]); } sd->second = cloned; } } if ( !externalAmplitudes().empty() ) { generator()->log() << "Initializing external amplitudes.\n" << flush; for ( set::tptr>::const_iterator ext = externalAmplitudes().begin(); ext != externalAmplitudes().end(); ++ext ) { if ( !(**ext).initializeExternal() ) { throw Exception() << "Failed to initialize amplitude '" << (**ext).name() << "'\n" << Exception::runerror; } } generator()->log() << "--------------------------------------------------------------------------------\n" << flush; } if ( !olpProcesses().empty() ) { generator()->log() << "Initializing one-loop provider(s).\n" << flush; map::tptr,map,int> > olps; for ( map::tptr,map,int> >::const_iterator oit = olpProcesses().begin(); oit != olpProcesses().end(); ++oit ) { olps[oit->first] = oit->second; } for ( map::tptr,map,int> >::const_iterator olpit = olps.begin(); olpit != olps.end(); ++olpit ) { if ( !olpit->first->startOLP(olpit->second) ) { throw Exception() << "MatchboxFactory: Failed to start OLP for amplitude '" << olpit->first->name() << "'\n" << Exception::runerror; } } generator()->log() << "--------------------------------------------------------------------------------\n" << flush; } generator()->log() << "Process setup finished.\n" << flush; ranSetup = true; } } void MatchboxFactory::SplittingChannel::print(ostream& os) const { os << "--- SplittingChannel setup -----------------------------------------------------\n"; os << " Born process "; const StandardXComb& bxc = *bornXComb; os << bxc.mePartonData()[0]->PDGName() << " " << bxc.mePartonData()[1]->PDGName() << " -> "; for ( cPDVector::const_iterator p = bxc.mePartonData().begin() + 2; p != bxc.mePartonData().end(); ++p ) { os << (**p).PDGName() << " "; } os << "\n"; os << " to real emission process "; const StandardXComb& rxc = *realXComb; os << rxc.mePartonData()[0]->PDGName() << " " << rxc.mePartonData()[1]->PDGName() << " -> "; for ( cPDVector::const_iterator p = rxc.mePartonData().begin() + 2; p != rxc.mePartonData().end(); ++p ) { os << (**p).PDGName() << " "; } os << "\n"; os << " with dipole:\n"; dipole->print(os); os << "--------------------------------------------------------------------------------\n"; os << flush; } list MatchboxFactory::getSplittingChannels(tStdXCombPtr xcptr) const { if ( xcptr->lastProjector() ) xcptr = xcptr->lastProjector(); const StandardXComb& xc = *xcptr; cPDVector proc = xc.mePartonData(); map::ptr> >::const_iterator splitEntries = splittingDipoles().find(proc); list res; if ( splitEntries == splittingDipoles().end() ) return res; const set::ptr>& splitDipoles = splitEntries->second; SplittingChannel channel; if ( !splitDipoles.empty() ) { Ptr::tptr bornME = const_ptr_cast::tptr>((**splitDipoles.begin()).underlyingBornME()); channel.bornXComb = bornME->makeXComb(xc.maxEnergy(),xc.particles(),xc.eventHandlerPtr(), const_ptr_cast(xc.subProcessHandler()), xc.pExtractor(),xc.CKKWHandler(), xc.partonBins(),xc.cuts(),xc.diagrams(),xc.mirror(), PartonPairVec()); } for ( set::ptr>::const_iterator sd = splitDipoles.begin(); sd != splitDipoles.end(); ++sd ) { channel.dipole = *sd; vector realXCombs = (**sd).makeRealXCombs(channel.bornXComb); for ( vector::const_iterator rxc = realXCombs.begin(); rxc != realXCombs.end(); ++rxc ) { channel.realXComb = *rxc; if ( showerApproximation()->needsTildeXCombs() ) { channel.tildeXCombs.clear(); assert(!channel.dipole->partnerDipoles().empty()); for ( vector::tptr>::const_iterator p = channel.dipole->partnerDipoles().begin(); p != channel.dipole->partnerDipoles().end(); ++p ) { StdXCombPtr txc = channel.dipole->makeBornXComb(channel.realXComb); if ( txc ) channel.tildeXCombs.push_back(txc); } } res.push_back(channel); } } if ( initVerbose() ) { generator()->log() << "--- MatchboxFactory splitting channels ----------------------------------------------\n"; const StandardXComb& bxc = *xcptr; generator()->log() << " hard process handled is: "; generator()->log() << bxc.mePartonData()[0]->PDGName() << " " << bxc.mePartonData()[1]->PDGName() << " -> "; for ( cPDVector::const_iterator p = bxc.mePartonData().begin() + 2; p != bxc.mePartonData().end(); ++p ) { generator()->log() << (**p).PDGName() << " "; } generator()->log() << "\n"; for ( list::const_iterator sp = res.begin(); sp != res.end(); ++sp ) { sp->print(generator()->log()); } generator()->log() << "-------------------------------------------------------------------------------------\n" << flush; } return res; } void MatchboxFactory::print(ostream& os) const { os << "--- MatchboxFactory setup -----------------------------------------------------------\n"; if ( !amplitudes().empty() ) { os << " generated Born matrix elements:\n"; for ( vector::ptr>::const_iterator m = bornMEs().begin(); m != bornMEs().end(); ++m ) { (**m).print(os); } os << flush; os << " generated real emission matrix elements:\n"; for ( vector::ptr>::const_iterator m = realEmissionMEs().begin(); m != realEmissionMEs().end(); ++m ) { (**m).print(os); } os << flush; } os << " generated Born+virtual matrix elements:\n"; for ( vector::ptr>::const_iterator bv = bornVirtualMEs().begin(); bv != bornVirtualMEs().end(); ++bv ) { (**bv).print(os); } os << " generated subtracted matrix elements:\n"; for ( vector::ptr>::const_iterator sub = subtractedMEs().begin(); sub != subtractedMEs().end(); ++sub ) { os << " '" << (**sub).name() << "'\n"; } os << "--------------------------------------------------------------------------------\n"; os << flush; } +void MatchboxFactory::summary(ostream& os) const { + os << "\n\n================================================================================\n" + << " Matchbox hard process summary\n" + << "================================================================================\n\n"; + + os << " Electro-weak parameter summary:\n" + << "--------------------------------------------------------------------------------\n\n"; + + os << " Electro-weak scheme : "; + switch ( SM().ewScheme() ) { + + case 0: os << "Default"; break; + case 1: os << "GMuScheme"; break; + case 2: os << "alphaMZScheme"; break; + case 3: os << "NoMass"; break; + case 4: os << "mW"; break; + case 5: os << "mZ"; break; + case 6: os << "Independent"; break; + case 7: os << "FeynRulesUFO"; break; + default: assert(false); + + } + + os << "\n"; + + os << " alphaEM is " + << (SM().ewScheme() == 0 && !theFixedQEDCouplings ? "running" : "fixed at alphaEM(m(Z))") << "\n"; + + if ( SM().ewScheme() == 0 && !theFixedQEDCouplings ) + os << " alphaEM is running at " << SM().alphaEMPtr()->nloops() + << " loops\n\n"; + else + os << "\n"; + + os << (SM().ewScheme() != 0 ? " Tree level relations " : " Best values ") + << "yield:\n\n" + << " m(Z)/GeV = " + << getParticleData(ParticleID::Z0)->hardProcessMass()/GeV + << "\n" + << " g(Z)/GeV = " + << getParticleData(ParticleID::Z0)->hardProcessWidth()/GeV + << "\n" + << " m(W)/GeV = " + << getParticleData(ParticleID::Wplus)->hardProcessMass()/GeV + << "\n" + << " g(W)/GeV = " + << getParticleData(ParticleID::Wplus)->hardProcessWidth()/GeV + << "\n" + << " m(H)/GeV = " + << getParticleData(ParticleID::h0)->hardProcessMass()/GeV + << "\n" + << " g(H)/GeV = " + << getParticleData(ParticleID::h0)->hardProcessWidth()/GeV + << "\n" + << " alphaEM(m(Z)) = " + << SM().alphaEMME(sqr(getParticleData(ParticleID::Z0)->hardProcessMass())) << "\n" + << " sin^2(theta) = " << SM().sin2ThetaW() + << "\n" + << " GeV^2 GF = " << GeV2*SM().fermiConstant() + << "\n\n"; + + os << " Quark masses and widths are:\n" + << "--------------------------------------------------------------------------------\n\n" + << " m(u)/GeV = " << getParticleData(ParticleID::u)->hardProcessMass()/GeV << "\n" + << " m(d)/GeV = " << getParticleData(ParticleID::d)->hardProcessMass()/GeV << "\n" + << " m(c)/GeV = " << getParticleData(ParticleID::c)->hardProcessMass()/GeV << "\n" + << " m(s)/GeV = " << getParticleData(ParticleID::s)->hardProcessMass()/GeV << "\n" + << " m(t)/GeV = " << getParticleData(ParticleID::t)->hardProcessMass()/GeV << "\n" + << " g(t)/GeV = " << getParticleData(ParticleID::t)->hardProcessWidth()/GeV << "\n" + << " m(b)/GeV = " << getParticleData(ParticleID::b)->hardProcessMass()/GeV << "\n\n"; + + os << " Lepton masses and widths are:\n" + << "--------------------------------------------------------------------------------\n\n" + << " m(n_e)/GeV = " << getParticleData(ParticleID::nu_e)->hardProcessMass()/GeV << "\n" + << " m(e)/GeV = " << getParticleData(ParticleID::eminus)->hardProcessMass()/GeV << "\n" + << " m(n_mu)/GeV = " << getParticleData(ParticleID::nu_mu)->hardProcessMass()/GeV << "\n" + << " m(mu)/GeV = " << getParticleData(ParticleID::muminus)->hardProcessMass()/GeV << "\n" + << " m(nu_tau)/GeV = " << getParticleData(ParticleID::nu_tau)->hardProcessMass()/GeV << "\n" + << " m(tau)/GeV = " << getParticleData(ParticleID::tauminus)->hardProcessMass()/GeV << "\n\n"; + + + os << " Strong coupling summary:\n" + << "--------------------------------------------------------------------------------\n\n"; + + os << " alphaS is"; + if ( !theFixedCouplings ) { + os << " running at " << SM().alphaSPtr()->nloops() + << " loops with\n" + << " alphaS(m(Z)) = " << SM().alphaSPtr()->value(sqr(getParticleData(ParticleID::Z0)->mass())) + << "\n\n"; + } else { + os << " fixed at " + << SM().alphaS() + << "\n\n"; + } + + if ( !theFixedCouplings ) { + os << " flavour thresholds are matched at\n"; + for ( long id = 1; id <= 6; ++id ) { + os << " m(" << id << ")/GeV = " + << (SM().alphaSPtr()->quarkMasses().empty() ? + getParticleData(id)->mass()/GeV : + SM().alphaSPtr()->quarkMasses()[id-1]/GeV) + << "\n"; + } + } + + os << "\n\n" << flush; + +} + + void MatchboxFactory::doinit() { theIsMatchboxRun() = true; if ( RunDirectories::empty() ) RunDirectories::pushRunId(generator()->runName()); setup(); if ( theShowerApproximation ) theShowerApproximation->init(); if ( initVerbose() && !ranSetup ) print(Repository::clog()); Ptr::tptr eh = dynamic_ptr_cast::tptr>(generator()->eventHandler()); assert(eh); + if ( initVerbose() && !ranSetup ) { + assert(standardModel()); + standardModel()->init(); + summary(Repository::clog()); + } SubProcessHandler::doinit(); } void MatchboxFactory::doinitrun() { theIsMatchboxRun() = true; if ( theShowerApproximation ) theShowerApproximation->initrun(); Ptr::tptr eh = dynamic_ptr_cast::tptr>(generator()->eventHandler()); assert(eh); SubProcessHandler::doinitrun(); } const string& MatchboxFactory::buildStorage() { return RunDirectories::buildStorage(); } const string& MatchboxFactory::runStorage() { return RunDirectories::runStorage(); } void MatchboxFactory::persistentOutput(PersistentOStream & os) const { os << theDiagramGenerator << theProcessData << theNLight << theNLightJetVec << theNHeavyJetVec << theNLightProtonVec << theOrderInAlphaS << theOrderInAlphaEW << theBornContributions << theVirtualContributions << theRealContributions << theIndependentVirtuals << theIndependentPKs << theSubProcessGroups << thePhasespace << theScaleChoice << theFactorizationScaleFactor << theRenormalizationScaleFactor << theFixedCouplings << theFixedQEDCouplings << theVetoScales << theAmplitudes << theBornMEs << theVirtuals << theRealEmissionMEs << theLoopInducedMEs << theBornVirtualMEs << theSubtractedMEs << theFiniteRealMEs << theVerbose<> theDiagramGenerator >> theProcessData >> theNLight >> theNLightJetVec >> theNHeavyJetVec >> theNLightProtonVec >> theOrderInAlphaS >> theOrderInAlphaEW >> theBornContributions >> theVirtualContributions >> theRealContributions >> theIndependentVirtuals >> theIndependentPKs >> theSubProcessGroups >> thePhasespace >> theScaleChoice >> theFactorizationScaleFactor >> theRenormalizationScaleFactor >> theFixedCouplings >> theFixedQEDCouplings >> theVetoScales >> theAmplitudes >> theBornMEs >> theVirtuals >> theRealEmissionMEs >> theLoopInducedMEs >> theBornVirtualMEs >> theSubtractedMEs >> theFiniteRealMEs >> theVerbose >> theDiagramWeightVerbose >> theDiagramWeightVerboseNBins >> theInitVerbose >> theSubtractionData >> theSubtractionPlotType >> theSubtractionScatterPlot >> thePoleData >> theParticleGroups >> processes >> loopInducedProcesses >> realEmissionProcesses >> theShowerApproximation >> theSplittingDipoles >> theRealEmissionScales >> theAllProcesses >> theOLPProcesses >> theExternalAmplitudes >> theSelectedAmplitudes >> theDeselectedAmplitudes >> theDipoleSet >> theReweighters >> thePreweighters >> theMECorrectionsOnly>> theLoopSimCorrections>>theHighestVirtualsize >> ranSetup >> theIncoming >> theFirstPerturbativePDF >> theSecondPerturbativePDF >> inProductionMode >> theSpinCorrelations >> theAlphaParameter >> theEnforceChargeConservation >> theEnforceColourConservation >> theEnforceLeptonNumberConservation >> theEnforceQuarkNumberConservation >> theLeptonFlavourDiagonal >> theQuarkFlavourDiagonal; } string MatchboxFactory::startParticleGroup(string name) { particleGroupName = StringUtils::stripws(name); particleGroup.clear(); return ""; } string MatchboxFactory::endParticleGroup(string) { if ( particleGroup.empty() ) throw Exception() << "MatchboxFactory: Empty particle group." << Exception::runerror; particleGroups()[particleGroupName] = particleGroup; particleGroup.clear(); return ""; } vector MatchboxFactory::parseProcess(string in) { vector process = StringUtils::split(in); if ( process.size() < 3 ) throw Exception() << "MatchboxFactory: Invalid process." << Exception::runerror; for ( vector::iterator p = process.begin(); p != process.end(); ++p ) { *p = StringUtils::stripws(*p); } vector pprocess; for ( vector::const_iterator p = process.begin(); p != process.end(); ++p ) { if ( *p == "->" ) continue; pprocess.push_back(*p); } return pprocess; } string MatchboxFactory::doProcess(string in) { processes.push_back(parseProcess(in)); return ""; } string MatchboxFactory::doLoopInducedProcess(string in) { loopInducedProcesses.push_back(parseProcess(in)); return ""; } string MatchboxFactory::doSingleRealProcess(string in) { realEmissionProcesses.push_back(parseProcess(in)); return ""; } struct SortPID { inline bool operator()(PDPtr a, PDPtr b) const { return a->id() < b->id(); } }; // // @TODO // // SP: After improving this for standard model process building this should // actually got into a separate process builder class or something along these // lines to have it better factored for use with BSM models. // // set MatchboxFactory:: makeSubProcesses(const vector& proc) const { if ( proc.empty() ) throw Exception() << "MatchboxFactory: No process specified." << Exception::runerror; vector groups; typedef map::const_iterator GroupIterator; for ( vector::const_iterator gr = proc.begin(); gr != proc.end(); ++gr ) { GroupIterator git = particleGroups().find(*gr); if ( git == particleGroups().end() ) { throw Exception() << "MatchboxFactory: Particle group '" << *gr << "' not defined." << Exception::runerror; } groups.push_back(git->second); } vector counts(groups.size(),0); PDVector proto(groups.size()); set allProcs; while ( true ) { for ( size_t k = 0; k < groups.size(); ++k ) proto[k] = groups[k][counts[k]]; int charge = 0; int colour = 0; int nleptons = 0; int nquarks = 0; int ncolour = 0; int nleptonsGen[4]; int nquarksGen[4]; for ( size_t i = 0; i < 4; ++i ) { nleptonsGen[i] = 0; nquarksGen[i] = 0; } for ( size_t k = 0; k < proto.size(); ++k ) { int sign = k > 1 ? 1 : -1; charge += sign * proto[k]->iCharge(); colour += sign * proto[k]->iColour(); if ( abs(proto[k]->id()) <= 8 ) { int generation = (abs(proto[k]->id()) - 1)/2; nquarks += sign * ( proto[k]->id() < 0 ? -1 : 1); nquarksGen[generation] += sign * ( proto[k]->id() < 0 ? -1 : 1); } if ( abs(proto[k]->id()) > 10 && abs(proto[k]->id()) <= 18 ) { int generation = (abs(proto[k]->id()) - 11)/2; nleptons += sign * ( proto[k]->id() < 0 ? -1 : 1); nleptonsGen[generation] += sign * ( proto[k]->id() < 0 ? -1 : 1); } if ( proto[k]->coloured() ) ++ncolour; } bool pass = true; if ( theEnforceChargeConservation ) pass &= (charge == 0); if ( theEnforceColourConservation ) pass &= (colour % 8 == 0) && (ncolour > 1); if ( theEnforceLeptonNumberConservation ) { pass &= (nleptons == 0); if ( theLeptonFlavourDiagonal ) { for ( size_t i = 0; i < 4; ++i ) pass &= (nleptonsGen[i] == 0); } } if ( theEnforceQuarkNumberConservation ) { pass &= (nquarks == 0); if ( theQuarkFlavourDiagonal ) { for ( size_t i = 0; i < 4; ++i ) pass &= (nquarksGen[i] == 0); } } if ( pass ) { for ( int i = 0; i < 2; ++i ) { if ( proto[i]->coloured() && proto[i]->hardProcessMass() != ZERO ) throw Exception() << "Inconsistent flavour scheme detected with massive incoming " << proto[i]->PDGName() << ". Check your setup." << Exception::runerror; } sort(proto.begin()+2,proto.end(),SortPID()); allProcs.insert(proto); } vector::reverse_iterator c = counts.rbegin(); vector::const_reverse_iterator g = groups.rbegin(); while ( c != counts.rend() ) { if ( ++(*c) == g->size() ) { *c = 0; ++c; ++g; } else { break; } } if ( c == counts.rend() ) break; } return allProcs; } void MatchboxFactory::Init() { static ClassDocumentation documentation ("MatchboxFactory", "NLO QCD corrections have been calculated " "using Matchbox \\cite{Platzer:2011bc}, \\cite{Matchbox:2015}", "%\\cite{Platzer:2011bc}\n" "\\bibitem{Platzer:2011bc}\n" "S.~Platzer and S.~Gieseke,\n" "``Dipole Showers and Automated NLO Matching in Herwig,''\n" "arXiv:1109.6256 [hep-ph].\n" "%%CITATION = ARXIV:1109.6256;%%\n" "%\\cite{Matchbox:2015}\n" "\\bibitem{Matchbox:2015}\n" "Herwig collaboration,\n" "``Precision LHC Event Generation with Herwig,''\n" "in preparation."); static Reference interfaceDiagramGenerator ("DiagramGenerator", "Set the diagram generator.", &MatchboxFactory::theDiagramGenerator, false, false, true, true, false); static Reference interfaceProcessData ("ProcessData", "Set the process data object to be used.", &MatchboxFactory::theProcessData, false, false, true, true, false); static Parameter interfaceOrderInAlphaS ("OrderInAlphaS", "The order in alpha_s to consider.", &MatchboxFactory::theOrderInAlphaS, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfaceOrderInAlphaEW ("OrderInAlphaEW", "The order in alpha_EW", &MatchboxFactory::theOrderInAlphaEW, 2, 0, 0, false, false, Interface::lowerlim); static Switch interfaceBornContributions ("BornContributions", "Switch on or off the Born contributions.", &MatchboxFactory::theBornContributions, true, false, false); static SwitchOption interfaceBornContributionsOn (interfaceBornContributions, "On", "Switch on Born contributions.", true); static SwitchOption interfaceBornContributionsOff (interfaceBornContributions, "Off", "Switch off Born contributions.", false); static Switch interfaceVirtualContributions ("VirtualContributions", "Switch on or off the virtual contributions.", &MatchboxFactory::theVirtualContributions, true, false, false); static SwitchOption interfaceVirtualContributionsOn (interfaceVirtualContributions, "On", "Switch on virtual contributions.", true); static SwitchOption interfaceVirtualContributionsOff (interfaceVirtualContributions, "Off", "Switch off virtual contributions.", false); static Switch interfaceRealContributions ("RealContributions", "Switch on or off the real contributions.", &MatchboxFactory::theRealContributions, true, false, false); static SwitchOption interfaceRealContributionsOn (interfaceRealContributions, "On", "Switch on real contributions.", true); static SwitchOption interfaceRealContributionsOff (interfaceRealContributions, "Off", "Switch off real contributions.", false); static Switch interfaceIndependentVirtuals ("IndependentVirtuals", "Switch on or off virtual contributions as separate subprocesses.", &MatchboxFactory::theIndependentVirtuals, true, false, false); static SwitchOption interfaceIndependentVirtualsOn (interfaceIndependentVirtuals, "On", "Switch on virtual contributions as separate subprocesses.", true); static SwitchOption interfaceIndependentVirtualsOff (interfaceIndependentVirtuals, "Off", "Switch off virtual contributions as separate subprocesses.", false); static Switch interfaceIndependentPKs ("IndependentPKOperators", "Switch on or off PK oeprators as separate subprocesses.", &MatchboxFactory::theIndependentPKs, true, false, false); static SwitchOption interfaceIndependentPKsOn (interfaceIndependentPKs, "On", "Switch on PK operators as separate subprocesses.", true); static SwitchOption interfaceIndependentPKsOff (interfaceIndependentPKs, "Off", "Switch off PK operators as separate subprocesses.", false); static Switch interfaceSubProcessGroups ("SubProcessGroups", "Switch on or off production of sub-process groups.", &MatchboxFactory::theSubProcessGroups, false, false, false); static SwitchOption interfaceSubProcessGroupsOn (interfaceSubProcessGroups, "On", "On", true); static SwitchOption interfaceSubProcessGroupsOff (interfaceSubProcessGroups, "Off", "Off", false); static Reference interfacePhasespace ("Phasespace", "Set the phasespace generator.", &MatchboxFactory::thePhasespace, false, false, true, true, false); static Reference interfaceScaleChoice ("ScaleChoice", "Set the scale choice object.", &MatchboxFactory::theScaleChoice, false, false, true, true, false); static Parameter interfaceFactorizationScaleFactor ("FactorizationScaleFactor", "The factorization scale factor.", &MatchboxFactory::theFactorizationScaleFactor, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceRenormalizationScaleFactor ("RenormalizationScaleFactor", "The renormalization scale factor.", &MatchboxFactory::theRenormalizationScaleFactor, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Switch interfaceFixedCouplings ("FixedCouplings", "Switch on or off fixed couplings.", &MatchboxFactory::theFixedCouplings, true, false, false); static SwitchOption interfaceFixedCouplingsOn (interfaceFixedCouplings, "On", "On", true); static SwitchOption interfaceFixedCouplingsOff (interfaceFixedCouplings, "Off", "Off", false); static Switch interfaceFixedQEDCouplings ("FixedQEDCouplings", "Switch on or off fixed QED couplings.", &MatchboxFactory::theFixedQEDCouplings, true, false, false); static SwitchOption interfaceFixedQEDCouplingsOn (interfaceFixedQEDCouplings, "On", "On", true); static SwitchOption interfaceFixedQEDCouplingsOff (interfaceFixedQEDCouplings, "Off", "Off", false); static Switch interfaceVetoScales ("VetoScales", "Switch on or setting veto scales.", &MatchboxFactory::theVetoScales, false, false, false); static SwitchOption interfaceVetoScalesOn (interfaceVetoScales, "On", "On", true); static SwitchOption interfaceVetoScalesOff (interfaceVetoScales, "Off", "Off", false); static RefVector interfaceAmplitudes ("Amplitudes", "The amplitude objects.", &MatchboxFactory::theAmplitudes, -1, false, false, true, true, false); static RefVector interfaceBornMEs ("BornMEs", "The Born matrix elements to be used", &MatchboxFactory::theBornMEs, -1, false, false, true, true, false); static RefVector interfaceVirtuals ("Virtuals", "The virtual corrections to include", &MatchboxFactory::theVirtuals, -1, false, false, true, true, false); static RefVector interfaceRealEmissionMEs ("RealEmissionMEs", "The RealEmission matrix elements to be used", &MatchboxFactory::theRealEmissionMEs, -1, false, false, true, true, false); static RefVector interfaceBornVirtuals ("BornVirtualMEs", "The generated Born/virtual contributions", &MatchboxFactory::theBornVirtualMEs, -1, false, true, true, true, false); static RefVector interfaceSubtractedMEs ("SubtractedMEs", "The generated subtracted real emission contributions", &MatchboxFactory::theSubtractedMEs, -1, false, true, true, true, false); static RefVector interfaceFiniteRealMEs ("FiniteRealMEs", "The generated finite real contributions", &MatchboxFactory::theFiniteRealMEs, -1, false, true, true, true, false); static Switch interfaceVerbose ("Verbose", "Print full infomation on each evaluated phase space point.", &MatchboxFactory::theVerbose, false, false, false); static SwitchOption interfaceVerboseOn (interfaceVerbose, "On", "On", true); static SwitchOption interfaceVerboseOff (interfaceVerbose, "Off", "Off", false); static Switch interfaceVerboseDia ("DiagramWeightVerbose", "Print full infomation on each evaluated phase space point.", &MatchboxFactory::theDiagramWeightVerbose, false, false, false); static SwitchOption interfaceVerboseDiaOn (interfaceVerboseDia, "On", "On", true); static SwitchOption interfaceVerboseDiaOff (interfaceVerboseDia, "Off", "Off", false); static Parameter interfaceVerboseDiaNbins ("DiagramWeightVerboseNBins", "No. of Bins for DiagramWeightVerbose Diagrams.", &MatchboxFactory::theDiagramWeightVerboseNBins, 200, 0, 0, false, false, Interface::lowerlim); static Switch interfaceInitVerbose ("InitVerbose", "Print setup information.", &MatchboxFactory::theInitVerbose, false, false, false); static SwitchOption interfaceInitVerboseOn (interfaceInitVerbose, "On", "On", true); static SwitchOption interfaceInitVerboseOff (interfaceInitVerbose, "Off", "Off", false); static Parameter interfaceSubtractionData ("SubtractionData", "Prefix for subtraction check data.", &MatchboxFactory::theSubtractionData, "", false, false); static Switch interfaceSubtractionPlotType ("SubtractionPlotType", "Switch for controlling what kind of plot is generated for checking the subtraction", &MatchboxFactory::theSubtractionPlotType, 1, false, false); static SwitchOption interfaceSubtractionPlotTypeLinearRatio (interfaceSubtractionPlotType, "LinRatio", "Switch on the linear plot of the ratio", 1); static SwitchOption interfaceSubtractionPlotTypeLogRelDiff (interfaceSubtractionPlotType, "LogRelDiff", "Switch on the logarithmic plot of the relative difference", 2); static Switch interfaceSubtractionScatterPlot ("SubtractionScatterPlot", "Switch for controlling whether subtraction data should be plotted for each phase space point individually", &MatchboxFactory::theSubtractionScatterPlot, false, false, false); static SwitchOption interfaceSubtractionScatterPlotOff (interfaceSubtractionScatterPlot, "Off", "Switch off the scatter plot", false); static SwitchOption interfaceSubtractionScatterPlotOn (interfaceSubtractionScatterPlot, "On", "Switch on the scatter plot", true); static Parameter interfacePoleData ("PoleData", "Prefix for subtraction check data.", &MatchboxFactory::thePoleData, "", false, false); static RefVector interfaceParticleGroup ("ParticleGroup", "The particle group just started.", &MatchboxFactory::particleGroup, -1, false, false, true, false, false); static Command interfaceStartParticleGroup ("StartParticleGroup", "Start a particle group.", &MatchboxFactory::startParticleGroup, false); static Command interfaceEndParticleGroup ("EndParticleGroup", "End a particle group.", &MatchboxFactory::endParticleGroup, false); static Command interfaceProcess ("Process", "Set the process(es) to consider.", &MatchboxFactory::doProcess, false); static Command interfaceLoopInducedProcess ("LoopInducedProcess", "Set the loop induced process(es) to consider.", &MatchboxFactory::doLoopInducedProcess, false); static Command interfaceSingleRealProcess ("SingleRealProcess", "Set the real emission process(es) to consider.", &MatchboxFactory::doSingleRealProcess, false); static Reference interfaceShowerApproximation ("ShowerApproximation", "Set the shower approximation to be considered.", &MatchboxFactory::theShowerApproximation, false, false, true, true, false); static Switch interfaceRealEmissionScales ("RealEmissionScales", "Switch on or off calculation of subtraction scales from real emission kinematics.", &MatchboxFactory::theRealEmissionScales, false, false, false); static SwitchOption interfaceRealEmissionScalesOn (interfaceRealEmissionScales, "On", "On", true); static SwitchOption interfaceRealEmissionScalesOff (interfaceRealEmissionScales, "Off", "Off", false); static Switch interfaceAllProcesses ("AllProcesses", "Consider all processes up to a maximum coupling order specified by the coupling order interfaces.", &MatchboxFactory::theAllProcesses, false, false, false); static SwitchOption interfaceAllProcessesYes (interfaceAllProcesses, "Yes", "Include all processes.", true); static SwitchOption interfaceAllProcessesNo (interfaceAllProcesses, "No", "Only consider processes matching the exact order in the couplings.", false); static RefVector interfaceSelectAmplitudes ("SelectAmplitudes", "The amplitude objects to be favoured in clashing responsibilities.", &MatchboxFactory::theSelectedAmplitudes, -1, false, false, true, true, false); static RefVector interfaceDeselectAmplitudes ("DeselectAmplitudes", "The amplitude objects to be disfavoured in clashing responsibilities.", &MatchboxFactory::theDeselectedAmplitudes, -1, false, false, true, true, false); static Switch interfaceDipoleSet ("DipoleSet", "The set of subtraction terms to be considered.", &MatchboxFactory::theDipoleSet, 0, false, false); static SwitchOption interfaceDipoleSetCataniSeymour (interfaceDipoleSet, "CataniSeymour", "Use default Catani-Seymour dipoles.", 0); static RefVector interfaceReweighters ("Reweighters", "Reweight objects for matrix elements.", &MatchboxFactory::theReweighters, -1, false, false, true, false, false); static RefVector interfacePreweighters ("Preweighters", "Preweight objects for matrix elements.", &MatchboxFactory::thePreweighters, -1, false, false, true, false, false); static Switch interfaceMECorrectionsOnly ("MECorrectionsOnly", "Prepare only ME corrections, but no NLO calculation.", &MatchboxFactory::theMECorrectionsOnly, false, false, false); static SwitchOption interfaceMECorrectionsOnlyYes (interfaceMECorrectionsOnly, "Yes", "Produce only ME corrections.", true); static SwitchOption interfaceMECorrectionsOnlyNo (interfaceMECorrectionsOnly, "No", "Produce full NLO.", false); static Switch interfaceLoopSimCorrections ("LoopSimCorrections", "Prepare LoopSim corrections.", &MatchboxFactory::theLoopSimCorrections, false, false, false); static SwitchOption interfaceLoopSimCorrectionsYes (interfaceLoopSimCorrections, "Yes", "Produce loopsim corrections.", true); static SwitchOption interfaceLoopSimCorrectionsNo (interfaceLoopSimCorrections, "No", "Produce full NLO.", false); static Switch interfaceFirstPerturbativePDF ("FirstPerturbativePDF", "", &MatchboxFactory::theFirstPerturbativePDF, true, false, false); static SwitchOption interfaceFirstPerturbativePDFYes (interfaceFirstPerturbativePDF, "Yes", "", true); static SwitchOption interfaceFirstPerturbativePDFNo (interfaceFirstPerturbativePDF, "No", "", false); static Switch interfaceSecondPerturbativePDF ("SecondPerturbativePDF", "", &MatchboxFactory::theSecondPerturbativePDF, true, false, false); static SwitchOption interfaceSecondPerturbativePDFYes (interfaceSecondPerturbativePDF, "Yes", "", true); static SwitchOption interfaceSecondPerturbativePDFNo (interfaceSecondPerturbativePDF, "No", "", false); static Command interfaceProductionMode ("ProductionMode", "Switch this factory to production mode.", &MatchboxFactory::doProductionMode, false); static Switch interfaceSpinCorrelations ("SpinCorrelations", "Fill information for the spin correlations, if possible.", &MatchboxFactory::theSpinCorrelations, false, false, false); static SwitchOption interfaceSpinCorrelationsYes (interfaceSpinCorrelations, "Yes", "", true); static SwitchOption interfaceSpinCorrelationsNo (interfaceSpinCorrelations, "No", "", false); static Parameter interfaceAlphaParameter ("AlphaParameter", "Nagy-AlphaParameter.", &MatchboxFactory::theAlphaParameter, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Switch interfaceEnforceChargeConservation ("EnforceChargeConservation", "Enforce charge conservation while generating the hard process.", &MatchboxFactory::theEnforceChargeConservation, true, false, false); static SwitchOption interfaceEnforceChargeConservationYes (interfaceEnforceChargeConservation, "Yes", "Enforce charge conservation.", true); static SwitchOption interfaceEnforceChargeConservationNo (interfaceEnforceChargeConservation, "No", "Do not enforce charge conservation.", false); static Switch interfaceEnforceColourConservation ("EnforceColourConservation", "Enforce colour conservation while generating the hard process.", &MatchboxFactory::theEnforceColourConservation, false, false, false); static SwitchOption interfaceEnforceColourConservationYes (interfaceEnforceColourConservation, "Yes", "Enforce colour conservation.", true); static SwitchOption interfaceEnforceColourConservationNo (interfaceEnforceColourConservation, "No", "Do not enforce colour conservation.", false); static Switch interfaceEnforceLeptonNumberConservation ("EnforceLeptonNumberConservation", "Enforce lepton number conservation while generating the hard process.", &MatchboxFactory::theEnforceLeptonNumberConservation, false, false, false); static SwitchOption interfaceEnforceLeptonNumberConservationYes (interfaceEnforceLeptonNumberConservation, "Yes", "Enforce lepton number conservation.", true); static SwitchOption interfaceEnforceLeptonNumberConservationNo (interfaceEnforceLeptonNumberConservation, "No", "Do not enforce lepton number conservation.", false); static Switch interfaceEnforceQuarkNumberConservation ("EnforceQuarkNumberConservation", "Enforce quark number conservation while generating the hard process.", &MatchboxFactory::theEnforceQuarkNumberConservation, false, false, false); static SwitchOption interfaceEnforceQuarkNumberConservationYes (interfaceEnforceQuarkNumberConservation, "Yes", "Enforce quark number conservation.", true); static SwitchOption interfaceEnforceQuarkNumberConservationNo (interfaceEnforceQuarkNumberConservation, "No", "Do not enforce quark number conservation.", false); static Switch interfaceLeptonFlavourDiagonal ("LeptonFlavourDiagonal", "Assume that lepton interactions are flavour diagonal while generating the hard process.", &MatchboxFactory::theLeptonFlavourDiagonal, false, false, false); static SwitchOption interfaceLeptonFlavourDiagonalYes (interfaceLeptonFlavourDiagonal, "Yes", "Assume that lepton interactions are flavour diagonal.", true); static SwitchOption interfaceLeptonFlavourDiagonalNo (interfaceLeptonFlavourDiagonal, "No", "Do not assume that lepton interactions are flavour diagonal.", false); static Switch interfaceQuarkFlavourDiagonal ("QuarkFlavourDiagonal", "Assume that quark interactions are flavour diagonal while generating the hard process.", &MatchboxFactory::theQuarkFlavourDiagonal, false, false, false); static SwitchOption interfaceQuarkFlavourDiagonalYes (interfaceQuarkFlavourDiagonal, "Yes", "Assume that quark interactions are flavour diagonal.", true); static SwitchOption interfaceQuarkFlavourDiagonalNo (interfaceQuarkFlavourDiagonal, "No", "Do not assume that quark interactions are flavour diagonal.", false); } // *** 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 describeHerwigMatchboxFactory("Herwig::MatchboxFactory", "Herwig.so"); diff --git a/MatrixElement/Matchbox/MatchboxFactory.h b/MatrixElement/Matchbox/MatchboxFactory.h --- a/MatrixElement/Matchbox/MatchboxFactory.h +++ b/MatrixElement/Matchbox/MatchboxFactory.h @@ -1,1285 +1,1292 @@ // -*- C++ -*- // // MatchboxFactory.h 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. // #ifndef HERWIG_MatchboxFactory_H #define HERWIG_MatchboxFactory_H // // This is the declaration of the MatchboxFactory class. // #include "ThePEG/Handlers/SubProcessHandler.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h" #include "Herwig/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h" #include "Herwig/MatrixElement/Matchbox/Utility/ProcessData.h" #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h" #include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief MatchboxFactory automatically sets up a NLO * QCD calculation carried out in dipole subtraction. * * @see \ref MatchboxFactoryInterfaces "The interfaces" * defined for MatchboxFactory. */ class MatchboxFactory: public SubProcessHandler { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ MatchboxFactory(); /** * The destructor. */ virtual ~MatchboxFactory(); //@} public: /** * Flag to indicate that at least one MatchboxFactory object is in action */ static bool isMatchboxRun() { return theIsMatchboxRun(); } /** @name Process and diagram information */ //@{ /** * Return the diagram generator. */ Ptr::tptr diagramGenerator() const { return theDiagramGenerator; } /** * Set the diagram generator. */ void diagramGenerator(Ptr::ptr dg) { theDiagramGenerator = dg; } /** * Return the process data. */ Ptr::tptr processData() const { return theProcessData; } /** * Set the process data. */ void processData(Ptr::ptr pd) { theProcessData = pd; } /** * Return the number of light flavours, this matrix * element is calculated for. */ unsigned int nLight() const { return theNLight; } /** * Set the number of light flavours, this matrix * element is calculated for. */ void nLight(unsigned int n) { theNLight = n; } /** * Return the vector that contains the PDG ids of * the light flavours, which are contained in the * jet particle group. */ vector nLightJetVec() const { return theNLightJetVec; } /** * Set the elements of the vector that contains the PDG * ids of the light flavours, which are contained in the * jet particle group. */ void nLightJetVec(int n) { theNLightJetVec.push_back(n); } /** * Return the vector that contains the PDG ids of * the heavy flavours, which are contained in the * jet particle group. */ vector nHeavyJetVec() const { return theNHeavyJetVec; } /** * Set the elements of the vector that contains the PDG * ids of the heavy flavours, which are contained in the * jet particle group. */ void nHeavyJetVec(int n) { theNHeavyJetVec.push_back(n); } /** * Return the vector that contains the PDG ids of * the light flavours, which are contained in the * proton particle group. */ vector nLightProtonVec() const { return theNLightProtonVec; } /** * Set the elements of the vector that contains the PDG * ids of the light flavours, which are contained in the * proton particle group. */ void nLightProtonVec(int n) { theNLightProtonVec.push_back(n); } /** * Return the order in \f$\alpha_S\f$. */ unsigned int orderInAlphaS() const { return theOrderInAlphaS; } /** * Set the order in \f$\alpha_S\f$. */ void orderInAlphaS(unsigned int o) { theOrderInAlphaS = o; } /** * Return the order in \f$\alpha_{EM}\f$. */ unsigned int orderInAlphaEW() const { return theOrderInAlphaEW; } /** * Set the order in \f$\alpha_{EM}\f$. */ void orderInAlphaEW(unsigned int o) { theOrderInAlphaEW = o; } /** * Return true, if all processes up to a maximum order are considered */ bool allProcesses() const { return theAllProcesses; } /** * Switch on/off inclusino off all processes up to a maximum order */ void setAllProcesses(bool on = true) { theAllProcesses = on; } /** * Return true, if Born contributions should be included. */ bool bornContributions() const { return theBornContributions; } /** * Switch on or off Born contributions */ void setBornContributions(bool on = true) { theBornContributions = on; } /** * Return true, if virtual contributions should be included. */ bool virtualContributions() const { return theVirtualContributions; } /** * Switch on or off virtual contributions */ void setVirtualContributions(bool on = true) { theVirtualContributions = on; } /** * Produce matrix element corrections, but no NLO */ bool meCorrectionsOnly() const { return theMECorrectionsOnly; } /** * Switch to produce matrix element corrections, but no NLO */ void setMECorrectionsOnly(bool on = true) { theMECorrectionsOnly = on; } /** * Produce matrix element corrections, with LoopSim NLO */ bool loopSimCorrections() const { return theLoopSimCorrections; } /** * Switch to produce matrix element corrections, with LoopSim NLO */ void setLoopSimCorrections(bool on = true) { theLoopSimCorrections = on; } /** * Return true, if subtracted real emission contributions should be included. */ bool realContributions() const { return theRealContributions; } /** * Switch on or off subtracted real emission contributions */ void setRealContributions(bool on = true) { theRealContributions = on; } /** * Return true, if virtual contributions should be treated as independent subprocesses */ bool independentVirtuals() const { return theIndependentVirtuals; } /** * Switch on/off virtual contributions should be treated as independent subprocesses */ void setIndependentVirtuals(bool on = true) { theIndependentVirtuals = on; } /** * Return true, if PK operator contributions should be treated as independent subprocesses */ bool independentPKs() const { return theIndependentPKs; } /** * Switch on/off PK operator contributions should be treated as independent subprocesses */ void setIndependentPKs(bool on = true) { theIndependentPKs = on; } /** * Return true, if SubProcessGroups should be * setup from this MEGroup. If not, a single SubProcess * is constructed from the data provided by the * head matrix element. */ bool subProcessGroups() const { return theSubProcessGroups; } /** * Switch on or off producing subprocess groups. */ void setSubProcessGroups(bool on = true) { theSubProcessGroups = on; } /** * Return true, if subtraction scales should be caluclated from real emission kinematics */ bool realEmissionScales() const { return theRealEmissionScales; } /** * Switch on/off that subtraction scales should be caluclated from real emission kinematics */ void setRealEmissionScales(bool on = true) { theRealEmissionScales = on; } /** * Set the shower approximation. */ void showerApproximation(Ptr::tptr app) { theShowerApproximation = app; } /** * Return the shower approximation. */ Ptr::tptr showerApproximation() const { return theShowerApproximation; } //@} /** @name Phasespace generation and scale choice */ //@{ /** * Return the phase space generator to be used. */ Ptr::tptr phasespace() const { return thePhasespace; } /** * Set the phase space generator to be used. */ void phasespace(Ptr::ptr ps) { thePhasespace = ps; } /** * Set the scale choice object */ void scaleChoice(Ptr::ptr sc) { theScaleChoice = sc; } /** * Return the scale choice object */ Ptr::tptr scaleChoice() const { return theScaleChoice; } /** * Get the factorization scale factor */ double factorizationScaleFactor() const { return theFactorizationScaleFactor; } /** * Set the factorization scale factor */ void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; } /** * Get the renormalization scale factor */ double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; } /** * Set the renormalization scale factor */ void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; } /** * Return true, if fixed couplings are used. */ bool fixedCouplings() const { return theFixedCouplings; } /** * Switch on fixed couplings. */ void setFixedCouplings(bool on = true) { theFixedCouplings = on; } /** * Return true, if fixed couplings are used. */ bool fixedQEDCouplings() const { return theFixedQEDCouplings; } /** * Switch on fixed couplings. */ void setFixedQEDCouplings(bool on = true) { theFixedQEDCouplings = on; } /** * Return true, if veto scales should be set * for the real emission */ bool vetoScales() const { return theVetoScales; } /** * Switch on setting veto scales */ void doVetoScales() { theVetoScales = true; } /** * Switch off setting veto scales */ void noVetoScales() { theVetoScales = true; } //@} /** @name Amplitudes and caching */ //@{ /** * Return the amplitudes to be considered */ const vector::ptr>& amplitudes() const { return theAmplitudes; } /** * Access the amplitudes to be considered */ vector::ptr>& amplitudes() { return theAmplitudes; } //@} /** @name Matrix element objects. */ //@{ /** * Return the Born matrix elements to be considered */ const vector::ptr>& bornMEs() const { return theBornMEs; } /** * Access the Born matrix elements to be considered */ vector::ptr>& bornMEs() { return theBornMEs; } /** * Return the loop induced matrix elements to be considered */ const vector::ptr>& loopInducedMEs() const { return theLoopInducedMEs; } /** * Access the loop induced matrix elements to be considered */ vector::ptr>& loopInducedMEs() { return theLoopInducedMEs; } /** * Return the processes to be ordered from an OLP */ const map::tptr, map,int> >& olpProcesses() const { return theOLPProcesses; } /** * Access the processes to be ordered from an OLP */ map::tptr, map,int> >& olpProcesses() { return theOLPProcesses; } /** * Order an OLP process and return its id */ int orderOLPProcess(const Process& p, Ptr::tptr amp, int type); /** * Return the amplitudes which need external initialization */ const set::tptr>& externalAmplitudes() const { return theExternalAmplitudes; } /** * Access the amplitudes which need external initialization */ set::tptr>& externalAmplitudes() { return theExternalAmplitudes; } /** * Return the virtual corrections to be considered */ const vector::ptr>& virtuals() const { return theVirtuals; } /** * Access the virtual corrections to be considered */ vector::ptr>& virtuals() { return theVirtuals; } /** * Return the produced NLO matrix elements */ const vector::ptr>& bornVirtualMEs() const { return theBornVirtualMEs; } /** * Access the produced NLO matrix elements */ vector::ptr>& bornVirtualMEs() { return theBornVirtualMEs; } /** * Return the real emission matrix elements to be considered */ const vector::ptr>& realEmissionMEs() const { return theRealEmissionMEs; } /** * Access the real emission matrix elements to be considered */ vector::ptr>& realEmissionMEs() { return theRealEmissionMEs; } /** * Return, which set of dipoles should be considered */ int dipoleSet() const { return theDipoleSet; } /** * Return, which set of dipoles should be considered */ void dipoleSet(int s) { theDipoleSet = s; } /** * Return the produced subtracted matrix elements */ const vector::ptr>& subtractedMEs() const { return theSubtractedMEs; } /** * Access the produced subtracted matrix elements */ vector::ptr>& subtractedMEs() { return theSubtractedMEs; } /** * Return the produced finite real emission matrix elements */ const vector::ptr>& finiteRealMEs() const { return theFiniteRealMEs; } /** * Access the produced finite real emission elements */ vector::ptr>& finiteRealMEs() { return theFiniteRealMEs; } /** * Return the map of Born processes to splitting dipoles */ const map::ptr> >& splittingDipoles() const { return theSplittingDipoles; } /** * Identify a splitting channel */ struct SplittingChannel { /** * The Born XComb */ StdXCombPtr bornXComb; /** * The real XComb */ StdXCombPtr realXComb; /** * The set of tilde XCombs to consider for the real xcomb */ vector tildeXCombs; /** * The dipole in charge of the splitting */ Ptr::ptr dipole; /** * Dump the setup */ void print(ostream&) const; }; /** * Generate all splitting channels for the Born process handled by * the given XComb */ list getSplittingChannels(tStdXCombPtr xc) const; /** * Return the reweight objects for matrix elements */ const vector& reweighters() const { return theReweighters; } /** * Access the reweight objects for matrix elements */ vector& reweighters() { return theReweighters; } /** * Return the preweight objects for matrix elements */ const vector& preweighters() const { return thePreweighters; } /** * Access the preweight objects for matrix elements */ vector& preweighters() { return thePreweighters; } //@} /** @name Setup the matrix elements */ //@{ /** * Return true if this object needs to be initialized before all * other objects (except those for which this function also returns * true). This default version always returns false, but subclasses * may override it to return true. */ virtual bool preInitialize() const { return true; } /** * Prepare a matrix element. */ void prepareME(Ptr::ptr); /** * Check consistency and switch to porduction mode. */ void productionMode(); /** * Setup everything */ virtual void setup(); /** * The highest multiplicity of legs having virtual contributions.(needed for madgraph) */ size_t highestVirt(){return theHighestVirtualsize;} //@} /** @name Diagnostic information */ //@{ /** * Return true, if verbose */ bool verbose() const { return theVerbose; } /** * Switch on diagnostic information. */ void setVerbose(bool on = true) { theVerbose = on; } /** * Return true, if diagram weight is verbose */ bool verboseDia() const { return theDiagramWeightVerbose; } /** * Number of bins for diagram weight verbosity */ int diagramWeightVerboseNBins() const {return theDiagramWeightVerboseNBins;} /** * Return true, if verbose while initializing */ bool initVerbose() const { return theInitVerbose || verbose(); } /** * Switch on diagnostic information while initializing */ void setInitVerbose(bool on = true) { theInitVerbose = on; } /** * Dump the setup */ void print(ostream&) const; /** * Return the subtraction data prefix. */ const string& subtractionData() const { return theSubtractionData; } /** * Set the subtraction data prefix. */ void subtractionData(const string& s) { theSubtractionData = s; } /** * Return the subtraction plot type. */ const int& subtractionPlotType() const { return theSubtractionPlotType; } /** * Set the subtraction plot type. */ void subtractionPlotType(const int& t) { theSubtractionPlotType = t; } /** * Return whether subtraction data should be plotted for all phase space points individually */ const bool& subtractionScatterPlot() const { return theSubtractionScatterPlot; } /** * Set whether subtraction data should be plotted for all phase space points individually */ void subtractionScatterPlot(const bool& s) { theSubtractionScatterPlot = s; } /** * Return the pole data prefix. */ const string& poleData() const { return thePoleData; } /** * Set the pole data prefix. */ void poleData(const string& s) { thePoleData = s; } /** * Return true, if cancellationn of epsilon poles should be checked. */ bool checkPoles() const { return poleData() != ""; } //@} /** @name Process generation */ //@{ /** * Return the particle groups. */ const map& particleGroups() const { return theParticleGroups; } /** * Access the particle groups. */ map& particleGroups() { return theParticleGroups; } /** * Return true, if the given particle is incoming */ bool isIncoming(cPDPtr p) const { return theIncoming.find(p->id()) != theIncoming.end(); } /** * Return true, if spin correlation information should be provided, if possible. */ bool spinCorrelations() const { return theSpinCorrelations; } /** * Indicate that spin correlation information should be provided, if possible. */ void setSpinCorrelations(bool yes) { theSpinCorrelations = yes; } //@} /** @name Truncated qtilde shower information */ //@{ /** * Return the subprocess of the real emission */ tSubProPtr hardTreeSubprocess() { return theHardtreeSubprocess; } /** * Set the subprocess of the real emission for use in calculating the shower hardtree */ void setHardTreeSubprocess(tSubProPtr hardTree) { theHardtreeSubprocess = hardTree; } /** * Return the born emitter */ int hardTreeEmitter() { return theHardtreeEmitter; } /** * Set the born emitter for use in calculating the shower hardtree */ void setHardTreeEmitter(int emitter) { theHardtreeEmitter = emitter; } /** * Return the born spectator */ int hardTreeSpectator() { return theHardtreeSpectator; } /** * Set the born spectator for use in calculating the shower hardtree */ void setHardTreeSpectator(int spectator) { theHardtreeSpectator = spectator; } //@} /** @name Data handling */ //@{ /** * Return (and possibly create) a directory to contain amplitude * information. */ const string& buildStorage(); /** * Return (and possibly create) a directory to contain integration grid * information. */ const string& runStorage(); /** * alpha of http://arxiv.org/pdf/hep-ph/0307268v2.pdf to restrict * dipole phase space */ double alphaParameter() const { return theAlphaParameter; } /** * set the alpha parameter (needed for massive PK-Operator) */ void setAlphaParameter(double a) { theAlphaParameter = a; } //@} public: + /** + * Print a summary of the parameters used + */ + void summary(ostream&) const; + +public: + /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); //@} private: /** * Flag to indicate that at least one MatchboxFactory object is in action */ static bool& theIsMatchboxRun(); /** * The diagram generator. */ Ptr::ptr theDiagramGenerator; /** * The process data object to be used */ Ptr::ptr theProcessData; /** * The number of light flavours, this matrix * element is calculated for. */ unsigned int theNLight; /** * Vector with the PDG ids of the light quark flavours, * which are contained in the jet particle group. */ vector theNLightJetVec; /** * Vector with the PDG ids of the heavy quark flavours, * which are contained in the jet particle group. */ vector theNHeavyJetVec; /** * Vector with the PDG ids of the light quark flavours, * which are contained in the proton particle group. */ vector theNLightProtonVec; /** * The order in \f$\alpha_S\f$. */ unsigned int theOrderInAlphaS; /** * The order in \f$\alpha_{EM}\f$. */ unsigned int theOrderInAlphaEW; /** * Switch on or off Born contributions */ bool theBornContributions; /** * Switch on or off virtual contributions */ bool theVirtualContributions; /** * Switch on or off subtracted real emission contributions should be included. */ bool theRealContributions; /** * True if virtual contributions should be treated as independent subprocesses */ bool theIndependentVirtuals; /** * True if PK operator contributions should be treated as independent subprocesses */ bool theIndependentPKs; /** * True, if SubProcessGroups should be * setup from this MEGroup. If not, a single SubProcess * is constructed from the data provided by the * head matrix element. */ bool theSubProcessGroups; /** * The phase space generator to be used. */ Ptr::ptr thePhasespace; /** * The scale choice object */ Ptr::ptr theScaleChoice; /** * The factorization scale factor. */ double theFactorizationScaleFactor; /** * The renormalization scale factor. */ double theRenormalizationScaleFactor; /** * Use non-running couplings. */ bool theFixedCouplings; /** * Use non-running couplings. */ bool theFixedQEDCouplings; /** * True, if veto scales should be set * for the real emission */ bool theVetoScales; /** * The amplitudes to be considered */ vector::ptr> theAmplitudes; /** * The Born matrix elements to be considered */ vector::ptr> theBornMEs; /** * The loop induced matrix elements to be considered */ vector::ptr> theLoopInducedMEs; /** * The virtual corrections to be considered */ vector::ptr> theVirtuals; /** * The real emission matrix elements to be considered */ vector::ptr> theRealEmissionMEs; /** * The produced NLO matrix elements */ vector::ptr> theBornVirtualMEs; /** * The produced subtracted matrix elements */ vector::ptr> theSubtractedMEs; /** * The produced finite real emission matrix elements */ vector::ptr> theFiniteRealMEs; /** * Which set of dipoles should be considered */ int theDipoleSet; /** * Switch on or off verbosity */ bool theVerbose; /** * Switch on or off diagram weight verbosity */ bool theDiagramWeightVerbose; /** * Number of bins for diagram weight verbosity */ int theDiagramWeightVerboseNBins; /** * True, if verbose while initializing */ bool theInitVerbose; /** * Prefix for subtraction data */ string theSubtractionData; /** * Set the type of plot that is to be generated for subtraction checking */ int theSubtractionPlotType; /** * Set whether subtraction data should be plotted for all phase space points individually */ bool theSubtractionScatterPlot; /** * Prefix for pole data. */ string thePoleData; /** * Command to limit the real emission process to be considered. */ string doSingleRealProcess(string); /** * The real emission process to be included; if empty, all possible * ones will be considered. */ vector > realEmissionProcesses; /** * Particle groups. */ map theParticleGroups; /** * Command to start a particle group. */ string startParticleGroup(string); /** * The name of the particle group currently edited. */ string particleGroupName; /** * The particle group currently edited. */ PDVector particleGroup; /** * Command to end a particle group. */ string endParticleGroup(string); /** * Parse a process description */ vector parseProcess(string); /** * Command to set the process. */ string doProcess(string); /** * Command to set the process. */ string doLoopInducedProcess(string); /** * The process to consider in terms of particle groups. */ vector > processes; /** * The loop induced process to consider in terms of particle groups. */ vector > loopInducedProcesses; /** * Generate subprocesses. */ set makeSubProcesses(const vector&) const; /** * Generate matrix element objects for the given process. */ vector::ptr> makeMEs(const vector&, unsigned int orderas, bool virt); /** * The shower approximation. */ Ptr::ptr theShowerApproximation; /** * The map of Born processes to splitting dipoles */ map::ptr> > theSplittingDipoles; /** * True, if subtraction scales should be caluclated from real emission kinematics */ bool theRealEmissionScales; /** * Consider all processes with order in couplings specifying the * maximum order. */ bool theAllProcesses; /** * The processes to be ordered from an OLP */ map::tptr,map,int> > theOLPProcesses; /** * Amplitudes which need external initialization */ set::tptr> theExternalAmplitudes; /** * Amplitudes to be selected on clashing responsibilities. */ vector::ptr> theSelectedAmplitudes; /** * Amplitudes to be deselected on clashing responsibilities. */ vector::ptr> theDeselectedAmplitudes; /** * Reweight objects for matrix elements */ vector theReweighters; /** * Preweight objects for matrix elements */ vector thePreweighters; /** * Produce matrix element corrections, but no NLO */ bool theMECorrectionsOnly; /** * The highest multiplicity of legs having virtual contributions.(needed for madgraph) */ int theHighestVirtualsize; /** * Produce matrix element corrections, with LoopSim NLO */ bool theLoopSimCorrections; /** * True, if the setup has already been run. */ bool ranSetup; /** * PDG ids of incoming particles */ set theIncoming; /** * True, if first incoming partons originate from perturbative PDF */ bool theFirstPerturbativePDF; /** * True, if second incoming partons originate from perturbative PDF */ bool theSecondPerturbativePDF; /** * True, if this Factory is in production mode. */ bool inProductionMode; /** * The real emission subprocess used when calculating the hardtree * in the truncated qtilde shower */ tSubProPtr theHardtreeSubprocess; /** * The born emitter used when calculating the hardtree in * the truncated shower */ int theHardtreeEmitter; /** * The born spectator used when calculating the hardtree in * the truncated shower */ int theHardtreeSpectator; /** * True, if spin correlation information should be provided, if possible. */ bool theSpinCorrelations; /** * The alpha parameter to be used for the dipole subtraction */ double theAlphaParameter; /** * Wether or not charge conservation should be enforced for the processes * constructed. */ bool theEnforceChargeConservation; /** * Wether or not colour conservation should be enforced for the processes * constructed. */ bool theEnforceColourConservation; /** * Wether or not lepton number conservation should be enforced for the processes * constructed. */ bool theEnforceLeptonNumberConservation; /** * Wether or not quark number conservation should be enforced for the processes * constructed. */ bool theEnforceQuarkNumberConservation; /** * Assume flavour diagonal lepton interactions */ bool theLeptonFlavourDiagonal; /** * Assume flavour diagonal quark interactions */ bool theQuarkFlavourDiagonal; /** * Command for production mode */ string doProductionMode(string) { productionMode(); return ""; } private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MatchboxFactory & operator=(const MatchboxFactory &); }; } #endif /* HERWIG_MatchboxFactory_H */ diff --git a/MatrixElement/Matchbox/Scales/MatchboxHtScale.cc b/MatrixElement/Matchbox/Scales/MatchboxHtScale.cc --- a/MatrixElement/Matchbox/Scales/MatchboxHtScale.cc +++ b/MatrixElement/Matchbox/Scales/MatchboxHtScale.cc @@ -1,164 +1,164 @@ // -*- C++ -*- // // MatchboxHtScale.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 MatchboxHtScale class. // #include "MatchboxHtScale.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; MatchboxHtScale::MatchboxHtScale() : theIncludeMT(false), theHTFactor(1.0), theMTFactor(1.0),theScalePtCut(15.*GeV) {} MatchboxHtScale::~MatchboxHtScale() {} IBPtr MatchboxHtScale::clone() const { return new_ptr(*this); } IBPtr MatchboxHtScale::fullclone() const { return new_ptr(*this); } Energy2 MatchboxHtScale::renormalizationScale() const { tcPDVector pd (mePartonData().begin() + 2, mePartonData().end()); vector p (meMomenta().begin() + 2, meMomenta().end()); tcPDPtr t1 = mePartonData()[0]; tcPDPtr t2 = mePartonData()[1]; tcCutsPtr cuts = lastCutsPtr(); theJetFinder->cluster(pd, p, cuts, t1, t2); initWeightFactors(pd,p,theJetFinder); // momentum of the non-jet system LorentzMomentum nonJetMomentum(ZERO,ZERO,ZERO,ZERO); // (weighted) pt of the jet systems Energy ptJetSum = ZERO; bool gotone = false; tcPDVector::const_iterator pdata = pd.begin(); vector::const_iterator mom = p.begin(); for ( ; mom != p.end(); ++pdata, ++mom ) { if ( theJetFinder->unresolvedMatcher()->check(**pdata)&& mom->perp()>theScalePtCut){ //abs(mom->rapidity()+(!lastXCombPtr()->head()?lastXCombPtr()->lastY():lastXCombPtr()->head()->lastY()))<5.01 gotone = true; ptJetSum += jetPtWeight(*mom)*mom->perp(); } else if ( theIncludeMT ) { nonJetMomentum += *mom; } } if ( !gotone && lastXCombPtr()->willPassCuts() ) throw Exception() << "MatchboxHtScale::renormalizationScale(): " << "No jets could be found. Check your setup." << "\nHint: The HT scale is defined with a PtMin cut on jets. (default:) " << "\n set /Herwig/MatrixElements/Matchbox/ScalesHTScale:JetPtCut 15.*GeV " << Exception::runerror; Energy mtNonJetSum = sqrt(nonJetMomentum.perp2() + nonJetMomentum.m2()); mtNonJetSum *= theMTFactor; ptJetSum *= theHTFactor; return sqr(ptJetSum + mtNonJetSum); } Energy2 MatchboxHtScale::factorizationScale() const { return renormalizationScale(); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void MatchboxHtScale::persistentOutput(PersistentOStream & os) const { os << theJetFinder << theIncludeMT << theHTFactor << theMTFactor << ounit(theScalePtCut,GeV); } void MatchboxHtScale::persistentInput(PersistentIStream & is, int) { is >> theJetFinder >> theIncludeMT >> theHTFactor >> theMTFactor >> iunit(theScalePtCut,GeV); } // *** 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 describeHerwigMatchboxHtScale("Herwig::MatchboxHtScale", "HwMatchboxScales.so"); void MatchboxHtScale::Init() { static ClassDocumentation documentation ("MatchboxHtScale implements scale choices related to transverse momenta."); static Reference interfaceJetFinder ("JetFinder", "A reference to the jet finder.", &MatchboxHtScale::theJetFinder, false, false, true, false, false); static Switch interfaceIncludeMT ("IncludeMT", "Include the transverse masses of the non-jet objects.", &MatchboxHtScale::theIncludeMT, false, false, false); static SwitchOption interfaceIncludeMTYes (interfaceIncludeMT, "Yes", "", true); static SwitchOption interfaceIncludeMTNo (interfaceIncludeMT, "No", "", false); static Parameter interfaceHTFactor ("HTFactor", "A factor to scale the HT contribution.", &MatchboxHtScale::theHTFactor, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceMTFactor ("MTFactor", "A factor to scale the MT contribution.", &MatchboxHtScale::theMTFactor, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceScalePtCut ("JetPtCut", "The Pt cut to define jets in the sum.", - &MatchboxHtScale::theScalePtCut, 15.*GeV, 0.*GeV, 0.*GeV, + &MatchboxHtScale::theScalePtCut, GeV, 15.*GeV, 0.*GeV, 0.*GeV, false, false, Interface::lowerlim); } diff --git a/MatrixElement/Matchbox/Utility/AmplitudeCache.h b/MatrixElement/Matchbox/Utility/AmplitudeCache.h --- a/MatrixElement/Matchbox/Utility/AmplitudeCache.h +++ b/MatrixElement/Matchbox/Utility/AmplitudeCache.h @@ -1,346 +1,346 @@ // -*- C++ -*- // // AmplitudeCache.h 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. // #ifndef HERWIG_AmplitudeCache_H #define HERWIG_AmplitudeCache_H #include "Herwig/MatrixElement/Matchbox/Utility/SpinorHelicity.h" #include "ThePEG/Config/algorithm.h" +#include namespace Herwig { using namespace ThePEG; +using std::array; namespace SpinorHelicity { /** * \ingroup Matchbox * \author Simon Platzer * * \brief Caching for amplitudes using spinor helicity techniques. * */ -template +template class AmplitudeCache { typedef map > AmplitudeCacheMap; typedef map > > CurrentCacheMap; /** + * Maximum N we can handle, SYM_N is storage size for a symmetric matrix of N * N elements + */ + enum { MAX_N = 7, SYM_N = MAX_N*(MAX_N+1)/2 }; + + /** * The number of points */ int theNPoints; /** * The energy scale to obtain dimensionless * quantities. */ mutable Energy theScale; /** * Masses indexed by partons */ - mutable vector theMasses; + mutable array theMasses; /** * Momenta indexed by partons */ - mutable vector theMomenta; + mutable array theMomenta; /** * Crossing signs indexed by partons */ - mutable vector theCrossingSigns; + mutable array theCrossingSigns; /** * Plus spinors indexed by partons */ - mutable vector thePlusSpinors; + mutable array thePlusSpinors; /** * Plus conjugate spinors indexed by partons */ - mutable vector thePlusConjugateSpinors; + mutable array thePlusConjugateSpinors; /** * Invariants indexed by partons */ - mutable vector > theInvariants; + mutable array theInvariants; /** * Flag products to be recalculated */ - mutable vector > getInvariant; + mutable array getInvariant; /** * Spinor products indexed by partons */ - mutable vector > thePlusProducts; + mutable array thePlusProducts; /** * Flag products to be recalculated */ - mutable vector > getPlusProduct; + mutable array getPlusProduct; /** * Spinor currents indexed by partons */ - mutable vector > > thePlusCurrents; + mutable array,SYM_N> thePlusCurrents; /** * Flag currents to be recalculated */ - mutable vector > getPlusCurrent; + mutable array getPlusCurrent; /** * Cache intermediate amplitudes by index */ mutable AmplitudeCacheMap theCachedAmplitudes; /** * The last query for a cached amplitude */ mutable typename AmplitudeCacheMap::iterator theLastAmplitude; /** * Cache intermediate currents by index */ mutable CurrentCacheMap theCachedCurrents; /** * The last query for a cached current */ mutable typename CurrentCacheMap::iterator theLastCurrent; /** + * Helper function to index symmetric arrays, assumes i <= j. + * Usual indexing function (N*i + j) corrected by triangle number for i-th row. + */ + inline size_t idx(size_t i, size_t j) const { + return MAX_N * i - (i + 1) * i / 2 + j; + } + + /** * Helper to reset flags */ struct boolResetter { - void operator()(vector::reference flag) const { - flag = true; - } void operator()(pair >& flag) const { flag.second.first = true; } void operator()(pair > >& flag) const { flag.second.first = true; } }; - /** - * Helper to reset flags - */ - struct boolVectorResetter { - void operator()(vector& flags) const { - for_each(flags.begin(),flags.end(),boolResetter()); - } - }; - public: /** * Constructor */ - AmplitudeCache() - : theNPoints(0) {} + AmplitudeCache() : theNPoints(0) {} /** * Prepare for n-point amplitude */ void nPoints(int n); /** * Return the number of points */ int nPoints() const { return theNPoints; } /** * Set the energy scale to obtain dimensionless * quantities and flag all quantities to be recalculated. */ void amplitudeScale(Energy s) const; /** * Set the momentum for the k'th parton * and its associated mass. */ void momentum(int k, const LorentzMomentum& p, bool getSpinors = true, Energy mass = ZERO) const; /** * Reset flags */ void reset() const; public: /** * Return the momentum for the k'th parton */ LorentzVector momentum(int k) const { return theMomenta[k]/theScale; } /** * Get the energy scale to obtain dimensionless * quantities and flag all quantities to be recalculated. */ Energy amplitudeScale() const { return theScale; } /** * Return the mass associated to the k'th parton */ double mass(int k) const { return theMasses[k]; } /** * Return the crossing sign for the * i'th parton */ int crossingSign(int i) const { return theCrossingSigns[i]; } /** * Return the crossing sign for the * i'th and j'th parton */ double crossingSign(int i, int j) const { return theCrossingSigns[i]*theCrossingSigns[j]; } /** * Return (ij) */ double invariant(int i, int j) const { - if ( i== j ) - return 0.; - if ( i > j ) - swap(i,j); - if ( getInvariant[i][j] ) { - getInvariant[i][j] = false; - theInvariants[i][j] = 2.*(momentum(i)*momentum(j)); + if ( i == j ) return 0.; + if ( i > j ) swap(i,j); + if ( getInvariant[idx(i,j)] ) { + getInvariant[idx(i,j)] = false; + theInvariants[idx(i,j)] = 2.*(momentum(i)*momentum(j)); } - return theInvariants[i][j]; + return theInvariants[idx(i,j)]; } /** * Return */ Complex plusProduct(int i, int j) const { if ( i== j ) return 0.; bool swapij = (i > j); if ( swapij ) swap(i,j); - if ( getPlusProduct[i][j] ) { - getPlusProduct[i][j] = false; - thePlusProducts[i][j] = + if ( getPlusProduct[idx(i,j)] ) { + getPlusProduct[idx(i,j)] = false; + thePlusProducts[idx(i,j)] = PlusSpinorProduct(thePlusConjugateSpinors[i], thePlusSpinors[j]).eval() / theScale; } - return swapij ? -thePlusProducts[i][j] : thePlusProducts[i][j]; + return swapij ? -thePlusProducts[idx(i,j)] : thePlusProducts[idx(i,j)]; } /** * Return [ij] */ Complex minusProduct(int i, int j) const { if ( i== j ) return 0.; return -crossingSign(i,j)*conj(plusProduct(i,j)); } /** * Return plusCurrent(int i, int j) const { bool swapij = (i > j); if ( swapij ) swap(i,j); - if ( getPlusCurrent[i][j] ) { - getPlusCurrent[i][j] = false; + if ( getPlusCurrent[idx(i,j)] ) { + getPlusCurrent[idx(i,j)] = false; if ( i != j ) { - thePlusCurrents[i][j] = + thePlusCurrents[idx(i,j)] = PlusSpinorCurrent(thePlusConjugateSpinors[i], MinusSpinor(theMomenta[j])).eval() / theScale; } else { - thePlusCurrents[i][j] = 2.*momentum(i); + thePlusCurrents[idx(i,j)] = 2.*momentum(i); } } - return swapij ? crossingSign(i,j)*thePlusCurrents[i][j].conjugate() : thePlusCurrents[i][j]; + return swapij ? crossingSign(i,j)*thePlusCurrents[idx(i,j)].conjugate() : thePlusCurrents[idx(i,j)]; } /** * Return [i|\gamma^\mu|j> */ LorentzVector minusCurrent(int i, int j) const { return plusCurrent(j,i); } public: /** * Return true, if the given amplitude * needs to be recalculated. */ bool getAmplitude(const AmplitudeKey& key) const { static Complex czero; if ( ( theLastAmplitude = theCachedAmplitudes.find(key) ) == theCachedAmplitudes.end() ) { theLastAmplitude = theCachedAmplitudes.insert(make_pair(key,make_pair(true,czero))).first; } return theLastAmplitude->second.first; } /** * Cache an amplitude */ void cacheAmplitude(Complex amp) const { theLastAmplitude->second = make_pair(false,amp); } /** * Return a cached amplitude */ const Complex& cachedAmplitude() const { return theLastAmplitude->second.second; } /** * Return true, if the given current * needs to be recalculated. */ bool getCurrent(const AmplitudeKey& key) const { static LorentzVector czero; if ( ( theLastCurrent = theCachedCurrents.find(key) ) == theCachedCurrents.end() ) { theLastCurrent = theCachedCurrents.insert(make_pair(key,make_pair(true,czero))).first; } return theLastCurrent->second.first; } /** * Cache an current */ void cacheCurrent(const LorentzVector& curr) const { theLastCurrent->second = make_pair(false,curr); } /** * Return a cached current */ const LorentzVector& cachedCurrent() const { return theLastCurrent->second.second; } }; } } #include "AmplitudeCache.tcc" #endif // HERWIG_AmplitudeCache_H diff --git a/MatrixElement/Matchbox/Utility/AmplitudeCache.tcc b/MatrixElement/Matchbox/Utility/AmplitudeCache.tcc --- a/MatrixElement/Matchbox/Utility/AmplitudeCache.tcc +++ b/MatrixElement/Matchbox/Utility/AmplitudeCache.tcc @@ -1,73 +1,60 @@ // -*- C++ -*- // // AmplitudeCache.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. // namespace Herwig { namespace SpinorHelicity { -template +template void AmplitudeCache::nPoints(int n) { + assert( n <= MAX_N ); + theNPoints = n; - theMasses.clear(); - theMomenta.clear(); - theCrossingSigns.clear(); - thePlusSpinors.clear(); - thePlusConjugateSpinors.clear(); - theInvariants.clear(); - thePlusProducts.clear(); - thePlusCurrents.clear(); - getInvariant.clear(); - getPlusProduct.clear(); - getPlusCurrent.clear(); - - theMasses.resize(n); - theMomenta.resize(n); - theCrossingSigns.resize(n); - thePlusSpinors.resize(n); - thePlusConjugateSpinors.resize(n); - theInvariants.resize(n,vector(n)); - thePlusProducts.resize(n,vector(n)); - thePlusCurrents.resize(n,vector >(n)); - getInvariant.resize(n,vector(n)); - getPlusProduct.resize(n,vector(n)); - getPlusCurrent.resize(n,vector(n)); + theMasses.fill({}); + theMomenta.fill({}); + theCrossingSigns.fill({}); + thePlusSpinors.fill(PlusSpinor()); + thePlusConjugateSpinors.fill(PlusConjugateSpinor()); + theInvariants.fill({}); + thePlusProducts.fill({}); + thePlusCurrents.fill({}); reset(); } -template +template void AmplitudeCache::amplitudeScale(Energy s) const { theScale = s; reset(); } -template +template void AmplitudeCache::momentum(int k, const LorentzMomentum& p, bool getSpinors, Energy mass) const { theMasses[k] = mass/theScale; theMomenta[k] = p; if ( getSpinors ) { theCrossingSigns[k] = p.t() > ZERO ? 1 : -1; thePlusSpinors[k] = PlusSpinor(p); thePlusConjugateSpinors[k] = PlusConjugateSpinor(p); } } -template +template void AmplitudeCache::reset() const { - for_each(getInvariant.begin(),getInvariant.end(),boolVectorResetter()); - for_each(getPlusProduct.begin(),getPlusProduct.end(),boolVectorResetter()); - for_each(getPlusCurrent.begin(),getPlusCurrent.end(),boolVectorResetter()); + getInvariant.fill(true); + getPlusProduct.fill(true); + getPlusCurrent.fill(true); for_each(theCachedAmplitudes.begin(),theCachedAmplitudes.end(),boolResetter()); for_each(theCachedCurrents.begin(),theCachedCurrents.end(),boolResetter()); } }} diff --git a/Sampling/BinSampler.cc b/Sampling/BinSampler.cc --- a/Sampling/BinSampler.cc +++ b/Sampling/BinSampler.cc @@ -1,717 +1,728 @@ // -*- C++ -*- // // BinSampler.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 BinSampler class. // #include "BinSampler.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 "ThePEG/Repository/Repository.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Handlers/StandardEventHandler.h" #include "ThePEG/Handlers/StandardXComb.h" #include #include "GeneralSampler.h" using namespace Herwig; BinSampler::BinSampler() : MultiIterationStatistics(), theBias(1.), theWeighted(false), theInitialPoints(1000000), theNIterations(1), theEnhancementFactor(1.0), theNonZeroInPresampling(false), theHalfPoints(false), theMaxNewMax(30), theReferenceWeight(1.0), theBin(-1), theInitialized(false), theRemapperPoints(0), theRemapChannelDimension(false), theLuminosityMapperBins(0), theGeneralMapperBins(0), theRemapperMinSelection(0.00001), theIntegrated(false), theRemappersFilled(false), - theHasGrids(false) {} + theHasGrids(false), + theKappa(1.){} BinSampler::~BinSampler() {} IBPtr BinSampler::clone() const { return new_ptr(*this); } IBPtr BinSampler::fullclone() const { return new_ptr(*this); } void BinSampler::sampler(Ptr::tptr s) { theSampler = s; } Ptr::tptr BinSampler::sampler() const { return theSampler; } string BinSampler::process() const { ostringstream os(""); const StandardEventHandler& eh = *theEventHandler; const StandardXComb& xc = *eh.xCombs()[theBin]; os << xc.matrixElement()->name() << " : "; os << xc.mePartonData()[0]->PDGName() << " " << xc.mePartonData()[1]->PDGName() << " -> "; for ( cPDVector::const_iterator pid = xc.mePartonData().begin() + 2; pid != xc.mePartonData().end(); ++pid ) os << (**pid).PDGName() << " "; return os.str(); } string BinSampler::shortprocess() const { ostringstream os(""); const StandardEventHandler& eh = *theEventHandler; const StandardXComb& xc = *eh.xCombs()[theBin]; os << xc.mePartonData()[0]->id() << " " << xc.mePartonData()[1]->id() << " : "; for ( cPDVector::const_iterator pid = xc.mePartonData().begin() + 2; pid != xc.mePartonData().end(); ++pid ) os << (**pid).id() << " "; return os.str(); } string BinSampler::id() const { ostringstream os(""); const StandardEventHandler& eh = *theEventHandler; const StandardXComb& xc = *eh.xCombs()[theBin]; string name = xc.matrixElement()->name(); string::size_type i = name.find_first_of("["); string nameFirst = name.substr(0,i); i = name.find_first_of("]"); string nameSecond = name.substr(i+1); os << nameFirst << nameSecond << ":"; for ( cPDVector::const_iterator pid = xc.mePartonData().begin(); pid != xc.mePartonData().end(); ++pid ) os << (**pid).id() << (pid != (--xc.mePartonData().end()) ? "," : ""); return os.str(); } double BinSampler::evaluate(vector p, bool remap) { double w = 1.0; if ( remap && !remappers.empty() ) { for ( size_t k = 0; k < p.size(); ++k ) { map::const_iterator r = remappers.find(k); if ( r != remappers.end() ) { pair f = r->second.generate(p[k]); p[k] = f.first; w /= f.second; } } } try { w *= eventHandler()->dSigDR(p) / nanobarn; } catch (Veto&) { w = 0.0; } catch (...) { throw; } if (randomNumberString()!="") for ( size_t k = 0; k < p.size(); ++k ) { RandomNumberHistograms[RandomNumberIndex(id(),k)].first.book(p[k],w); RandomNumberHistograms[RandomNumberIndex(id(),k)].second+=w; } return w; } double BinSampler::generate() { double w = 1.; for ( size_t k = 0; k < lastPoint().size(); ++k ) { lastPoint()[k] = UseRandom::rnd(); } try { w = evaluate(lastPoint()); } catch (Veto&) { w = 0.0; } catch (...) { throw; } if ( !weighted() && initialized() ) { - double p = min(abs(w),referenceWeight())/referenceWeight(); + double p = min(abs(w),kappa()*referenceWeight())/(kappa()*referenceWeight()); double sign = w >= 0. ? 1. : -1.; if ( p < 1 && UseRandom::rnd() > p ) w = 0.; else - w = sign*max(abs(w),referenceWeight()); + w = sign*max(abs(w),referenceWeight()*kappa()); } select(w); if ( w != 0.0 ) accept(); + assert(kappa()==1.||sampler()->almostUnweighted()); return w; } void BinSampler::fillRemappers(bool progress) { if ( remappers.empty() ) return; unsigned long nanPoints = 0; boost::progress_display* progressBar = 0; if ( progress ) { Repository::clog() << "warming up " << process(); progressBar = new boost::progress_display(theRemapperPoints,Repository::clog()); } unsigned long countzero =0; for ( unsigned long k = 0; k < theRemapperPoints; ++k,++countzero ) { if (countzero>=theRemapperPoints)break; double w = 1.; for ( size_t j = 0; j < lastPoint().size(); ++j ) { lastPoint()[j] = UseRandom::rnd(); } try { w = evaluate(lastPoint(),false); } catch (Veto&) { w = 0.0; } catch (...) { throw; } if ( isnan(w) || isinf(w) ) ++nanPoints; if ( theNonZeroInPresampling && w==0. ){ k--; continue; } if ( w != 0.0 ) { countzero=0; for ( map::iterator r = remappers.begin(); r != remappers.end(); ++r ) r->second.fill(lastPoint()[r->first],w); } if ( progressBar ) ++(*progressBar); } if ( progressBar ) { delete progressBar; } if ( nanPoints ) { Repository::clog() << "Warning: " << nanPoints << " out of " << theRemapperPoints << " points with nan or inf " << "weight encountered while filling remappers.\n" << flush; } } void BinSampler::saveIntegrationData() const { XML::Element stats = MultiIterationStatistics::toXML(); stats.appendAttribute("process",id()); sampler()->grids().append(stats); } void BinSampler::readIntegrationData() { if ( theIntegrated ) return; bool haveStats = false; list::iterator sit = sampler()->grids().children().begin(); for ( ; sit != sampler()->grids().children().end(); ++sit ) { if ( sit->type() != XML::ElementTypes::Element ) continue; if ( sit->name() != "MultiIterationStatistics" ) continue; string proc; sit->getFromAttribute("process",proc); if ( proc == id() ) { haveStats = true; break; } } if ( haveStats ) { MultiIterationStatistics::fromXML(*sit); sampler()->grids().erase(sit); theIntegrated = true; } else { throw Exception() << "\n--------------------------------------------------------------------------------\n\n" << "Expected integration data.\n\n" << "* When using the build setup make sure the integrate command has been run.\n\n" << "* Check the [EventGenerator].log file for further information.\n\n" << "* Make sure that the Herwig folder can be found and that it contains a HerwigGrids.xml file.\n\n" << "* If you have split the integration jobs, make sure that each integration job was finished.\n" << " Afterwards delete the global HerwigGrids.xml file in the Herwig subfolder\n" << " to automatically create an updated version of the global HerwigGrids.xml file.\n\n" << "--------------------------------------------------------------------------------\n" << Exception::abortnow; } } void BinSampler::saveRemappers() const { if ( remappers.empty() ) return; XML::Element maps(XML::ElementTypes::Element,"Remappers"); maps.appendAttribute("process",id()); for ( map::const_iterator r = remappers.begin(); r != remappers.end(); ++r ) { XML::Element rmap = r->second.toXML(); rmap.appendAttribute("dimension",r->first); maps.append(rmap); } sampler()->grids().append(maps); } void BinSampler::setupRemappers(bool progress) { if ( !theRemapperPoints ) return; if ( theRemappersFilled ) return; lastPoint().resize(dimension()); bool haveGrid = false; list::iterator git = sampler()->grids().children().begin(); for ( ; git != sampler()->grids().children().end(); ++git ) { if ( git->type() != XML::ElementTypes::Element ) continue; if ( git->name() != "Remappers" ) continue; string proc; git->getFromAttribute("process",proc); if ( proc == id() ) { haveGrid = true; break; } } if ( haveGrid ) { for ( list::iterator cit = git->children().begin(); cit != git->children().end(); ++cit ) { if ( cit->type() != XML::ElementTypes::Element ) continue; if ( cit->name() != "Remapper" ) continue; size_t dimension = 0; cit->getFromAttribute("dimension",dimension); remappers[dimension].fromXML(*cit); } sampler()->grids().erase(git); } if ( !haveGrid ) { const StandardEventHandler& eh = *eventHandler(); const StandardXComb& xc = *eh.xCombs()[bin()]; const pair& pdims = xc.partonDimensions(); set remapped; if ( theRemapChannelDimension && xc.diagrams().size() > 1 && dimension() > pdims.first + pdims.second ) { remappers[pdims.first] = Remapper(xc.diagrams().size(),theRemapperMinSelection,false); remapped.insert(pdims.first); } if ( theLuminosityMapperBins > 1 && dimension() >= pdims.first + pdims.second ) { for ( int n = 0; n < pdims.first; ++n ) { remappers[n] = Remapper(theLuminosityMapperBins,theRemapperMinSelection,true); remapped.insert(n); } for ( int n = dimension() - pdims.second; n < dimension(); ++n ) { remappers[n] = Remapper(theLuminosityMapperBins,theRemapperMinSelection,true); remapped.insert(n); } } if ( theGeneralMapperBins > 1 ) { for ( int n = 0; n < dimension(); n++ ) { if ( remapped.find(n) == remapped.end() ) { remappers[n] = Remapper(theGeneralMapperBins,theRemapperMinSelection,true); remapped.insert(n); } } } fillRemappers(progress); for ( map::iterator r = remappers.begin(); r != remappers.end(); ++r ) { r->second.finalize(); } } theRemappersFilled = true; } void BinSampler::runIteration(unsigned long points, bool progress) { boost::progress_display* progressBar = 0; if ( progress ) { Repository::clog() << "integrating " << process() << " , iteration " << (iterations().size() + 1); progressBar = new boost::progress_display(points,Repository::clog()); } double w=0.; double maxweight=0; int numlastmax=0; unsigned long countzero =0; int newmax=0; for ( unsigned long k = 0; k < points; ++k,++countzero ) { if (countzero>=points)break; w=abs(generate()); if(theNonZeroInPresampling && w==0.0){ k--; continue; } if (w!=0.0) countzero =0; numlastmax++; if (theHalfPoints&&maxweighttheMaxNewMax){ throw Exception() << "\n--------------------------------------------------------------------------------\n\n" << "To many new Maxima.\n\n" << "* With the option:\n\n" << "* set Sampler:BinSampler:HalfPoints Yes\n\n" << "* for every new maximum weight found until the half of the persampling points\n" << "* the counter is set to zero. We count the number of new maxima.\n" << "* You have reached: "<grids().children().empty() ) { nIterations(1); } if ( !integrated() ) { unsigned long points = initialPoints(); for ( unsigned long k = 0; k < nIterations(); ++k ) { runIteration(points,progress); if ( k < nIterations() - 1 ) { points = (unsigned long)(points*enhancementFactor()); adapt(); nextIteration(); } } } isInitialized(); } void BinSampler::finalize(bool){ if (theRandomNumbers!="") for ( map >:: const_iterator b = RandomNumberHistograms.begin(); b != RandomNumberHistograms.end(); ++b ) { b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second); } } BinSampler::RandomNumberHistogram:: RandomNumberHistogram(double low, double up, unsigned int nbins) : lower(low) { nbins = nbins + 1; double c = up / (nbins-1.); for ( unsigned int k = 1; k < nbins; ++k ) { bins[low+c*k] = 0.; binsw1[low+c*k] = 0.; } } void BinSampler::RandomNumberHistogram:: dump(const std::string& folder,const std::string& prefix, const std::string& process, const int NR) const { ostringstream fname(""); std::string prefix2; std::string prefix3=prefix; std::remove_copy(prefix.begin(), prefix.end(), std::back_inserter(prefix2), '.'); prefix3=prefix2;prefix2.clear(); std::remove_copy(prefix3.begin(), prefix3.end(), std::back_inserter(prefix2), ':'); prefix3=prefix2;prefix2.clear(); std::remove_copy(prefix3.begin(), prefix3.end(), std::back_inserter(prefix2), ','); fname << "RN-"<< NR ; ofstream out((folder+"/"+prefix2+fname.str()+".dat").c_str()); double sumofweights=0.; for ( map::const_iterator b = bins.begin();b != bins.end(); ++b ) sumofweights+=b->second; double sumofweights2=0.; for ( map::const_iterator b = binsw1.begin();b != binsw1.end(); ++b ) sumofweights2+=b->second; map::const_iterator b2 = binsw1.begin(); if ( sumofweights == 0 ) { cerr << "Not enough statistic accumulated for " << process << " skipping random number diagnostic.\n" << flush; return; } for ( map::const_iterator b = bins.begin(); b != bins.end(); ++b, ++b2) { out << " " << b->first << " " << b->second/sumofweights*100. << " " << b2->second/sumofweights2*100. << "\n" << flush; } double xmin = -0.01; double xmax = 1.01; ofstream gpout((folder+"/"+prefix2+fname.str()+".gp").c_str()); gpout << "set terminal epslatex color solid\n" << "set output '" << prefix2+fname.str() << "-plot.tex'\n" << "set xrange [" << xmin << ":" << xmax << "]\n"; gpout << "set xlabel 'rn "<> theBias >> theWeighted >> theInitialPoints >> theNIterations >> theEnhancementFactor >> theNonZeroInPresampling >> theHalfPoints >> theMaxNewMax >> theReferenceWeight >> theBin >> theInitialized >> theLastPoint >> theEventHandler >> theSampler >> theRandomNumbers >> theRemapperPoints >> theRemapChannelDimension - >> theLuminosityMapperBins >> theGeneralMapperBins; + >> theLuminosityMapperBins >> theGeneralMapperBins >> theKappa; } // *** 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 describeHerwigBinSampler("Herwig::BinSampler", "HwSampling.so"); void BinSampler::Init() { static ClassDocumentation documentation ("BinSampler samples XCombs bins. This default implementation performs flat MC integration."); static Parameter interfaceInitialPoints ("InitialPoints", "The number of points to use for initial integration.", &BinSampler::theInitialPoints, 1000000, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceNIterations ("NIterations", "The number of iterations to perform initially.", &BinSampler::theNIterations, 1, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceEnhancementFactor ("EnhancementFactor", "The enhancement factor for the number of points in the next iteration.", &BinSampler::theEnhancementFactor, 2.0, 1.0, 0, false, false, Interface::lowerlim); static Switch interfaceNonZeroInPresampling ("NonZeroInPresampling", "Switch on to count only non zero weights in presampling.", &BinSampler::theNonZeroInPresampling, true, false, false); static SwitchOption interfaceNonZeroInPresamplingYes (interfaceNonZeroInPresampling, "Yes", "", true); static SwitchOption interfaceNonZeroInPresamplingNo (interfaceNonZeroInPresampling, "No", "", false); static Switch interfaceHalfPoints ("HalfPoints", "Switch on to reset the counter of points if new maximumis was found in the first 1/2 points.", &BinSampler::theHalfPoints, true, false, false); static SwitchOption interfaceHalfPointsYes (interfaceHalfPoints, "Yes", "", true); static SwitchOption interfaceHalfPointsNo (interfaceHalfPoints, "No", "", false); static Parameter interfaceMaxNewMax ("MaxNewMax", "The maximum number of allowed new maxima in combination with the HalfPoints option.", &BinSampler::theMaxNewMax, 30, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceRandomNumbers ("RandomNumbers", "Prefix for distributions of the random numbers.", &BinSampler::theRandomNumbers, "", false, false); static Parameter interfaceRemapperPoints ("RemapperPoints", "The number of points to be used for filling remappers.", &BinSampler::theRemapperPoints, 10000, 0, 0, false, false, Interface::lowerlim); static Switch interfaceRemapChannelDimension ("RemapChannelDimension", "Switch on remapping of the channel dimension.", &BinSampler::theRemapChannelDimension, true, false, false); static SwitchOption interfaceRemapChannelDimensionYes (interfaceRemapChannelDimension, "Yes", "", true); static SwitchOption interfaceRemapChannelDimensionNo (interfaceRemapChannelDimension, "No", "", false); static Parameter interfaceLuminosityMapperBins ("LuminosityMapperBins", "The number of bins to be used for remapping parton luminosities.", &BinSampler::theLuminosityMapperBins, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfaceGeneralMapperBins ("GeneralMapperBins", "The number of bins to be used for remapping other phase space dimensions.", &BinSampler::theGeneralMapperBins, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfaceRemapperMinSelection ("RemapperMinSelection", "The minimum bin selection probability for remappers.", &BinSampler::theRemapperMinSelection, 0.00001, 0.0, 1.0, false, false, Interface::limited); + + static Parameter interfaceKappa + ("Kappa", + "In AllmostUnweighted mode unweight to Kappa ReferenceWeight.", + &BinSampler::theKappa, 1., 0.000001, 1.0, + false, false, Interface::limited); + + + } diff --git a/Sampling/BinSampler.h b/Sampling/BinSampler.h --- a/Sampling/BinSampler.h +++ b/Sampling/BinSampler.h @@ -1,567 +1,587 @@ // -*- C++ -*- // // BinSampler.h 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. // #ifndef Herwig_BinSampler_H #define Herwig_BinSampler_H // // This is the declaration of the BinSampler class. // #include "ThePEG/Handlers/StandardEventHandler.h" #include "ThePEG/Utilities/Exception.h" #include "ThePEG/Repository/UseRandom.h" #include "MultiIterationStatistics.h" #include "Remapper.h" namespace Herwig { using namespace ThePEG; class GeneralSampler; /** * \ingroup Matchbox * \author Simon Platzer * * \brief BinSampler samples XCombs bins. This default implementation * performs flat MC integration. * * @see \ref BinSamplerInterfaces "The interfaces" * defined for BinSampler. */ class BinSampler: public Herwig::MultiIterationStatistics { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ BinSampler(); /** * The destructor. */ virtual ~BinSampler(); //@} public: /** * Clone this object. */ Ptr::ptr cloneMe() const { return dynamic_ptr_cast::ptr>(clone()); } public: /** * Evaluate the cross section */ double evaluate(vector p, bool remap = true); /** * Return the bias with which this sampler is selected. The sampler * needs to divide out this bias in its weight calculation. */ double bias() const { return theBias; } /** * Set the bias with which this sampler is selected. */ void bias(double b) { theBias = b; } /** * Set the event handler */ void eventHandler(tStdEHPtr eh) { theEventHandler = eh; } /** * Return the event handler */ tStdEHPtr eventHandler() const { return theEventHandler; } /** * Set the containing sampler */ void sampler(Ptr::tptr); /** * Get the containing sampler */ Ptr::tptr sampler() const; /** * Return the bin */ int bin() const { return theBin; } /** * Set the bin */ void bin(int b) { theBin = b; } /** * Return a string describing the process handled by this sampler. */ string process() const; /** * Return a short string describing the process handled by this sampler. */ string shortprocess() const; /** * Return a string identifying the process handled by this sampler. */ string id() const; /** * Return the last generated point. */ const vector& lastPoint() const { return theLastPoint; } /** * Access the last generated point. */ vector& lastPoint() { return theLastPoint; } /** * Return the reference weight to be used */ double referenceWeight() const { return theReferenceWeight; } /** * Set the reference weight to be used */ void referenceWeight(double w) { theReferenceWeight = w; } /** * Return true, if this sampler can provide unweighted events; if * the proposal density is not an overestimate, weights larger than * one can be generated, the handling of these points being subject * to the GeneralSampler class. */ virtual bool canUnweight() const { return true; } /** * Return true, if this sampler adapts on the fly while generating * events. Cross sections in the GeneralSampler class are calculated * from adding up the cross sections quoted by individual samplers. */ virtual bool adaptsOnTheFly() const { return false; } /** * If this sampler features a compensation algorithm, return true if * more events need to be generated to finish the compensation. */ virtual bool compensating() const { return false; } /** * Return true, if weighted events should be generated */ bool weighted() const { return theWeighted; } /** * Indicate that weighted events should be generated */ void doWeighted(bool yes = true) { theWeighted = yes; } /** * Exception to be thrown if cross section information should be updated. */ struct NextIteration {}; /** * Generate the next point and return its weight; store the point in * lastPoint(). */ virtual double generate(); /** * Fill and finalize the remappers present */ void fillRemappers(bool progress); /** * Write remappers to grid file */ void saveRemappers() const; /** * Write integration data to grid files */ void saveIntegrationData() const; /** * Save grid data */ virtual void saveGrid() const {} /** * Read integration data from grid files */ void readIntegrationData(); /** * Read remappers from grid file */ void setupRemappers(bool progress); /** * Run a single iteration of n points, optionally printing a * progress bar to cout. Calls generate n times. */ void runIteration(unsigned long n, bool progress); /** * Adapt this sampler after an iteration has been run */ virtual void adapt() {} /** * Initialize this bin sampler. This default version calls runIteration. */ virtual void initialize(bool progress); /** * Return true, if this sampler has already been initialized. */ bool initialized() const { return theInitialized; } /** * Indicate that this sampler has already been initialized. */ void isInitialized() { theInitialized = true; } /** * Return true, if integration has already been performed */ bool integrated() const { return theIntegrated; } /** * Return true, if remappers have been set up */ bool remappersFilled() const { return theRemappersFilled; } /** * Return true, if this sampler has already read grid data. */ bool hasGrids() const { return theHasGrids; } /** * Indicate that this sampler has already read grid data. */ void didReadGrids() { theHasGrids = true; } /** * Finalize this sampler. */ virtual void finalize(bool); /** * Return the total integrated cross section determined from the * Monte Carlo sampling so far. */ virtual CrossSection integratedXSec() const { return averageWeight()*nanobarn; } /** * Return the error on the total integrated cross section determined * from the Monte Carlo sampling so far. */ virtual CrossSection integratedXSecErr() const { return sqrt(abs(averageWeightVariance()))*nanobarn; } /** * Define the key for the collinear subtraction data. */ struct RandomNumberHistogram { /** * The lower bound */ double lower; /** * The bins, indexed by upper bound. */ map bins; map binsw1; /** * Constructor */ RandomNumberHistogram(double low = 0.0, double up = 1., unsigned int nbins = 20); /** * Book an event. */ void book(double inv, double weight) { map::iterator b = bins.upper_bound(inv); if ( b == bins.end() ) return; b->second = b->second+weight; map::iterator b2 = binsw1.upper_bound(inv); if ( b2 == binsw1.end() ) return; b2->second = b2->second+1.; } /** * Write to file given name and invariant. */ void dump(const std::string& folder,const std::string& prefix, const std::string& process,const int NR)const; }; typedef pair RandomNumberIndex; map > RandomNumberHistograms; public: /** * Return the dimension. */ int dimension() const { return theEventHandler->nDim(bin()); } /** * Return the number of points to be used for initial integration. */ unsigned long initialPoints() const { return theInitialPoints; } /** * Set the number of points to be used for initial integration. */ void initialPoints(unsigned long n) { theInitialPoints = n; } /** * Return the number of iterations to be considered for initialization. */ size_t nIterations() const { return theNIterations; } /** * Set the number of iterations to be considered for initialization. */ void nIterations(size_t n) { theNIterations = n; } /** * Set the factor to enhance the number of points for the next * iteration. */ void enhancementFactor(double f) { theEnhancementFactor = f; } /** * Return the factor to enhance the number of points for the next * iteration. */ double enhancementFactor() const { return theEnhancementFactor; } /** * Return the folder for the random number plots. */ string randomNumberString() const {return theRandomNumbers;} + /** + * In the AlmostUnweighted mode we do not need to unweight + * the events to the reference weight. + * Kappa reduces effectivly the reference weight. + * This can be used for processes, where unweighting + * is hardly feasable. + */ + double kappa() const {return theKappa;} + public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The bias with which this sampler is selected. */ double theBias; /** * True, if weighted events should be generated */ bool theWeighted; /** * The number of points to use for initial integration. */ unsigned long theInitialPoints; /** * The number of iterations to be considered for initialization. */ size_t theNIterations; /** * Factor to enhance the number of points for the next iteration. */ double theEnhancementFactor; /** * Switch to count only non zero weights in presampling. */ bool theNonZeroInPresampling; /** * Switch to require that we get half of the points * in each iteration below the maximum weight of the iteration. */ bool theHalfPoints; /** * The maximum number of allowed new maxima, * in combination with HalfPoints, in order to prevent unstable * processes. */ int theMaxNewMax; /** * The reference weight to be used */ double theReferenceWeight; /** * The bin to be sampled. */ int theBin; /** * Wether or not this sampler has already been initialized. */ bool theInitialized; /** * The last generated point. */ vector theLastPoint; /** * The event handler to be used. */ tStdEHPtr theEventHandler; /** * The containing sampler */ Ptr::tptr theSampler; /** * Folder for the random number plots. */ string theRandomNumbers; /** * Remapper objects indexed by dimension */ map remappers; /** * The number of points to be used for initial filling of the remappers */ unsigned long theRemapperPoints; /** * True if channels should get a remapper */ bool theRemapChannelDimension; /** * The number of bins to be used for luminosity dimensions */ unsigned long theLuminosityMapperBins; /** * The number of bins to be used for any other dimension */ unsigned long theGeneralMapperBins; /** * The minimum selection probability for remapper bins */ double theRemapperMinSelection; /** * True, if integration has already be performed */ bool theIntegrated; /** * True, if remappers have been set up */ bool theRemappersFilled; /** * True, if this sampler has already read grid data. */ bool theHasGrids; + + + /** + * In the AlmostUnweighted mode we do not need to unweight + * the events to the reference weight. + * Kappa reduces effectivly the reference weight. + * This can be used for processes, where unweighting + * is hardly feasable. + */ + double theKappa; + private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ BinSampler & operator=(const BinSampler &); }; } #endif /* Herwig_BinSampler_H */ diff --git a/Sampling/CellGrids/CellGridSampler.cc b/Sampling/CellGrids/CellGridSampler.cc --- a/Sampling/CellGrids/CellGridSampler.cc +++ b/Sampling/CellGrids/CellGridSampler.cc @@ -1,344 +1,345 @@ // -*- C++ -*- // // CellGridSampler.cpp 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 CellGridSampler class. // #include "CellGridSampler.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Repository/Repository.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Handlers/StandardEventHandler.h" #include "ThePEG/Handlers/StandardXComb.h" #include #include "CellGridSampler.h" #include "Herwig/Sampling/GeneralSampler.h" using namespace Herwig; using namespace ExSample; CellGridSampler::CellGridSampler() : BinSampler(), SimpleCellGrid(), theExplorationPoints(1000), theExplorationSteps(8), theGain(0.3), theEpsilon(0.01), theMinimumSelection(0.0001), theLuminositySplits(0), theChannelSplits(0), theAllChannelSplits(false), theUnweightCells(true) {} CellGridSampler::~CellGridSampler() {} IBPtr CellGridSampler::clone() const { return new_ptr(*this); } IBPtr CellGridSampler::fullclone() const { return new_ptr(*this); } double CellGridSampler::generate() { UseRandom rnd; double w = SimpleCellGrid::sample(rnd,*this,lastPoint(), !weighted() && initialized() && theUnweightCells, !initialized()); if ( !weighted() && initialized() ) { - double p = min(abs(w),referenceWeight())/referenceWeight(); + double p = min(abs(w),kappa()*referenceWeight())/(kappa()*referenceWeight()); double sign = w >= 0. ? 1. : -1.; if ( p < 1 && UseRandom::rnd() > p ) w = 0.; else - w = sign*max(abs(w),referenceWeight()); + w = sign*max(abs(w),referenceWeight()*kappa()); } select(w); if ( w != 0.0 ) accept(); + assert(kappa()==1.||sampler()->almostUnweighted()); return w; } void CellGridSampler::adapt() { UseRandom rnd; set newCells; SimpleCellGrid::adapt(theGain,theEpsilon,newCells); SimpleCellGrid::explore(theExplorationPoints,rnd,*this,newCells,Repository::clog()); SimpleCellGrid::setWeights(); SimpleCellGrid::updateIntegral(); SimpleCellGrid::minimumSelection(theMinimumSelection); } void CellGridSampler::saveGrid() const { XML::Element grid = SimpleCellGrid::toXML(); grid.appendAttribute("process",id()); sampler()->grids().append(grid); } void CellGridSampler::initialize(bool progress) { bool haveGrid = false; list::iterator git = sampler()->grids().children().begin(); for ( ; git != sampler()->grids().children().end(); ++git ) { if ( git->type() != XML::ElementTypes::Element ) continue; if ( git->name() != "CellGrid" ) continue; string proc; git->getFromAttribute("process",proc); if ( proc == id() ) { haveGrid = true; break; } } if ( haveGrid ) { SimpleCellGrid::fromXML(*git); sampler()->grids().erase(git); didReadGrids(); } lastPoint().resize(dimension()); if (randomNumberString()!="") for(size_t i=0;i(dimension(),0.0),vector(dimension(),1.0)); SimpleCellGrid::weightInformation().resize(dimension()); UseRandom rnd; boost::progress_display* progressBar = 0; if ( progress ) { Repository::clog() << "exploring " << process(); progressBar = new boost::progress_display(theExplorationSteps,Repository::clog()); } std::set newCells; if ( pre_adaption_splits().empty() && (theLuminositySplits || theChannelSplits || theAllChannelSplits) ) { const StandardEventHandler& eh = *eventHandler(); const StandardXComb& xc = *eh.xCombs()[bin()]; the_pre_adaption_splits.resize(dimension(),0); const pair& pdims = xc.partonDimensions(); if ( theLuminositySplits && dimension() >= pdims.first + pdims.second ) { for ( int n = 0; n < pdims.first; ++n ) the_pre_adaption_splits[n] = theLuminositySplits; for ( int n = dimension() - pdims.second; n < dimension(); ++n ) the_pre_adaption_splits[n] = theLuminositySplits; } if ( theChannelSplits && xc.diagrams().size() && dimension() > pdims.first + pdims.second ) { the_pre_adaption_splits[pdims.first] = theChannelSplits; } if ( theAllChannelSplits && xc.diagrams().size() > 1 && dimension() > pdims.first + pdims.second ) { the_pre_adaption_splits[pdims.first] = xc.diagrams().size() - 1; } } for(int splitdim=0; splitdimgrids().append(grid); if (randomNumberString()!="") for ( map >:: const_iterator b = RandomNumberHistograms.begin(); b != RandomNumberHistograms.end(); ++b ) { b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second); } } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void CellGridSampler::persistentOutput(PersistentOStream & os) const { os << theExplorationPoints << theExplorationSteps << theGain << theEpsilon << theMinimumSelection << the_pre_adaption_splits << theLuminositySplits << theChannelSplits << theAllChannelSplits << theUnweightCells; } void CellGridSampler::persistentInput(PersistentIStream & is, int) { is >> theExplorationPoints >> theExplorationSteps >> theGain >> theEpsilon >> theMinimumSelection >> the_pre_adaption_splits >> theLuminositySplits >> theChannelSplits >> theAllChannelSplits >> theUnweightCells; } // *** 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 describeHerwigCellGridSampler("Herwig::CellGridSampler", "HwSampling.so"); void CellGridSampler::Init() { static ClassDocumentation documentation ("CellGridSampler samples XCombs bins using CellGrids."); static Parameter interfaceExplorationPoints ("ExplorationPoints", "The number of points to use for cell exploration.", &CellGridSampler::theExplorationPoints, 1000, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceExplorationSteps ("ExplorationSteps", "The number of exploration steps to perform.", &CellGridSampler::theExplorationSteps, 8, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceGain ("Gain", "The gain factor used for adaption.", &CellGridSampler::theGain, 0.3, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceEpsilon ("Epsilon", "The efficieny threshold used for adaption.", &CellGridSampler::theEpsilon, 0.01, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceMinimumSelection ("MinimumSelection", "The minimum cell selection probability.", &CellGridSampler::theMinimumSelection, 0.0001, 0.0, 1.0, false, false, Interface::limited); static ParVector interfacethe_pre_adaption_splits ("preadaptionsplit", "The splittings for each dimension befor adaption.", &CellGridSampler::the_pre_adaption_splits, 1., -1, 0.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceLuminositySplits ("LuminositySplits", "", &CellGridSampler::theLuminositySplits, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfaceChannelSplits ("ChannelSplits", "", &CellGridSampler::theChannelSplits, 0, 0, 0, false, false, Interface::lowerlim); static Switch interfaceAllChannelSplits ("AllChannelSplits", "", &CellGridSampler::theAllChannelSplits, false, false, false); static SwitchOption interfaceAllChannelSplitsOn (interfaceAllChannelSplits, "On", "", true); static SwitchOption interfaceAllChannelSplitsOff (interfaceAllChannelSplits, "Off", "", false); static Switch interfaceUnweightCells ("UnweightCells", "", &CellGridSampler::theUnweightCells, true, false, false); static SwitchOption interfaceUnweightCellsYes (interfaceUnweightCells, "Yes", "", true); static SwitchOption interfaceUnweightCellsNo (interfaceUnweightCells, "No", "", false); } diff --git a/Sampling/GeneralSampler.cc b/Sampling/GeneralSampler.cc --- a/Sampling/GeneralSampler.cc +++ b/Sampling/GeneralSampler.cc @@ -1,1029 +1,1036 @@ // -*- C++ -*- // // GeneralSampler.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 GeneralSampler class. // #include "GeneralSampler.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Repository/Repository.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Utilities/LoopGuard.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Handlers/StandardEventHandler.h" #include "ThePEG/Handlers/StandardXComb.h" #include "Herwig/Utilities/RunDirectories.h" #include "Herwig/Utilities/XML/ElementIO.h" #include #include #include #include using namespace Herwig; GeneralSampler::GeneralSampler() : theVerbose(false), theIntegratedXSec(ZERO), theIntegratedXSecErr(ZERO), theUpdateAfter(1), crossSectionCalls(0), gotCrossSections(false), theSumWeights(0.), theSumWeights2(0.), theAttempts(0), theAccepts(0), theMaxWeight(0.0), theAddUpSamplers(false), theGlobalMaximumWeight(true), theFlatSubprocesses(false), isSampling(false), theMinSelection(0.01), runCombinationData(false), theAlmostUnweighted(false), maximumExceeds(0), maximumExceededBy(0.), correctWeights(0.),theMaxEnhancement(1.05), didReadGrids(false), theParallelIntegration(false), theIntegratePerJob(0), theIntegrationJobs(0), theIntegrationJobsCreated(0), justAfterIntegrate(false), theWriteGridsOnFinish(false) {} GeneralSampler::~GeneralSampler() {} IBPtr GeneralSampler::clone() const { return new_ptr(*this); } IBPtr GeneralSampler::fullclone() const { return new_ptr(*this); } double sign(double x) { return x >= 0. ? 1. : -1.; } void GeneralSampler::initialize() { if ( theParallelIntegration && runLevel() == ReadMode ) throw Exception() << "\n--------------------------------------------------------------------------------\n\n" << "Parallel integration is only supported in the build/integrate/run mode\n\n" << "--------------------------------------------------------------------------------\n" << Exception::abortnow; if ( runLevel() == ReadMode || runLevel() == IntegrationMode ) { assert(theSamplers.empty()); if ( !theGrids.children().empty() ) Repository::clog() << "--------------------------------------------------------------------------------\n\n" << "Using an existing grid. Please consider re-running the grid adaption\n" << "when there have been significant changes to parameters, cuts, etc.\n\n" << "--------------------------------------------------------------------------------\n" << flush; } if ( theParallelIntegration ) { if ( !theIntegratePerJob && !theIntegrationJobs ) throw Exception() << "Please specify the number of subprocesses per integration job or the " << "number of integration jobs to be created." << Exception::abortnow; if ( theIntegrationJobs ) { unsigned int nintegrate = eventHandler()->nBins()/theIntegrationJobs; if ( eventHandler()->nBins() % theIntegrationJobs != 0 ) ++nintegrate; theIntegratePerJob = max(theIntegratePerJob,nintegrate); } unsigned int jobCount = 0; ofstream* jobList = 0; generator()->log() << "--------------------------------------------------------------------------------\n" << "preparing integration jobs ...\n" << flush; vector randomized; vector pickfrom; for ( int b = 0; b < eventHandler()->nBins(); ++b ) pickfrom.push_back(b); //set check; while ( !pickfrom.empty() ) { size_t idx = UseRandom::irnd(pickfrom.size()); randomized.push_back(pickfrom[idx]); pickfrom.erase(pickfrom.begin() + idx); } int b = 0; for ( vector::const_iterator bx = randomized.begin(); bx != randomized.end(); ++bx, ++b ) { if ( b == 0 || b % theIntegratePerJob == 0 ) { if ( jobList ) { jobList->close(); delete jobList; jobList = 0; } ostringstream name; string prefix = RunDirectories::buildStorage(); if ( prefix.empty() ) prefix = "./"; else if ( *prefix.rbegin() != '/' ) prefix += "/"; name << prefix << "integrationJob" << jobCount; ++jobCount; string fname = name.str(); jobList = new ofstream(fname.c_str()); if ( !*jobList ) { delete jobList; throw Exception() << "Failed to write integration job list" << Exception::abortnow; } } *jobList << *bx << " "; } theIntegrationJobsCreated = jobCount; generator()->log() << "--------------------------------------------------------------------------------\n\n" << "Wrote " << jobCount << " integration jobs\n" << "Please submit integration jobs with the\nintegrate --jobid=x\ncommand for job ids " << "from 0 to " << (jobCount-1) << "\n\n" + << "e.g.:\n\n" + << " for i in $(seq 0 "<< (jobCount-1) <<");do Herwig integrate --jobid=$i "<runName()<<".run & done\n\n" << "--------------------------------------------------------------------------------\n" << flush; if ( jobList ) { jobList->close(); delete jobList; jobList = 0; } theParallelIntegration = false; return; } if ( runLevel() == BuildMode ) return; if ( !samplers().empty() ) return; if ( binSampler()->adaptsOnTheFly() ) { if ( !theAddUpSamplers ) { Repository::clog() << "Warning: On-the-fly adapting samplers require cross section calculation from " << "adding up individual samplers. The AddUpSamplers flag will be switched on."; } theAddUpSamplers = true; } if ( !weighted() && !binSampler()->canUnweight() ) throw Exception() << "Unweighted events requested from weighted bin sampler object."; if ( theFlatSubprocesses && !theGlobalMaximumWeight ) { Repository::clog() << "Warning: Can only use a global maximum weight when selecting subprocesses " << "uniformly. The GlobalMaximumWeight flag will be switched on."; theGlobalMaximumWeight = true; } set binsToIntegrate; if ( integrationList() != "" ) { string prefix = RunDirectories::buildStorage(); if ( prefix.empty() ) prefix = "./"; else if ( *prefix.rbegin() != '/' ) prefix += "/"; string fname = prefix + integrationList(); ifstream jobList(fname.c_str()); if ( jobList ) { int b = 0; while ( jobList >> b ) binsToIntegrate.insert(b); } else { Repository::clog() << "Job list '" << integrationList() << "' not found.\n" << "Assuming empty integration job\n" << flush; return; } } if ( binsToIntegrate.empty() ) { for ( int b = 0; b < eventHandler()->nBins(); ++b ) binsToIntegrate.insert(b); } boost::progress_display* progressBar = 0; if ( !theVerbose && !justAfterIntegrate ) { Repository::clog() << "integrating subprocesses"; progressBar = new boost::progress_display(binsToIntegrate.size(),Repository::clog()); } - + int count=0; for ( set::const_iterator bit = binsToIntegrate.begin(); bit != binsToIntegrate.end(); ++bit ) { + count++; + if(theVerbose&& + (runLevel() == ReadMode || + runLevel() == IntegrationMode)) + cout<<"\nIntegrate "<< count <<" of "<::ptr s = theBinSampler->cloneMe(); s->eventHandler(eventHandler()); s->sampler(this); s->bin(*bit); lastSampler(s); s->doWeighted(eventHandler()->weighted()); s->setupRemappers(theVerbose); if ( justAfterIntegrate ) s->readIntegrationData(); s->initialize(theVerbose); samplers()[*bit] = s; if ( !theVerbose && !justAfterIntegrate ) ++(*progressBar); if ( s->nanPoints() && theVerbose ) { Repository::clog() << "warning: " << s->nanPoints() << " of " << s->allPoints() << " points with nan or inf weight.\n" << flush; } } if ( progressBar ) { delete progressBar; progressBar = 0; } if ( runLevel() == IntegrationMode ) { theGrids = XML::Element(XML::ElementTypes::Element,"Grids"); for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { s->second->saveGrid(); s->second->saveRemappers(); s->second->saveIntegrationData(); } writeGrids(); return; } if ( theVerbose ) { bool oldAdd = theAddUpSamplers; theAddUpSamplers = true; try { Repository::clog() << "estimated total cross section is ( " << integratedXSec()/nanobarn << " +/- " << integratedXSecErr()/nanobarn << " ) nb\n" << flush; } catch (...) { theAddUpSamplers = oldAdd; throw; } theAddUpSamplers = oldAdd; } updateSamplers(); if ( samplers().empty() ) { throw Exception() << "No processes with non-zero cross section present." << Exception::abortnow; } if ( !justAfterIntegrate ) { theGrids = XML::Element(XML::ElementTypes::Element,"Grids"); for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { s->second->saveGrid(); s->second->saveRemappers(); } writeGrids(); } } double GeneralSampler::generate() { long excptTries = 0; gotCrossSections = false; lastSampler(samplers().upper_bound(UseRandom::rnd())->second); double weight = 0.; while ( true ) { try { weight = 1.0; double p = lastSampler()->referenceWeight()/lastSampler()->bias()/theMaxWeight; if ( weighted() ) weight *= p; else if ( p < UseRandom::rnd() ){ weight = 0.0; // The lastSampler was picked according to the bias of the process. --excptTries; } if ( weight != 0.0 ) weight *= lastSampler()->generate()/lastSampler()->referenceWeight(); } catch(BinSampler::NextIteration) { updateSamplers(); lastSampler(samplers().upper_bound(UseRandom::rnd())->second); if ( ++excptTries == eventHandler()->maxLoop() ) break; continue; } catch (...) { throw; } if ( isnan(lastSampler()->lastWeight()) || isinf(lastSampler()->lastWeight()) ) { lastSampler() = samplers().upper_bound(UseRandom::rnd())->second; if ( ++excptTries == eventHandler()->maxLoop() ) break; continue; } theAttempts += 1; if ( abs(weight) == 0.0 ) { lastSampler(samplers().upper_bound(UseRandom::rnd())->second); if ( ++excptTries == eventHandler()->maxLoop() ) break; continue; } if ( !eventHandler()->weighted() && !theAlmostUnweighted ) { if ( abs(weight) > 1. ) { ++maximumExceeds; maximumExceededBy += abs(weight)-1.; } correctWeights+=weight; if ( weight > 0.0 ) weight = 1.; else weight = -1.; } break; } theAccepts += 1; if ( excptTries == eventHandler()->maxLoop() ) throw Exception() << "GeneralSampler::generate() : Maximum number of tries to re-run event " << "selection reached. Aborting now." << Exception::runerror; lastPoint() = lastSampler()->lastPoint(); lastSampler()->accept(); theSumWeights += weight; theSumWeights2 += sqr(weight); return weight; } void GeneralSampler::rejectLast() { if ( !lastSampler() ) return; double w = 0.0; if ( weighted() ) w = lastSampler()->lastWeight()/lastSampler()->bias()/theMaxWeight; else w = lastSampler()->lastWeight()/lastSampler()->referenceWeight(); lastSampler()->reject(); theSumWeights -= w; theSumWeights2 -= sqr(w); theAttempts -= 1; theAccepts -= 1; } void GeneralSampler::updateSamplers() { map::ptr> checkedSamplers; for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { if ( s->second->averageAbsWeight() == 0.0 ) { generator()->log() << "Warning: no phase space points with non-zero cross section\n" << "could be obtained for the process: " << s->second->process() << "\n" << "This process will not be considered. Try increasing InitialPoints.\n" << flush; if ( s->second->nanPoints() ) { generator()->log() << "Warning: " << s->second->nanPoints() << " of " << s->second->allPoints() << " points with nan or inf weight\n" << "in " << s->second->process() << "\n" << flush; } continue; } checkedSamplers.insert(*s); } theSamplers = checkedSamplers; if ( samplers().empty() ) return; double allMax = 0.0; double sumbias = 0.; for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { double bias = 1.; if ( !theFlatSubprocesses ) bias *= s->second->averageAbsWeight(); s->second->bias(bias); sumbias += bias; allMax = max(allMax,s->second->maxWeight()*theMaxEnhancement); } double nsumbias = 0.0; bool needAdjust = false; for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { needAdjust |= s->second->bias()/sumbias < theMinSelection; s->second->bias(max(s->second->bias()/sumbias,theMinSelection)); nsumbias += s->second->bias(); } if ( nsumbias == 0.0 ) { samplers().clear(); return; } if ( needAdjust ) { for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { s->second->bias(s->second->bias()/nsumbias); } } theMaxWeight = 0.0; for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { double wref = theGlobalMaximumWeight ? allMax : s->second->maxWeight()*theMaxEnhancement; s->second->referenceWeight(wref); theMaxWeight = max(theMaxWeight,wref/s->second->bias()); if ( (isSampling && s->second == lastSampler()) || !isSampling ) s->second->nextIteration(); } map::ptr> newSamplers; double current = 0.; for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { if ( s->second->bias() == 0.0 ) continue; current += s->second->bias(); newSamplers[current] = s->second; } samplers() = newSamplers; } void GeneralSampler::currentCrossSections() const { if ( !theAddUpSamplers ) { double n = attempts(); if ( n > 1 ) { theIntegratedXSec = sumWeights()*maxXSec()/attempts(); double sw = sumWeights(); double sw2 = sumWeights2(); theIntegratedXSecErr = maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1)); } else { theIntegratedXSec = ZERO; theIntegratedXSecErr = ZERO; } return; } if ( gotCrossSections ) return; if ( crossSectionCalls > 0 ) { if ( ++crossSectionCalls == theUpdateAfter ) { crossSectionCalls = 0; } else return; } ++crossSectionCalls; gotCrossSections = true; theIntegratedXSec = ZERO; double var = 0.0; for ( map::ptr>::const_iterator s = samplers().begin(); s != samplers().end(); ++s ) { theIntegratedXSec += s->second->integratedXSec(); var += sqr(s->second->integratedXSecErr()/nanobarn); } theIntegratedXSecErr = sqrt(var)*nanobarn; } void GeneralSampler::prepare() { readGrids(); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void GeneralSampler::doinit() { if ( RunDirectories::empty() ) RunDirectories::pushRunId(generator()->runName()); if ( integratePerJob() || integrationJobs() ) { theParallelIntegration = true; theIntegratePerJob = integratePerJob(); theIntegrationJobs = integrationJobs(); } readGrids(); if ( theGrids.children().empty() && runLevel() == RunMode ) generator()->log() << "\n--------------------------------------------------------------------------------\n\n" << "Warning: No grid file could be found at the start of this run.\n\n" << "* For a read/run setup intented to be used with --setupfile please consider\n" << " using the build/integrate/run setup.\n" << "* For a build/integrate/run setup to be used with --setupfile please ensure\n" << " that the same setupfile is provided to both, the integrate and run steps.\n\n" << "--------------------------------------------------------------------------------\n" << flush; if ( samplers().empty() && runLevel() == RunMode ) justAfterIntegrate = true; SamplerBase::doinit(); } void GeneralSampler::dofinish() { set compensating; for ( map::ptr>::const_iterator s = samplers().begin(); s != samplers().end(); ++s ) { if ( s->second->compensating() ) { compensating.insert(s->second->process()); } if ( s->second->nanPoints() ) { generator()->log() << "warning: " << s->second->nanPoints() << " of " << s->second->allPoints() << " points with nan or inf weight\n" << "in " << s->second->process() << "\n" << flush; } s->second->finalize(theVerbose); } if ( theVerbose ) { if ( !compensating.empty() ) { generator()->log() << "warning: sampling for the following processes is still compensating:\n"; for ( set::const_iterator c = compensating.begin(); c != compensating.end(); ++c ) generator()->log() << *c << "\n"; } generator()->log() << "final integrated cross section is ( " << integratedXSec()/nanobarn << " +/- " << integratedXSecErr()/nanobarn << " ) nb\n" << flush; } if ( !compensating.empty() ) { generator()->log() << "Warning: Some samplers are still in compensating mode.\n" << flush; } if ( maximumExceeds != 0 ) { //generator()->log() << maximumExceeds << " of " << theAttempts // << " attempted points exceeded the guessed maximum weight\n" // << "with an average relative deviation of " // << maximumExceededBy/maximumExceeds << "\n\n" << flush; generator()->log() <<"\n\n\nNote: In this run "<log() <<"This corresponds to a cross section difference between:\n" <<" UnitWeights: "<< theMaxWeight*theSumWeights/theAttempts<<"nb\n" <<" AlmostUnweighted: "<< theMaxWeight*correctWeights/theAttempts<< "nb\n"<< " use 'set Sampler:AlmostUnweighted On' to switch to non-unit weights.\n\n"; generator()->log() <<"The maximum weight determined in the read/integrate step has been enhanced by \n"<< " set /Herwig/Samplers/Sampler:MaxEnhancement "<< theMaxEnhancement<< ".\nIf the rate of excessions ("<<(double)maximumExceeds*100/(double)theAccepts<< "%) or the change of the cross section is large,\nyou can try to:\n\n"<< "Enhance the number of points used in the read/integrate step\n"<< " set /Herwig/Samplers/Sampler:BinSampler:InitialPoints ...\n\n"<< "and/or enhance the reference weight found in the read/integrate step\n"<< " set /Herwig/Samplers/Sampler:MaxEnhancement 1.x\n\n"<< "If this does not help (and your process is well defined by cuts)\n"<< "don't hesitate to contact herwig@projects.hepforge.org.\n\n"; } if ( runCombinationData ) { string dataName = RunDirectories::runStorage(); if ( dataName.empty() ) dataName = "./"; else if ( *dataName.rbegin() != '/' ) dataName += "/"; dataName += "HerwigSampling.dat"; ofstream data(dataName.c_str()); double runXSec = theMaxWeight*theSumWeights/theAttempts; double runXSecErr = sqr(theMaxWeight)*(1./theAttempts)*(1./(theAttempts-1.))* abs(theSumWeights2 - sqr(theSumWeights)/theAttempts); data << setprecision(17); data << "CrossSectionCombined " << (integratedXSec()/nanobarn) << " +/- " << (integratedXSecErr()/nanobarn) << "\n" << "CrossSectionRun " << runXSec << " +/- " << sqrt(runXSecErr) << "\n" << "PointsAttempted " << theAttempts << "\n" << "PointsAccepted " << theAccepts << "\n" << "SumWeights " << theSumWeights*theMaxWeight << "\n" << "SumWeights2 " << theSumWeights2*sqr(theMaxWeight) << "\n" << flush; } theGrids = XML::Element(XML::ElementTypes::Element,"Grids"); for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { s->second->saveGrid(); s->second->saveRemappers(); if ( justAfterIntegrate ) s->second->saveIntegrationData(); } if ( theWriteGridsOnFinish ) writeGrids(); SamplerBase::dofinish(); } void GeneralSampler::doinitrun() { readGrids(); if ( theGrids.children().empty() && !didReadGrids ) generator()->log() << "\n--------------------------------------------------------------------------------\n\n" << "Warning:No grid file could be found at the start of this run.\n\n" << "* For a read/run setup intented to be used with --setupfile please consider\n" << " using the build/integrate/run setup.\n" << "* For a build/integrate/run setup to be used with --setupfile please ensure\n" << " that the same setupfile is provided to both, the integrate and run steps.\n\n" << "--------------------------------------------------------------------------------\n" << flush; if ( samplers().empty() ) { justAfterIntegrate = true; if ( !hasSetupFile() ) initialize(); } else { for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) { s->second->setupRemappers(theVerbose); if ( justAfterIntegrate ) s->second->readIntegrationData(); s->second->initialize(theVerbose); } } isSampling = true; SamplerBase::doinitrun(); } void GeneralSampler::rebind(const TranslationMap & trans) { for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) s->second = trans.translate(s->second); SamplerBase::rebind(trans); } IVector GeneralSampler::getReferences() { IVector ret = SamplerBase::getReferences(); for ( map::ptr>::iterator s = samplers().begin(); s != samplers().end(); ++s ) ret.push_back(s->second); return ret; } void GeneralSampler::writeGrids() const { if ( theGrids.children().empty() ) return; string dataName = RunDirectories::runStorage(); if ( dataName.empty() ) dataName = "./"; else if ( *dataName.rbegin() != '/' ) dataName += "/"; dataName += "HerwigGrids.xml"; ofstream out(dataName.c_str()); XML::ElementIO::put(theGrids,out); } void GeneralSampler::readGrids() { // return if grids were already read if ( didReadGrids ) return; // check for global HerwigGrids.xml file or combine integration jobs to a global HerwigGrids.xml file // Show messages of integration job combination only in the first run (if no global HerwigGrids.xml file is found in one of the directories) // or in case of an error // Check if a globalHerwigGridsFileFound was found and keep messages in a stringstream buffer beforehand bool globalHerwigGridsFileFound = false; bool integrationJobCombinationSuccessful = true; std::stringstream messageBuffer; RunDirectories directories; while ( directories && !didReadGrids ) { string dataName = directories.nextRunStorage(); if ( dataName.empty() ) dataName = "./"; else if ( *dataName.rbegin() != '/' ) dataName += "/"; string directoryName = dataName; dataName += "HerwigGrids.xml"; ifstream in(dataName.c_str()); if ( in ) { theGrids = XML::ElementIO::get(in); didReadGrids = true; // Set to true if in any of the directories a global HerwigGrid.xml file was found globalHerwigGridsFileFound = true; } else { // Check if integrationJob was split and try to merge single integrationJobs together // integrationJobsCreated() == 0 indicates that parallel integration has not been // requested, while the parallel integration parameters may well yield a single job if(integrationJobsCreated() >= 1 && runLevel() == RunMode) { messageBuffer << "\n\n* Global HerwigGrids.xml file does not exist yet" << "\n and integration jobs were split into " << integrationJobsCreated() << " integration jobs." << "\n Trying to combine single integration jobs to a global HerwigGrids.xml file" << "\n using the following directory " << directoryName << "."; theGrids = XML::Element(XML::ElementTypes::Element,"Grids"); integrationJobCombinationSuccessful = true; for(unsigned int currentProcessedIntegrationJobNum = 0; currentProcessedIntegrationJobNum < integrationJobsCreated(); ++currentProcessedIntegrationJobNum) { ostringstream currentProcessedIntegrationJob; currentProcessedIntegrationJob << directoryName << "integrationJob" << currentProcessedIntegrationJobNum << "/HerwigGrids.xml"; if(boost::filesystem::exists(boost::filesystem::path(currentProcessedIntegrationJob.str()))) { ifstream localGridFileIN(currentProcessedIntegrationJob.str().c_str()); if(localGridFileIN) { theGrids = theGrids + XML::ElementIO::get(localGridFileIN); messageBuffer << "\n* Added integration job " << currentProcessedIntegrationJobNum << " to global HerwigGrids.xml file."; } else { integrationJobCombinationSuccessful = false; messageBuffer << "\n* Could not open/add integration job " << currentProcessedIntegrationJobNum << " to global HerwigGrids.xml file."; } } else { integrationJobCombinationSuccessful = false; messageBuffer << "\n* Could not find integration job " << currentProcessedIntegrationJob.str(); } } if(integrationJobCombinationSuccessful) { string globalGridFile = directoryName + "HerwigGrids.xml"; ofstream globalGridFileOF(globalGridFile.c_str()); XML::ElementIO::put(theGrids,globalGridFileOF); messageBuffer << "\n* Global HerwigGrids.xml file was created, the integration jobs 0 to " << integrationJobsCreated()-1 << " were combined." << "\n* If previous warnings in regards to the HerwigGrids.xml file occured, these can be safely ignored." << "\n* Note: This message will occur only in the first run and will be suppressed in further runs.\n" << flush; didReadGrids = true; } else { messageBuffer << "\n* Global HerwigGrids.xml file could not be created due to failed combination of integration jobs." << "\n Please check the above-mentioned missing/failed integration jobs which are needed for the combination." << "\n* Note: It can be that the HerwigGrids.xml file is searched and can be found in further directories." << "\n In this case you can ignore this warning message.\n" << flush; } } } } // Show messages if global HerwigGrids.xml file was not found or first combination run if (!globalHerwigGridsFileFound && (theVerbose || !integrationJobCombinationSuccessful)) BaseRepository::cout() << messageBuffer.str() << "\n" << flush; if ( !didReadGrids ) theGrids = XML::Element(XML::ElementTypes::Element,"Grids"); } void GeneralSampler::persistentOutput(PersistentOStream & os) const { os << theVerbose << theBinSampler << theSamplers << theLastSampler << theUpdateAfter << crossSectionCalls << gotCrossSections << ounit(theIntegratedXSec,nanobarn) << ounit(theIntegratedXSecErr,nanobarn) << theSumWeights << theSumWeights2 << theAttempts << theAccepts << theMaxWeight << theAddUpSamplers << theGlobalMaximumWeight << theFlatSubprocesses << isSampling << theMinSelection << runCombinationData << theAlmostUnweighted << maximumExceeds << maximumExceededBy << correctWeights << theMaxEnhancement << theParallelIntegration << theIntegratePerJob << theIntegrationJobs << theIntegrationJobsCreated << theWriteGridsOnFinish; } void GeneralSampler::persistentInput(PersistentIStream & is, int) { is >> theVerbose >> theBinSampler >> theSamplers >> theLastSampler >> theUpdateAfter >> crossSectionCalls >> gotCrossSections >> iunit(theIntegratedXSec,nanobarn) >> iunit(theIntegratedXSecErr,nanobarn) >> theSumWeights >> theSumWeights2 >> theAttempts >> theAccepts >> theMaxWeight >> theAddUpSamplers >> theGlobalMaximumWeight >> theFlatSubprocesses >> isSampling >> theMinSelection >> runCombinationData >> theAlmostUnweighted >> maximumExceeds >> maximumExceededBy >> correctWeights >> theMaxEnhancement >> theParallelIntegration >> theIntegratePerJob >> theIntegrationJobs >> theIntegrationJobsCreated >> theWriteGridsOnFinish; } // *** 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 describeHerwigGeneralSampler("Herwig::GeneralSampler", "HwSampling.so"); void GeneralSampler::Init() { static ClassDocumentation documentation ("A GeneralSampler class"); static Reference interfaceBinSampler ("BinSampler", "The bin sampler to be used.", &GeneralSampler::theBinSampler, false, false, true, false, false); static Parameter interfaceUpdateAfter ("UpdateAfter", "Update cross sections every number of events.", &GeneralSampler::theUpdateAfter, 1, 1, 0, false, false, Interface::lowerlim); static Switch interfaceVerbose ("Verbose", "", &GeneralSampler::theVerbose, false, false, false); static SwitchOption interfaceVerboseOn (interfaceVerbose, "On", "", true); static SwitchOption interfaceVerboseOff (interfaceVerbose, "Off", "", false); static Switch interfaceAddUpSamplers ("AddUpSamplers", "Calculate cross sections from adding up individual samplers.", &GeneralSampler::theAddUpSamplers, false, false, false); static SwitchOption interfaceAddUpSamplersOn (interfaceAddUpSamplers, "On", "", true); static SwitchOption interfaceAddUpSamplersOff (interfaceAddUpSamplers, "Off", "", false); static Switch interfaceGlobalMaximumWeight ("GlobalMaximumWeight", "Use a global maximum weight instead of partial unweighting.", &GeneralSampler::theGlobalMaximumWeight, true, false, false); static SwitchOption interfaceGlobalMaximumWeightOn (interfaceGlobalMaximumWeight, "On", "", true); static SwitchOption interfaceGlobalMaximumWeightOff (interfaceGlobalMaximumWeight, "Off", "", false); static Parameter interfaceMaxEnhancement ("MaxEnhancement", "Enhance the maximum reference weight found in the read step.", &GeneralSampler::theMaxEnhancement, 1.1, 1.0, 1.5, false, false, Interface::limited); static Switch interfaceFlatSubprocesses ("FlatSubprocesses", "[debug] Perform a flat subprocess selection.", &GeneralSampler::theFlatSubprocesses, false, false, false); static SwitchOption interfaceFlatSubprocessesOn (interfaceFlatSubprocesses, "On", "", true); static SwitchOption interfaceFlatSubprocessesOff (interfaceFlatSubprocesses, "Off", "", false); static Parameter interfaceMinSelection ("MinSelection", "A minimum subprocess selection probability.", &GeneralSampler::theMinSelection, 0.01, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceRunCombinationData ("RunCombinationData", "", &GeneralSampler::runCombinationData, false, false, false); static SwitchOption interfaceRunCombinationDataOn (interfaceRunCombinationData, "On", "", true); static SwitchOption interfaceRunCombinationDataOff (interfaceRunCombinationData, "Off", "", false); static Switch interfaceAlmostUnweighted ("AlmostUnweighted", "", &GeneralSampler::theAlmostUnweighted, false, false, false); static SwitchOption interfaceAlmostUnweightedOn (interfaceAlmostUnweighted, "On", "", true); static SwitchOption interfaceAlmostUnweightedOff (interfaceAlmostUnweighted, "Off", "", false); static Switch interfaceParallelIntegration ("ParallelIntegration", "Prepare parallel jobs for integration.", &GeneralSampler::theParallelIntegration, false, false, false); static SwitchOption interfaceParallelIntegrationYes (interfaceParallelIntegration, "Yes", "", true); static SwitchOption interfaceParallelIntegrationNo (interfaceParallelIntegration, "No", "", false); static Parameter interfaceIntegratePerJob ("IntegratePerJob", "The number of subprocesses to integrate per job.", &GeneralSampler::theIntegratePerJob, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfaceIntegrationJobs ("IntegrationJobs", "The maximum number of integration jobs to create.", &GeneralSampler::theIntegrationJobs, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfaceIntegrationJobsCreated ("IntegrationJobsCreated", "The number of integration jobs which were actually created.", &GeneralSampler::theIntegrationJobsCreated, 1, 1, 0, false, false, Interface::lowerlim); static Switch interfaceWriteGridsOnFinish ("WriteGridsOnFinish", "Write grids on finishing a run.", &GeneralSampler::theWriteGridsOnFinish, false, false, false); static SwitchOption interfaceWriteGridsOnFinishYes (interfaceWriteGridsOnFinish, "Yes", "", true); static SwitchOption interfaceWriteGridsOnFinishNo (interfaceWriteGridsOnFinish, "No", "", false); } diff --git a/Sampling/GeneralSampler.h b/Sampling/GeneralSampler.h --- a/Sampling/GeneralSampler.h +++ b/Sampling/GeneralSampler.h @@ -1,489 +1,496 @@ // -*- C++ -*- // // GeneralSampler.h 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. // #ifndef Herwig_GeneralSampler_H #define Herwig_GeneralSampler_H // // This is the declaration of the GeneralSampler class. // #include "ThePEG/Handlers/SamplerBase.h" #include "BinSampler.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief A GeneralSampler class * * @see \ref GeneralSamplerInterfaces "The interfaces" * defined for GeneralSampler. */ class GeneralSampler: public SamplerBase { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ GeneralSampler(); /** * The destructor. */ virtual ~GeneralSampler(); //@} public: /** @name Virtual functions from SamplerBase. */ //@{ /** * Initialize the the sampler, possibly doing presampling of the * phase space. */ virtual void initialize(); /** * Generarate a new phase space point and return a weight associated * with it. This weight should preferably be 1. */ virtual double generate(); /** * Reject the last chosen phase space point. */ virtual void rejectLast(); /** * If the sampler is able to sample several different functions * separately, this function should return the last chosen * function. This default version always returns 0. */ virtual int lastBin() const { return lastSampler() ? lastSampler()->bin() : 0; } /** * Return the total integrated cross section determined from the * Monte Carlo sampling so far. */ virtual CrossSection integratedXSec() const { currentCrossSections(); return theIntegratedXSec; } /** * Return the error on the total integrated cross section determined * from the Monte Carlo sampling so far. */ virtual CrossSection integratedXSecErr() const { currentCrossSections(); return theIntegratedXSecErr; } /** * Return the overestimated integrated cross section. */ virtual CrossSection maxXSec() const { if ( theAddUpSamplers ) return SamplerBase::maxXSec(); return theMaxWeight*nanobarn; } /** * Return the sum of the weights returned by generate() so far (of * the events that were not rejeted). */ virtual double sumWeights() const { return theSumWeights; } /** * Return the sum of the weights squaredreturned by generate() so far (of * the events that were not rejeted). */ virtual double sumWeights2() const { return theSumWeights2; } /** * Return the number of attempts */ virtual double attempts() const { if ( theAddUpSamplers ) return SamplerBase::attempts(); return theAttempts; } /** * Return the number of accepts */ double accepts() const { return theAccepts; } //@} /** * Return the samplers */ const map::ptr>& samplers() const { return theSamplers; } /** * Return the bin sampler */ Ptr::ptr binSampler() const { return theBinSampler; } /** * Return the last selected bin sampler */ Ptr::tptr lastSampler() const { return theLastSampler; } /** * True if we should do weighted events */ bool weighted() const { return eventHandler()->weighted(); } + + /** + * True if the sampler runs in Allmostunweighted mode. + */ + + bool almostUnweighted() const { return theAlmostUnweighted; } + public: /** * Return the XML element containing the grids */ const XML::Element& grids() const { return theGrids; } /** * Access the XML element containing the grids */ XML::Element& grids() { return theGrids; } /** * Write out grids */ void writeGrids() const; /** * Read in grids */ void readGrids(); /** * Return the number of integration jobs which were actually created. */ unsigned int integrationJobsCreated() { return theIntegrationJobsCreated; } /** * An external hook to prepare the sampler for generating events, e.g. by * combining grid files from parallel integration runs. */ virtual void prepare(); protected: /** * Access the samplers */ map::ptr>& samplers() { return theSamplers; } /** * Set the last selected bin sampler */ void lastSampler(Ptr::tptr s) { theLastSampler = s; } /** * Calculate cross sections from samplers at current state. */ void currentCrossSections() const; /** * Update the sampler selection */ void updateSamplers(); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). protected: /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); /** * Rebind pointer to other Interfaced objects. Called in the setup phase * after all objects used in an EventGenerator has been cloned so that * the pointers will refer to the cloned objects afterwards. * @param trans a TranslationMap relating the original objects to * their respective clones. * @throws RebindException if no cloned object was found for a given * pointer. */ virtual void rebind(const TranslationMap & trans); /** * Return a vector of all pointers to Interfaced objects used in this * object. * @return a vector of pointers. */ virtual IVector getReferences(); private: /** * Whether or not additional information should be printed to cout. */ bool theVerbose; /** * The XML element containing the grids */ XML::Element theGrids; /** * The bin sampler to use. */ Ptr::ptr theBinSampler; /** * The selector map for the bin samplers. */ map::ptr> theSamplers; /** * The last selected bin sampler. */ Ptr::tptr theLastSampler; /** * The integrated cross section */ mutable CrossSection theIntegratedXSec; /** * The integrated cross section error */ mutable CrossSection theIntegratedXSecErr; /** * The number of events after which cross sections should truly be * updated. This is used to prevent exhaustive combination of * statistics when HepMC events are written out. */ size_t theUpdateAfter; /** * The number of calls to currentCrossSections since the last * update. */ mutable size_t crossSectionCalls; /** * True, if currentCrossSections has been called since the last call * to generate. */ mutable bool gotCrossSections; /** * The sum of weights */ double theSumWeights; /** * The sum of weights squared */ double theSumWeights2; /** * The number of attempts */ double theAttempts; /** * The number of accepts */ double theAccepts; /** * The maximum weight encountered */ double theMaxWeight; /** * True, if cross sections are to be combined from each sampler * individually */ bool theAddUpSamplers; /** * True, if the global maximum weight should be used as * reference. If not, the maximum weights of individual samplers are * used, and selection probabilities fro the samplers are adjusted * accordingly. */ bool theGlobalMaximumWeight; /** * True, if subprocesses should be selected flat. This is a debug * flag, cross section information and distributions will not be * correct. */ bool theFlatSubprocesses; /** * True, if we are generating events. */ bool isSampling; /** * A minimum selection probability for each sampler */ double theMinSelection; /** * True, if information for combining unnormalized runs should be * printed out */ bool runCombinationData; /** * True, if we should perform an almost unweighted sampling */ bool theAlmostUnweighted; /** * Number of points which exceeded the maximum */ unsigned long maximumExceeds; /** * The average relative deviation from the maximum weight */ double maximumExceededBy; /** * The correct cross section as one would exspect with * almostUnweighted. */ double correctWeights; /** * Enhancement factor to the maximum weight. * This is to get less maximumExceeds. */ double theMaxEnhancement; /** * True, if grids have already been read. */ bool didReadGrids; /** * True, if parallel subprocess integration should be enabled */ bool theParallelIntegration; /** * The number of subprocesses to integrate per job */ unsigned int theIntegratePerJob; /** * The maximum number of integration jobs to be created */ unsigned int theIntegrationJobs; /** * The number of integration jobs which were actually created */ unsigned int theIntegrationJobsCreated; /** * Indicate that initialization is only reading a grid. */ bool justAfterIntegrate; /** * True, if grids should be written at the end of a run */ bool theWriteGridsOnFinish; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ GeneralSampler & operator=(const GeneralSampler &); }; } #endif /* Herwig_GeneralSampler_H */ diff --git a/Sampling/MonacoSampler.cc b/Sampling/MonacoSampler.cc --- a/Sampling/MonacoSampler.cc +++ b/Sampling/MonacoSampler.cc @@ -1,398 +1,399 @@ // -*- C++ -*- // // MonacoSampler.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 MonacoSampler class. // #include "MonacoSampler.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 "ThePEG/Repository/Repository.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Handlers/StandardEventHandler.h" #include "ThePEG/Handlers/StandardXComb.h" #include #include "MonacoSampler.h" #include "Herwig/Sampling/GeneralSampler.h" using namespace Herwig; MonacoSampler::MonacoSampler() : BinSampler(), theAlpha(0.875), theGridDivisions(48), theIterationPoints(0) {} MonacoSampler::~MonacoSampler() {} IBPtr MonacoSampler::clone() const { return new_ptr(*this); } IBPtr MonacoSampler::fullclone() const { return new_ptr(*this); } double MonacoSampler::generate() { double w = 1.; // cout<<"\npoint: "; std::valarray upperb(dimension()); for ( int k = 0; k < dimension(); ++k ) { double div = (1 - UseRandom::rnd()) * theGridDivisions; upperb[k] = static_cast(div); double gupper, glower; if ( upperb[k] <= 0 ) { upperb[k] = 0; glower = 0.; gupper = theGrid(k,0); } else if (upperb[k] >= static_cast(theGridDivisions)) { upperb[k] = theGridDivisions-1; glower = theGrid(k,theGridDivisions-2); gupper = theGrid(k,theGridDivisions-1); } else { glower = theGrid(k,upperb[k]-1); gupper = theGrid(k,upperb[k]); } double gdiff = gupper - glower; lastPoint()[k] = glower + (div-upperb[k])*gdiff; w *= gdiff * theGridDivisions; } // cout<dSigDR(lastPoint()) / nanobarn; } catch (Veto&) { w = 0.0; } catch (...) { throw; } // only store numbers double wgt = w; if ( isnan(wgt) || isinf(wgt) ) wgt = 0; // save results for later grid optimization theIterationPoints++; for ( int k = 0; k < dimension(); ++k ) { theGridData(k,upperb[k]) += wgt*wgt; } if (randomNumberString()!="") for ( size_t k = 0; k < lastPoint().size(); ++k ) { RandomNumberHistograms[RandomNumberIndex(id(),k)].first.book(lastPoint()[k],wgt); RandomNumberHistograms[RandomNumberIndex(id(),k)].second+=wgt; } if ( !weighted() && initialized() ) { - double p = min(abs(w),referenceWeight())/referenceWeight(); + double p = min(abs(w),kappa()*referenceWeight())/(kappa()*referenceWeight()); double sign = w >= 0. ? 1. : -1.; if ( p < 1 && UseRandom::rnd() > p ) w = 0.; else - w = sign*max(abs(w),referenceWeight()); + w = sign*max(abs(w),kappa()*referenceWeight()); } select(w); + assert(kappa()==1.||sampler()->almostUnweighted()); if ( w != 0.0 ) accept(); return w; } void MonacoSampler::saveGrid() const { XML::Element grid = toXML(); grid.appendAttribute("process",id()); sampler()->grids().append(grid); } void MonacoSampler::initialize(bool progress) { //read in grid bool haveGrid = false; list::iterator git = sampler()->grids().children().begin(); for ( ; git != sampler()->grids().children().end(); ++git ) { if ( git->type() != XML::ElementTypes::Element ) continue; if ( git->name() != "Monaco" ) continue; string proc; git->getFromAttribute("process",proc); if ( proc == id() ) { haveGrid = true; break; } } if ( haveGrid ) { fromXML(*git); sampler()->grids().erase(git); didReadGrids(); } else { // flat grid theGrid.resize(dimension(),theGridDivisions); for (int k = 0; k < dimension(); k++) for (size_t l = 0; l < theGridDivisions; l++) theGrid(k,l) = (l+1)/static_cast(theGridDivisions); theGridData = boost::numeric::ublas::zero_matrix(dimension(),theGridDivisions); theIterationPoints = 0; } lastPoint().resize(dimension()); if (randomNumberString()!="") for(size_t i=0;igrids().children().empty() ) { // nIterations(1); // } unsigned long points = initialPoints(); for ( unsigned long k = 0; k < nIterations(); ++k ) { runIteration(points,progress); if ( k < nIterations() - 1 ) { points = (unsigned long)(points*enhancementFactor()); adapt(); nextIteration(); } } adapt(); didReadGrids(); isInitialized(); } void MonacoSampler::adapt() { int dim = dimension(); // refine grid std::valarray gridcumul(dim); for (int k=0; k ri(theGridDivisions); for (size_t l=0; l= 0) && (gridcumul[k] != 0)) { theGridData(k,l) = max( 1.0e-30, theGridData(k,l) ); double gpart = gridcumul[k] / theGridData(k,l); ri[l] = pow( (gpart - 1.0) / (gpart * log( gpart )), theAlpha); } else { ri[l] = pow( 1. / log( 1e30 ), theAlpha); } rc += ri[l]; } rc /= theGridDivisions; double gridold = 0, gridnew = 0.; double deltar = 0.; unsigned int m = 0; std::valarray theGridRowNew(theGridDivisions); for (size_t l = 0; l < theGridDivisions; ++l) { deltar += ri[l]; gridold = gridnew; gridnew = theGrid(k,l); for (; deltar > rc; m++) { deltar -= rc; theGridRowNew[m] = gridnew - (gridnew - gridold) * deltar / ri[l]; } } for (size_t l = 0; l < theGridDivisions-1; ++l) { theGrid(k,l) = theGridRowNew[l]; } theGrid(k,theGridDivisions-1) = 1.0; } theGridData = boost::numeric::ublas::zero_matrix(dimension(),theGridDivisions); theIterationPoints = 0; } void MonacoSampler::finalize(bool) { // save grid adapt(); XML::Element grid = MonacoSampler::toXML(); grid.appendAttribute("process",id()); sampler()->grids().append(grid); if (randomNumberString()!="") for ( map >:: const_iterator b = RandomNumberHistograms.begin(); b != RandomNumberHistograms.end(); ++b ) { b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second); } } void MonacoSampler::fromXML(const XML::Element& grid) { int dim = 0; grid.getFromAttribute("Dimension",dim); if ( dim != dimension() ) { throw std::runtime_error("[MonacoSampler] Number of dimensions in grid file does not match expectation."); } size_t griddivisions = 0; grid.getFromAttribute("GridDivisions",griddivisions); boost::numeric::ublas::matrix tmpgrid(dim,griddivisions); pair,list::iterator>::const_iterator,multimap,list::iterator>::const_iterator> cit; cit = grid.findAll(XML::ElementTypes::Element,"GridVector"); if ( cit.first->second == grid.children().end() ) throw std::runtime_error("[MonacoSampler] Expected a GridVector element."); for (multimap,list::iterator>::const_iterator iit=cit.first; iit!=cit.second; ++iit) { const XML::Element& gridvector = *iit->second; int k = 0; gridvector.getFromAttribute("Index",k); if ( k >= dim ) { throw std::runtime_error("[MonacoSampler] Index of grid dimension larger than grid size."); } else { list::const_iterator git; git = gridvector.findFirst(XML::ElementTypes::ParsedCharacterData,""); if ( git == gridvector.children().end() ) throw std::runtime_error("[MonacoSampler] Expected grid data."); istringstream bdata(git->content()); for ( size_t l = 0; l < griddivisions; ++l ) { bdata >> tmpgrid(k,l); } } } // store back into main variable // if griddivisions do not match, rebin preserving bin density theGrid.resize(dim,theGridDivisions); theIterationPoints = 0; double divratio = griddivisions / static_cast(theGridDivisions); for (int k = 0; k < dim; k++) { double xold = 0, xnew = 0, deltar = 0; size_t l = 0; for (size_t m = 0; m < griddivisions; m++) { deltar += 1; xold = xnew; xnew = tmpgrid(k,m); for (; deltar > divratio; l++) { deltar -= divratio; theGrid(k,l) = xnew - (xnew - xold) * deltar; } } theGrid(k,theGridDivisions-1) = 1.0; } theGridData = boost::numeric::ublas::zero_matrix(dimension(),theGridDivisions); } XML::Element MonacoSampler::toXML() const { XML::Element grid(XML::ElementTypes::Element,"Monaco"); grid.appendAttribute("Dimension",dimension()); grid.appendAttribute("GridDivisions",theGridDivisions); for ( int k = 0; k < dimension(); ++k ) { XML::Element gridvector(XML::ElementTypes::Element,"GridVector"); gridvector.appendAttribute("Index",k); ostringstream bdata; bdata << setprecision(17); for ( size_t l = 0; l < theGridDivisions; ++l ) bdata << theGrid(k,l) << " "; XML::Element belem(XML::ElementTypes::ParsedCharacterData,bdata.str()); gridvector.append(belem); grid.append(gridvector); } return grid; } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void MonacoSampler::persistentOutput(PersistentOStream & os) const { BinSampler::put(os); os << theAlpha << theGridDivisions; } void MonacoSampler::persistentInput(PersistentIStream & is, int) { BinSampler::get(is); is >> theAlpha >> theGridDivisions; } // *** 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 describeHerwigMonacoSampler("Herwig::MonacoSampler", "HwSampling.so"); void MonacoSampler::Init() { static ClassDocumentation documentation ("MonacoSampler samples XCombs bins. This implementation performs weighted MC integration using Monaco, an adapted Vegas algorithm."); static Parameter interfaceAlpha ("Alpha", "Rate of grid modification (0 for no modification).", &MonacoSampler::theAlpha, 0.875, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceGridDivisions ("GridDivisions", "The number of divisions per grid dimension.", &MonacoSampler::theGridDivisions, 48, 1, 0, false, false, Interface::lowerlim); } diff --git a/Shower/Dipole/AlphaS/lo_alpha_s.h b/Shower/Dipole/AlphaS/lo_alpha_s.h --- a/Shower/Dipole/AlphaS/lo_alpha_s.h +++ b/Shower/Dipole/AlphaS/lo_alpha_s.h @@ -1,164 +1,167 @@ // -*- C++ -*- // couplings/lo_alpha_s.h is part of matchbox // (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de #ifndef matchbox_couplings_lo_alpha_s_h #define matchbox_couplings_lo_alpha_s_h #include "alpha_s.h" namespace matchbox { using namespace ThePEG; /** * LO running alpha_s * * @see \ref lo_alpha_sInterfaces "The interfaces" * defined for lo_alpha_s. */ class lo_alpha_s : public alpha_s { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ lo_alpha_s(); /** * The destructor. */ virtual ~lo_alpha_s(); //@} public: /// return alpha_s as function of scale, QCD scale /// and number of active flavours virtual double operator () (Energy2 scale, Energy2 lambda2, unsigned int nf) const; + /// return the number of loops which determine this running + virtual unsigned int nloops () const { return 1; } + public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @name os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @name is the persistent input stream read from. * @name version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual inline void doinit() throw(InitException) { freezing_scale_ *= scale_factor(); alpha_s::doinit(); } //@} protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * The static object used to initialize the description of this class. * Indicates that this is an abstract class with persistent data. */ static ClassDescription initlo_alpha_s; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ lo_alpha_s & operator=(const lo_alpha_s &); private: Energy freezing_scale_; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of lo_alpha_s. */ template <> struct BaseClassTrait { /** Typedef of the first base class of lo_alpha_s. */ typedef matchbox::alpha_s NthBase; }; /** This template specialization informs ThePEG about the name of * the lo_alpha_s class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "matchbox::lo_alpha_s"; } /** * The name of a file containing the dynamic library where the class * lo_alpha_s is implemented. It may also include several, space-separated, * libraries if the class lo_alpha_s depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShowerAlphaS.so"; } }; /** @endcond */ } #endif /* matchbox_couplings_lo_alpha_s_h */ diff --git a/Shower/Dipole/AlphaS/nlo_alpha_s.h b/Shower/Dipole/AlphaS/nlo_alpha_s.h --- a/Shower/Dipole/AlphaS/nlo_alpha_s.h +++ b/Shower/Dipole/AlphaS/nlo_alpha_s.h @@ -1,194 +1,197 @@ // -*- C++ -*- // couplings/nlo_alpha_s.h is part of matchbox // (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de #ifndef matchbox_couplings_nlo_alpha_s_h #define matchbox_couplings_nlo_alpha_s_h #include "alpha_s.h" namespace matchbox { using namespace ThePEG; /** * NLO running alpha_s * * @see \ref nlo_alpha_sInterfaces "The interfaces" * defined for nlo_alpha_s. */ class nlo_alpha_s : public alpha_s { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ nlo_alpha_s(); /** * The destructor. */ virtual ~nlo_alpha_s(); //@} public: /// return alpha_s as function of scale, QCD scale /// and number of active flavours virtual double operator () (Energy2 scale, Energy2 lambda2, unsigned int nf) const; + /// return the number of loops which determine this running + virtual unsigned int nloops () const { return 2; } + public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @name os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @name is the persistent input stream read from. * @name version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual inline void doinit() throw(InitException) { freezing_scale_ *= scale_factor(); alpha_s::doinit(); } //@} protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * The static object used to initialize the description of this class. * Indicates that this is an abstract class with persistent data. */ static ClassDescription initnlo_alpha_s; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ nlo_alpha_s & operator=(const nlo_alpha_s &); private: struct rg_solution { inline double operator () (double alpha) { double beta0 = (33.-2.*nf)/(12.*Constants::pi); double beta1 = (153.-19.*nf)/(24.*sqr(Constants::pi)); return ((1./alpha)+(beta1/beta0)*log(alpha/(beta0+beta1*alpha))- beta0*slog); } double slog; unsigned int nf; }; Energy freezing_scale_; bool exact_evaluation_; static rg_solution& rg () { static rg_solution rg_; return rg_; } static gsl::bisection_root_solver& rg_solver () { static gsl::bisection_root_solver rg_solver_(rg()); return rg_solver_; } bool two_largeq_terms_; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of nlo_alpha_s. */ template <> struct BaseClassTrait { /** Typedef of the first base class of nlo_alpha_s. */ typedef matchbox::alpha_s NthBase; }; /** This template specialization informs ThePEG about the name of * the nlo_alpha_s class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "matchbox::nlo_alpha_s"; } /** * The name of a file containing the dynamic library where the class * nlo_alpha_s is implemented. It may also include several, space-separated, * libraries if the class nlo_alpha_s depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShowerAlphaS.so"; } }; /** @endcond */ } #endif /* matchbox_couplings_nlo_alpha_s_h */ diff --git a/Utilities/Makefile.am b/Utilities/Makefile.am --- a/Utilities/Makefile.am +++ b/Utilities/Makefile.am @@ -1,50 +1,51 @@ SUBDIRS = XML Statistics noinst_LTLIBRARIES = libHwUtils.la pkglib_LTLIBRARIES = libHwRunDirectories.la libHwUtils_la_SOURCES = \ EnumParticles.h \ Interpolator.tcc Interpolator.h \ Kinematics.cc Kinematics.h \ Maths.h Maths.cc \ StandardSelectors.cc StandardSelectors.h\ Histogram.cc Histogram.fh Histogram.h \ GaussianIntegrator.cc GaussianIntegrator.h \ GaussianIntegrator.tcc \ Statistic.h HerwigStrategy.cc HerwigStrategy.h \ GSLIntegrator.h GSLIntegrator.tcc \ -GSLBisection.h GSLBisection.tcc GSLHelper.h +GSLBisection.h GSLBisection.tcc GSLHelper.h \ +expm-1.h nodist_libHwUtils_la_SOURCES = hgstamp.inc BUILT_SOURCES = hgstamp.inc CLEANFILES = hgstamp.inc AUTOMAKE_OPTIONS = -Wno-portability HGVERSION := $(shell hg -R $(top_srcdir) parents --template '"Herwig {node|short} ({branch})"' 2> /dev/null || echo \"$(PACKAGE_STRING)\" || true ) .PHONY: update_hgstamp hgstamp.inc: update_hgstamp @[ -f $@ ] || touch $@ @echo '$(HGVERSION)' | cmp -s $@ - || echo '$(HGVERSION)' > $@ libHwUtils_la_LIBADD = \ XML/libHwXML.la \ Statistics/libHwStatistics.la libHwRunDirectories_la_SOURCES = \ RunDirectories.h RunDirectories.cc libHwRunDirectories_la_LDFLAGS = $(AM_LDFLAGS) -version-info 1:0:0 check_PROGRAMS = utilities_test utilities_test_SOURCES = \ tests/utilitiesTestsMain.cc \ tests/utilitiesTestsGlobalFixture.h \ tests/utilitiesTestsKinematics.h \ tests/utilitiesTestMaths.h \ tests/utilitiesTestsStatistic.h utilities_test_LDADD = $(BOOST_UNIT_TEST_FRAMEWORK_LIBS) $(BOOST_FILESYSTEM_LIBS) $(BOOST_SYSTEM_LIBS) $(THEPEGLIB) -ldl libHwUtils.la utilities_test_LDFLAGS = $(AM_LDFLAGS) -export-dynamic $(BOOST_UNIT_TEST_FRAMEWORK_LDFLAGS) $(BOOST_SYSTEM_LDFLAGS) $(BOOST_FILESYSTEM_LDFLAGS) $(THEPEGLDFLAGS) utilities_test_CPPFLAGS = $(AM_CPPFLAGS) $(BOOST_CPPFLAGS) -DHERWIG_PKGDATADIR="\"$(pkgdatadir)\"" -DHERWIG_PKGLIBDIR="\"$(pkglibdir)\"" -DTHEPEG_PKGLIBDIR="\"$(THEPEGLIBPATH)\"" TESTS = utilities_test diff --git a/Utilities/expm-1.h b/Utilities/expm-1.h new file mode 100644 --- /dev/null +++ b/Utilities/expm-1.h @@ -0,0 +1,147 @@ +// +// Copyright (c) 2007 +// Tsai, Dung-Bang +// National Taiwan University, Department of Physics +// +// E-Mail : dbtsai (at) gmail.com +// Begine : 2007/11/20 +// Last modify : 2007/11/22 +// Version : v0.1 +// +// EXPGM_PAD computes the matrix exponential exp(H) for general matrixs, +// including complex and real matrixs using the irreducible (p,p) degree +// rational Pade approximation to the exponential +// exp(z) = r(z)=(+/-)( I+2*(Q(z)/P(z))). +// +// Usage : +// +// U = expm_pad(H) +// U = expm_pad(H, p) +// +// where p is internally set to 6 (recommended and gererally satisfactory). +// +// See also MATLAB supplied functions, EXPM and EXPM1. +// +// Reference : +// EXPOKIT, Software Package for Computing Matrix Exponentials. +// ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 +// +// Permission to use, copy, modify, distribute and sell this software +// and its documentation for any purpose is hereby granted without fee, +// provided that the above copyright notice appear in all copies and +// that both that copyright notice and this permission notice appear +// in supporting documentation. The authors make no representations +// about the suitability of this software for any purpose. +// It is provided "as is" without express or implied warranty. +// + +#ifndef _BOOST_UBLAS_EXPM_ +#define _BOOST_UBLAS_EXPM_ +#include +#include +#include +#include + +namespace boost { namespace numeric { namespace ublas { + +template MATRIX expm_pad(const MATRIX &H, const int p = 6) { + typedef typename MATRIX::value_type value_type; + typedef typename MATRIX::size_type size_type; + typedef double real_value_type; // Correct me. Need to modify. + assert(H.size1() == H.size2()); + const size_type n = H.size1(); + const identity_matrix I(n); + matrix U(n,n),H2(n,n),P(n,n),Q(n,n); + real_value_type norm = 0.0; + + // Calcuate Pade coefficients (1-based instead of 0-based as in the c vector) + vector c(p+2); + c(1)=1; + for(size_type i = 1; i <= p; ++i) + c(i+1) = c(i) * ((p + 1.0 - i)/(i * (2.0 * p + 1 - i))); + // Calcuate the infinty norm of H, which is defined as the largest row sum of a matrix + for(size_type i=0; i(H(j,i)); // Correct me, if H is complex, can I use that abs? + norm = std::max(norm, temp); + } + if (norm == 0.0) + { + boost::throw_exception(boost::numeric::ublas::bad_argument()); + std::cerr<<"Error! Null input in the routine EXPM_PAD.\n"; + exit(0); + } + // Scaling, seek s such that || H*2^(-s) || < 1/2, and set scale = 2^(-s) + int s = 0; + real_value_type scale = 1.0; + if(norm > 0.5) { + s = std::max(0, static_cast((log(norm) / log(2.0) + 2.0))); + scale /= static_cast(std::pow(2.0, s)); + U.assign(scale * H); // Here U is used as temp value due to that H is const + } + else + U.assign(H); + // Horner evaluation of the irreducible fraction, see the following ref above. + // Initialise P (numerator) and Q (denominator) + H2.assign( prod(U, U) ); + Q.assign( c(p+1)*I ); + P.assign( c(p)*I ); + size_type odd = 1; + for( size_type k = p - 1; k > 0; --k) + { + if( odd == 1) + { + Q = ( prod(Q, H2) + c(k) * I ); + } + else + { + P = ( prod(P, H2) + c(k) * I ); + } + odd = 1 - odd; + } + if( odd == 1) + { + Q = ( prod(Q, U) ); + Q -= P ; + //U.assign( -(I + 2*(Q\P))); + } + else + { + P = (prod(P, U)); + Q -= P; + //U.assign( I + 2*(Q\P)); + } + // In origine expokit package, they use lapack ZGESV to obtain inverse matrix, + // and in that ZGESV routine, it uses LU decomposition for obtaing inverse matrix. + // Since in ublas, there is no matrix inversion template, I simply use the build-in + // LU decompostion package in ublas, and back substitute by myself. + // + //////////////// Implement Matrix Inversion /////////////////////// + permutation_matrix pm(n); + int res = lu_factorize(Q, pm); + if( res != 0) + { + std::cerr << "Error in the matrix inversion in template expm_pad.\n"; + exit(0); + } + H2 = I; // H2 is not needed anymore, so it is temporary used as identity matrix for substituting. + + lu_substitute(Q, pm, H2); + if( odd == 1) + U.assign( -(I + 2.0 * prod(H2, P))); + else + U.assign( I + 2.0 * prod(H2, P)); + // Squaring + for(size_t i = 0; i < s; ++i) + { + U = (prod(U,U)); + } + return U; + } + +}}} + + +#endif diff --git a/src/defaults/MatchboxDefaults.in.in b/src/defaults/MatchboxDefaults.in.in --- a/src/defaults/MatchboxDefaults.in.in +++ b/src/defaults/MatchboxDefaults.in.in @@ -1,781 +1,781 @@ # -*- ThePEG-repository -*- ################################################################################ # # Default setup for Matchbox matrix element generation. # You do not need to make any change in here; processes of # interest can be chosen in the standard input files. # ################################################################################ ################################################################################ # Load libraries ################################################################################ library JetCuts.so library FastJetFinder.so library HwMatchboxScales.so library HwMatchboxCuts.so library HwSampling.so library HwColorFull.so library HwMatchboxBuiltin.so ################################################################################ # Integration/sampling ################################################################################ mkdir /Herwig/Samplers cd /Herwig/Samplers create Herwig::BinSampler FlatBinSampler set FlatBinSampler:InitialPoints 1000 set FlatBinSampler:UseAllIterations No create Herwig::CellGridSampler CellGridSampler set CellGridSampler:InitialPoints 10000 set CellGridSampler:ExplorationPoints 500 set CellGridSampler:ExplorationSteps 4 set CellGridSampler:Gain 0.3 set CellGridSampler:Epsilon 1.0 set CellGridSampler:MinimumSelection 0.000001 set CellGridSampler:NIterations 1 set CellGridSampler:EnhancementFactor 1 set CellGridSampler:UseAllIterations No set CellGridSampler:RemapperPoints 50000 set CellGridSampler:RemapperMinSelection 0.00001 set CellGridSampler:RemapChannelDimension Yes set CellGridSampler:LuminosityMapperBins 20 set CellGridSampler:GeneralMapperBins 0 set CellGridSampler:HalfPoints No set CellGridSampler:MaxNewMax 30 set CellGridSampler:NonZeroInPresampling Yes create Herwig::MonacoSampler MonacoSampler set MonacoSampler:InitialPoints 15000 set MonacoSampler:NIterations 4 set MonacoSampler:EnhancementFactor 1.2 set MonacoSampler:UseAllIterations No set MonacoSampler:RemapChannelDimension No set MonacoSampler:LuminosityMapperBins 0 set MonacoSampler:HalfPoints No set MonacoSampler:MaxNewMax 30 set MonacoSampler:NonZeroInPresampling Yes create Herwig::GeneralSampler Sampler set Sampler:UpdateAfter 1000 set Sampler:BinSampler CellGridSampler set Sampler:AddUpSamplers Off set Sampler:GlobalMaximumWeight Off set Sampler:FlatSubprocesses Off set Sampler:MinSelection 0.000001 set Sampler:AlmostUnweighted Off set Sampler:RunCombinationData Off set Sampler:WriteGridsOnFinish No set Sampler:MaxEnhancement 1.1 ################################################################################ # Setup the factory object ################################################################################ mkdir /Herwig/MatrixElements/Matchbox cd /Herwig/MatrixElements/Matchbox create Herwig::MatchboxFactory Factory do Factory:StartParticleGroup p insert Factory:ParticleGroup 0 /Herwig/Particles/b insert Factory:ParticleGroup 0 /Herwig/Particles/bbar insert Factory:ParticleGroup 0 /Herwig/Particles/c insert Factory:ParticleGroup 0 /Herwig/Particles/cbar insert Factory:ParticleGroup 0 /Herwig/Particles/s insert Factory:ParticleGroup 0 /Herwig/Particles/sbar insert Factory:ParticleGroup 0 /Herwig/Particles/d insert Factory:ParticleGroup 0 /Herwig/Particles/dbar insert Factory:ParticleGroup 0 /Herwig/Particles/u insert Factory:ParticleGroup 0 /Herwig/Particles/ubar insert Factory:ParticleGroup 0 /Herwig/Particles/g do Factory:EndParticleGroup do Factory:StartParticleGroup pbar insert Factory:ParticleGroup 0 /Herwig/Particles/b insert Factory:ParticleGroup 0 /Herwig/Particles/bbar insert Factory:ParticleGroup 0 /Herwig/Particles/c insert Factory:ParticleGroup 0 /Herwig/Particles/cbar insert Factory:ParticleGroup 0 /Herwig/Particles/s insert Factory:ParticleGroup 0 /Herwig/Particles/sbar insert Factory:ParticleGroup 0 /Herwig/Particles/d insert Factory:ParticleGroup 0 /Herwig/Particles/dbar insert Factory:ParticleGroup 0 /Herwig/Particles/u insert Factory:ParticleGroup 0 /Herwig/Particles/ubar insert Factory:ParticleGroup 0 /Herwig/Particles/g do Factory:EndParticleGroup do Factory:StartParticleGroup j insert Factory:ParticleGroup 0 /Herwig/Particles/b insert Factory:ParticleGroup 0 /Herwig/Particles/bbar insert Factory:ParticleGroup 0 /Herwig/Particles/c insert Factory:ParticleGroup 0 /Herwig/Particles/cbar insert Factory:ParticleGroup 0 /Herwig/Particles/s insert Factory:ParticleGroup 0 /Herwig/Particles/sbar insert Factory:ParticleGroup 0 /Herwig/Particles/d insert Factory:ParticleGroup 0 /Herwig/Particles/dbar insert Factory:ParticleGroup 0 /Herwig/Particles/u insert Factory:ParticleGroup 0 /Herwig/Particles/ubar insert Factory:ParticleGroup 0 /Herwig/Particles/g do Factory:EndParticleGroup do Factory:StartParticleGroup u insert Factory:ParticleGroup 0 /Herwig/Particles/u do Factory:EndParticleGroup do Factory:StartParticleGroup ubar insert Factory:ParticleGroup 0 /Herwig/Particles/ubar do Factory:EndParticleGroup do Factory:StartParticleGroup d insert Factory:ParticleGroup 0 /Herwig/Particles/d do Factory:EndParticleGroup do Factory:StartParticleGroup dbar insert Factory:ParticleGroup 0 /Herwig/Particles/dbar do Factory:EndParticleGroup do Factory:StartParticleGroup s insert Factory:ParticleGroup 0 /Herwig/Particles/s do Factory:EndParticleGroup do Factory:StartParticleGroup sbar insert Factory:ParticleGroup 0 /Herwig/Particles/sbar do Factory:EndParticleGroup do Factory:StartParticleGroup c insert Factory:ParticleGroup 0 /Herwig/Particles/c do Factory:EndParticleGroup do Factory:StartParticleGroup cbar insert Factory:ParticleGroup 0 /Herwig/Particles/cbar do Factory:EndParticleGroup do Factory:StartParticleGroup b insert Factory:ParticleGroup 0 /Herwig/Particles/b do Factory:EndParticleGroup do Factory:StartParticleGroup bbar insert Factory:ParticleGroup 0 /Herwig/Particles/bbar do Factory:EndParticleGroup do Factory:StartParticleGroup t insert Factory:ParticleGroup 0 /Herwig/Particles/t do Factory:EndParticleGroup do Factory:StartParticleGroup tbar insert Factory:ParticleGroup 0 /Herwig/Particles/tbar do Factory:EndParticleGroup do Factory:StartParticleGroup g insert Factory:ParticleGroup 0 /Herwig/Particles/g do Factory:EndParticleGroup do Factory:StartParticleGroup gamma insert Factory:ParticleGroup 0 /Herwig/Particles/gamma do Factory:EndParticleGroup do Factory:StartParticleGroup h0 insert Factory:ParticleGroup 0 /Herwig/Particles/h0 do Factory:EndParticleGroup do Factory:StartParticleGroup W+ insert Factory:ParticleGroup 0 /Herwig/Particles/W+ do Factory:EndParticleGroup do Factory:StartParticleGroup W- insert Factory:ParticleGroup 0 /Herwig/Particles/W- do Factory:EndParticleGroup do Factory:StartParticleGroup Z0 insert Factory:ParticleGroup 0 /Herwig/Particles/Z0 do Factory:EndParticleGroup do Factory:StartParticleGroup e+ insert Factory:ParticleGroup 0 /Herwig/Particles/e+ do Factory:EndParticleGroup do Factory:StartParticleGroup e- insert Factory:ParticleGroup 0 /Herwig/Particles/e- do Factory:EndParticleGroup do Factory:StartParticleGroup mu+ insert Factory:ParticleGroup 0 /Herwig/Particles/mu+ do Factory:EndParticleGroup do Factory:StartParticleGroup mu- insert Factory:ParticleGroup 0 /Herwig/Particles/mu- do Factory:EndParticleGroup do Factory:StartParticleGroup tau+ insert Factory:ParticleGroup 0 /Herwig/Particles/tau+ do Factory:EndParticleGroup do Factory:StartParticleGroup tau- insert Factory:ParticleGroup 0 /Herwig/Particles/tau- do Factory:EndParticleGroup do Factory:StartParticleGroup nu_e insert Factory:ParticleGroup 0 /Herwig/Particles/nu_e do Factory:EndParticleGroup do Factory:StartParticleGroup nu_mu insert Factory:ParticleGroup 0 /Herwig/Particles/nu_mu do Factory:EndParticleGroup do Factory:StartParticleGroup nu_tau insert Factory:ParticleGroup 0 /Herwig/Particles/nu_tau do Factory:EndParticleGroup do Factory:StartParticleGroup nu_ebar insert Factory:ParticleGroup 0 /Herwig/Particles/nu_ebar do Factory:EndParticleGroup do Factory:StartParticleGroup nu_mubar insert Factory:ParticleGroup 0 /Herwig/Particles/nu_mubar do Factory:EndParticleGroup do Factory:StartParticleGroup nu_taubar insert Factory:ParticleGroup 0 /Herwig/Particles/nu_taubar do Factory:EndParticleGroup do Factory:StartParticleGroup l insert Factory:ParticleGroup 0 /Herwig/Particles/e+ insert Factory:ParticleGroup 0 /Herwig/Particles/mu+ insert Factory:ParticleGroup 0 /Herwig/Particles/e- insert Factory:ParticleGroup 0 /Herwig/Particles/mu- do Factory:EndParticleGroup do Factory:StartParticleGroup nu insert Factory:ParticleGroup 0 /Herwig/Particles/nu_e insert Factory:ParticleGroup 0 /Herwig/Particles/nu_mu insert Factory:ParticleGroup 0 /Herwig/Particles/nu_ebar insert Factory:ParticleGroup 0 /Herwig/Particles/nu_mubar do Factory:EndParticleGroup do Factory:StartParticleGroup l+ insert Factory:ParticleGroup 0 /Herwig/Particles/e+ insert Factory:ParticleGroup 0 /Herwig/Particles/mu+ do Factory:EndParticleGroup do Factory:StartParticleGroup l- insert Factory:ParticleGroup 0 /Herwig/Particles/e- insert Factory:ParticleGroup 0 /Herwig/Particles/mu- do Factory:EndParticleGroup ################################################################################ # Default settings for hard process widths ################################################################################ set /Herwig/Particles/mu+:HardProcessWidth 0*GeV set /Herwig/Particles/mu-:HardProcessWidth 0*GeV set /Herwig/Particles/tau+:HardProcessWidth 0*GeV set /Herwig/Particles/tau-:HardProcessWidth 0*GeV ################################################################################ # Setup amplitudes ################################################################################ cd /Herwig/MatrixElements/Matchbox mkdir Amplitudes cd Amplitudes create ColorFull::TraceBasis TraceBasis create Herwig::MatchboxHybridAmplitude GenericProcesses @LOAD_MADGRAPH@ HwMatchboxMadGraph.so @CREATE_MADGRAPH@ Herwig::MadGraphAmplitude MadGraph @SET_MADGRAPH@ MadGraph:ColourBasis TraceBasis @LOAD_GOSAM@ HwMatchboxGoSam.so @CREATE_GOSAM@ Herwig::GoSamAmplitude GoSam @LOAD_NJET@ HwMatchboxNJet.so @CREATE_NJET@ Herwig::NJetsAmplitude NJet @DO_NJET@ NJet:Massless 5 @DO_NJET@ NJet:Massless -5 @LOAD_OPENLOOPS@ HwMatchboxOpenLoops.so @CREATE_OPENLOOPS@ Herwig::OpenLoopsAmplitude OpenLoops @LOAD_VBFNLO@ HwMatchboxVBFNLO.so @CREATE_VBFNLO@ Herwig::VBFNLOAmplitude VBFNLO mkdir Builtin cd Builtin create Herwig::SimpleColourBasis SimpleColourBasis create Herwig::SimpleColourBasis2 SimpleColourBasis2 create Herwig::MatchboxAmplitudellbarqqbar Amplitudellbarqqbar set Amplitudellbarqqbar:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudellbarqqbarg Amplitudellbarqqbarg set Amplitudellbarqqbarg:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudellbarqqbargg Amplitudellbarqqbargg set Amplitudellbarqqbargg:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudellbarqqbarqqbar Amplitudellbarqqbarqqbar set Amplitudellbarqqbarqqbar:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudelnuqqbar Amplitudelnuqqbar set Amplitudelnuqqbar:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudelnuqqbarg Amplitudelnuqqbarg set Amplitudelnuqqbarg:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudelnuqqbargg Amplitudelnuqqbargg set Amplitudelnuqqbargg:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudelnuqqbarqqbar Amplitudelnuqqbarqqbar set Amplitudelnuqqbarqqbar:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudehgg Amplitudehgg set Amplitudehgg:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudehggg Amplitudehggg set Amplitudehggg:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudehqqbarg Amplitudehqqbarg set Amplitudehqqbarg:ColourBasis SimpleColourBasis create Herwig::MatchboxAmplitudeqqbarttbar Amplitudeqqbarttbar set Amplitudeqqbarttbar:ColourBasis SimpleColourBasis2 create Herwig::MatchboxAmplitudeqqbarttbarg Amplitudeqqbarttbarg set Amplitudeqqbarttbarg:ColourBasis SimpleColourBasis2 create Herwig::MatchboxAmplitudeggttbar Amplitudeggttbar set Amplitudeggttbar:ColourBasis SimpleColourBasis2 create Herwig::MatchboxAmplitudeggttbarg Amplitudeggttbarg set Amplitudeggttbarg:ColourBasis SimpleColourBasis2 insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbar insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbarg insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbargg insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbarqqbar insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudelnuqqbar insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudelnuqqbarg insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudelnuqqbargg insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudelnuqqbarqqbar insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehgg insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehggg insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudehqqbarg insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudeqqbarttbar insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudeqqbarttbarg insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudeggttbar insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudeggttbarg ################################################################################ # Setup phasespace generators ################################################################################ cd /Herwig/MatrixElements/Matchbox mkdir Phasespace cd Phasespace create Herwig::PhasespaceCouplings PhasespaceCouplings create Herwig::MatchboxRambo Rambo set Rambo:CouplingData PhasespaceCouplings create Herwig::FlatInvertiblePhasespace InvertiblePhasespace set InvertiblePhasespace:CouplingData PhasespaceCouplings create Herwig::FlatInvertibleLabframePhasespace InvertibleLabframePhasespace set InvertibleLabframePhasespace:CouplingData PhasespaceCouplings set InvertibleLabframePhasespace:LogSHat False create Herwig::TreePhasespaceChannels TreePhasespaceChannels create Herwig::TreePhasespace TreePhasespace set TreePhasespace:ChannelMap TreePhasespaceChannels set TreePhasespace:M0 0.0001*GeV -set TreePhasespace:MC 0.000001*GeV +set TreePhasespace:MC 0.00005*GeV set TreePhasespace:CouplingData PhasespaceCouplings do TreePhasespace:SetPhysicalCoupling 21 -1 1 0.059 do TreePhasespace:SetPhysicalCoupling 21 -2 2 0.059 do TreePhasespace:SetPhysicalCoupling 21 -3 3 0.059 do TreePhasespace:SetPhysicalCoupling 21 -4 4 0.059 do TreePhasespace:SetPhysicalCoupling 21 -5 5 0.059 do TreePhasespace:SetPhysicalCoupling 21 -6 6 0.059 do TreePhasespace:SetPhysicalCoupling 21 1 -1 0.059 do TreePhasespace:SetPhysicalCoupling 21 2 -2 0.059 do TreePhasespace:SetPhysicalCoupling 21 3 -3 0.059 do TreePhasespace:SetPhysicalCoupling 21 4 -4 0.059 do TreePhasespace:SetPhysicalCoupling 21 5 -5 0.059 do TreePhasespace:SetPhysicalCoupling 21 6 -6 0.059 do TreePhasespace:SetPhysicalCoupling 1 21 1 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 2 21 2 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 3 21 3 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 4 21 4 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 5 21 5 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 6 21 6 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -1 21 -1 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -2 21 -2 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -3 21 -3 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -4 21 -4 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -5 21 -5 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -6 21 -6 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 1 1 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 2 2 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 3 3 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 4 4 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 5 5 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling 6 6 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -1 -1 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -2 -2 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -3 -3 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -4 -4 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -5 -5 21 0.15733333333333333333 do TreePhasespace:SetPhysicalCoupling -6 -6 21 0.15733333333333333333 do TreePhasespace:SetCoupling 25 -1 1 0 do TreePhasespace:SetCoupling 25 -2 2 0 do TreePhasespace:SetCoupling 25 -3 3 0.00000001184279069851 do TreePhasespace:SetCoupling 25 -4 4 0.00000205034465001885 do TreePhasespace:SetCoupling 25 -5 5 0.00002314757096085280 do TreePhasespace:SetCoupling 25 -6 6 0.03982017320025470767 do TreePhasespace:SetCoupling 25 -11 11 0.00000000000034264835 do TreePhasespace:SetCoupling 25 -12 12 0 do TreePhasespace:SetCoupling 25 -13 13 0.00000001464912263400 do TreePhasespace:SetCoupling 25 -14 14 0 do TreePhasespace:SetCoupling 25 -15 15 0.00000414359033108195 do TreePhasespace:SetCoupling 25 -16 16 0 do TreePhasespace:SetCoupling 22 -1 1 0.00083932358497608365 do TreePhasespace:SetCoupling 22 -2 2 0.00335729433990433461 do TreePhasespace:SetCoupling 22 -3 3 0.00083932358497608365 do TreePhasespace:SetCoupling 22 -4 4 0.00335729433990433461 do TreePhasespace:SetCoupling 22 -5 5 0.00083932358497608365 do TreePhasespace:SetCoupling 22 -6 6 0.00335729433990433461 do TreePhasespace:SetCoupling 22 -11 11 0.00755391226478475287 do TreePhasespace:SetCoupling 22 -13 13 0.00755391226478475287 do TreePhasespace:SetCoupling 22 -15 15 0.00755391226478475287 do TreePhasespace:SetCoupling 24 -2 1 0.01652748072644379386 do TreePhasespace:SetCoupling 24 -4 1 0.00382028458188709739 do TreePhasespace:SetCoupling 24 -6 1 0.00014707756360995175 do TreePhasespace:SetCoupling 24 -2 3 0.00382265953677814621 do TreePhasespace:SetCoupling 24 -4 3 0.01651340063673257587 do TreePhasespace:SetCoupling 24 -6 3 0.00068534412570265868 do TreePhasespace:SetCoupling 24 -2 5 0.00005954351191129535 do TreePhasespace:SetCoupling 24 -4 5 0.00069891529650865192 do TreePhasespace:SetCoupling 24 -6 5 0.01694947628265615369 do TreePhasespace:SetCoupling 24 -12 11 0.01696396350749155147 do TreePhasespace:SetCoupling 24 -14 13 0.01696396350749155147 do TreePhasespace:SetCoupling 24 -16 15 0.01696396350749155147 do TreePhasespace:SetCoupling -24 2 -1 0.01652748072644379386 do TreePhasespace:SetCoupling -24 4 -1 0.00382028458188709739 do TreePhasespace:SetCoupling -24 6 -1 0.00014707756360995175 do TreePhasespace:SetCoupling -24 2 -3 0.00382265953677814621 do TreePhasespace:SetCoupling -24 4 -3 0.01651340063673257587 do TreePhasespace:SetCoupling -24 6 -3 0.00068534412570265868 do TreePhasespace:SetCoupling -24 2 -5 0.00005954351191129535 do TreePhasespace:SetCoupling -24 4 -5 0.00069891529650865192 do TreePhasespace:SetCoupling -24 6 -5 0.01694947628265615369 do TreePhasespace:SetCoupling -24 12 -11 0.01696396350749155147 do TreePhasespace:SetCoupling -24 14 -13 0.01696396350749155147 do TreePhasespace:SetCoupling -24 16 -15 0.01696396350749155147 do TreePhasespace:SetCoupling 23 -1 1 0.00407649129960709158 do TreePhasespace:SetCoupling 23 -2 2 0.00317809816318353030 do TreePhasespace:SetCoupling 23 -3 3 0.00407649129960709158 do TreePhasespace:SetCoupling 23 -4 4 0.00317809816318353030 do TreePhasespace:SetCoupling 23 -5 5 0.00407649129960709158 do TreePhasespace:SetCoupling 23 -6 6 0.00317809816318353030 do TreePhasespace:SetCoupling 23 -11 11 0.00276049468148072129 do TreePhasespace:SetCoupling 23 -12 12 0.00545567409075140513 do TreePhasespace:SetCoupling 23 -13 13 0.00276049468148072129 do TreePhasespace:SetCoupling 23 -14 14 0.00545567409075140513 do TreePhasespace:SetCoupling 23 -15 15 0.00276049468148072129 do TreePhasespace:SetCoupling 23 -16 16 0.00545567409075140513 do TreePhasespace:SetCoupling 21 21 21 0.354 do TreePhasespace:SetCoupling 25 21 21 0.00000000016160437564 do TreePhasespace:SetCoupling 25 25 25 0.18719783125611995353 do TreePhasespace:SetCoupling 25 22 22 0.00000000006295673620 do TreePhasespace:SetCoupling 25 24 -24 219.30463760755686425818 do TreePhasespace:SetCoupling 25 23 23 362.91922658249853887524 do TreePhasespace:SetCoupling 22 24 -24 0.00755391226478475287 do TreePhasespace:SetCoupling 23 24 -24 0.02637401475019835008 @CREATE_VBFNLO@ Herwig::VBFNLOPhasespace VBFNLOPhasespace @SET_VBFNLO@ VBFNLOPhasespace:CouplingData PhasespaceCouplings set /Herwig/MatrixElements/Matchbox/Factory:Phasespace TreePhasespace ################################################################################ # Setup utilities for matching ################################################################################ cd /Herwig/MatrixElements/Matchbox create Herwig::HardScaleProfile HardScaleProfile create Herwig::MEMatching MEMatching set MEMatching:RestrictPhasespace On set MEMatching:HardScaleProfile /Herwig/MatrixElements/Matchbox/HardScaleProfile set MEMatching:BornScaleInSubtraction BornScale set MEMatching:RealEmissionScaleInSubtraction RealScale set MEMatching:EmissionScaleInSubtraction RealScale set MEMatching:BornScaleInSplitting ShowerScale set MEMatching:RealEmissionScaleInSplitting ShowerScale set MEMatching:EmissionScaleInSplitting ShowerScale set MEMatching:TruncatedShower Yes set MEMatching:MaxPtIsMuF Yes set MEMatching:FFPtCut 1.0*GeV set MEMatching:FIPtCut 1.0*GeV set MEMatching:IIPtCut 1.0*GeV set MEMatching:SafeCut 0.*GeV create Herwig::ShowerApproximationGenerator MECorrectionHandler set MECorrectionHandler:ShowerApproximation MEMatching set MECorrectionHandler:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/InvertiblePhasespace set MECorrectionHandler:PresamplingPoints 50000 set MECorrectionHandler:FreezeGrid 100000 create Herwig::DipoleMatching DipoleMatching HwDipoleMatching.so # set in DipoleShowerDefaults.in as not available at this point # set DipoleMatching:ShowerHandler /Herwig/DipoleShower/DipoleShowerHandler set DipoleMatching:BornScaleInSubtraction BornScale set DipoleMatching:RealEmissionScaleInSubtraction BornScale set DipoleMatching:EmissionScaleInSubtraction BornScale set DipoleMatching:FFPtCut 1.0*GeV set DipoleMatching:FIPtCut 1.0*GeV set DipoleMatching:IIPtCut 1.0*GeV set DipoleMatching:SafeCut 4.*GeV create Herwig::QTildeMatching QTildeMatching HwQTildeMatching.so set QTildeMatching:ShowerHandler /Herwig/Shower/ShowerHandler set QTildeMatching:BornScaleInSubtraction BornScale set QTildeMatching:RealEmissionScaleInSubtraction BornScale set QTildeMatching:EmissionScaleInSubtraction BornScale set QTildeMatching:QTildeFinder /Herwig/Shower/PartnerFinder set QTildeMatching:SafeCut 4.*GeV # just a dummy, since SudakovCommonn can't be used # it's only used to get the value of the kinCutoffScale set QTildeMatching:QTildeSudakov /Herwig/Shower/QtoQGSudakov ################################################################################ # Setup utilities for process generation ################################################################################ cd /Herwig/MatrixElements/Matchbox mkdir Utility cd Utility create Herwig::Tree2toNGenerator DiagramGenerator insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFGVertex insert DiagramGenerator:Vertices 0 /Herwig/Vertices/GGGVertex insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFPVertex insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFZVertex cp /Herwig/Vertices/FFWVertex /Herwig/Vertices/FFWMatchboxVertex insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFWMatchboxVertex insert DiagramGenerator:Vertices 0 /Herwig/Vertices/WWHVertex insert DiagramGenerator:Vertices 0 /Herwig/Vertices/WWWVertex insert DiagramGenerator:Vertices 0 /Herwig/Vertices/HGGVertex insert DiagramGenerator:Vertices 0 /Herwig/Vertices/HHHVertex cp /Herwig/Vertices/FFHVertex /Herwig/Vertices/TTHVertex set /Herwig/Vertices/TTHVertex:Fermion 6 insert DiagramGenerator:Vertices 0 /Herwig/Vertices/TTHVertex cp /Herwig/Vertices/FFHVertex /Herwig/Vertices/BBHVertex set /Herwig/Vertices/BBHVertex:Fermion 5 create Herwig::ProcessData ProcessData set /Herwig/MatrixElements/Matchbox/Factory:DiagramGenerator DiagramGenerator set /Herwig/MatrixElements/Matchbox/Factory:ProcessData ProcessData ################################################################################ # Setup jet cuts ################################################################################ cd /Herwig/Cuts create Herwig::MatchboxFactoryMatcher MatchboxJetMatcher set MatchboxJetMatcher:Group j create ThePEG::FastJetFinder JetFinder set JetFinder:UnresolvedMatcher MatchboxJetMatcher set JetFinder:Variant AntiKt set JetFinder:RecombinationScheme E set JetFinder:Mode Inclusive set JetFinder:ConeRadius 0.7 create ThePEG::JetRegion FirstJet set FirstJet:PtMin 20.*GeV do FirstJet:YRange -5.0 5.0 set FirstJet:Fuzzy Yes set FirstJet:EnergyCutWidth 4.0*GeV set FirstJet:RapidityCutWidth 0.4 insert FirstJet:Accepts[0] 1 create ThePEG::JetRegion SecondJet set SecondJet:PtMin 20.*GeV do SecondJet:YRange -5.0 5.0 set SecondJet:Fuzzy Yes set SecondJet:EnergyCutWidth 4.0*GeV set SecondJet:RapidityCutWidth 0.4 insert SecondJet:Accepts[0] 2 create ThePEG::JetRegion ThirdJet set ThirdJet:PtMin 20.*GeV do ThirdJet:YRange -5.0 5.0 set ThirdJet:Fuzzy Yes set ThirdJet:EnergyCutWidth 4.0*GeV set ThirdJet:RapidityCutWidth 0.4 insert ThirdJet:Accepts[0] 3 create ThePEG::JetRegion FourthJet set FourthJet:PtMin 20.*GeV do FourthJet:YRange -5.0 5.0 set FourthJet:Fuzzy Yes set FourthJet:EnergyCutWidth 4.0*GeV set FourthJet:RapidityCutWidth 0.4 insert FourthJet:Accepts[0] 4 create ThePEG::FuzzyTheta FuzzyTheta set FuzzyTheta:EnergyWidth 4.0*GeV set FuzzyTheta:RapidityWidth 0.4 set FuzzyTheta:AngularWidth 0.4 create ThePEG::NJetsCut NJetsCut set NJetsCut:UnresolvedMatcher MatchboxJetMatcher set NJetsCut:NJetsMin 2 create ThePEG::JetCuts JetCuts set JetCuts:UnresolvedMatcher MatchboxJetMatcher set JetCuts:Ordering OrderPt create Herwig::IdentifiedParticleCut IdentifiedParticleCut cp IdentifiedParticleCut LeptonCut set LeptonCut:Matcher /Herwig/Matchers/Lepton cp IdentifiedParticleCut ChargedLeptonCut set ChargedLeptonCut:Matcher /Herwig/Matchers/ChargedLepton cp IdentifiedParticleCut BottomQuarkCut set BottomQuarkCut:Matcher /Herwig/Matchers/Bottom cp IdentifiedParticleCut TopQuarkCut set TopQuarkCut:Matcher /Herwig/Matchers/Top cp IdentifiedParticleCut WBosonCut set WBosonCut:Matcher /Herwig/Matchers/WBoson cp IdentifiedParticleCut ZBosonCut set ZBosonCut:Matcher /Herwig/Matchers/ZBoson cp IdentifiedParticleCut HiggsBosonCut set HiggsBosonCut:Matcher /Herwig/Matchers/HiggsBoson cp IdentifiedParticleCut PhotonCut set PhotonCut:Matcher /Herwig/Matchers/Photon create Herwig::FrixionePhotonSeparationCut PhotonIsolationCut set PhotonIsolationCut:UnresolvedMatcher MatchboxJetMatcher create Herwig::MatchboxDeltaRCut MatchboxDeltaRCut cp MatchboxDeltaRCut LeptonDeltaRCut set LeptonDeltaRCut:FirstMatcher /Herwig/Matchers/Lepton set LeptonDeltaRCut:SecondMatcher /Herwig/Matchers/Lepton cp MatchboxDeltaRCut ChargedLeptonDeltaRCut set ChargedLeptonDeltaRCut:FirstMatcher /Herwig/Matchers/ChargedLepton set ChargedLeptonDeltaRCut:SecondMatcher /Herwig/Matchers/ChargedLepton create Herwig::InvariantMassCut InvariantMassCut cp InvariantMassCut LeptonPairMassCut set LeptonPairMassCut:FirstMatcher /Herwig/Matchers/Lepton set LeptonPairMassCut:SecondMatcher /Herwig/Matchers/Lepton cp InvariantMassCut ChargedLeptonPairMassCut set ChargedLeptonPairMassCut:FirstMatcher /Herwig/Matchers/ChargedLepton set ChargedLeptonPairMassCut:SecondMatcher /Herwig/Matchers/ChargedLepton create Herwig::MissingPtCut MissingPtCut set MissingPtCut:Matcher /Herwig/Matchers/Neutrino ################################################################################ # Setup scale choices ################################################################################ cd /Herwig/MatrixElements/Matchbox mkdir Scales cd Scales create Herwig::MatchboxScaleChoice SHatScale cp SHatScale FixedScale set FixedScale:FixedScale 100.*GeV create Herwig::MatchboxPtScale MaxJetPtScale set MaxJetPtScale:JetFinder /Herwig/Cuts/JetFinder create Herwig::MatchboxLeptonMassScale LeptonPairMassScale create Herwig::MatchboxLeptonPtScale LeptonPairPtScale create Herwig::MatchboxHtScale HTScale create Herwig::MatchboxTopMassScale TopPairMassScale create Herwig::MatchboxTopMTScale TopPairMTScale set HTScale:JetFinder /Herwig/Cuts/JetFinder set HTScale:IncludeMT No set HTScale:JetPtCut 15.*GeV cp HTScale HTPrimeScale set HTPrimeScale:IncludeMT Yes set HTPrimeScale:JetPtCut 15.*GeV cp LeptonPairMassScale LeptonQ2Scale set /Herwig/MatrixElements/Matchbox/Factory:ScaleChoice LeptonPairMassScale ################################################################################ # Factories for different colliders # only provided for backwards compatibility; refer to Matchbox/*.in input file # snippets for generic handling ################################################################################ cd /Herwig/MatrixElements/Matchbox cp Factory EEFactory set EEFactory:PartonExtractor /Herwig/Partons/EEExtractor set EEFactory:Cuts /Herwig/Cuts/EECuts set EEFactory:FirstPerturbativePDF No set EEFactory:SecondPerturbativePDF No cp Factory DISFactory set DISFactory:PartonExtractor /Herwig/Partons/DISExtractor set DISFactory:Cuts /Herwig/Cuts/DISCuts set DISFactory:FirstPerturbativePDF No set DISFactory:SecondPerturbativePDF Yes cp Factory PPFactory set PPFactory:PartonExtractor /Herwig/Partons/QCDExtractor set PPFactory:Cuts /Herwig/Cuts/QCDCuts set PPFactory:FirstPerturbativePDF Yes set PPFactory:SecondPerturbativePDF Yes cd /