Page MenuHomeHEPForge
Files F7879623

No OneTemporary
Actions

View Options

diff --git a/DipoleShower/DipoleShowerHandler.cc b/DipoleShower/DipoleShowerHandler.cc
--- a/DipoleShower/DipoleShowerHandler.cc
+++ b/DipoleShower/DipoleShowerHandler.cc
@@ -1,1126 +1,1098 @@
// -*- C++ -*-
//
// DipoleShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 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 DipoleShowerHandler class.
//
#include <config.h>
#include "DipoleShowerHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
// include theses to have complete types
#include "Herwig/Shower/Base/Evolver.h"
#include "Herwig/Shower/Base/ShowerParticle.h"
#include "Herwig/PDF/MPIPDF.h"
#include "Herwig/PDF/MinBiasPDF.h"
#include "Herwig/Shower/Base/ShowerTree.h"
#include "Herwig/Shower/Base/KinematicsReconstructor.h"
#include "Herwig/Shower/Base/PartnerFinder.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include "Herwig/DipoleShower/Utility/DipolePartonSplitter.h"
#include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
using namespace Herwig;
bool DipoleShowerHandler::firstWarn = true;
DipoleShowerHandler::DipoleShowerHandler()
: ShowerHandler(), chainOrderVetoScales(true),
nEmissions(0), discardNoEmissions(false), firstMCatNLOEmission(false),
doFSR(true), doISR(true), realignmentScheme(0),
verbosity(0), printEvent(0), nTries(0),
didRadiate(false), didRealign(false),
theRenormalizationScaleFreeze(1.*GeV),
theFactorizationScaleFreeze(2.*GeV),
isMCatNLOSEvent(false),
isMCatNLOHEvent(false), theDoCompensate(false),
theFreezeGrid(500000), maxPt(ZERO),
muPt(ZERO) {}
DipoleShowerHandler::~DipoleShowerHandler() {}
IBPtr DipoleShowerHandler::clone() const {
return new_ptr(*this);
}
IBPtr DipoleShowerHandler::fullclone() const {
return new_ptr(*this);
}
tPPair DipoleShowerHandler::cascade(tSubProPtr sub, XCPtr,
Energy optHardPt, Energy optCutoff) {
useMe();
prepareCascade(sub);
- currentWeights.clear();
-
if ( !doFSR && ! doISR )
return sub->incoming();
eventRecord().clear();
eventRecord().prepare(sub,dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr()),pdfs());
if ( eventRecord().outgoing().empty() && !doISR )
return sub->incoming();
if ( !eventRecord().incoming().first->coloured() &&
!eventRecord().incoming().second->coloured() &&
!doFSR )
return sub->incoming();
nTries = 0;
while ( true ) {
try {
didRadiate = false;
didRealign = false;
isMCatNLOSEvent = false;
isMCatNLOHEvent = false;
if ( eventRecord().xcombPtr() ) {
Ptr<SubtractedME>::tptr subme =
dynamic_ptr_cast<Ptr<SubtractedME>::tptr>(eventRecord().xcombPtr()->matrixElement());
Ptr<MatchboxMEBase>::tptr me =
dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(eventRecord().xcombPtr()->matrixElement());
Ptr<SubtractionDipole>::tptr dipme =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(eventRecord().xcombPtr()->matrixElement());
if ( subme ) {
if ( subme->showerApproximation() ) {
// don't do this for POWHEG-type corrections
if ( !subme->showerApproximation()->needsSplittingGenerator() ) {
theShowerApproximation = subme->showerApproximation();
if ( subme->realShowerSubtraction() )
isMCatNLOHEvent = true;
else if ( subme->virtualShowerSubtraction() )
isMCatNLOSEvent = true;
}
}
} else if ( me ) {
if ( me->factory()->showerApproximation() ) {
if ( !me->factory()->showerApproximation()->needsSplittingGenerator() ) {
theShowerApproximation = me->factory()->showerApproximation();
isMCatNLOSEvent = true;
}
}
}
string error = "Inconsistent hard emission set-up in DipoleShowerHandler::cascade. ";
if (evolver()->hardEmissionMode()==1 || evolver()->hardEmissionMode()==3 )
throw Exception() << error
<< "Cannot generate POWHEG corrections "
<< "for particle decays using DipoleShowerHandler. "
<< "Check value of Evolver:HardEmissionMode."
<< Exception::runerror;
if ( ( isMCatNLOSEvent || isMCatNLOHEvent ) && evolver()->hardEmissionMode()==2)
throw Exception() << error
<< "Cannot generate POWHEG matching with MC@NLO shower approximation. "
<< "Add 'set Evolver:HardEmissionMode 0' to input file."
<< Exception::runerror;
if (me && me->factory()->showerApproximation()){
if(me->factory()->showerApproximation()->needsTruncatedShower())
throw Exception() << error
<< "No truncated shower needed with DipoleShowerHandler. Add "
<< "'set MEMatching:TruncatedShower No' to input file."
<< Exception::runerror;
if (!( isMCatNLOSEvent || isMCatNLOHEvent ) &&
evolver()->hardEmissionMode()==0 && firstWarn){
firstWarn=false;
throw Exception() << error
<< "Evolver:HardEmissionMode will be set to 'MatchboxPOWHEG'"
<< Exception::warning;
}
}
else if (subme && subme->factory()->showerApproximation()){
if(subme->factory()->showerApproximation()->needsTruncatedShower())
throw Exception() << error
<< "No truncated shower needed with DipoleShowerHandler. Add "
<< "'set MEMatching:TruncatedShower No' to input file."
<< Exception::runerror;
if (!( isMCatNLOSEvent || isMCatNLOHEvent ) &&
evolver()->hardEmissionMode()==0 && firstWarn){
firstWarn=false;
throw Exception() << error
<< "Evolver:HardEmissionMode will be set to 'MatchboxPOWHEG'"
<< Exception::warning;
}
}
else if (dipme && evolver()->hardEmissionMode() == 0 && firstWarn){
firstWarn=false;
throw Exception() << error
<< "Evolver:HardEmissionMode will be set to 'MatchboxPOWHEG'"
<< Exception::warning;
}
else if (!dipme && evolver()->hardEmissionMode()==2 &&
ShowerHandler::currentHandler()->firstInteraction())
throw Exception() << error
<< "POWHEG matching requested for LO events. Include "
<< "'set Factory:ShowerApproximation MEMatching' in input file."
<< Exception::runerror;
}
hardScales(lastXCombPtr()->lastCentralScale());
if ( verbosity > 1 ) {
generator()->log() << "DipoleShowerHandler starting off:\n";
eventRecord().debugLastEvent(generator()->log());
generator()->log() << flush;
}
unsigned int nEmitted = 0;
if ( firstMCatNLOEmission ) {
if ( !isMCatNLOHEvent )
nEmissions = 1;
else
nEmissions = 0;
}
if ( !firstMCatNLOEmission ) {
doCascade(nEmitted,optHardPt,optCutoff);
if ( discardNoEmissions ) {
if ( !didRadiate )
throw Veto();
if ( nEmissions )
if ( nEmissions < nEmitted )
throw Veto();
}
} else {
if ( nEmissions == 1 )
doCascade(nEmitted,optHardPt,optCutoff);
}
if ( intrinsicPtGenerator ) {
if ( eventRecord().incoming().first->coloured() &&
eventRecord().incoming().second->coloured() ) {
SpinOneLorentzRotation rot =
intrinsicPtGenerator->kick(eventRecord().incoming(),
eventRecord().intermediates());
eventRecord().transform(rot);
}
}
didRealign = realign();
constituentReshuffle();
break;
} catch (RedoShower&) {
if ( ++nTries > maxtry() )
throw ShowerTriesVeto(maxtry());
eventRecord().clear();
eventRecord().prepare(sub,dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr()),pdfs());
continue;
} catch (...) {
throw;
}
}
- tEventPtr event = eventHandler()->currentEvent();
-
- if ( event->optionalWeights().empty() ) {
- event->optionalWeights() = currentWeights;
- } else {
- map<string,double> combineMap;
- for ( map<string,double>::const_iterator w =
- event->optionalWeights().begin();
- w != event->optionalWeights().end(); ++w ) {
- for ( map<string,double>::const_iterator v = currentWeights.begin();
- v != currentWeights.end(); ++v ) {
- combineMap[w->first + "-" + v->first] = w->second*v->second;
- }
- }
- event->optionalWeights() = combineMap;
- }
-
return eventRecord().fillEventRecord(newStep(),firstInteraction(),didRealign);
}
void DipoleShowerHandler::constituentReshuffle() {
if ( constituentReshuffler ) {
constituentReshuffler->reshuffle(eventRecord().outgoing(),
eventRecord().incoming(),
eventRecord().intermediates());
}
}
void DipoleShowerHandler::hardScales(Energy2 muf) {
maxPt = generator()->maximumCMEnergy();
if ( restrictPhasespace() ) {
if ( !hardScaleIsMuF() || !firstInteraction() ) {
if ( !eventRecord().outgoing().empty() ) {
for ( PList::const_iterator p = eventRecord().outgoing().begin();
p != eventRecord().outgoing().end(); ++p )
maxPt = min(maxPt,(**p).momentum().mt());
} else {
assert(!eventRecord().hard().empty());
Lorentz5Momentum phard(ZERO,ZERO,ZERO,ZERO);
for ( PList::const_iterator p = eventRecord().hard().begin();
p != eventRecord().hard().end(); ++p )
phard += (**p).momentum();
Energy mhard = phard.m();
maxPt = mhard;
}
maxPt *= hardScaleFactor();
} else {
maxPt = hardScaleFactor()*sqrt(muf);
}
muPt = maxPt;
} else {
muPt = hardScaleFactor()*sqrt(muf);
}
for ( list<DipoleChain>::iterator ch = eventRecord().chains().begin();
ch != eventRecord().chains().end(); ++ch ) {
Energy minVetoScale = -1.*GeV;
for ( list<Dipole>::iterator dip = ch->dipoles().begin();
dip != ch->dipoles().end(); ++dip ) {
// max scale per config
Energy maxFirst = 0.0*GeV;
Energy maxSecond = 0.0*GeV;
for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k =
kernels.begin(); k != kernels.end(); ++k ) {
pair<bool,bool> conf = make_pair(true,false);
if ( (**k).canHandle(dip->index(conf)) ) {
Energy scale =
evolutionOrdering()->hardScale(dip->emitter(conf),dip->spectator(conf),
dip->emitterX(conf),dip->spectatorX(conf),
**k,dip->index(conf));
maxFirst = max(maxFirst,scale);
}
conf = make_pair(false,true);
if ( (**k).canHandle(dip->index(conf)) ) {
Energy scale =
evolutionOrdering()->hardScale(dip->emitter(conf),dip->spectator(conf),
dip->emitterX(conf),dip->spectatorX(conf),
**k,dip->index(conf));
maxSecond = max(maxSecond,scale);
}
}
if ( dip->leftParticle()->vetoScale() >= ZERO ) {
maxFirst = min(maxFirst,sqrt(dip->leftParticle()->vetoScale()));
if ( minVetoScale >= ZERO )
minVetoScale = min(minVetoScale,sqrt(dip->leftParticle()->vetoScale()));
else
minVetoScale = sqrt(dip->leftParticle()->vetoScale());
}
if ( dip->rightParticle()->vetoScale() >= ZERO ) {
maxSecond = min(maxSecond,sqrt(dip->rightParticle()->vetoScale()));
if ( minVetoScale >= ZERO )
minVetoScale = min(minVetoScale,sqrt(dip->rightParticle()->vetoScale()));
else
minVetoScale = sqrt(dip->rightParticle()->vetoScale());
}
maxFirst = min(maxPt,maxFirst);
dip->emitterScale(make_pair(true,false),maxFirst);
maxSecond = min(maxPt,maxSecond);
dip->emitterScale(make_pair(false,true),maxSecond);
}
if ( !evolutionOrdering()->independentDipoles() &&
chainOrderVetoScales &&
minVetoScale >= ZERO ) {
for ( list<Dipole>::iterator dip = ch->dipoles().begin();
dip != ch->dipoles().end(); ++dip ) {
dip->leftScale(min(dip->leftScale(),minVetoScale));
dip->rightScale(min(dip->rightScale(),minVetoScale));
}
}
}
}
Energy DipoleShowerHandler::getWinner(DipoleSplittingInfo& winner,
const Dipole& dip,
pair<bool,bool> conf,
Energy optHardPt,
Energy optCutoff) {
return
getWinner(winner,dip.index(conf),
dip.emitterX(conf),dip.spectatorX(conf),
conf,dip.emitter(conf),dip.spectator(conf),
dip.emitterScale(conf),optHardPt,optCutoff);
}
Energy DipoleShowerHandler::getWinner(SubleadingSplittingInfo& winner,
Energy optHardPt,
Energy optCutoff) {
return
getWinner(winner,winner.index(),
winner.emitterX(),winner.spectatorX(),
winner.configuration(),
winner.emitter(),winner.spectator(),
winner.startScale(),optHardPt,optCutoff);
}
Energy DipoleShowerHandler::getWinner(DipoleSplittingInfo& winner,
const DipoleIndex& index,
double emitterX, double spectatorX,
pair<bool,bool> conf,
tPPtr emitter, tPPtr spectator,
Energy startScale,
Energy optHardPt,
Energy optCutoff) {
if ( !index.initialStateEmitter() &&
!doFSR ) {
winner.didStopEvolving();
return 0.0*GeV;
}
if ( index.initialStateEmitter() &&
!doISR ) {
winner.didStopEvolving();
return 0.0*GeV;
}
DipoleSplittingInfo candidate;
candidate.index(index);
candidate.configuration(conf);
candidate.emitterX(emitterX);
candidate.spectatorX(spectatorX);
if ( generators().find(candidate.index()) == generators().end() )
getGenerators(candidate.index());
//
// NOTE -- needs proper fixing at some point
//
// For some very strange reason, equal_range gives back
// key ranges it hasn't been asked for. This particularly
// happens e.g. for FI dipoles of the same kind, but different
// PDF (hard vs MPI PDF). I can't see a reason for this,
// as DipoleIndex properly implements comparison for equality
// and (lexicographic) ordering; for the time being, we
// use equal_range, extented by an explicit check for wether
// the key is indeed what we wanted. See line after (*) comment
// below.
//
pair<GeneratorMap::iterator,GeneratorMap::iterator> gens
= generators().equal_range(candidate.index());
Energy winnerScale = 0.0*GeV;
GeneratorMap::iterator winnerGen = generators().end();
for ( GeneratorMap::iterator gen = gens.first; gen != gens.second; ++gen ) {
// (*) see NOTE above
if ( !(gen->first == candidate.index()) )
continue;
if ( startScale <= gen->second->splittingKinematics()->IRCutoff() )
continue;
Energy dScale =
gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),
spectator->momentum());
// in very exceptional cases happening in DIS
if ( isnan(dScale/GeV ) )
throw RedoShower();
candidate.scale(dScale);
candidate.continuesEvolving();
Energy hardScale = evolutionOrdering()->maxPt(startScale,candidate,*(gen->second->splittingKernel()));
Energy maxPossible =
gen->second->splittingKinematics()->ptMax(candidate.scale(),
candidate.emitterX(), candidate.spectatorX(),
candidate.index(),
*gen->second->splittingKernel());
Energy ircutoff =
optCutoff < gen->second->splittingKinematics()->IRCutoff() ?
gen->second->splittingKinematics()->IRCutoff() :
optCutoff;
if ( maxPossible <= ircutoff ) {
continue;
}
if ( maxPossible >= hardScale )
candidate.hardPt(hardScale);
else {
hardScale = maxPossible;
candidate.hardPt(maxPossible);
}
- gen->second->generate(candidate,currentWeights,optHardPt,optCutoff);
+ gen->second->generate(candidate,currentWeights_,optHardPt,optCutoff);
Energy nextScale = evolutionOrdering()->evolutionScale(gen->second->lastSplitting(),*(gen->second->splittingKernel()));
if ( nextScale > winnerScale ) {
winner.fill(candidate);
gen->second->completeSplitting(winner);
winnerGen = gen;
winnerScale = nextScale;
}
}
if ( winnerGen == generators().end() ) {
winner.didStopEvolving();
return 0.0*GeV;
}
if ( winner.stoppedEvolving() )
return 0.0*GeV;
return winnerScale;
}
void DipoleShowerHandler::doCascade(unsigned int& emDone,
Energy optHardPt,
Energy optCutoff) {
if ( nEmissions )
if ( emDone == nEmissions )
return;
DipoleSplittingInfo winner;
DipoleSplittingInfo dipoleWinner;
while ( eventRecord().haveChain() ) {
if ( verbosity > 2 ) {
generator()->log() << "DipoleShowerHandler selecting splittings for the chain:\n"
<< eventRecord().currentChain() << flush;
}
list<Dipole>::iterator winnerDip = eventRecord().currentChain().dipoles().end();
Energy winnerScale = 0.0*GeV;
Energy nextLeftScale = 0.0*GeV;
Energy nextRightScale = 0.0*GeV;
for ( list<Dipole>::iterator dip = eventRecord().currentChain().dipoles().begin();
dip != eventRecord().currentChain().dipoles().end(); ++dip ) {
nextLeftScale = getWinner(dipoleWinner,*dip,make_pair(true,false),optHardPt,optCutoff);
if ( nextLeftScale > winnerScale ) {
winnerScale = nextLeftScale;
winner = dipoleWinner;
winnerDip = dip;
}
nextRightScale = getWinner(dipoleWinner,*dip,make_pair(false,true),optHardPt,optCutoff);
if ( nextRightScale > winnerScale ) {
winnerScale = nextRightScale;
winner = dipoleWinner;
winnerDip = dip;
}
if ( evolutionOrdering()->independentDipoles() ) {
Energy dipScale = max(nextLeftScale,nextRightScale);
if ( dip->leftScale() > dipScale )
dip->leftScale(dipScale);
if ( dip->rightScale() > dipScale )
dip->rightScale(dipScale);
}
}
if ( verbosity > 1 ) {
if ( winnerDip != eventRecord().currentChain().dipoles().end() )
generator()->log() << "DipoleShowerHandler selected the splitting:\n"
<< winner << " for the dipole\n"
<< (*winnerDip) << flush;
else
generator()->log() << "DipoleShowerHandler could not select a splitting above the IR cutoff\n"
<< flush;
}
// pop the chain if no dipole did radiate
if ( winnerDip == eventRecord().currentChain().dipoles().end() ) {
if ( theEventReweight ) {
double w = theEventReweight->weightNoEmission(eventRecord().incoming(),
eventRecord().outgoing(),
eventRecord().hard(),theGlobalAlphaS);
Ptr<StandardEventHandler>::tptr eh =
dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(eventHandler());
assert(eh);
eh->reweight(w);
}
eventRecord().popChain();
continue;
}
// otherwise perform the splitting
didRadiate = true;
isMCatNLOSEvent = false;
isMCatNLOHEvent = false;
pair<list<Dipole>::iterator,list<Dipole>::iterator> children;
DipoleChain* firstChain = 0;
DipoleChain* secondChain = 0;
eventRecord().split(winnerDip,winner,children,firstChain,secondChain);
assert(firstChain && secondChain);
evolutionOrdering()->setEvolutionScale(winnerScale,winner,*firstChain,children);
if ( !secondChain->dipoles().empty() )
evolutionOrdering()->setEvolutionScale(winnerScale,winner,*secondChain,children);
if ( verbosity > 1 ) {
generator()->log() << "DipoleShowerHandler did split the last selected dipole into:\n"
<< (*children.first) << (*children.second) << flush;
}
if ( verbosity > 2 ) {
generator()->log() << "After splitting the last selected dipole, "
<< "DipoleShowerHandler encountered the following chains:\n"
<< (*firstChain) << (*secondChain) << flush;
}
if ( theEventReweight ) {
double w = theEventReweight->weight(eventRecord().incoming(),
eventRecord().outgoing(),
eventRecord().hard(),theGlobalAlphaS);
Ptr<StandardEventHandler>::tptr eh =
dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(eventHandler());
assert(eh);
eh->reweight(w);
}
if ( nEmissions )
if ( ++emDone == nEmissions )
return;
}
}
bool DipoleShowerHandler::realign() {
if ( !didRadiate && !intrinsicPtGenerator )
return false;
if ( eventRecord().incoming().first->coloured() ||
eventRecord().incoming().second->coloured() ) {
if ( eventRecord().incoming().first->momentum().perp2()/GeV2 < 1e-10 &&
eventRecord().incoming().second->momentum().perp2()/GeV2 < 1e-10 )
return false;
pair<Lorentz5Momentum,Lorentz5Momentum> inMomenta
(eventRecord().incoming().first->momentum(),
eventRecord().incoming().second->momentum());
SpinOneLorentzRotation transform((inMomenta.first+inMomenta.second).findBoostToCM());
Axis dir = (transform * inMomenta.first).vect().unit();
Axis rot (-dir.y(),dir.x(),0);
double theta = dir.theta();
if ( lastParticles().first->momentum().z() < ZERO )
theta = -theta;
transform.rotate(-theta,rot);
inMomenta.first = transform*inMomenta.first;
inMomenta.second = transform*inMomenta.second;
assert(inMomenta.first.z() > ZERO &&
inMomenta.second.z() < ZERO);
Energy2 sHat =
(eventRecord().incoming().first->momentum() +
eventRecord().incoming().second->momentum()).m2();
pair<Energy,Energy> masses(eventRecord().incoming().first->mass(),
eventRecord().incoming().second->mass());
pair<Energy,Energy> qs;
if ( !eventRecord().incoming().first->coloured() ) {
assert(masses.second == ZERO);
qs.first = eventRecord().incoming().first->momentum().z();
qs.second = (sHat-sqr(masses.first))/(2.*(qs.first+sqrt(sqr(masses.first)+sqr(qs.first))));
} else if ( !eventRecord().incoming().second->coloured() ) {
assert(masses.first == ZERO);
qs.second = eventRecord().incoming().second->momentum().z();
qs.first = (sHat-sqr(masses.second))/(2.*(qs.second+sqrt(sqr(masses.second)+sqr(qs.second))));
} else {
assert(masses.first == ZERO && masses.second == ZERO);
if ( realignmentScheme == 0 ) {
double yX = eventRecord().pX().rapidity();
double yInt = (transform*eventRecord().pX()).rapidity();
double dy = yX-yInt;
qs.first = (sqrt(sHat)/2.)*exp(dy);
qs.second = (sqrt(sHat)/2.)*exp(-dy);
} else if ( realignmentScheme == 1 ) {
Energy sS = sqrt((lastParticles().first->momentum() +
lastParticles().second->momentum()).m2());
qs.first = eventRecord().fractions().first * sS / 2.;
qs.second = eventRecord().fractions().second * sS / 2.;
}
}
double beta =
(qs.first-qs.second) /
( sqrt(sqr(masses.first)+sqr(qs.first)) +
sqrt(sqr(masses.second)+sqr(qs.second)) );
transform.boostZ(beta);
Lorentz5Momentum tmp;
if ( eventRecord().incoming().first->coloured() ) {
tmp = eventRecord().incoming().first->momentum();
tmp = transform * tmp;
eventRecord().incoming().first->set5Momentum(tmp);
}
if ( eventRecord().incoming().second->coloured() ) {
tmp = eventRecord().incoming().second->momentum();
tmp = transform * tmp;
eventRecord().incoming().second->set5Momentum(tmp);
}
eventRecord().transform(transform);
return true;
}
return false;
}
void DipoleShowerHandler::resetAlphaS(Ptr<AlphaSBase>::tptr as) {
for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k = kernels.begin();
k != kernels.end(); ++k ) {
(**k).alphaS(as);
(**k).renormalizationScaleFreeze(theRenormalizationScaleFreeze);
(**k).factorizationScaleFreeze(theFactorizationScaleFreeze);
}
// clear the generators to be rebuild
// actually, there shouldn't be any generators
// when this happens.
generators().clear();
}
void DipoleShowerHandler::resetReweight(Ptr<DipoleSplittingReweight>::tptr rw) {
for ( GeneratorMap::iterator k = generators().begin();
k != generators().end(); ++k )
k->second->splittingReweight(rw);
}
void DipoleShowerHandler::getGenerators(const DipoleIndex& ind,
Ptr<DipoleSplittingReweight>::tptr rw) {
bool gotone = false;
for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k =
kernels.begin(); k != kernels.end(); ++k ) {
if ( (**k).canHandle(ind) ) {
if ( verbosity > 0 ) {
generator()->log() << "DipoleShowerHandler encountered the dipole configuration\n"
<< ind << " in event number "
<< eventHandler()->currentEvent()->number()
<< "\nwhich can be handled by the splitting kernel '"
<< (**k).name() << "'.\n" << flush;
}
gotone = true;
Ptr<DipoleSplittingGenerator>::ptr nGenerator =
new_ptr(DipoleSplittingGenerator());
nGenerator->doCompensate(theDoCompensate);
nGenerator->splittingKernel(*k);
nGenerator->splittingKernel()->renormalizationScaleFactor(renormalizationScaleFactor());
nGenerator->splittingKernel()->factorizationScaleFactor(factorizationScaleFactor());
nGenerator->splittingKernel()->freezeGrid(theFreezeGrid);
GeneratorMap::const_iterator equivalent = generators().end();
for ( GeneratorMap::const_iterator eq = generators().begin();
eq != generators().end(); ++eq ) {
if ( !eq->second->wrapping() )
if ( (**k).canHandleEquivalent(ind,*(eq->second->splittingKernel()),eq->first) ) {
equivalent = eq;
if ( verbosity > 0 ) {
generator()->log() << "The dipole configuration "
<< ind
<< " can equivalently be handled by the existing\n"
<< "generator for configuration "
<< eq->first << " using the kernel '"
<< eq->second->splittingKernel()->name()
<< "'\n" << flush;
}
break;
}
}
if ( equivalent != generators().end() ) {
nGenerator->wrap(equivalent->second);
}
DipoleSplittingInfo dummy;
dummy.index(ind);
nGenerator->splittingReweight(rw);
nGenerator->prepare(dummy);
generators().insert(make_pair(ind,nGenerator));
}
}
if ( !gotone ) {
generator()->logWarning(Exception()
<< "DipoleShowerHandler could not "
<< "find a splitting kernel which is able "
<< "to handle splittings off the dipole "
<< ind << ".\n"
<< "Please check the input files."
<< Exception::warning);
}
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleShowerHandler::doinit() {
ShowerHandler::doinit();
if ( theGlobalAlphaS )
resetAlphaS(theGlobalAlphaS);
- for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k = kernels.begin();
- k != kernels.end(); ++k ) {
- if ( !renormalizationScaleVariations().empty() )
- addToVariations(renormalizationScaleVariations(),
- (**k).renormalizationScaleVariations());
- if ( !factorizationScaleVariations().empty() )
- addToVariations(factorizationScaleVariations(),
- (**k).factorizationScaleVariations());
- }
}
void DipoleShowerHandler::dofinish() {
ShowerHandler::dofinish();
}
void DipoleShowerHandler::doinitrun() {
ShowerHandler::doinitrun();
}
void DipoleShowerHandler::persistentOutput(PersistentOStream & os) const {
os << kernels << theEvolutionOrdering
<< constituentReshuffler << intrinsicPtGenerator
<< theGlobalAlphaS << chainOrderVetoScales
<< nEmissions << discardNoEmissions << firstMCatNLOEmission << doFSR << doISR
<< realignmentScheme << verbosity << printEvent
<< ounit(theRenormalizationScaleFreeze,GeV)
<< ounit(theFactorizationScaleFreeze,GeV)
<< isMCatNLOSEvent << isMCatNLOHEvent << theShowerApproximation
<< theDoCompensate << theFreezeGrid
<< theEventReweight << ounit(maxPt,GeV)
<< ounit(muPt,GeV);
}
void DipoleShowerHandler::persistentInput(PersistentIStream & is, int) {
is >> kernels >> theEvolutionOrdering
>> constituentReshuffler >> intrinsicPtGenerator
>> theGlobalAlphaS >> chainOrderVetoScales
>> nEmissions >> discardNoEmissions >> firstMCatNLOEmission >> doFSR >> doISR
>> realignmentScheme >> verbosity >> printEvent
>> iunit(theRenormalizationScaleFreeze,GeV)
>> iunit(theFactorizationScaleFreeze,GeV)
>> isMCatNLOSEvent >> isMCatNLOHEvent >> theShowerApproximation
>> theDoCompensate >> theFreezeGrid
>> theEventReweight >> iunit(maxPt,GeV)
>> iunit(muPt,GeV);
}
ClassDescription<DipoleShowerHandler> DipoleShowerHandler::initDipoleShowerHandler;
// Definition of the static class description member.
void DipoleShowerHandler::Init() {
static ClassDocumentation<DipoleShowerHandler> documentation
("The DipoleShowerHandler class manages the showering using "
"the dipole shower algorithm.",
"The shower evolution was performed using the algorithm described in "
"\\cite{Platzer:2009jq} and \\cite{Platzer:2011bc}.",
"%\\cite{Platzer:2009jq}\n"
"\\bibitem{Platzer:2009jq}\n"
"S.~Platzer and S.~Gieseke,\n"
"``Coherent Parton Showers with Local Recoils,''\n"
" JHEP {\\bf 1101}, 024 (2011)\n"
"arXiv:0909.5593 [hep-ph].\n"
"%%CITATION = ARXIV:0909.5593;%%\n"
"%\\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;%%");
static RefVector<DipoleShowerHandler,DipoleSplittingKernel> interfaceKernels
("Kernels",
"Set the splitting kernels to be used by the dipole shower.",
&DipoleShowerHandler::kernels, -1, false, false, true, false, false);
static Reference<DipoleShowerHandler,DipoleEvolutionOrdering> interfaceEvolutionOrdering
("EvolutionOrdering",
"Set the evolution ordering to be used.",
&DipoleShowerHandler::theEvolutionOrdering, false, false, true, false, false);
static Reference<DipoleShowerHandler,ConstituentReshuffler> interfaceConstituentReshuffler
("ConstituentReshuffler",
"The object to be used to reshuffle partons to their constitutent mass shells.",
&DipoleShowerHandler::constituentReshuffler, false, false, true, true, false);
static Reference<DipoleShowerHandler,IntrinsicPtGenerator> interfaceIntrinsicPtGenerator
("IntrinsicPtGenerator",
"Set the object in charge to generate intrinsic pt for incoming partons.",
&DipoleShowerHandler::intrinsicPtGenerator, false, false, true, true, false);
static Reference<DipoleShowerHandler,AlphaSBase> interfaceGlobalAlphaS
("GlobalAlphaS",
"Set a global strong coupling for all splitting kernels.",
&DipoleShowerHandler::theGlobalAlphaS, false, false, true, true, false);
static Switch<DipoleShowerHandler,bool> interfaceDoFSR
("DoFSR",
"Switch on or off final state radiation.",
&DipoleShowerHandler::doFSR, true, false, false);
static SwitchOption interfaceDoFSROn
(interfaceDoFSR,
"On",
"Switch on final state radiation.",
true);
static SwitchOption interfaceDoFSROff
(interfaceDoFSR,
"Off",
"Switch off final state radiation.",
false);
static Switch<DipoleShowerHandler,bool> interfaceDoISR
("DoISR",
"Switch on or off initial state radiation.",
&DipoleShowerHandler::doISR, true, false, false);
static SwitchOption interfaceDoISROn
(interfaceDoISR,
"On",
"Switch on initial state radiation.",
true);
static SwitchOption interfaceDoISROff
(interfaceDoISR,
"Off",
"Switch off initial state radiation.",
false);
static Switch<DipoleShowerHandler,int> interfaceRealignmentScheme
("RealignmentScheme",
"The realignment scheme to use.",
&DipoleShowerHandler::realignmentScheme, 0, false, false);
static SwitchOption interfaceRealignmentSchemePreserveRapidity
(interfaceRealignmentScheme,
"PreserveRapidity",
"Preserve the rapidity of non-coloured outgoing system.",
0);
static SwitchOption interfaceRealignmentSchemeEvolutionFractions
(interfaceRealignmentScheme,
"EvolutionFractions",
"Use momentum fractions as generated by the evolution.",
1);
static SwitchOption interfaceRealignmentSchemeCollisionFrame
(interfaceRealignmentScheme,
"CollisionFrame",
"Determine realignment from collision frame.",
2);
static Switch<DipoleShowerHandler,bool> interfaceChainOrderVetoScales
("ChainOrderVetoScales",
"[experimental] Switch on or off the chain ordering for veto scales.",
&DipoleShowerHandler::chainOrderVetoScales, true, false, false);
static SwitchOption interfaceChainOrderVetoScalesOn
(interfaceChainOrderVetoScales,
"On",
"Switch on chain ordering for veto scales.",
true);
static SwitchOption interfaceChainOrderVetoScalesOff
(interfaceChainOrderVetoScales,
"Off",
"Switch off chain ordering for veto scales.",
false);
interfaceChainOrderVetoScales.rank(-1);
static Parameter<DipoleShowerHandler,unsigned int> interfaceNEmissions
("NEmissions",
"[debug option] Limit the number of emissions to be generated. Zero does not limit the number of emissions.",
&DipoleShowerHandler::nEmissions, 0, 0, 0,
false, false, Interface::lowerlim);
interfaceNEmissions.rank(-1);
static Switch<DipoleShowerHandler,bool> interfaceDiscardNoEmissions
("DiscardNoEmissions",
"[debug option] Discard events without radiation.",
&DipoleShowerHandler::discardNoEmissions, false, false, false);
static SwitchOption interfaceDiscardNoEmissionsOn
(interfaceDiscardNoEmissions,
"On",
"Discard events without radiation.",
true);
static SwitchOption interfaceDiscardNoEmissionsOff
(interfaceDiscardNoEmissions,
"Off",
"Do not discard events without radiation.",
false);
interfaceDiscardNoEmissions.rank(-1);
static Switch<DipoleShowerHandler,bool> interfaceFirstMCatNLOEmission
("FirstMCatNLOEmission",
"[debug option] Only perform the first MC@NLO emission.",
&DipoleShowerHandler::firstMCatNLOEmission, false, false, false);
static SwitchOption interfaceFirstMCatNLOEmissionOn
(interfaceFirstMCatNLOEmission,
"On",
"",
true);
static SwitchOption interfaceFirstMCatNLOEmissionOff
(interfaceFirstMCatNLOEmission,
"Off",
"",
false);
interfaceFirstMCatNLOEmission.rank(-1);
static Parameter<DipoleShowerHandler,int> interfaceVerbosity
("Verbosity",
"[debug option] Set the level of debug information provided.",
&DipoleShowerHandler::verbosity, 0, 0, 0,
false, false, Interface::lowerlim);
interfaceVerbosity.rank(-1);
static Parameter<DipoleShowerHandler,int> interfacePrintEvent
("PrintEvent",
"[debug option] The number of events for which debugging information should be provided.",
&DipoleShowerHandler::printEvent, 0, 0, 0,
false, false, Interface::lowerlim);
interfacePrintEvent.rank(-1);
static Parameter<DipoleShowerHandler,Energy> interfaceRenormalizationScaleFreeze
("RenormalizationScaleFreeze",
"The freezing scale for the renormalization scale.",
&DipoleShowerHandler::theRenormalizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<DipoleShowerHandler,Energy> interfaceFactorizationScaleFreeze
("FactorizationScaleFreeze",
"The freezing scale for the factorization scale.",
&DipoleShowerHandler::theFactorizationScaleFreeze, GeV, 2.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Switch<DipoleShowerHandler,bool> interfaceDoCompensate
("DoCompensate",
"",
&DipoleShowerHandler::theDoCompensate, false, false, false);
static SwitchOption interfaceDoCompensateYes
(interfaceDoCompensate,
"Yes",
"",
true);
static SwitchOption interfaceDoCompensateNo
(interfaceDoCompensate,
"No",
"",
false);
static Parameter<DipoleShowerHandler,unsigned long> interfaceFreezeGrid
("FreezeGrid",
"",
&DipoleShowerHandler::theFreezeGrid, 500000, 1, 0,
false, false, Interface::lowerlim);
static Reference<DipoleShowerHandler,DipoleEventReweight> interfaceEventReweight
("EventReweight",
"",
&DipoleShowerHandler::theEventReweight, false, false, true, true, false);
}
diff --git a/DipoleShower/DipoleShowerHandler.h b/DipoleShower/DipoleShowerHandler.h
--- a/DipoleShower/DipoleShowerHandler.h
+++ b/DipoleShower/DipoleShowerHandler.h
@@ -1,509 +1,504 @@
// -*- C++ -*-
//
// DipoleShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 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_DipoleShowerHandler_H
#define HERWIG_DipoleShowerHandler_H
//
// This is the declaration of the DipoleShowerHandler class.
//
#include "Herwig/Shower/ShowerHandler.h"
#include "Herwig/DipoleShower/Base/DipoleSplittingInfo.h"
#include "Herwig/DipoleShower/Base/DipoleSplittingReweight.h"
#include "Herwig/DipoleShower/Kernels/DipoleSplittingKernel.h"
#include "Herwig/DipoleShower/Base/DipoleSplittingGenerator.h"
#include "Herwig/DipoleShower/Base/DipoleEventRecord.h"
#include "Herwig/DipoleShower/Base/DipoleEvolutionOrdering.h"
#include "Herwig/DipoleShower/Base/DipoleEventReweight.h"
#include "Herwig/DipoleShower/Utility/ConstituentReshuffler.h"
#include "Herwig/DipoleShower/Utility/IntrinsicPtGenerator.h"
#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer
*
* \brief The DipoleShowerHandler class manages the showering using
* the dipole shower algorithm.
*
* @see \ref DipoleShowerHandlerInterfaces "The interfaces"
* defined for DipoleShowerHandler.
*/
class DipoleShowerHandler: public ShowerHandler {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DipoleShowerHandler();
/**
* The destructor.
*/
virtual ~DipoleShowerHandler();
//@}
public:
/**
* Indicate a problem in the shower.
*/
struct RedoShower {};
/**
* Insert an additional splitting kernel.
*/
void addSplitting(Ptr<DipoleSplittingKernel>::ptr sp) {
kernels.push_back(sp);
}
/**
* Reset the alpha_s for all splitting kernels.
*/
void resetAlphaS(Ptr<AlphaSBase>::tptr);
/**
* Reset the splitting reweight for all splitting kernels.
*/
void resetReweight(Ptr<DipoleSplittingReweight>::tptr);
/**
* Return true, if the shower handler can generate a truncated
* shower for POWHEG style events generated using Matchbox
*/
virtual bool canHandleMatchboxTrunc() const { return false; }
/**
* Return true, if this cascade handler will perform reshuffling from hard
* process masses.
*/
virtual bool isReshuffling() const { return false; }
/**
* Return the relevant hard scale to be used in the profile scales
*/
virtual Energy hardScale() const {
return muPt;
}
protected:
typedef multimap<DipoleIndex,Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap;
/**
* The main method which manages the showering of a subprocess.
*/
virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb) {
return cascade(sub,xcomb,ZERO,ZERO);
}
/**
* The main method which manages the showering of a subprocess.
*/
tPPair cascade(tSubProPtr sub, XCPtr xcomb,
Energy optHardPt, Energy optCutoff);
/**
* Build splitting generators for the given
* dipole index.
*/
void getGenerators(const DipoleIndex&,
Ptr<DipoleSplittingReweight>::tptr rw =
Ptr<DipoleSplittingReweight>::tptr());
/**
* Setup the hard scales.
*/
void hardScales(Energy2 scale);
/**
* Return the evolution ordering
*/
Ptr<DipoleEvolutionOrdering>::tptr evolutionOrdering() const { return theEvolutionOrdering; }
/**
* Reshuffle to constituent mass shells
*/
void constituentReshuffle();
/**
* Access the generator map
*/
GeneratorMap& generators() { return theGenerators; }
/**
* Access the event record
*/
DipoleEventRecord& eventRecord() { return theEventRecord; }
/**
* Return the event record
*/
const DipoleEventRecord& eventRecord() const { return theEventRecord; }
/**
* Return the splitting kernels.
*/
const vector<Ptr<DipoleSplittingKernel>::ptr>& splittingKernels() const {
return kernels;
}
/**
* Realign the event such as to have the incoming partons along thre
* beam axes.
*/
bool realign();
private:
/**
* Perform the cascade.
*/
void doCascade(unsigned int& emDone,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
/**
* Get the winning splitting for the
* given dipole and configuration.
*/
Energy getWinner(DipoleSplittingInfo& winner,
const Dipole& dip,
pair<bool,bool> conf,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
/**
* Get the winning splitting for the
* given dipole and configuration.
*/
Energy getWinner(SubleadingSplittingInfo& winner,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
/**
* Get the winning splitting for the
* given dipole and configuration.
*/
Energy getWinner(DipoleSplittingInfo& winner,
const DipoleIndex& index,
double emitterX, double spectatorX,
pair<bool,bool> conf,
tPPtr emitter, tPPtr spectator,
Energy startScale,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
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:
/** @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();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
private:
/**
* The splitting kernels to be used.
*/
vector<Ptr<DipoleSplittingKernel>::ptr> kernels;
/**
* The evolution ordering considered
*/
Ptr<DipoleEvolutionOrdering>::ptr theEvolutionOrdering;
/**
* The ConstituentReshuffler to be used
*/
Ptr<ConstituentReshuffler>::ptr constituentReshuffler;
/**
* The intrinsic pt generator to be used.
*/
Ptr<IntrinsicPtGenerator>::ptr intrinsicPtGenerator;
/**
* A global alpha_s to be used for all splitting kernels.
*/
Ptr<AlphaSBase>::ptr theGlobalAlphaS;
/**
* Apply chain ordering to events from matrix
* element corrections.
*/
bool chainOrderVetoScales;
/**
* Limit the number of emissions.
* Limit applied if > 0.
*/
unsigned int nEmissions;
/**
* Discard events which did not radiate.
*/
bool discardNoEmissions;
/**
* Perform the first MC@NLO emission only.
*/
bool firstMCatNLOEmission;
/**
* Switch on or off final state radiation.
*/
bool doFSR;
/**
* Switch on or off initial state radiation.
*/
bool doISR;
/**
* The realignment scheme
*/
int realignmentScheme;
private:
/**
* The verbosity level.
* 0 - print no info
* 1 - print diagnostic information on setting up
* splitting generators etc.
* 2 - print detailed event information for up to
* printEvent events.
* 3 - print dipole chains after each splitting.
*/
int verbosity;
/**
* See verbosity.
*/
int printEvent;
private:
/**
* The splitting generators indexed by the dipole
* indices they can work on.
*/
GeneratorMap theGenerators;
/**
* The evnt record used.
*/
DipoleEventRecord theEventRecord;
/**
* The number of shoer tries so far.
*/
unsigned int nTries;
/**
* Whether or not we did radiate anything
*/
bool didRadiate;
/**
* Whether or not we did realign the event
*/
bool didRealign;
private:
/**
* A freezing value for the renormalization scale
*/
Energy theRenormalizationScaleFreeze;
/**
* A freezing value for the factorization scale
*/
Energy theFactorizationScaleFreeze;
/**
* True, if we are showering on a MC@NLO S event
*/
bool isMCatNLOSEvent;
/**
* True, if we are showering on a MC@NLO H event
*/
bool isMCatNLOHEvent;
/**
* The matching subtraction, if appropriate
*/
Ptr<ShowerApproximation>::tptr theShowerApproximation;
/**
* True, if sampler should apply compensation
*/
bool theDoCompensate;
/**
* Return the number of accepted points after which the grid should
* be frozen
*/
unsigned long theFreezeGrid;
/**
* A pointer to the dipole event reweight object
*/
Ptr<DipoleEventReweight>::ptr theEventReweight;
/**
* True if no warnings have been issued yet
*/
static bool firstWarn;
/**
* The shower starting scale for the last event encountered
*/
Energy maxPt;
/**
* The shower hard scale for the last event encountered
*/
Energy muPt;
- /**
- * The shower variation weights
- */
- map<string,double> currentWeights;
-
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<DipoleShowerHandler> initDipoleShowerHandler;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleShowerHandler & operator=(const DipoleShowerHandler &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleShowerHandler. */
template <>
struct BaseClassTrait<Herwig::DipoleShowerHandler,1> {
/** Typedef of the first base class of DipoleShowerHandler. */
typedef Herwig::ShowerHandler NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleShowerHandler class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleShowerHandler>
: public ClassTraitsBase<Herwig::DipoleShowerHandler> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleShowerHandler"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleShowerHandler is implemented. It may also include several, space-separated,
* libraries if the class DipoleShowerHandler depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_DipoleShowerHandler_H */
diff --git a/DipoleShower/Kernels/DipoleSplittingKernel.cc b/DipoleShower/Kernels/DipoleSplittingKernel.cc
--- a/DipoleShower/Kernels/DipoleSplittingKernel.cc
+++ b/DipoleShower/Kernels/DipoleSplittingKernel.cc
@@ -1,259 +1,251 @@
// -*- C++ -*-
//
// DipoleSplittingKernel.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 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 DipoleSplittingKernel class.
//
#include "DipoleSplittingKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/ShowerHandler.h"
using namespace Herwig;
DipoleSplittingKernel::DipoleSplittingKernel()
: HandlerBase(), theScreeningScale(0.0*GeV),
thePresamplingPoints(2000), theMaxtry(100000),
theFreezeGrid(500000),
theStrictLargeN(false),
theFactorizationScaleFactor(1.0),
theRenormalizationScaleFactor(1.0),
theRenormalizationScaleFreeze(1.*GeV),
theFactorizationScaleFreeze(1.*GeV),
theVirtualitySplittingScale(false) {}
DipoleSplittingKernel::~DipoleSplittingKernel() {}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleSplittingKernel::persistentOutput(PersistentOStream & os) const {
os << theAlphaS << ounit(theScreeningScale,GeV) << theSplittingKinematics << thePDFRatio
<< thePresamplingPoints << theMaxtry << theFreezeGrid
<< theFlavour << theMCCheck << theStrictLargeN
<< theFactorizationScaleFactor
<< theRenormalizationScaleFactor
<< ounit(theRenormalizationScaleFreeze,GeV)
<< ounit(theFactorizationScaleFreeze,GeV)
- << theVirtualitySplittingScale
- << theRenormalizationScaleVariations
- << theFactorizationScaleVariations;
+ << theVirtualitySplittingScale;
}
void DipoleSplittingKernel::persistentInput(PersistentIStream & is, int) {
is >> theAlphaS >> iunit(theScreeningScale,GeV) >> theSplittingKinematics >> thePDFRatio
>> thePresamplingPoints >> theMaxtry >> theFreezeGrid
>> theFlavour >> theMCCheck >> theStrictLargeN
>> theFactorizationScaleFactor
>> theRenormalizationScaleFactor
>> iunit(theRenormalizationScaleFreeze,GeV)
>> iunit(theFactorizationScaleFreeze,GeV)
- >> theVirtualitySplittingScale
- >> theRenormalizationScaleVariations
- >> theFactorizationScaleVariations;
+ >> theVirtualitySplittingScale;
}
double DipoleSplittingKernel::alphaPDF(const DipoleSplittingInfo& split,
- Energy optScale) const {
+ Energy optScale,
+ double rScaleFactor,
+ double fScaleFactor) const {
+
+ // TODO cache / precalculate the variations
Energy pt = optScale == ZERO ? split.lastPt() : optScale;
Energy2 scale = ZERO;
if ( !virtualitySplittingScale() ) {
scale = sqr(pt) + sqr(theScreeningScale);
} else {
scale = sqr(splittingKinematics()->QFromPt(pt,split)) + sqr(theScreeningScale);
}
- Energy2 rScale = sqr(theRenormalizationScaleFactor)*scale;
+ Energy2 rScale = sqr(theRenormalizationScaleFactor*rScaleFactor)*scale;
rScale = rScale > sqr(renormalizationScaleFreeze()) ? rScale : sqr(renormalizationScaleFreeze());
- Energy2 fScale = sqr(theFactorizationScaleFactor)*scale;
+ Energy2 fScale = sqr(theFactorizationScaleFactor*fScaleFactor)*scale;
fScale = fScale > sqr(factorizationScaleFreeze()) ? fScale : sqr(factorizationScaleFreeze());
double ret = alphaS()->value(rScale) / (2.*Constants::pi);
if ( split.index().initialStateEmitter() ) {
assert(pdfRatio());
ret *=
split.lastEmitterZ() *
(*pdfRatio())(split.index().emitterPDF(), fScale,
split.index().emitterData(),split.emitterData(),
split.emitterX(),split.lastEmitterZ());
}
if ( split.index().initialStateSpectator() ) {
assert(pdfRatio());
ret *=
split.lastSpectatorZ() *
(*pdfRatio())(split.index().spectatorPDF(), fScale,
split.index().spectatorData(),split.spectatorData(),
split.spectatorX(),split.lastSpectatorZ());
}
if ( ret < 0. )
ret = 0.;
return ret;
}
-void DipoleSplittingKernel::accept(const DipoleSplittingInfo&,
+void DipoleSplittingKernel::accept(const DipoleSplittingInfo& split,
double, double,
- map<string,double>&) const {
+ map<string,double>& vars) const {
}
-void DipoleSplittingKernel::veto(const DipoleSplittingInfo&,
- double, double,
- map<string,double>&) const {
-}
-
-string DipoleSplittingKernel::doRenormalizationScaleVariations(string in) {
- return ShowerHandler::parseVariations(in,theRenormalizationScaleVariations);
-}
-
-string DipoleSplittingKernel::doFactorizationScaleVariations(string in) {
- return ShowerHandler::parseVariations(in,theFactorizationScaleVariations);
+void DipoleSplittingKernel::veto(const DipoleSplittingInfo& split,
+ double p, double r,
+ map<string,double>& var) const {
}
AbstractClassDescription<DipoleSplittingKernel> DipoleSplittingKernel::initDipoleSplittingKernel;
// Definition of the static class description member.
void DipoleSplittingKernel::Init() {
static ClassDocumentation<DipoleSplittingKernel> documentation
("DipoleSplittingKernel is the base class for all kernels "
"used within the dipole shower.");
static Reference<DipoleSplittingKernel,AlphaSBase> interfaceAlphaS
("AlphaS",
"The strong coupling to be used by this splitting kernel.",
&DipoleSplittingKernel::theAlphaS, false, false, true, true, false);
static Parameter<DipoleSplittingKernel,Energy> interfaceScreeningScale
("ScreeningScale",
"A colour screening scale",
&DipoleSplittingKernel::theScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Reference<DipoleSplittingKernel,DipoleSplittingKinematics> interfaceSplittingKinematics
("SplittingKinematics",
"The splitting kinematics to be used by this splitting kernel.",
&DipoleSplittingKernel::theSplittingKinematics, false, false, true, false, false);
static Reference<DipoleSplittingKernel,PDFRatio> interfacePDFRatio
("PDFRatio",
"Set the optional PDF ratio object to evaluate this kernel",
&DipoleSplittingKernel::thePDFRatio, false, false, true, true, false);
static Parameter<DipoleSplittingKernel,unsigned long> interfacePresamplingPoints
("PresamplingPoints",
"The number of points used to presample this kernel.",
&DipoleSplittingKernel::thePresamplingPoints, 2000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKernel,unsigned long> interfaceMaxtry
("Maxtry",
"The maximum number of attempts to generate a splitting.",
&DipoleSplittingKernel::theMaxtry, 10000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKernel,unsigned long> interfaceFreezeGrid
("FreezeGrid",
"",
&DipoleSplittingKernel::theFreezeGrid, 500000, 1, 0,
false, false, Interface::lowerlim);
static Reference<DipoleSplittingKernel,ParticleData> interfaceFlavour
("Flavour",
"Set the flavour to be produced if ambiguous.",
&DipoleSplittingKernel::theFlavour, false, false, true, true, false);
static Reference<DipoleSplittingKernel,DipoleMCCheck> interfaceMCCheck
("MCCheck",
"[debug option] MCCheck",
&DipoleSplittingKernel::theMCCheck, false, false, true, true, false);
interfaceMCCheck.rank(-1);
static Switch<DipoleSplittingKernel,bool> interfaceStrictLargeN
("StrictLargeN",
"Work in a strict large-N limit.",
&DipoleSplittingKernel::theStrictLargeN, false, false, false);
static SwitchOption interfaceStrictLargeNOn
(interfaceStrictLargeN,
"On",
"Replace C_F -> C_A/2 where present",
true);
static SwitchOption interfaceStrictLargeNOff
(interfaceStrictLargeN,
"Off",
"Keep C_F=4/3",
false);
interfaceStrictLargeN.rank(-2);
static Parameter<DipoleSplittingKernel,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&DipoleSplittingKernel::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
interfaceFactorizationScaleFactor.rank(-2);
static Parameter<DipoleSplittingKernel,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&DipoleSplittingKernel::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
interfaceRenormalizationScaleFactor.rank(-2);
static Parameter<DipoleSplittingKernel,Energy> interfaceRenormalizationScaleFreeze
("RenormalizationScaleFreeze",
"The freezing scale for the renormalization scale.",
&DipoleSplittingKernel::theRenormalizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKernel,Energy> interfaceFactorizationScaleFreeze
("FactorizationScaleFreeze",
"The freezing scale for the factorization scale.",
&DipoleSplittingKernel::theFactorizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Switch<DipoleSplittingKernel,bool> interfaceVirtualitySplittingScale
("VirtualitySplittingScale",
"Use the virtuality as the splitting scale.",
&DipoleSplittingKernel::theVirtualitySplittingScale, false, false, false);
static SwitchOption interfaceVirtualitySplittingScaleYes
(interfaceVirtualitySplittingScale,
"Yes",
"Use vrituality.",
true);
static SwitchOption interfaceVirtualitySplittingScaleNo
(interfaceVirtualitySplittingScale,
"No",
"Use transverse momentum.",
false);
}
diff --git a/DipoleShower/Kernels/DipoleSplittingKernel.h b/DipoleShower/Kernels/DipoleSplittingKernel.h
--- a/DipoleShower/Kernels/DipoleSplittingKernel.h
+++ b/DipoleShower/Kernels/DipoleSplittingKernel.h
@@ -1,507 +1,453 @@
// -*- C++ -*-
//
// DipoleSplittingKernel.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 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_DipoleSplittingKernel_H
#define HERWIG_DipoleSplittingKernel_H
//
// This is the declaration of the DipoleSplittingKernel class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/StandardModel/AlphaSBase.h"
#include "ThePEG/PDF/PDF.h"
#include "Herwig/DipoleShower/Utility/PDFRatio.h"
#include "Herwig/DipoleShower/Base/DipoleSplittingInfo.h"
#include "Herwig/DipoleShower/Kinematics/DipoleSplittingKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer
*
* \brief DipoleSplittingKernel is the base class for all kernels
* used within the dipole shower.
*
* @see \ref DipoleSplittingKernelInterfaces "The interfaces"
* defined for DipoleSplittingKernel.
*/
class DipoleSplittingKernel: public HandlerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DipoleSplittingKernel();
/**
* The destructor.
*/
virtual ~DipoleSplittingKernel();
//@}
public:
/**
* Return the alpha_s to be used
*/
Ptr<AlphaSBase>::tptr alphaS() const { return theAlphaS; }
/**
* Set the alpha_s to be used
*/
void alphaS(Ptr<AlphaSBase>::tptr ap) { theAlphaS = ap; }
/**
* Return the splitting kinematics object
*/
Ptr<DipoleSplittingKinematics>::tptr splittingKinematics() const { return theSplittingKinematics; }
/**
* Return the mc check object
*/
Ptr<DipoleMCCheck>::ptr mcCheck() const { return theMCCheck; }
/**
* Set the splitting kinematics object
*/
void splittingKinematics(Ptr<DipoleSplittingKinematics>::tptr sp) { theSplittingKinematics = sp; }
/**
* Return the PDFRatio object
*/
Ptr<PDFRatio>::tptr pdfRatio() const { return thePDFRatio; }
/**
* Set the PDFRatio object
*/
void pdfRatio(Ptr<PDFRatio>::tptr sp) { thePDFRatio = sp; }
/**
* Return the number of additional parameter
* random variables needed to evaluate this kernel
* except the momentum fractions of incoming partons.
* These will be accessible through the
* lastSplittingParameters() container of the splitting
* info object.
*/
virtual int nDimAdditional() const { return 0; }
/**
* Set the freezing value for the renormalization scale
*/
void renormalizationScaleFreeze(Energy s) { theRenormalizationScaleFreeze = s; }
/**
* Set the freezing value for the factorization scale
*/
void factorizationScaleFreeze(Energy s) { theFactorizationScaleFreeze = s; }
/**
* Get the freezing value for the renormalization scale
*/
Energy renormalizationScaleFreeze() const { return theRenormalizationScaleFreeze; }
/**
* Get the freezing value for the factorization scale
*/
Energy factorizationScaleFreeze() const { return theFactorizationScaleFreeze; }
public:
/**
* Return true, if this splitting kernel
* applies to the given dipole index.
*/
virtual bool canHandle(const DipoleIndex&) const = 0;
/**
* Return true, if this splitting kernel is
* the same for the given index a, as the given
* splitting kernel for index b.
*/
virtual bool canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const = 0;
/**
* Return the emitter data after splitting, given
* a dipole index.
*/
virtual tcPDPtr emitter(const DipoleIndex&) const = 0;
/**
* Return the emission data after splitting, given
* a dipole index.
*/
virtual tcPDPtr emission(const DipoleIndex&) const = 0;
/**
* Return the spectator data after splitting, given
* a dipole index.
*/
virtual tcPDPtr spectator(const DipoleIndex&) const = 0;
/**
* Return the flavour produced, if this cannot
* be determined from the dipole.
*/
PDPtr flavour() const { return theFlavour; }
/**
* Return true, if this splitting kernel is supposed to work in a
* strict large-N limit, i.e. replacing C_F by C_A/2
*/
bool strictLargeN() const { return theStrictLargeN; }
public:
/**
* Inform this splitting kernel, that it is being
* presampled until a call to stopPresampling
*/
virtual void startPresampling(const DipoleIndex&) {}
/**
* Inform this splitting kernel, that it is not being
* presampled until a call to startPresampling
*/
virtual void stopPresampling(const DipoleIndex&) {}
/**
* Return the number of points to presample this
* splitting generator.
*/
unsigned long presamplingPoints() const { return thePresamplingPoints; }
/**
* Return the maximum number of trials
* to generate a splitting.
*/
unsigned long maxtry() const { return theMaxtry; }
/**
* Return the number of accepted points after which the grid should
* be frozen
*/
unsigned long freezeGrid() const { return theFreezeGrid; }
/**
* Set the number of accepted points after which the grid should
* be frozen
*/
void freezeGrid(unsigned long n) { theFreezeGrid = n; }
/**
* Evaluate this splitting kernel for the given
* dipole splitting.
*/
virtual double evaluate(const DipoleSplittingInfo&) const = 0;
/**
* Update the variations vector at the given splitting using the indicated
* kernel and overestimate values.
*/
virtual void accept(const DipoleSplittingInfo&, double, double, map<string,double>&) const;
/**
* Update the variations vector at the given splitting using the indicated
* kernel and overestimate values.
*/
virtual void veto(const DipoleSplittingInfo&, double, double, map<string,double>&) const;
/**
* Return true, if this kernel is capable of
* delivering an overestimate to the kernel, and
* of inverting the integral over the overestimate
* w.r.t. the phasepsace provided by the given
* DipoleSplittingInfo object.
*/
virtual bool haveOverestimate(const DipoleSplittingInfo&) const { return false; }
/**
* Return the overestimate to this splitting kernel
* for the given dipole splitting.
*/
virtual double overestimate(const DipoleSplittingInfo&) const { return -1.; }
/**
* Invert the integral over the overestimate
* w.r.t. the phasepsace provided by the given
* DipoleSplittingInfo object to equal
* the given value.
*/
virtual double invertOverestimateIntegral(const DipoleSplittingInfo&, double) const { return -1.; }
public:
/**
* 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; }
- /**
- * Vary the factorization scale by the indicated values and using the
- * indicated variation labels.
- */
- map<string,double>& renormalizationScaleVariations() {
- return theRenormalizationScaleVariations;
- }
-
- /**
- * Vary the factorization scale by the indicated values and using the
- * indicated variation labels.
- */
- const map<string,double>& renormalizationScaleVariations() const {
- return theRenormalizationScaleVariations;
- }
-
- /**
- * Vary the renormalization scale by the indicated values and using the
- * indicated variation labels.
- */
- map<string,double>& factorizationScaleVariations() {
- return theFactorizationScaleVariations;
- }
-
- /**
- * Vary the renormalization scale by the indicated values and using the
- * indicated variation labels.
- */
- const map<string,double>& factorizationScaleVariations() const {
- return theFactorizationScaleVariations;
- }
-
protected:
/**
* Return the common factor of (alphas/2pi)*(pdf ratio)
*/
double alphaPDF(const DipoleSplittingInfo&,
- Energy optScale = ZERO) const;
+ Energy optScale = ZERO,
+ double rScaleFactor = 1.0,
+ double fScaleFactor = 1.0) const;
/**
* Return true, if the virtuality of the splitting should be used as the
* argument of alphas rather than the pt
*/
bool virtualitySplittingScale() const { return theVirtualitySplittingScale; }
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();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The alpha_s to be used.
*/
Ptr<AlphaSBase>::ptr theAlphaS;
/**
* An optional 'colour screening' scale
* for alternative intrinsic pt generation.
*/
Energy theScreeningScale;
/**
* The splitting kinematics to be used.
*/
Ptr<DipoleSplittingKinematics>::ptr theSplittingKinematics;
/**
* An optional PDF ratio object to be used
* when evaluating this kernel.
*/
Ptr<PDFRatio>::ptr thePDFRatio;
/**
* The number of points to presample this
* splitting generator.
*/
unsigned long thePresamplingPoints;
/**
* The maximum number of trials
* to generate a splitting.
*/
unsigned long theMaxtry;
/**
* Return the number of accepted points after which the grid should
* be frozen
*/
unsigned long theFreezeGrid;
/**
* The flavour produced, if this cannot
* be determined from the dipole.
*/
PDPtr theFlavour;
/**
* Pointer to a check histogram object
*/
Ptr<DipoleMCCheck>::ptr theMCCheck;
/**
* True, if this splitting kernel is supposed to work in a
* strict large-N limit, i.e. replacing C_F by C_A/2
*/
bool theStrictLargeN;
/**
* The factorization scale factor.
*/
double theFactorizationScaleFactor;
/**
* The renormalization scale factor.
*/
double theRenormalizationScaleFactor;
/**
* A freezing value for the renormalization scale
*/
Energy theRenormalizationScaleFreeze;
/**
* A freezing value for the factorization scale
*/
Energy theFactorizationScaleFreeze;
/**
* True, if the virtuality of the splitting should be used as the
* argument of alphas rather than the pt
*/
bool theVirtualitySplittingScale;
- /**
- * Vary the renormalization scale by the indicated values and using the
- * indicated variation labels.
- */
- map<string,double> theRenormalizationScaleVariations;
-
- /**
- * Command to vary the renormalization scale by the indicated values and using
- * the indicated variation labels.
- */
- string doRenormalizationScaleVariations(string);
-
- /**
- * Vary the factorization scale by the indicated values and using the
- * indicated variation labels.
- */
- map<string,double> theFactorizationScaleVariations;
-
- /**
- * Command to vary the factorization scale by the indicated values and using
- * the indicated variation labels.
- */
- string doFactorizationScaleVariations(string);
-
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is an abstract class with persistent data.
*/
static AbstractClassDescription<DipoleSplittingKernel> initDipoleSplittingKernel;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleSplittingKernel & operator=(const DipoleSplittingKernel &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleSplittingKernel. */
template <>
struct BaseClassTrait<Herwig::DipoleSplittingKernel,1> {
/** Typedef of the first base class of DipoleSplittingKernel. */
typedef HandlerBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleSplittingKernel class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleSplittingKernel>
: public ClassTraitsBase<Herwig::DipoleSplittingKernel> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleSplittingKernel"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleSplittingKernel is implemented. It may also include several, space-separated,
* libraries if the class DipoleSplittingKernel depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_DipoleSplittingKernel_H */
diff --git a/Shower/ShowerHandler.cc b/Shower/ShowerHandler.cc
--- a/Shower/ShowerHandler.cc
+++ b/Shower/ShowerHandler.cc
@@ -1,845 +1,860 @@
// -*- C++ -*-
//
// ShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ShowerHandler class.
//
#include "ShowerHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Interface/Command.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/PDF/PartonBinInstance.h"
#include "Herwig/PDT/StandardMatchers.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Utilities/Throw.h"
#include "Herwig/Shower/Base/Evolver.h"
#include "Herwig/Shower/Base/ShowerParticle.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig/Utilities/EnumParticles.h"
#include "Herwig/PDF/MPIPDF.h"
#include "Herwig/PDF/MinBiasPDF.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "Herwig/Shower/Base/ShowerTree.h"
#include "Herwig/Shower/Base/HardTree.h"
#include "Herwig/Shower/Base/KinematicsReconstructor.h"
#include "Herwig/Shower/Base/PartnerFinder.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include <cassert>
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
DescribeClass<ShowerHandler,CascadeHandler>
describeShowerHandler ("Herwig::ShowerHandler","HwShower.so");
ShowerHandler::~ShowerHandler() {}
ShowerHandler * ShowerHandler::currentHandler_ = 0;
void ShowerHandler::doinit() {
CascadeHandler::doinit();
// copy particles to decay before showering from input vector to the
// set used in the simulation
if ( particlesDecayInShower_.empty() )
particlesDecayInShower_.insert(inputparticlesDecayInShower_.begin(),
inputparticlesDecayInShower_.end());
ShowerTree::_decayInShower = particlesDecayInShower_;
ShowerTree::_vmin2 = vMin_;
ShowerTree::_spaceTime = includeSpaceTime_;
}
IBPtr ShowerHandler::clone() const {
return new_ptr(*this);
}
IBPtr ShowerHandler::fullclone() const {
return new_ptr(*this);
}
ShowerHandler::ShowerHandler() :
pdfFreezingScale_(2.5*GeV),
maxtry_(10),maxtryMPI_(10),maxtryDP_(10),
includeSpaceTime_(false), vMin_(0.1*GeV2), subProcess_(),
factorizationScaleFactor_(1.0),
renormalizationScaleFactor_(1.0),
hardScaleFactor_(1.0), scaleFactorOption_(0),
restrictPhasespace_(true), maxPtIsMuF_(false),
splitHardProcess_(true) {
inputparticlesDecayInShower_.push_back( 6 ); // top
inputparticlesDecayInShower_.push_back( 23 ); // Z0
inputparticlesDecayInShower_.push_back( 24 ); // W+/-
inputparticlesDecayInShower_.push_back( 25 ); // h0
}
void ShowerHandler::doinitrun(){
CascadeHandler::doinitrun();
//can't use isMPIOn here, because the EventHandler is not set at that stage
if(MPIHandler_){
MPIHandler_->initialize();
if(MPIHandler_->softInt())
remDec_->initSoftInteractions(MPIHandler_->Ptmin(), MPIHandler_->beta());
}
ShowerTree::_decayInShower = particlesDecayInShower_;
ShowerTree::_vmin2 = vMin_;
ShowerTree::_spaceTime = includeSpaceTime_;
}
void ShowerHandler::dofinish(){
CascadeHandler::dofinish();
if(MPIHandler_) MPIHandler_->finalize();
}
void ShowerHandler::persistentOutput(PersistentOStream & os) const {
os << evolver_ << remDec_ << ounit(pdfFreezingScale_,GeV) << maxtry_
<< maxtryMPI_ << maxtryDP_ << inputparticlesDecayInShower_
<< particlesDecayInShower_ << MPIHandler_ << PDFA_ << PDFB_
<< PDFARemnant_ << PDFBRemnant_
<< includeSpaceTime_ << ounit(vMin_,GeV2)
<< factorizationScaleFactor_ << renormalizationScaleFactor_
<< hardScaleFactor_ << scaleFactorOption_
<< restrictPhasespace_ << maxPtIsMuF_ << hardScaleProfile_
- << splitHardProcess_
- << renormalizationScaleVariations_
- << factorizationScaleVariations_;
+ << splitHardProcess_ << showerVariations_;
}
void ShowerHandler::persistentInput(PersistentIStream & is, int) {
is >> evolver_ >> remDec_ >> iunit(pdfFreezingScale_,GeV) >> maxtry_
>> maxtryMPI_ >> maxtryDP_ >> inputparticlesDecayInShower_
>> particlesDecayInShower_ >> MPIHandler_ >> PDFA_ >> PDFB_
>> PDFARemnant_ >> PDFBRemnant_
>> includeSpaceTime_ >> iunit(vMin_,GeV2)
>> factorizationScaleFactor_ >> renormalizationScaleFactor_
>> hardScaleFactor_ >> scaleFactorOption_
>> restrictPhasespace_ >> maxPtIsMuF_ >> hardScaleProfile_
- >> splitHardProcess_
- >> renormalizationScaleVariations_
- >> factorizationScaleVariations_;
+ >> splitHardProcess_ >> showerVariations_;
}
void ShowerHandler::Init() {
static ClassDocumentation<ShowerHandler> documentation
("Main driver class for the showering.",
"The Shower evolution was performed using an algorithm described in "
"\\cite{Marchesini:1983bm,Marchesini:1987cf,Gieseke:2003rz,Bahr:2008pv}.",
"%\\cite{Marchesini:1983bm}\n"
"\\bibitem{Marchesini:1983bm}\n"
" G.~Marchesini and B.~R.~Webber,\n"
" ``Simulation Of QCD Jets Including Soft Gluon Interference,''\n"
" Nucl.\\ Phys.\\ B {\\bf 238}, 1 (1984).\n"
" %%CITATION = NUPHA,B238,1;%%\n"
"%\\cite{Marchesini:1987cf}\n"
"\\bibitem{Marchesini:1987cf}\n"
" G.~Marchesini and B.~R.~Webber,\n"
" ``Monte Carlo Simulation of General Hard Processes with Coherent QCD\n"
" Radiation,''\n"
" Nucl.\\ Phys.\\ B {\\bf 310}, 461 (1988).\n"
" %%CITATION = NUPHA,B310,461;%%\n"
"%\\cite{Gieseke:2003rz}\n"
"\\bibitem{Gieseke:2003rz}\n"
" S.~Gieseke, P.~Stephens and B.~Webber,\n"
" ``New formalism for QCD parton showers,''\n"
" JHEP {\\bf 0312}, 045 (2003)\n"
" [arXiv:hep-ph/0310083].\n"
" %%CITATION = JHEPA,0312,045;%%\n"
);
static Reference<ShowerHandler,Evolver>
interfaceEvolver("Evolver",
"A reference to the Evolver object",
&Herwig::ShowerHandler::evolver_,
false, false, true, false);
static Reference<ShowerHandler,HwRemDecayer>
interfaceRemDecayer("RemDecayer",
"A reference to the Remnant Decayer object",
&Herwig::ShowerHandler::remDec_,
false, false, true, false);
static Parameter<ShowerHandler,Energy> interfacePDFFreezingScale
("PDFFreezingScale",
"The PDF freezing scale",
&ShowerHandler::pdfFreezingScale_, GeV, 2.5*GeV, 2.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<ShowerHandler,unsigned int> interfaceMaxTry
("MaxTry",
"The maximum number of attempts for the main showering loop",
&ShowerHandler::maxtry_, 10, 1, 100,
false, false, Interface::limited);
static Parameter<ShowerHandler,unsigned int> interfaceMaxTryMPI
("MaxTryMPI",
"The maximum number of regeneration attempts for an additional scattering",
&ShowerHandler::maxtryMPI_, 10, 0, 100,
false, false, Interface::limited);
static Parameter<ShowerHandler,unsigned int> interfaceMaxTryDP
("MaxTryDP",
"The maximum number of regeneration attempts for an additional hard scattering",
&ShowerHandler::maxtryDP_, 10, 0, 100,
false, false, Interface::limited);
static ParVector<ShowerHandler,long> interfaceDecayInShower
("DecayInShower",
"PDG codes of the particles to be decayed in the shower",
&ShowerHandler::inputparticlesDecayInShower_, -1, 0l, -10000000l, 10000000l,
false, false, Interface::limited);
static Reference<ShowerHandler,UEBase> interfaceMPIHandler
("MPIHandler",
"The object that administers all additional scatterings.",
&ShowerHandler::MPIHandler_, false, false, true, true);
static Reference<ShowerHandler,PDFBase> interfacePDFA
("PDFA",
"The PDF for beam particle A. Overrides the particle's own PDF setting."
"By default used for both the shower and forced splitting in the remnant",
&ShowerHandler::PDFA_, false, false, true, true, false);
static Reference<ShowerHandler,PDFBase> interfacePDFB
("PDFB",
"The PDF for beam particle B. Overrides the particle's own PDF setting."
"By default used for both the shower and forced splitting in the remnant",
&ShowerHandler::PDFB_, false, false, true, true, false);
static Reference<ShowerHandler,PDFBase> interfacePDFARemnant
("PDFARemnant",
"The PDF for beam particle A used to generate forced splittings of the remnant."
" This overrides both the particle's own PDF setting and the value set by PDFA if used.",
&ShowerHandler::PDFARemnant_, false, false, true, true, false);
static Reference<ShowerHandler,PDFBase> interfacePDFBRemnant
("PDFBRemnant",
"The PDF for beam particle B used to generate forced splittings of the remnant."
" This overrides both the particle's own PDF setting and the value set by PDFB if used.",
&ShowerHandler::PDFBRemnant_, false, false, true, true, false);
static Switch<ShowerHandler,bool> interfaceIncludeSpaceTime
("IncludeSpaceTime",
"Whether to include the model for the calculation of space-time distances",
&ShowerHandler::includeSpaceTime_, false, false, false);
static SwitchOption interfaceIncludeSpaceTimeYes
(interfaceIncludeSpaceTime,
"Yes",
"Include the model",
true);
static SwitchOption interfaceIncludeSpaceTimeNo
(interfaceIncludeSpaceTime,
"No",
"Only include the displacement from the particle-s lifetime for decaying particles",
false);
static Parameter<ShowerHandler,Energy2> interfaceMinimumVirtuality
("MinimumVirtuality",
"The minimum virtuality for the space-time model",
&ShowerHandler::vMin_, GeV2, 0.1*GeV2, 0.0*GeV2, 1000.0*GeV2,
false, false, Interface::limited);
static Parameter<ShowerHandler,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&ShowerHandler::factorizationScaleFactor_, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<ShowerHandler,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&ShowerHandler::renormalizationScaleFactor_, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<ShowerHandler,double> interfaceHardScaleFactor
("HardScaleFactor",
"The hard scale factor.",
&ShowerHandler::hardScaleFactor_, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Switch<ShowerHandler,int> interfaceScaleFactorOption
("ScaleFactorOption",
"Where to apply scale factors.",
&ShowerHandler::scaleFactorOption_, 0, false, false);
static SwitchOption interfaceScaleFactorOptionAll
(interfaceScaleFactorOption,
"All",
"Apply in first and secondary scatterings.",
0);
static SwitchOption interfaceScaleFactorOptionHard
(interfaceScaleFactorOption,
"Hard",
"Only apply for the hard process.",
1);
static SwitchOption interfaceScaleFactorOptionSecondary
(interfaceScaleFactorOption,
"Secondary",
"Only apply to secondary scatterings.",
2);
static Reference<ShowerHandler,HardScaleProfile> interfaceHardScaleProfile
("HardScaleProfile",
"The hard scale profile to use.",
&ShowerHandler::hardScaleProfile_, false, false, true, true, false);
static Switch<ShowerHandler,bool> interfaceMaxPtIsMuF
("MaxPtIsMuF",
"",
&ShowerHandler::maxPtIsMuF_, false, false, false);
static SwitchOption interfaceMaxPtIsMuFYes
(interfaceMaxPtIsMuF,
"Yes",
"",
true);
static SwitchOption interfaceMaxPtIsMuFNo
(interfaceMaxPtIsMuF,
"No",
"",
false);
static Switch<ShowerHandler,bool> interfaceRestrictPhasespace
("RestrictPhasespace",
"Switch on or off phasespace restrictions",
&ShowerHandler::restrictPhasespace_, true, false, false);
static SwitchOption interfaceRestrictPhasespaceOn
(interfaceRestrictPhasespace,
"On",
"Perform phasespace restrictions",
true);
static SwitchOption interfaceRestrictPhasespaceOff
(interfaceRestrictPhasespace,
"Off",
"Do not perform phasespace restrictions",
false);
static Switch<ShowerHandler,bool> interfaceSplitHardProcess
("SplitHardProcess",
"Whether or not to try and split the hard process into production and decay processes",
&ShowerHandler::splitHardProcess_, true, false, false);
static SwitchOption interfaceSplitHardProcessYes
(interfaceSplitHardProcess,
"Yes",
"Split the hard process",
true);
static SwitchOption interfaceSplitHardProcessNo
(interfaceSplitHardProcess,
"No",
"Don't split the hard process",
false);
+
+ static Command<ShowerHandler> interfaceAddVariation
+ ("AddVariation",
+ "Add a shower variation.",
+ &ShowerHandler::doAddVariation, false);
}
Energy ShowerHandler::hardScale() const {
return evolver_->hardScale();
}
void ShowerHandler::cascade() {
tcPDFPtr first = firstPDF().pdf();
tcPDFPtr second = secondPDF().pdf();
if ( PDFA_ ) first = PDFA_;
if ( PDFB_ ) second = PDFB_;
resetPDFs(make_pair(first,second));
if( ! rempdfs_.first)
rempdfs_.first = PDFARemnant_ ? PDFPtr(PDFARemnant_) : const_ptr_cast<PDFPtr>(first);
if( ! rempdfs_.second)
rempdfs_.second = PDFBRemnant_ ? PDFPtr(PDFBRemnant_) : const_ptr_cast<PDFPtr>(second);
// get the incoming partons
tPPair incomingPartons =
eventHandler()->currentCollision()->primarySubProcess()->incoming();
// and the parton bins
PBIPair incomingBins =
make_pair(lastExtractor()->partonBinInstance(incomingPartons.first),
lastExtractor()->partonBinInstance(incomingPartons.second));
// and the incoming hadrons
tPPair incomingHadrons =
eventHandler()->currentCollision()->incoming();
remDec_->setHadronContent(incomingHadrons);
// check if incoming hadron == incoming parton
// and get the incoming hadron if exists or parton otherwise
incoming_ = make_pair(incomingBins.first ?
incomingBins.first ->particle() : incomingPartons.first,
incomingBins.second ?
incomingBins.second->particle() : incomingPartons.second);
// check the collision is of the beam particles
// and if not boost collision to the right frame
// i.e. the hadron-hadron CMF of the collision
bool btotal(false);
LorentzRotation rtotal;
if(incoming_.first != incomingHadrons.first ||
incoming_.second != incomingHadrons.second ) {
btotal = true;
boostCollision(false);
}
// set the current ShowerHandler
currentHandler_ = this;
// first shower the hard process
useMe();
try {
SubProPtr sub = eventHandler()->currentCollision()->primarySubProcess();
incomingPartons = cascade(sub,lastXCombPtr());
}
catch(ShowerTriesVeto &veto){
throw Exception() << "Failed to generate the shower after "
<< veto.tries
<< " attempts in ShowerHandler::cascade()"
<< Exception::eventerror;
}
if(showerHardProcessVeto()) throw Veto();
// if a non-hadron collision return (both incoming non-hadronic)
if( ( !incomingBins.first||
!isResolvedHadron(incomingBins.first ->particle()))&&
( !incomingBins.second||
!isResolvedHadron(incomingBins.second->particle()))) {
// boost back to lab if needed
if(btotal) boostCollision(true);
// unset the current ShowerHandler
currentHandler_ = 0;
return;
}
// get the remnants for hadronic collision
pair<tRemPPtr,tRemPPtr> remnants(getRemnants(incomingBins));
// set the starting scale of the forced splitting to the PDF freezing scale
remDec_->initialize(remnants, incoming_, *currentStep(), pdfFreezingScale());
// do the first forcedSplitting
try {
remDec_->doSplit(incomingPartons, make_pair(rempdfs_.first,rempdfs_.second), true);
}
catch (ExtraScatterVeto) {
throw Exception() << "Remnant extraction failed in "
<< "ShowerHandler::cascade() from primary interaction"
<< Exception::eventerror;
}
// if no MPI return
if( !isMPIOn() ) {
remDec_->finalize();
// boost back to lab if needed
if(btotal) boostCollision(true);
// unset the current ShowerHandler
currentHandler_ = 0;
return;
}
// generate the multiple scatters use modified pdf's now:
setMPIPDFs();
// additional "hard" processes
unsigned int tries(0);
// This is the loop over additional hard scatters (most of the time
// only one, but who knows...)
for(unsigned int i=1; i <= getMPIHandler()->additionalHardProcs(); i++){
//counter for regeneration
unsigned int multSecond = 0;
// generate the additional scatters
while( multSecond < getMPIHandler()->multiplicity(i) ) {
// generate the hard scatter
tStdXCombPtr lastXC = getMPIHandler()->generate(i);
SubProPtr sub = lastXC->construct();
// add to the Step
newStep()->addSubProcess(sub);
// increment the counters
tries++;
multSecond++;
if(tries == maxtryDP_)
throw Exception() << "Failed to establish the requested number "
<< "of additional hard processes. If this error "
<< "occurs often, your selection of additional "
<< "scatter is probably unphysical"
<< Exception::eventerror;
// Generate the shower. If not possible veto the event
try {
incomingPartons = cascade(sub,lastXC);
}
catch(ShowerTriesVeto &veto){
throw Exception() << "Failed to generate the shower of "
<< "a secondary hard process after "
<< veto.tries
<< " attempts in Evolver::showerHardProcess()"
<< Exception::eventerror;
}
try {
// do the forcedSplitting
remDec_->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false);
}
catch(ExtraScatterVeto){
//remove all particles associated with the subprocess
newStep()->removeParticle(incomingPartons.first);
newStep()->removeParticle(incomingPartons.second);
//remove the subprocess from the list
newStep()->removeSubProcess(sub);
//regenerate the scattering
multSecond--;
continue;
}
// connect with the remnants but don't set Remnant colour,
// because that causes problems due to the multiple colour lines.
if ( !remnants.first ->extract(incomingPartons.first , false) ||
!remnants.second->extract(incomingPartons.second, false) )
throw Exception() << "Remnant extraction failed in "
<< "ShowerHandler::cascade() for additional scatter"
<< Exception::runerror;
}
}
// the underlying event processes
unsigned int ptveto(1), veto(0);
unsigned int max(getMPIHandler()->multiplicity());
for(unsigned int i=0; i<max; i++) {
// check how often this scattering has been regenerated
if(veto > maxtryMPI_) break;
//generate PSpoint
tStdXCombPtr lastXC = getMPIHandler()->generate();
SubProPtr sub = lastXC->construct();
//If Algorithm=1 additional scatters of the signal type
// with pt > ptmin have to be vetoed
//with probability 1/(m+1), where m is the number of occurances in this event
if( getMPIHandler()->Algorithm() == 1 ){
//get the pT
Energy pt = sub->outgoing().front()->momentum().perp();
if(pt > getMPIHandler()->PtForVeto() && UseRandom::rnd() < 1./(ptveto+1) ){
ptveto++;
i--;
continue;
}
}
// add to the SubProcess to the step
newStep()->addSubProcess(sub);
// Run the Shower. If not possible veto the scattering
try {
incomingPartons = cascade(sub,lastXC);
}
// discard this extra scattering, but try the next one
catch(ShowerTriesVeto) {
newStep()->removeSubProcess(sub);
//regenerate the scattering
veto++;
i--;
continue;
}
try{
//do the forcedSplitting
remDec_->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false);
}
catch (ExtraScatterVeto) {
//remove all particles associated with the subprocess
newStep()->removeParticle(incomingPartons.first);
newStep()->removeParticle(incomingPartons.second);
//remove the subprocess from the list
newStep()->removeSubProcess(sub);
//regenerate the scattering
veto++;
i--;
continue;
}
//connect with the remnants but don't set Remnant colour,
//because that causes problems due to the multiple colour lines.
if ( !remnants.first ->extract(incomingPartons.first , false) ||
!remnants.second->extract(incomingPartons.second, false) )
throw Exception() << "Remnant extraction failed in "
<< "ShowerHandler::cascade() for MPI hard scattering"
<< Exception::runerror;
//reset veto counter
veto = 0;
}
// finalize the remnants
remDec_->finalize(getMPIHandler()->colourDisrupt(),
getMPIHandler()->softMultiplicity());
// boost back to lab if needed
if(btotal) boostCollision(true);
// unset the current ShowerHandler
currentHandler_ = 0;
getMPIHandler()->clean();
}
void ShowerHandler::fillEventRecord() {
// create a new step
StepPtr pstep = newStep();
assert(!done_.empty());
assert(done_[0]->isHard());
// insert the steps
for(unsigned int ix=0;ix<done_.size();++ix) {
done_[ix]->fillEventRecord(pstep,
evolver_->isISRadiationON(),
evolver_->isFSRadiationON());
}
}
void ShowerHandler::prepareCascade(tSubProPtr sub) {
current_ = currentStep();
subProcess_ = sub;
+ for ( map<string,double>::iterator w = currentWeights_.begin();
+ w != currentWeights_.end(); ++w ) {
+ w->second = 1.0;
+ }
}
+void ShowerHandler::combineWeights() {
+ tEventPtr event = eventHandler()->currentEvent();
+ for ( map<string,double>::const_iterator w =
+ currentWeights_.begin(); w != currentWeights_.end(); ++w ) {
+ map<string,double>::iterator ew = event->optionalWeights().find(w->first);
+ if ( ew != event->optionalWeights().end() )
+ ew->second *= w->second;
+ else
+ event->optionalWeights()[w->first] = w->second;
+ }
+}
+
+void ShowerHandler::ShowerVariation::fromInFile(const string&) {
+ // need agree on a syntax to steer this
+ assert(false && "implementation missing");
+}
+
+void ShowerHandler::ShowerVariation::put(PersistentOStream& os) const {
+ os << renormalizationScaleFactor << factorizationScaleFactor
+ << firstInteraction << secondaryInteractions;
+}
+
+void ShowerHandler::ShowerVariation::get(PersistentIStream& is) {
+ is >> renormalizationScaleFactor >> factorizationScaleFactor
+ >> firstInteraction >> secondaryInteractions;
+}
+
+string ShowerHandler::doAddVariation(string) {
+ // need agree on a syntax to steer this
+ assert(false && "implementation missing");
+}
+
tPPair ShowerHandler::cascade(tSubProPtr sub,
XCPtr xcomb) {
prepareCascade(sub);
// set the scale variation factors; needs to go after prepareCascade
// to trigger possible different variations for hard and secondary
// scatters
evolver()->renormalizationScaleFactor(renormalizationScaleFactor());
evolver()->factorizationScaleFactor(factorizationScaleFactor());
evolver()->restrictPhasespace(restrictPhasespace());
evolver()->hardScaleIsMuF(hardScaleIsMuF());
// start of the try block for the whole showering process
unsigned int countFailures=0;
while (countFailures<maxtry_) {
try {
decay_.clear();
done_.clear();
ShowerTree::constructTrees(currentSubProcess(),hard_,decay_,
firstInteraction() ? tagged() :
tPVector(currentSubProcess()->outgoing().begin(),
currentSubProcess()->outgoing().end()),
splitHardProcess_);
// if no hard process
if(!hard_) throw Exception() << "Shower starting with a decay"
<< "is not implemented"
<< Exception::runerror;
// perform the shower for the hard process
evolver_->showerHardProcess(hard_,xcomb);
done_.push_back(hard_);
hard_->updateAfterShower(decay_);
// if no decaying particles to shower break out of the loop
if(decay_.empty()) break;
// shower the decay products
while(!decay_.empty()) {
// find particle whose production process has been showered
ShowerDecayMap::iterator dit = decay_.begin();
while(!dit->second->parent()->hasShowered() && dit!=decay_.end()) ++dit;
assert(dit!=decay_.end());
// get the particle
ShowerTreePtr decayingTree = dit->second;
// remove it from the multimap
decay_.erase(dit);
// make sure the particle has been decayed
decayingTree->decay(decay_);
// now shower the decay
evolver_->showerDecay(decayingTree);
done_.push_back(decayingTree);
decayingTree->updateAfterShower(decay_);
}
// suceeded break out of the loop
break;
}
catch (KinematicsReconstructionVeto) {
++countFailures;
}
}
// if loop exited because of too many tries, throw event away
if (countFailures >= maxtry_) {
hard_=ShowerTreePtr();
decay_.clear();
done_.clear();
throw Exception() << "Too many tries for main while loop "
<< "in ShowerHandler::cascade()."
<< Exception::eventerror;
}
//enter the particles in the event record
fillEventRecord();
// clear storage
hard_=ShowerTreePtr();
decay_.clear();
done_.clear();
// non hadronic case return
if (!isResolvedHadron(incoming_.first ) &&
!isResolvedHadron(incoming_.second) )
return incoming_;
// remake the remnants (needs to be after the colours are sorted
// out in the insertion into the event record)
if ( firstInteraction() ) return remakeRemnant(sub->incoming());
//Return the new pair of incoming partons. remakeRemnant is not
//necessary here, because the secondary interactions are not yet
//connected to the remnants.
return make_pair(findFirstParton(sub->incoming().first ),
findFirstParton(sub->incoming().second));
}
ShowerHandler::RemPair
ShowerHandler::getRemnants(PBIPair incomingBins) {
RemPair remnants;
// first beam particle
if(incomingBins.first&&!incomingBins.first->remnants().empty()) {
remnants.first =
dynamic_ptr_cast<tRemPPtr>(incomingBins.first->remnants()[0] );
if(remnants.first) {
ParticleVector children=remnants.first->children();
for(unsigned int ix=0;ix<children.size();++ix) {
if(children[ix]->dataPtr()==remnants.first->dataPtr())
remnants.first = dynamic_ptr_cast<RemPPtr>(children[ix]);
}
//remove existing colour lines from the remnants
if(remnants.first->colourLine())
remnants.first->colourLine()->removeColoured(remnants.first);
if(remnants.first->antiColourLine())
remnants.first->antiColourLine()->removeAntiColoured(remnants.first);
}
}
// seconnd beam particle
if(incomingBins.second&&!incomingBins. second->remnants().empty()) {
remnants.second =
dynamic_ptr_cast<tRemPPtr>(incomingBins.second->remnants()[0] );
if(remnants.second) {
ParticleVector children=remnants.second->children();
for(unsigned int ix=0;ix<children.size();++ix) {
if(children[ix]->dataPtr()==remnants.second->dataPtr())
remnants.second = dynamic_ptr_cast<RemPPtr>(children[ix]);
}
//remove existing colour lines from the remnants
if(remnants.second->colourLine())
remnants.second->colourLine()->removeColoured(remnants.second);
if(remnants.second->antiColourLine())
remnants.second->antiColourLine()->removeAntiColoured(remnants.second);
}
}
assert(remnants.first || remnants.second);
return remnants;
}
tPPair ShowerHandler::remakeRemnant(tPPair oldp){
// get the parton extractor
PartonExtractor & pex = *lastExtractor();
// get the new partons
tPPair newp = make_pair(findFirstParton(oldp.first ),
findFirstParton(oldp.second));
// if the same do nothing
if(newp == oldp) return oldp;
// Creates the new remnants and returns the new PartonBinInstances
// ATTENTION Broken here for very strange configuration
PBIPair newbins = pex.newRemnants(oldp, newp, newStep());
newStep()->addIntermediate(newp.first);
newStep()->addIntermediate(newp.second);
// return the new partons
return newp;
}
PPtr ShowerHandler::findFirstParton(tPPtr seed) const{
if(seed->parents().empty()) return seed;
tPPtr parent = seed->parents()[0];
//if no parent there this is a loose end which will
//be connected to the remnant soon.
if(!parent || parent == incoming_.first ||
parent == incoming_.second ) return seed;
else return findFirstParton(parent);
}
namespace {
void addChildren(tPPtr in,set<tPPtr> & particles) {
particles.insert(in);
for(unsigned int ix=0;ix<in->children().size();++ix)
addChildren(in->children()[ix],particles);
}
}
void ShowerHandler::boostCollision(bool boost) {
// calculate boost from lab to rest
if(!boost) {
Lorentz5Momentum ptotal=incoming_.first ->momentum()+incoming_.second->momentum();
boost_ = LorentzRotation(-ptotal.boostVector());
Axis axis((boost_*incoming_.first ->momentum()).vect().unit());
if(axis.perp2()>0.) {
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
boost_.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
}
}
// first call performs the boost and second inverse
// get the particles to be boosted
set<tPPtr> particles;
addChildren(incoming_.first,particles);
addChildren(incoming_.second,particles);
// apply the boost
for(set<tPPtr>::const_iterator cit=particles.begin();
cit!=particles.end();++cit) {
(*cit)->transform(boost_);
}
if(!boost) boost_.invert();
}
void ShowerHandler::setMPIPDFs() {
if ( !mpipdfs_.first ) {
// first have to check for MinBiasPDF
tcMinBiasPDFPtr first = dynamic_ptr_cast<tcMinBiasPDFPtr>(firstPDF().pdf());
if(first)
mpipdfs_.first = new_ptr(MPIPDF(first->originalPDF()));
else
mpipdfs_.first = new_ptr(MPIPDF(firstPDF().pdf()));
}
if ( !mpipdfs_.second ) {
tcMinBiasPDFPtr second = dynamic_ptr_cast<tcMinBiasPDFPtr>(secondPDF().pdf());
if(second)
mpipdfs_.second = new_ptr(MPIPDF(second->originalPDF()));
else
mpipdfs_.second = new_ptr(MPIPDF(secondPDF().pdf()));
}
if( !remmpipdfs_.first ) {
tcMinBiasPDFPtr first = dynamic_ptr_cast<tcMinBiasPDFPtr>(rempdfs_.first);
if(first)
remmpipdfs_.first = new_ptr(MPIPDF(first->originalPDF()));
else
remmpipdfs_.first = new_ptr(MPIPDF(rempdfs_.first));
}
if( !remmpipdfs_.second ) {
tcMinBiasPDFPtr second = dynamic_ptr_cast<tcMinBiasPDFPtr>(rempdfs_.second);
if(second)
remmpipdfs_.second = new_ptr(MPIPDF(second->originalPDF()));
else
remmpipdfs_.second = new_ptr(MPIPDF(rempdfs_.second));
}
// reset the PDFs stored in the base class
resetPDFs(mpipdfs_);
}
bool ShowerHandler::isResolvedHadron(tPPtr particle) {
if(!HadronMatcher::Check(particle->data())) return false;
for(unsigned int ix=0;ix<particle->children().size();++ix) {
if(particle->children()[ix]->id()==ParticleID::Remnant) return true;
}
return false;
}
HardTreePtr ShowerHandler::generateCKKW(ShowerTreePtr ) const {
return HardTreePtr();
}
-string ShowerHandler::parseVariations(string in, map<string,double>& vars) {
- if ( in.empty() )
- return "nothing specified";
- istringstream read(in);
- string id; double v;
- while ( read ) {
- read >> id;
- if ( !read )
- return "no value specified for tag '" + id + "'";
- read >> v;
- addToVariations(id, v, vars);
- }
- return "";
-}
-
-string ShowerHandler::doRenormalizationScaleVariations(string in) {
- return parseVariations(in,renormalizationScaleVariations_);
-}
-
-string ShowerHandler::doFactorizationScaleVariations(string in) {
- return parseVariations(in,factorizationScaleVariations_);
-}
-
diff --git a/Shower/ShowerHandler.h b/Shower/ShowerHandler.h
--- a/Shower/ShowerHandler.h
+++ b/Shower/ShowerHandler.h
@@ -1,659 +1,676 @@
// -*- C++ -*-
//
// ShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_ShowerHandler_H
#define HERWIG_ShowerHandler_H
//
// This is the declaration of the ShowerHandler class.
//
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "Herwig/Shower/UEBase.h"
#include "Herwig/Shower/Base/Evolver.fh"
#include "Herwig/Shower/Base/ShowerParticle.fh"
#include "Herwig/Shower/Base/ShowerTree.fh"
#include "Herwig/Shower/Base/HardTree.fh"
#include "Herwig/PDF/HwRemDecayer.fh"
#include "ThePEG/EventRecord/RemnantParticle.fh"
#include "ShowerHandler.fh"
#include "Herwig/MatrixElement/Matchbox/Matching/HardScaleProfile.h"
namespace Herwig {
/**
* Typedef for the ShowerTree for the decays
*/
typedef multimap<Energy,ShowerTreePtr,std::greater<Energy> > ShowerDecayMap;
using namespace ThePEG;
/** \ingroup Shower
*
* This class is the main driver of the shower: it is responsible for
* the proper handling of all other specific collaborating classes
* and for the storing of the produced particles in the event record.
*
* @see \ref ShowerHandlerInterfaces "The interfaces"
*
* @see ThePEG::CascadeHandler
* @see MPIHandler
* @see HwRemDecayer
*/
class ShowerHandler: public CascadeHandler {
public:
/** Typedef for a pair of ThePEG::RemnantParticle pointers. */
typedef pair<tRemPPtr, tRemPPtr> RemPair;
/**
* The default constructor.
*/
ShowerHandler();
/**
* Destructor
*/
virtual ~ShowerHandler();
public:
/**
* The main method which manages the multiple interactions and starts
* the shower by calling cascade(sub, lastXC).
*/
virtual void cascade();
/**
* Hook to allow vetoing of event after showering hard sub-process
* as in e.g. MLM merging.
*/
virtual bool showerHardProcessVeto() { return false; };
/**
* Return true, if this cascade handler will perform reshuffling from hard
* process masses.
*/
virtual bool isReshuffling() const { return true; }
public:
/**@name Methods related to PDF freezing */
//@{
/**
* Get the PDF freezing scale
*/
Energy pdfFreezingScale() const {return pdfFreezingScale_;}
//@}
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();
public:
/** @name Functions to access information. */
//@{
/**
* Return true if currently the primary subprocess is showered.
*/
bool firstInteraction() const {
return ( subProcess_ ==
eventHandler()->currentCollision()->primarySubProcess() );
}
/**
* Return the currently used SubProcess.
*/
tSubProPtr currentSubProcess() const {
assert(subProcess_);
return subProcess_;
}
/**
* Return true if multiple parton interactions are switched on
* and can be used for this beam setup.
*/
bool isMPIOn() const {
return MPIHandler_ && MPIHandler_->beamOK();
}
/**
* Return the remnant decayer.
*/
tHwRemDecPtr remnantDecayer() const { return remDec_; }
//@}
/**
* Access to the Evolver
*/
tEvolverPtr evolver() const {return evolver_;}
/**
* Generate hard emissions for CKKW etc
*/
virtual HardTreePtr generateCKKW(ShowerTreePtr tree) const;
/**
* Return true, if the shower handler can generate a truncated
* shower for POWHEG style events generated using Matchbox
*/
virtual bool canHandleMatchboxTrunc() const { return false; }
/**
* The factorization scale factor.
*/
double factorizationScaleFactor() const {
if ( scaleFactorOption_ == 0 || !subProcess_ )
return factorizationScaleFactor_;
if ( scaleFactorOption_ == 1 )
return firstInteraction() ? factorizationScaleFactor_ : 1.0;
if ( scaleFactorOption_ == 2 )
return !firstInteraction() ? factorizationScaleFactor_ : 1.0;
return 1.0;
}
/**
* The renormalization scale factor.
*/
double renormalizationScaleFactor() const {
if ( scaleFactorOption_ == 0 || !subProcess_ )
return renormalizationScaleFactor_;
if ( scaleFactorOption_ == 1 )
return firstInteraction() ? renormalizationScaleFactor_ : 1.0;
if ( scaleFactorOption_ == 2 )
return !firstInteraction() ? renormalizationScaleFactor_ : 1.0;
return 1.0;
}
/**
* The scale factor for the hard scale
*/
double hardScaleFactor() const {
if ( scaleFactorOption_ == 0 || !subProcess_ )
return hardScaleFactor_;
if ( scaleFactorOption_ == 1 )
return firstInteraction() ? hardScaleFactor_ : 1.0;
if ( scaleFactorOption_ == 2 )
return !firstInteraction() ? hardScaleFactor_ : 1.0;
return 1.0;
}
/**
* The option on when to apply the scale factors
*/
int scaleFactorOption() const { return scaleFactorOption_; }
/**
* Return true, if the phase space restrictions of the dipole shower should
* be applied.
*/
bool restrictPhasespace() const { return restrictPhasespace_; }
/**
* Return profile scales
*/
Ptr<HardScaleProfile>::tptr profileScales() const { return hardScaleProfile_; }
/**
* Return the relevant hard scale to be used in the profile scales
*/
virtual Energy hardScale() const;
/**
* Return true if maximum pt should be deduced from the factorization scale
*/
bool hardScaleIsMuF() const { return maxPtIsMuF_; }
/**
- * Parse a variation statement
+ * A struct identifying a shower variation
*/
- static string parseVariations(string, map<string,double>&);
+ struct ShowerVariation {
+
+ /**
+ * Vary the renormalization scale by the given factor.
+ */
+ double renormalizationScaleFactor;
+
+ /**
+ * Vary the factorization scale by the given factor.
+ */
+ double factorizationScaleFactor;
+
+ /**
+ * Apply the variation to the first interaction
+ */
+ bool firstInteraction;
+
+ /**
+ * Apply the variation to the secondary interactions
+ */
+ bool secondaryInteractions;
+
+ /**
+ * Default constructor
+ */
+ ShowerVariation()
+ : renormalizationScaleFactor(1.0),
+ factorizationScaleFactor(1.0),
+ firstInteraction(true),
+ secondaryInteractions(false) {}
+
+ /**
+ * Parse from in file command
+ */
+ void fromInFile(const string&);
+
+ /**
+ * Put to persistent stream
+ */
+ void put(PersistentOStream& os) const;
+
+ /**
+ * Get from persistent stream
+ */
+ void get(PersistentIStream& is);
+
+ };
/**
- * Add additional variations.
+ * Access the shower variations
*/
- static void addToVariations(string in, double v, map<string,double>& vars) {
- vars[in] = v;
+ map<string,ShowerVariation>& showerVariations() {
+ return showerVariations_;
}
/**
- * Add additional variations.
+ * Return the shower variations
*/
- static void addToVariations(const map<string,double>& vin,
- map<string,double>& vars) {
- for ( map<string,double>::const_iterator w = vin.begin();
- w != vin.end(); ++w )
- addToVariations(w->first,w->second,vars);
+ const map<string,ShowerVariation>& showerVariations() const {
+ return showerVariations_;
}
- /**
- * Vary the factorization scale by the indicated values and using the
- * indicated variation labels.
- */
- map<string,double>& renormalizationScaleVariations() {
- return renormalizationScaleVariations_;
- }
+protected:
/**
- * Vary the factorization scale by the indicated values and using the
- * indicated variation labels.
+ * The shower variation weights
*/
- const map<string,double>& renormalizationScaleVariations() const {
- return renormalizationScaleVariations_;
- }
+ map<string,double> currentWeights_;
/**
- * Vary the renormalization scale by the indicated values and using the
- * indicated variation labels.
+ * Combine the variation weights which have been encountered
*/
- map<string,double>& factorizationScaleVariations() {
- return factorizationScaleVariations_;
- }
-
- /**
- * Vary the renormalization scale by the indicated values and using the
- * indicated variation labels.
- */
- const map<string,double>& factorizationScaleVariations() const {
- return factorizationScaleVariations_;
- }
+ void combineWeights();
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:
/**
* Prepare to shower the given subprocess
*/
void prepareCascade(tSubProPtr sub);
/**
* The main method which manages the showering of a subprocess.
*/
virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb);
/**
* Return the maximum number of attempts for showering
* a given subprocess.
*/
unsigned int maxtry() const { return maxtry_; }
/**
* At the end of the Showering, transform ShowerParticle objects
* into ThePEG particles and fill the event record with them.
* Notice that the parent/child relationships and the
* transformation from ShowerColourLine objects into ThePEG
* ColourLine ones must be properly handled.
*/
void fillEventRecord();
/**
* Find the parton extracted from the incoming particle after ISR
*/
PPtr findFirstParton(tPPtr seed) const;
/**
* Fix Remnant connections after ISR
*/
tPPair remakeRemnant(tPPair oldp);
/**
* Get the remnants from the ThePEG::PartonBinInstance es and
* do some checks.
*/
RemPair getRemnants(PBIPair incbins);
/**
* Make the remnant after the shower
*/
void makeRemnants();
/**
* Reset the PDF's after the hard collision has been showered
*/
void setMPIPDFs();
/**
* Boost all the particles in the collision so that the collision always occurs
* in the rest frame with the incoming particles along the z axis
*/
void boostCollision(bool boost);
/**
* Is a beam particle where hadronic structure is resolved
*/
bool isResolvedHadron(tPPtr);
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();
/**
* Called at the end of the run phase.
*/
virtual void dofinish();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerHandler & operator=(const ShowerHandler &);
private:
/**
* Access function for the MPIHandler, it should only be called after
* checking with isMPIOn.
*/
tUEBasePtr getMPIHandler() const {
assert(MPIHandler_);
return MPIHandler_;
}
private:
/**
* a MPIHandler to administer the creation of several (semihard)
* partonic interactions.
*/
UEBasePtr MPIHandler_;
/**
* Pointer to the evolver
*/
EvolverPtr evolver_;
/**
* Pointer to the HwRemDecayer
*/
HwRemDecPtr remDec_;
/**
* The PDF for beam particle A. Overrides the particle's own PDF setting.
*/
PDFPtr PDFA_;
/**
* The PDF for beam particle B. Overrides the particle's own PDF setting.
*/
PDFPtr PDFB_;
/**
* The PDF for beam particle A for remnant splitting. Overrides the particle's own PDF setting.
*/
PDFPtr PDFARemnant_;
/**
* The PDF for beam particle B for remnant splitting. Overrides the particle's own PDF setting.
*/
PDFPtr PDFBRemnant_;
/**
* The PDF freezing scale
*/
Energy pdfFreezingScale_;
/**
* Maximum number of attempts for the
* main showering loop
*/
unsigned int maxtry_;
/**
* Maximum number of attempts for the regeneration of an additional
* scattering, before the number of scatters is reduced.
*/
unsigned int maxtryMPI_;
/**
* Maximum number of attempts for the regeneration of an additional
* hard scattering, before this event is vetoed.
*/
unsigned int maxtryDP_;
/**
* PDG codes of the particles which decay during showering
* this is fast storage for use during running
*/
set<long> particlesDecayInShower_;
/**
* PDG codes of the particles which decay during showering
* this is a vector that is interfaced so they can be changed
*/
vector<long> inputparticlesDecayInShower_;
/**
* Whether or not to include spa-cetime distances in the shower
*/
bool includeSpaceTime_;
/**
* The minimum virtuality for the space-time model
*/
Energy2 vMin_;
/**
* The ShowerTree for the hard process
*/
ShowerTreePtr hard_;
/**
* The incoming beam particles for the current collision
*/
tPPair incoming_;
/**
* The ShowerTree for the decays
*/
ShowerDecayMap decay_;
/**
* The ShowerTrees for which the initial shower
*/
vector<ShowerTreePtr> done_;
/**
* Const pointer to the current step
*/
tcStepPtr current_;
/**
* Const pointer to the currently handeled ThePEG::SubProcess
*/
tSubProPtr subProcess_;
/**
* pointer to "this", the current ShowerHandler.
*/
static ShowerHandler * currentHandler_;
/**
* Boost to get back to the lab
*/
LorentzRotation boost_;
/**
* The MPI PDF's to be used for secondary scatters.
*/
pair <PDFPtr, PDFPtr> mpipdfs_;
/**
* The MPI PDF's to be used for secondary scatters.
*/
pair <PDFPtr, PDFPtr> rempdfs_;
/**
* The MPI PDF's to be used for secondary scatters.
*/
pair <PDFPtr, PDFPtr> remmpipdfs_;
/**
* The factorization scale factor.
*/
double factorizationScaleFactor_;
/**
* The renormalization scale factor.
*/
double renormalizationScaleFactor_;
/**
* The scale factor for the hard scale
*/
double hardScaleFactor_;
/**
* The option on when to apply the scale factors
*/
int scaleFactorOption_;
/**
* True, if the phase space restrictions of the dipole shower should
* be applied.
*/
bool restrictPhasespace_;
/**
* True if maximum pt should be deduced from the factorization scale
*/
bool maxPtIsMuF_;
/**
* The profile scales
*/
Ptr<HardScaleProfile>::ptr hardScaleProfile_;
/**
* Whether or not to split into hard and decay trees
*/
bool splitHardProcess_;
/**
- * Vary the renormalization scale by the indicated values and using the
- * indicated variation labels.
+ * The shower variations
*/
- map<string,double> renormalizationScaleVariations_;
+ map<string,ShowerVariation> showerVariations_;
/**
- * Command to vary the renormalization scale by the indicated values and using
- * the indicated variation labels.
+ * Command to add a shower variation
*/
- string doRenormalizationScaleVariations(string);
-
- /**
- * Vary the factorization scale by the indicated values and using the
- * indicated variation labels.
- */
- map<string,double> factorizationScaleVariations_;
-
- /**
- * Command to vary the factorization scale by the indicated values and using
- * the indicated variation labels.
- */
- string doFactorizationScaleVariations(string);
+ string doAddVariation(string);
public:
/**
* struct that is used to catch exceptions which are thrown
* due to energy conservation issues of additional scatters
*/
struct ExtraScatterVeto {};
/**
* struct that is used to catch exceptions which are thrown
* due to fact that the Shower has been invoked more than
* a defined threshold on a certain configuration
*/
struct ShowerTriesVeto {
/** variable to store the number of attempts */
const int tries;
/** constructor */
ShowerTriesVeto(int t) : tries(t) {}
};
/**
* pointer to "this", the current ShowerHandler.
*/
static const ShowerHandler * currentHandler() {
assert(currentHandler_);
return currentHandler_;
}
protected:
/**
* Set the current handler
*/
void setCurrentHandler() {
currentHandler_ = this;
}
};
+inline PersistentOStream& operator<<(PersistentOStream& os, const ShowerHandler::ShowerVariation& var) {
+ var.put(os); return os;
+}
+
+inline PersistentIStream& operator>>(PersistentIStream& is, ShowerHandler::ShowerVariation& var) {
+ var.get(is); return is;
+}
+
}
#endif /* HERWIG_ShowerHandler_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:39 PM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3799145
Default Alt Text
(123 KB)

Event Timeline