Page MenuHomeHEPForge
Files F8308853

No OneTemporary
Actions

View Options

diff --git a/Shower/Dipole/DipoleShowerHandler.cc b/Shower/Dipole/DipoleShowerHandler.cc
--- a/Shower/Dipole/DipoleShowerHandler.cc
+++ b/Shower/Dipole/DipoleShowerHandler.cc
@@ -1,1299 +1,1302 @@
// -*- 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 "Herwig/PDF/MPIPDF.h"
#include "Herwig/PDF/MinBiasPDF.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include "Herwig/Shower/Dipole/Utility/DipolePartonSplitter.h"
#include "Herwig/MatrixElement/Matchbox/Base/MergerBase.fh"
#include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include <queue>
using namespace Herwig;
bool DipoleShowerHandler::firstWarn = true;
DipoleShowerHandler::DipoleShowerHandler()
: ShowerHandler(), chainOrderVetoScales(true),
nEmissions(0), discardNoEmissions(false), firstMCatNLOEmission(false),
thePowhegDecayEmission(true),
realignmentScheme(0),
verbosity(0), printEvent(0), nTries(0),
didRadiate(false), didRealign(false),
theRenormalizationScaleFreeze(1.*GeV),
theFactorizationScaleFreeze(2.*GeV), theDoCompensate(false),
theFreezeGrid(500000), theDetuning(1.0),
maxPt(ZERO), muPt(ZERO) {}
DipoleShowerHandler::~DipoleShowerHandler() {}
IBPtr DipoleShowerHandler::clone() const {
return new_ptr(*this);
}
IBPtr DipoleShowerHandler::fullclone() const {
return new_ptr(*this);
}
void DipoleShowerHandler::cascade(tPVector finalstate) {
assert(false&&"Dipoleshower as second shower not implemented");
return;}
tPPair DipoleShowerHandler::cascade(tSubProPtr sub, XCombPtr,
Energy optHardPt, Energy optCutoff) {
useMe();
prepareCascade(sub);
resetWeights();
if ( !doFSR() && ! doISR() )
return sub->incoming();
eventRecord().clear();
eventRecord().prepare(sub, dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr()), pdfs(),
ShowerHandler::currentHandler()->generator()->currentEvent()->incoming());
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;
if ( eventRecord().truncatedShower() ) {
throw Exception() << "Inconsistent hard emission set-up in DipoleShowerHandler::cascade. "
<< "No truncated shower needed with DipoleShowerHandler. Add "
<< "'set MEMatching:TruncatedShower No' to input file."
<< Exception::runerror;
}
hardScales(lastXCombPtr()->lastShowerScale());
if ( verbosity > 1 ) {
generator()->log() << "DipoleShowerHandler starting off:\n";
eventRecord().debugLastEvent(generator()->log());
generator()->log() << flush;
}
unsigned int nEmitted = 0;
if ( firstMCatNLOEmission ) {
if ( !eventRecord().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();
// Decay and shower any particles that require decaying
while ( !eventRecord().decays().empty() ) {
map<PPtr,PerturbativeProcessPtr>::iterator decayIt = eventRecord().decays().begin();
// find the decay to do, one with greatest width and parent showered
while(find(eventRecord().outgoing().begin(),eventRecord().outgoing().end(),decayIt->first)==
eventRecord().outgoing().end() &&
find(eventRecord().hard().begin(),eventRecord().hard().end(),decayIt->first)==
eventRecord().hard().end()) ++decayIt;
assert(decayIt!=eventRecord().decays().end());
PPtr incoming = decayIt->first;
eventRecord().currentDecay(decayIt->second);
// Use this to record if an emission actually happens
bool powhegEmission = thePowhegDecayEmission;
// Decay the particle / sort out its pert proc
Energy showerScale = eventRecord().decay(incoming, powhegEmission);
// Following the decay, the bool powheg emission is updated to indicate whether or not an emission occurred
if ( powhegEmission )
nEmitted += 1;
// Check that there is only one particle incoming to the decay
assert(eventRecord().currentDecay()->incoming().size()==1);
// Prepare the event record for the showering of the decay
bool needToShower = eventRecord().prepareDecay(eventRecord().currentDecay());
// Only need to shower if we have coloured outgoing particles
if ( needToShower ) {
// The decays currently considered produce a maximum of 2 chains (with powheg emission)
// so all dipole should have the same scale as returned by the decay function.
assert( eventRecord().chains().size() <= 2 );
for ( list<DipoleChain>::iterator ch = eventRecord().chains().begin();
ch != eventRecord().chains().end(); ++ch) {
for ( list<Dipole>::iterator dip = ch->dipoles().begin();
dip != ch->dipoles().end(); ++dip ) {
assert ( showerScale > ZERO );
dip->leftScale( showerScale );
dip->rightScale( showerScale );
}
}
// Perform the cascade
doCascade(nEmitted,optHardPt,optCutoff,true);
// Do the constituent mass shell reshuffling
decayConstituentReshuffle(eventRecord().currentDecay(), true);
}
// Update the decays, adding any decays and updating momenta
eventRecord().updateDecays(eventRecord().currentDecay());
eventRecord().decays().erase(decayIt->first);
}
break;
}
catch (RedoShower&) {
resetWeights();
if ( ++nTries > maxtry() )
throw ShowerTriesVeto(maxtry());
eventRecord().clear();
eventRecord().prepare(sub,dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr()),pdfs(),ShowerHandler::currentHandler()->generator()->currentEvent()->incoming());
continue;
} catch (...) {
throw;
}
}
tPPair res=eventRecord().fillEventRecord(newStep(),firstInteraction(),didRealign);
setDidRunCascade(true);
return res;
}
// Reshuffle the outgoing partons from the hard process onto their constituent mass shells
void DipoleShowerHandler::constituentReshuffle() {
if ( constituentReshuffler ) {
if ( eventRecord().decays().empty() ) {
constituentReshuffler->reshuffle(eventRecord().outgoing(),
eventRecord().incoming(),
eventRecord().intermediates());
return;
}
// if(eventRecord().decays().empty()) return;
else {
PList decaying;
for(map<PPtr,PerturbativeProcessPtr>::iterator it=eventRecord().decays().begin();
it!=eventRecord().decays().end();++it) {
decaying.push_back(it->first);
}
constituentReshuffler->hardProcDecayReshuffle( decaying,
eventRecord().outgoing(),
eventRecord().hard(),
eventRecord().incoming(),
eventRecord().intermediates());
}
}
// After reshuffling the hard process, the decays need to be updated
// as this is not done in reshuffle
vector<pair<PPtr,PerturbativeProcessPtr> > decays;
for(map<PPtr,PerturbativeProcessPtr>::iterator it=eventRecord().decays().begin();
it!=eventRecord().decays().end();++it) {
decays.push_back(make_pair(it->first,it->second));
}
for(unsigned int ix=0;ix<decays.size();++ix) {
PPtr unstable = decays[ix].first;
PList::iterator pos = find(eventRecord().intermediates().begin(),eventRecord().intermediates().end(),decays[ix].first);
// Update the PPtr in theDecays
if(pos!=eventRecord().intermediates().end()) {
unstable = *pos;
while(!unstable->children().empty()) {
unstable = unstable->children()[0];
}
eventRecord().decays().erase(decays[ix].first);
eventRecord().decays()[unstable] = decays[ix].second;
// Update the momenta of any other particles in the decay chain
// (for externally provided events)
if ( !(eventRecord().decays()[unstable]->outgoing().empty()) )
eventRecord().updateDecayChainMom( unstable , eventRecord().decays()[unstable]);
}
else {
if ( !(eventRecord().decays()[unstable]->outgoing().empty()) ) {
// Update the momenta of any other particles in the decay chain
// (for externally provided events)
// Note this needs to be done for all decaying particles in the
// outgoing/hard regardless of whether that particle radiated
// or was involved in the reshuffling, this is due to the
// transformation performed for IILightKinematics.
if ( (find(eventRecord().outgoing().begin(),
eventRecord().outgoing().end(), unstable) != eventRecord().outgoing().end())
|| (find(eventRecord().hard().begin(),
eventRecord().hard().end(), unstable) != eventRecord().hard().end()) )
eventRecord().updateDecayChainMom( unstable , eventRecord().decays()[unstable]);
}
}
}
eventRecord().currentDecay(PerturbativeProcessPtr());
}
// Reshuffle outgoing partons from a decay process onto their constituent mass shells
void DipoleShowerHandler::decayConstituentReshuffle(PerturbativeProcessPtr decayProc, const bool test) {
if ( !test ) {
// decayReshuffle updates both the event record and the decay perturbative process
if ( constituentReshuffler ) {
constituentReshuffler->decayReshuffle(decayProc,
eventRecord().outgoing(),
eventRecord().hard(),
eventRecord().intermediates());
}
return;
}
if ( test ) {
// Test this function by comparing the
// invariant mass of the outgoing decay
// systems before and after reshuffling
Lorentz5Momentum testOutMomBefore (ZERO,ZERO,ZERO,ZERO);
Energy testInvMassBefore = ZERO;
for ( vector<pair<PPtr,PerturbativeProcessPtr> >::iterator testDecayOutItBefore = decayProc->outgoing().begin();
testDecayOutItBefore!= decayProc->outgoing().end(); ++testDecayOutItBefore ) {
testOutMomBefore += testDecayOutItBefore->first->momentum();
}
testInvMassBefore = testOutMomBefore.m();
// decayReshuffle updates both the event record and the decay perturbative process
if ( constituentReshuffler ) {
constituentReshuffler->decayReshuffle(decayProc,
eventRecord().outgoing(),
eventRecord().hard(),
eventRecord().intermediates());
}
Lorentz5Momentum testOutMomAfter (ZERO,ZERO,ZERO,ZERO);
Energy testInvMassAfter = ZERO;
for ( vector<pair<PPtr,PerturbativeProcessPtr> >::iterator testDecayOutItAfter = decayProc->outgoing().begin();
testDecayOutItAfter!= decayProc->outgoing().end(); ++testDecayOutItAfter ) {
testOutMomAfter += testDecayOutItAfter->first->momentum();
}
testInvMassAfter = testOutMomAfter.m();
Energy incomingMass = decayProc->incoming()[0].first->momentum().m();
assert( abs(testInvMassBefore-incomingMass)/GeV < 1e-5 );
assert( abs(testInvMassBefore-testInvMassAfter)/GeV < 1e-5);
}
}
// Sets the scale of each particle in the dipole chains by finding the smallest of several upper bound energy scales: the CMEnergy of the event, the transverse mass of outgoing particles, the hardScale (maxPT or maxQ) calculated for each dipole (in both configurations) and the veto scale for each particle
void DipoleShowerHandler::hardScales(Energy2 muf) {
// Initalise maximum pt as max CMEnergy of the event
maxPt = generator()->maximumCMEnergy();
if ( restrictPhasespace() ) {
// First interaction == hard collision (i.e. not a MPI collision)
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());
}
//Look at any non-coloured outgoing particles in the current subprocess
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 ) {
// Note that minVetoScale is a value for each DipoleChain, not each dipole
// It will contain the minimum veto scale from all of the dipoles in the chain
Energy minVetoScale = -1.*GeV;
for ( list<Dipole>::iterator dip = ch->dipoles().begin();
dip != ch->dipoles().end(); ++dip ) {
// max scale per config
Energy maxFirst = ZERO;
Energy maxSecond = ZERO;
// Loop over the kernels for the given dipole.
// For each dipole configuration, calculate ptMax (or QMax if virtuality ordering) for each kernel and find the maximum
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)) ) {
// Look in DipoleChainOrdering for this
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);
}
}
// Find the maximum value from comparing the maxScale found from maxPt and the vetoScale of the particle
if ( dip->leftParticle()->vetoScale() >= ZERO ) {
maxFirst = min(maxFirst,sqrt(dip->leftParticle()->vetoScale()));
// minVetoScale is a value for each DipoleChain, not each dipole
// It contains the minimum veto scale for all the dipoles in the entire DipoleChain
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());
}
// Set the emitterScale for both members of each dipole
maxFirst = min(maxPt,maxFirst);
dip->emitterScale(make_pair(true,false),maxFirst);
maxSecond = min(maxPt,maxSecond);
dip->emitterScale(make_pair(false,true),maxSecond);
}
// if the smallest veto scale (i.e. from all of the dipoles) is smaller than the scale calculated for a particular particle in a particular dipole, replace the scale with the veto scale
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;
}
// Currently do not split IF dipoles so
// don't evaluate them in order to avoid
// exceptions in the log
if ( index.incomingDecayEmitter() ) {
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(),theSplittingReweight);
//
// 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.
//
// SW - Update 04/01/2016: Note - This caused a bug for me as I did not
// include equality checks on the decay booleans in the == definition
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 ( std::isnan( double(dScale/MeV) ) )
throw RedoShower();
candidate.scale(dScale);
// Calculate the mass of the recoil system
// for decay dipoles
if (candidate.index().incomingDecayEmitter() || candidate.index().incomingDecaySpectator() ) {
Energy recoilMass = gen->second->splittingKinematics()->recoilMassKin(emitter->momentum(),
spectator->momentum());
candidate.recoilMass(recoilMass);
}
candidate.continuesEvolving();
Energy hardScale = evolutionOrdering()->maxPt(startScale,candidate,*(gen->second->splittingKernel()));
Energy maxPossible =
gen->second->splittingKinematics()->ptMax(candidate.scale(),
candidate.emitterX(), candidate.spectatorX(),
candidate,
*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);
Energy nextScale = evolutionOrdering()->evolutionScale(gen->second->lastSplitting(),*(gen->second->splittingKernel()));
if ( nextScale > winnerScale ) {
winner.fill(candidate);
gen->second->completeSplitting(winner);
winnerGen = gen;
winnerScale = nextScale;
}
reweight(reweight() * gen->second->splittingWeight());
}
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,
const bool decay) {
if ( nEmissions )
if ( emDone == nEmissions )
return;
DipoleSplittingInfo winner;
DipoleSplittingInfo dipoleWinner;
+
+
+ size_t currentMultiplicity=eventHandler()->currentStep()->particles().size();
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() ) {
eventRecord().popChain();
if ( theEventReweight && eventRecord().chains().empty() )
if ( (theEventReweight->firstInteraction() && firstInteraction()) ||
(theEventReweight->secondaryInteractions() && !firstInteraction()) ) {
double w = theEventReweight->weightCascade(eventRecord().incoming(),
eventRecord().outgoing(),
eventRecord().hard(),theGlobalAlphaS);
reweight(reweight()*w);
}
continue;
}
// otherwise perform the splitting
-
+
if (theMergingHelper&&eventHandler()->currentCollision()) {
- if (theMergingHelper->maxLegs()>eventRecord().outgoing().size()+eventRecord().hard().size()+2)
+ if (theMergingHelper->maxLegs()>currentMultiplicity)
if (theMergingHelper->mergingScale()<winnerScale){
bool transparent=true;
if (transparent) {
pair<list<Dipole>::iterator,list<Dipole>::iterator> tmpchildren;
DipoleSplittingInfo tmpwinner=winner;
DipoleChain* tmpfirstChain = 0;
DipoleChain* tmpsecondChain = 0;
PVector New=eventRecord().tmpsplit(winnerDip,tmpwinner,tmpchildren,tmpfirstChain,tmpsecondChain);
/*cout<<"\n ";
for (PVector::iterator p=New.begin(); p!=New.end(); p++) {
cout<<"\n "<<(*p)->momentum()/GeV<<flush;
}
*/
if (theMergingHelper->matrixElementRegion(New,winnerScale,theMergingHelper->mergingScale())) {
optHardPt=winnerScale;
continue;
}
}else{
optHardPt=winnerScale;
continue;
}
}
}
didRadiate = true;
eventRecord().isMCatNLOSEvent(false);
eventRecord().isMCatNLOHEvent(false);
pair<list<Dipole>::iterator,list<Dipole>::iterator> children;
DipoleChain* firstChain = 0;
DipoleChain* secondChain = 0;
// Note: the dipoles are updated in eventRecord().split(....) after the splitting,
// hence the entire cascade is handled in doCascade
// The dipole scales are updated in dip->split(....)
if ( decay )
winner.isDecayProc( true );
eventRecord().split(winnerDip,winner,children,firstChain,secondChain);
-
+ currentMultiplicity++;
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 )
if ( (theEventReweight->firstInteraction() && firstInteraction()) ||
(theEventReweight->secondaryInteractions() && !firstInteraction()) ) {
double w = theEventReweight->weight(eventRecord().incoming(),
eventRecord().outgoing(),
eventRecord().hard(),theGlobalAlphaS);
reweight(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 ) {
if ( !(**k).alphaS() )
(**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);
if ( renormalizationScaleFactor() != 1. )
nGenerator->splittingKernel()->renormalizationScaleFactor(renormalizationScaleFactor());
if ( factorizationScaleFactor() != 1. )
nGenerator->splittingKernel()->factorizationScaleFactor(factorizationScaleFactor());
if ( !nGenerator->splittingReweight() )
nGenerator->splittingReweight(rw);
nGenerator->splittingKernel()->freezeGrid(theFreezeGrid);
nGenerator->splittingKernel()->detuning(theDetuning);
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->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);
}
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
<< thePowhegDecayEmission
<< realignmentScheme << verbosity << printEvent
<< ounit(theRenormalizationScaleFreeze,GeV)
<< ounit(theFactorizationScaleFreeze,GeV)
<< theShowerApproximation
<< theDoCompensate << theFreezeGrid << theDetuning
<< theEventReweight << theSplittingReweight << ounit(maxPt,GeV)
<< ounit(muPt,GeV)<< theMergingHelper;
}
void DipoleShowerHandler::persistentInput(PersistentIStream & is, int) {
is >> kernels >> theEvolutionOrdering
>> constituentReshuffler >> intrinsicPtGenerator
>> theGlobalAlphaS >> chainOrderVetoScales
>> nEmissions >> discardNoEmissions >> firstMCatNLOEmission
>> thePowhegDecayEmission
>> realignmentScheme >> verbosity >> printEvent
>> iunit(theRenormalizationScaleFreeze,GeV)
>> iunit(theFactorizationScaleFreeze,GeV)
>> theShowerApproximation
>> theDoCompensate >> theFreezeGrid >> theDetuning
>> theEventReweight >> theSplittingReweight >> iunit(maxPt,GeV)
>> iunit(muPt,GeV)>>theMergingHelper;
}
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,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 Parameter<DipoleShowerHandler,double> interfaceDetuning
("Detuning",
"A value to detune the overestimate kernel.",
&DipoleShowerHandler::theDetuning, 1.0, 1.0, 0,
false, false, Interface::lowerlim);
static Reference<DipoleShowerHandler,DipoleEventReweight> interfaceEventReweight
("EventReweight",
"",
&DipoleShowerHandler::theEventReweight, false, false, true, true, false);
static Reference<DipoleShowerHandler,DipoleSplittingReweight> interfaceSplittingReweight
("SplittingReweight",
"Set the splitting reweight.",
&DipoleShowerHandler::theSplittingReweight, false, false, true, true, false);
static Switch<DipoleShowerHandler, bool> interfacePowhegDecayEmission
("PowhegDecayEmission",
"Use Powheg style emission for the decays",
&DipoleShowerHandler::thePowhegDecayEmission, true, false, false);
static SwitchOption interfacePowhegDecayEmissionYes
(interfacePowhegDecayEmission,"Yes","Powheg decay emission on", true);
static SwitchOption interfacePowhegDecayEmissionNo
(interfacePowhegDecayEmission,"No","Powheg decay emission off", false);
static Reference<DipoleShowerHandler,MergerBase> interfaceMergingHelper
("MergingHelper",
"",
&DipoleShowerHandler::theMergingHelper, false, false, true, true, false);
}
diff --git a/Shower/Dipole/Merging/Merger.cc b/Shower/Dipole/Merging/Merger.cc
--- a/Shower/Dipole/Merging/Merger.cc
+++ b/Shower/Dipole/Merging/Merger.cc
@@ -1,1556 +1,1557 @@
// -*- C++ -*-
//
// Merger.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 Merger class.
//
#include "Merger.h"
#include "Node.h"
#include "MFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "Herwig/Shower/Dipole/DipoleShowerHandler.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Config/Constants.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "fastjet/ClusterSequence.hh"
using namespace Herwig;
Merger::Merger()
: MergerBase() {
minusL=false;
Unlopsweights=true;
theCMWScheme=true;
theSmearing=0.;
theDipoleShowerHandler=
Ptr<DipoleShowerHandler>::ptr();
FFLTK=new_ptr(FFLightTildeKinematics());
FILTK=new_ptr(FILightTildeKinematics());
IFLTK=new_ptr(IFLightTildeKinematics());
IILTK=new_ptr(IILightTildeKinematics());
FFMTK=new_ptr(FFMassiveTildeKinematics());
FIMTK=new_ptr(FIMassiveTildeKinematics());
IFMTK=new_ptr(IFMassiveTildeKinematics());
isUnitarized=true;
isNLOUnitarized=true;
defMERegionByJetAlg=false;
theOpenInitialStateZ=false;
theChooseHistory=8;
xiRenME=1.;
xiFacME=1.;
xiRenSh=1.;
xiFacSh=1.;
xiQSh=1.;
theGamma=1.;
ee_ycut=1000000;
pp_dcut=1000000;
therenormscale=0.*GeV;
theIRSafePT=1.*GeV;
theMergePt=3.*GeV;
theCentralMergePt=3.*GeV;
}
Merger::~Merger() {}
IBPtr Merger::clone() const {
return new_ptr(*this);
}
IBPtr Merger::fullclone() const {
return new_ptr(*this);
}
CrossSection Merger::MergingDSigDRBornGamma(NPtr Node){
double weightCL=0.;
weight=1.;
if (!Node->children().empty()) {
PVector particles;
for (size_t i=0;i<Node->nodeME()->mePartonData().size();i++){
Ptr<ThePEG::Particle>::ptr p =Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]);
particles.push_back(p);
}
if(!matrixElementRegion(particles, Node->children()[0]->dipol()->lastPt(), theMergePt))weight*= 0.;
}
NPtr Born= Node-> getHistory(true,xiQSh);
NPtr CLNode;
NPtr BornCL;
if( !Node->children().empty()){
if (UseRandom::rnd()<.5){
CLNode= Node->randomChild();
bool inhist=CLNode->isInHistoryOf(Born);
weight*=inhist?(-2.*int(Node->children().size())):0.;
projected=true;Node->nodeME()->projectorStage(1);
weightCL=2.*int(Node->children().size());
BornCL=CLNode-> getHistory(false,xiQSh);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (treefactory()->onlymulti()!=-1&&
treefactory()->onlymulti()!=int(Node->legsize()-Node->nodeME()->projectorStage()))
return ZERO;
if(!Node->allAbove(mergePt()-0.1*GeV))weight=0.;
if(CLNode&&!CLNode->allAbove(mergePt()-0.1*GeV))weightCL=0.;
if (!Born->xcomb()->willPassCuts()){
if (Born->parent()) {
//cout<<"\n parent not pass";
}
return ZERO;
}
CrossSection res=ZERO;
bool maxMulti=Node->legsize() == int(maxLegsLO());
if(weight!=0.){
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node);
if (!fillProjector(projectedscale))weight=0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.&&weightCL==0.)return ZERO;
res+=weight*TreedSigDR(startscale,Node,(!maxMulti&&!projected)?theGamma:1.);
}
if(CLNode&&theGamma!=1.){
Energy startscale=CKKW_StartScale(BornCL);
Energy projectedscale=startscale;
fillHistory( startscale, BornCL, CLNode);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weightCL*=history.back().weight*alphaReweight()*pdfReweight();
CrossSection diff=ZERO;
Node->nodeME()->factory()->setAlphaParameter(1.);
diff-=weightCL*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME));
Node->nodeME()->factory()->setAlphaParameter(theGamma);
string alp=(CLNode->dipol()->aboveAlpha()?"Above":"Below");
diff+=weightCL*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME));
Node->nodeME()->factory()->setAlphaParameter(1.);
res+=diff;
}
return res;
}
CrossSection Merger::MergingDSigDRBornStandard(NPtr Node){
if (!Node->children().empty()) {
PVector particles;
for (size_t i=0;i<Node->nodeME()->mePartonData().size();i++){
Ptr<ThePEG::Particle>::ptr p =Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]);
particles.push_back(p);
}
if(!matrixElementRegion(particles, Node->children()[0]->dipol()->lastPt(), theMergePt))return ZERO;
}
NPtr Born= Node-> getHistory(true,xiQSh);
if( Born!= Node){
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;Node->nodeME()->projectorStage(1);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (treefactory()->onlymulti()!=-1&&
treefactory()->onlymulti()!=int(Node->xcomb()->mePartonData().size()-Node->nodeME()->projectorStage()))
return ZERO;
assert(Node->allAbove(mergePt()-0.1*GeV));
if (!Born->xcomb()->willPassCuts()){
return ZERO;
}
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
double w1=history.back().weight;
double w2=alphaReweight();
double w3=pdfReweight();
if(std::isnan(w1)){cout<<"\nhistory weight is nan";w1=0.;};
if(std::isnan(w2)){cout<<"\nalphaReweight weight is nan";w1=0.;};
if(std::isnan(w3)){cout<<"\npdfReweight weight is nan";w1=0.;};
weight*=w1*w2*w3;
if(weight==0.)return ZERO;
bool maxMulti=Node->xcomb()->meMomenta().size() == maxLegsLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
CrossSection w4=TreedSigDR(startscale,Node,1.);
if(std::isnan(double(w4/nanobarn))){cout<<"\nTreedSigDR is nan";w1=0.;};
return weight*w4;
}
CrossSection Merger::MergingDSigDRVirtualStandard(NPtr Node ){
NPtr Born= Node-> getHistory(true,xiQSh);
if( Born!= Node){
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;Node->nodeME()->projectorStage(1);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node);
//if (history.size()==1&&Node->children().size()!=0) {
// cout<<"\n1-->"<<startscale/GeV<<" "<<weight;
//}
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
bool maxMulti=Node->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
CrossSection matrixElement=LoopdSigDR(startscale,Node);
CrossSection Bornweight=Node->nodeME()->dSigHatDRB();
CrossSection unlopsweight =(-sumpdfReweightUnlops()
-sumalphaReweightUnlops()
-sumfillHistoryUnlops())
*Bornweight
*SM().alphaS()/(2.*ThePEG::Constants::pi);
return weight*
as(startscale*xiRenSh)/SM().alphaS()*
(matrixElement+unlopsweight);
}
CrossSection Merger::MergingDSigDRRealStandard(NPtr Node){
bool allAbove=Node->allAbove(mergePt());
if(!Node->allAbove(theIRSafePT))return ZERO;
if (allAbove)return MergingDSigDRRealAllAbove(Node);
if (UseRandom::rnd()<.5)
return 2.*MergingDSigDRRealBelowSubReal( Node );
return 2.*MergingDSigDRRealBelowSubInt( Node);
}
CrossSection Merger::MergingDSigDRRealAllAbove(NPtr Node){
NPtr Born= Node-> getHistory(true,xiQSh);
NPtr CLNode= Node->randomChild();
bool inhist=CLNode->isInHistoryOf(Born);
Born=CLNode-> getHistory(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(2);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(1);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(1);
}
if (!CLNode->allAbove(mergePt()))return ZERO;
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, CLNode);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
bool maxMulti=CLNode->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
CrossSection res= weight*as(startscale*xiRenSh)/SM().alphaS()*
(double)Node->children().size()*
((inhist?TreedSigDR(startscale,Node):ZERO)
+CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME)));
// cout<<"\nCLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn "<<CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
return res;
}
CrossSection Merger::MergingDSigDRRealBelowSubReal(NPtr Node){
NPtrVec children=Node->children();
Selector<NPtr> HistNodeSel;
Energy minScale=generator()->maximumCMEnergy();
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
if ((*child)->dipol()->lastPt()<minScale) {
minScale=(*child)->dipol()->lastPt();
HistNodeSel.insert(1.,*child);
}
}
NPtr HistNode=HistNodeSel.select(UseRandom::rnd());
NPtr Born=HistNode-> getHistory(false,xiQSh);
if( Born!= HistNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(1);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(0);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(0);
}
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, HistNode);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
bool maxMulti=HistNode->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
CrossSection sumPS=ZERO;
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
if ((*child)->allAbove(mergePt()))
sumPS-=(*child)->calcPs(startscale*xiFacME);
}
CrossSection me=TreedSigDR(startscale,Node);
//cout<<"\nSubReal "<<Node->miniPt()/GeV<<" "<<me<<" "<<sumPS<<" "<<me/sumPS;
return weight*as(startscale*xiRenSh)/SM().alphaS()*
(me-sumPS);
}
CrossSection Merger::MergingDSigDRRealBelowSubInt(NPtr Node){
NPtr CLNode= Node->randomChild();
NPtr Born=CLNode-> getHistory(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(2);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(1);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(1);
}
if (!CLNode->allAbove(mergePt()))return ZERO;
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, CLNode);
if (!fillProjector(projectedscale))return ZERO;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
bool maxMulti=CLNode->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
pair<CrossSection,CrossSection> DipAndPs=
CLNode->calcDipandPS(startscale*xiFacME);
return weight*as(startscale*xiRenSh)/SM().alphaS()*
(double)Node->children().size()*(DipAndPs.second-DipAndPs.first);
}
CrossSection Merger::TreedSigDR(Energy startscale,NPtr Node,double diffAlpha){
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
theCalculateInNode=false;
CrossSection res=DeepHead->nodeME()->dSigHatDRB();
if (diffAlpha!=1.) {
res+=DeepHead->nodeME()->dSigHatDRAlphaDiff(diffAlpha);
}
renormscale(0.0*GeV);
theCalculateInNode=true;
return res;
}
CrossSection Merger::LoopdSigDR(Energy startscale,NPtr Node){
// The deephead should be calculated here.
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
+ DeepHead->nodeME()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
theCalculateInNode=false;
DeepHead->nodeME()->doOneLoopNoBorn();
CrossSection res=DeepHead->nodeME()->dSigHatDRV()+DeepHead->nodeME()->dSigHatDRI();
DeepHead->nodeME()->noOneLoopNoBorn();
renormscale(0.0*GeV);
theCalculateInNode=true;
return res;
}
bool Merger::fillProjector(Energy& prerunning){
// in the shower handler the scale is multiplied by xiQSh
// so here we need to devide this
double xiQShfactor=history.begin()->node->legsize()==N0()?xiQSh:1.;
if(history.begin()->node->deepHead()->nodeME()->projectorStage() == 0){
prerunning=(history.size()==1?1.:(1./xiQShfactor))*history.back().scale;
return true;
}
for (History::iterator it=history.begin();it!=history.end();it++){
if (projectorStage((*it).node)&&history.begin()->node->deepHead()->nodeME()->projectorStage() != 0){
history.begin()->node->deepHead()->xcomb()->lastProjector((*it).node->xcomb());
prerunning=(it==history.begin()?1.:(1./xiQShfactor))*(*it).scale;
return true;
}
}
return false;
}
double Merger::pdfReweight(){
double res=1.;
double facfactor=(history[0].node->legsize()==N0()?xiFacME:xiFacSh);
for(int side=0;side!=2;side++){
if(history[0].node->xcomb()->mePartonData()[side]->coloured()){
History::iterator it=history.begin();
for (;it+1!=history.end();it++){
res *= pdfratio((*it).node, facfactor*(*it).scale,xiFacSh*((*it).node->dipol()->lastPt()), side);
facfactor=xiFacSh;
}
res*=pdfratio(history.back().node,facfactor*history.back().scale ,history[0].scale*xiFacME, side);
}
}
return res;
}
double Merger::alphaReweight(){
double res=1.;
Energy Q_R=(history[0].node->legsize()==N0()?xiRenME:xiRenSh)*history[0].scale;
res *= pow(as(Q_R) / SM().alphaS(), history[0].node->nodeME()->orderInAlphaS());
res *= pow(SM().alphaEMME(history[0].node->nodeME()->factory()->scaleChoice()->renormalizationScaleQED())/ SM().alphaEMMZ(), history[0].node->nodeME()->orderInAlphaEW());
if (!(history[0].node->children().empty())){
res *=pow((theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*Nf(Q_R))*as(Q_R))/2./Constants::pi):1.),int(history[0].node->legsize()-N0()));
}
for (History::iterator it=history.begin();(it+1)!=history.end();it++){
if ((*it).node->parent()){
Energy q_i=xiRenSh* (*it).node->dipol()->lastPt();
res *= as(q_i)/ SM().alphaS()
*(theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*Nf(q_i))*as(q_i))/2./Constants::pi):1.);
}
}
return res;
}
void Merger::fillHistory(Energy scale, NPtr Begin, NPtr EndNode){
history.clear();
double sudakov0_n=1.;
history.push_back(HistoryStep(Begin,sudakov0_n,scale));
double xiQShfactor=history.begin()->node->legsize()==N0()?xiQSh:1.;
scale*=xiQShfactor;
if (Begin->parent()||!isUnitarized) {
while (Begin->parent() && (Begin != EndNode)) {
if (!dosudakov(Begin,scale, Begin->dipol()->lastPt(),sudakov0_n)){
history.push_back(HistoryStep(Begin->parent(),0.,scale));
}
scale=Begin->dipol()->lastPt();
history.push_back(HistoryStep(Begin->parent(),sudakov0_n,Begin->dipol()->lastPt()));
Begin = Begin->parent();
}
Energy notunirunning=scale;
if (!isUnitarized&&N()+N0() > int(Begin->deepHead()->legsize())) {
if (!dosudakov(Begin,notunirunning,mergePt(),sudakov0_n)){
history.back().weight=0.;
}else{
history.back().weight=sudakov0_n;
}
}
}
if( history.size()==1)scale/=xiQSh;
}
double Merger::sumpdfReweightUnlops(){
double res=0.;
Energy beam1Scale=history[0].scale*(history[0].node->legsize()==N0()?xiFacME:xiFacSh);
Energy beam2Scale=history[0].scale*(history[0].node->legsize()==N0()?xiFacME:xiFacSh);
History::iterator it=history.begin();
for (;it!=history.end();it++){
History::iterator ittmp=it;
ittmp++;
if((*it).node->xcomb()->mePartonData()[0]->coloured()&&(*it).node->nodeME()->lastX1()>0.99)return 0.;
if((*it).node->xcomb()->mePartonData()[1]->coloured()&&(*it).node->nodeME()->lastX2()>0.99)return 0.;
if (ittmp!=history.end()){
if((*it).node->nodeME()->lastX1()<0.00001)return 0.;
if((*it).node->nodeME()->lastX2()<0.00001)return 0.;
if ((*it).node->dipol()->bornEmitter() == 0 ){
res +=pdfUnlops((*it).node->nodeME()->lastParticles().first->dataPtr(),
(*it).node->nodeME()->lastPartons().first->dataPtr(),
(*it).node->xcomb()->partonBins().first->pdf(),
beam1Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX1(),
Nf(history[0].scale),
history[0].scale);
beam1Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornEmitter() == 1 ){
res +=pdfUnlops((*it).node->nodeME()->lastParticles().second->dataPtr(),
(*it).node->nodeME()->lastPartons().second->dataPtr(),
(*it).node->xcomb()->partonBins().second->pdf(),
beam2Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX2(),
Nf(history[0].scale),
history[0].scale);
beam2Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornSpectator() == 0 &&(*it).node->dipol()->bornEmitter() >1){//
res +=pdfUnlops((*it).node->nodeME()->lastParticles().first->dataPtr(),
(*it).node->nodeME()->lastPartons().first->dataPtr(),
(*it).node->xcomb()->partonBins().first->pdf(),
beam1Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX1(),
Nf(history[0].scale),
history[0].scale);
//pdfratio((*it).node, beam1Scale, sqrt((*it).node->dipol()->lastPt()), 1);
beam1Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornSpectator() == 1 &&(*it).node->dipol()->bornEmitter() >1){//
res +=pdfUnlops((*it).node->nodeME()->lastParticles().second->dataPtr(),
(*it).node->nodeME()->lastPartons().second->dataPtr(),
(*it).node->xcomb()->partonBins().second->pdf(),
beam2Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX2(),
Nf(history[0].scale),
history[0].scale);
//pdfratio((*it).node, beam2Scale , sqrt((*it).node->dipol()->lastPt()), 2);
beam2Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
}
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[0]->coloured()){
res +=pdfUnlops(history.back().node->nodeME()->lastParticles().first->dataPtr(),
history.back().node->nodeME()->lastPartons().first->dataPtr(),
history.back().node->xcomb()->partonBins().first->pdf(),
beam1Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX1(),
Nf(history[0].scale),
history[0].scale);
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[1]->coloured()) {
res +=pdfUnlops(history.back().node->nodeME()->lastParticles().second->dataPtr(),
history.back().node->nodeME()->lastPartons().second->dataPtr(),
history.back().node->xcomb()->partonBins().second->pdf(),
beam2Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX2(),
Nf(history[0].scale),
history[0].scale);
//history[0].node->deepHead()->nodeME()->pdf2(sqr(beam2Scale))/history[0].node->deepHead()->nodeME()->pdf2(sqr(history[0].scale));
}
return res;
}
double Merger::pdfUnlops(tcPDPtr particle,tcPDPtr parton,tcPDFPtr pdf,Energy running, Energy next,double x,int nlp,Energy fixedScale) {
//copied from PKOperator
double res=0.;
int number=10;//*int((max(running/next,next/running)));
for (int nr=0;nr<number;nr++){
double restmp=0;
using namespace RandomHelpers;
double r = UseRandom::rnd();
double eps = 1e-3;
pair<double,double> zw =
generate((piecewise(),
flat(0.0,x),
match(inverse(0.0,x,1.0) + inverse(1.0+eps,x,1.0))),r);
double z = zw.first;
double mapz = zw.second;
double PDFxparton=pdf->xfx(particle,parton,sqr(fixedScale),x)/x;
double CA = SM().Nc();
double CF = (SM().Nc()*SM().Nc()-1.0)/(2.*SM().Nc());
if (abs(parton->id()) < 7) {
double PDFxByzgluon=pdf->xfx(particle,getParticleData(ParticleID::g),sqr(fixedScale),x/z)*z/x;
double PDFxByzparton=pdf->xfx(particle,parton,sqr(fixedScale),x/z)*z/x;
assert(abs(parton->id()) < 7);
restmp += CF*(3./2.+2.*log(1.-x)) * PDFxparton;
if ( z > x ) {
restmp+= 0.5 * ( sqr(z) + sqr(1.-z) ) * PDFxByzgluon / z;
restmp += CF*2.*(PDFxByzparton - z*PDFxparton)/(z*(1.-z));
restmp -= CF*PDFxByzparton * (1.+z)/z;
}
}else{
assert(parton->id() == ParticleID::g);
double PDFxByzgluon=pdf->xfx(particle,getParticleData(ParticleID::g),sqr(fixedScale),x/z)*z/x;
if ( z > x ){
double factor = CF * ( 1. + sqr(1.-z) ) / sqr(z);
for ( int f = -nlp; f <= nlp; ++f ) {
if ( f == 0 )
continue;
restmp += pdf->xfx(particle,getParticleData(f),sqr(fixedScale),x/z)*z/x*factor;
}
}
restmp += ( (11./6.) * CA - (1./3.)*Nf(history[0].scale) + 2.*CA*log(1.-x) ) *PDFxparton;
if ( z > x ) {
restmp += 2. * CA * ( PDFxByzgluon - z*PDFxparton ) / (z*(1.-z));
restmp += 2.* CA *( (1.-z)/z - 1. + z*(1.-z) ) * PDFxByzgluon / z;
}
}
if (PDFxparton<1e-8)restmp= 0.;
res+=1*restmp*log(sqr(running/next))/PDFxparton*mapz;
}
return res/number;
}
double Merger::sumalphaReweightUnlops(){
double res=0.;
if (!(history[0].node->children().empty())){
res +=alphasUnlops(history[0].scale*xiRenSh,
history[0].scale*xiRenME)*int(history[0].node->legsize()-N0());
assert(int(history[0].node->legsize()-N0()>0));
}
// dsig is calculated with fixed alpha_s
for (History::iterator it=history.begin();(it+1)!=history.end();it++){
assert((*it).node->parent());
res +=alphasUnlops((*it).node->dipol()->lastPt()*xiRenSh ,history[0].scale);
}
return res;
}
double Merger::sumfillHistoryUnlops(){
double res=0.;
double xiQShfactor=history.begin()->node->legsize()==N0()?xiQSh:1.;
for (History::iterator it = history.begin(); (it+1) != history.end();it++){
doUNLOPS((*it).node,(it == history.begin()?xiQShfactor:1.)*(*it).scale , (*it).node->dipol()->lastPt() , history[0].scale, res);
}
return res;
}
Ptr<MFactory>::ptr Merger::treefactory(){return theTreeFactory;}
void Merger::doinit(){
assert(DSH()->hardScaleIsMuF());
}
CrossSection Merger::MergingDSigDR() {
history.clear();
if (theFirstNodeMap.find(theCurrentME)==theFirstNodeMap.end()) {
cout<<"\nnot in map:"<<theFirstNodeMap.size();
}
NPtr Node = theFirstNodeMap[theCurrentME]; assert(Node);
xiRenME=theCurrentME->renormalizationScaleFactor();
xiFacME=theCurrentME->factorizationScaleFactor();
xiRenSh=DSH()->renormalizationScaleFactor();
xiFacSh=DSH()->factorizationScaleFactor();
xiQSh=DSH()->hardScaleFactor();
DSH()->setCurrentHandler();
DSH()->eventHandler(generator()->eventHandler());
assert(DSH()->hardScaleIsMuF());
CrossSection res=ZERO;
if(Node->deepHead()->subtractedReal()){
res=MergingDSigDRRealStandard(Node);
theCurrentMaxLegs=maxLegsNLO();
}else if(Node->deepHead()->virtualContribution()){
res=MergingDSigDRVirtualStandard(Node);
theCurrentMaxLegs=maxLegsNLO();
}else if(theGamma!=1.){
res=MergingDSigDRBornGamma(Node);
theCurrentMaxLegs=maxLegsLO();
}else{
res=MergingDSigDRBornStandard(Node);
theCurrentMaxLegs=maxLegsLO();
}
//cout<<"\nrunning "<<Node->runningPt()/GeV<<" "<<Node->legsize();
theCurrentME->lastXCombPtr()->lastCentralScale(sqr(Node->runningPt()));
theCurrentME->lastXCombPtr()->lastShowerScale(sqr(Node->runningPt()));
if(theCurrentME->lastXCombPtr()->lastProjector()){
theCurrentME->lastXCombPtr()->lastProjector()->lastCentralScale(sqr(Node->runningPt()));
theCurrentME->lastXCombPtr()->lastProjector()->lastShowerScale(sqr(Node->runningPt()));
}
renormscale(0.0*GeV);
if (res == ZERO){
history.clear();
return res;
}
cleanup(Node);
DSH()->eventRecord().clear();
theCurrentME->lastXCombPtr()->subProcess(SubProPtr());
history.clear();
if(std::isnan(double(res/nanobarn))|| !std::isfinite(double(res/nanobarn))){
cout<<"Warning merger weight is "<<res/nanobarn<<" -> setting to 0";
return ZERO;
}
return res;
}
//----------------------------------Reviewed--------------------------------------//
void Merger::CKKW_PrepareSudakov(NPtr Born,Energy running){
//cleanup(Born);
tSubProPtr sub = Born->xcomb()->construct();
DSH()->resetPDFs(make_pair(Born->xcomb()->partonBins().first->pdf(),
Born->xcomb()->partonBins().second->pdf()),
Born->xcomb()->partonBins());
DSH()->setCurrentHandler();
DSH()->currentHandler()->generator()->currentEventHandler(Born->deepHead()->xcomb()->eventHandlerPtr());
DSH()->currentHandler()->remnantDecayer()->setHadronContent(Born->deepHead()->xcomb()->lastParticles());
DSH()->eventRecord().clear();
DSH()->eventRecord().slimprepare(sub, dynamic_ptr_cast<tStdXCombPtr>(Born->xcomb()), DSH()->pdfs(), Born->deepHead()->xcomb()->lastParticles());
DSH()->hardScales(sqr(running));
}
Energy Merger::CKKW_StartScale(NPtr Born){
Energy res=generator()->maximumCMEnergy();;
if(!Born->children().empty()){
for (int i=0;i<Born->legsize();i++){
if (!Born->nodeME()->mePartonData()[i]->coloured())continue;
for (int j=i+1;j<Born->legsize();j++){
if (i==j||!Born->nodeME()->mePartonData()[j]->coloured())continue;
res= min(res,sqrt(2.*Born->nodeME()->lastMEMomenta()[i]*Born->nodeME()->lastMEMomenta()[j]));
}
}
}else{
Born->nodeME()->factory()->scaleChoice()->setXComb(Born->xcomb());
res= sqrt(Born->nodeME()->factory()->scaleChoice()->renormalizationScale());
}
Born->nodeME()->factory()->scaleChoice()->setXComb(Born->xcomb());
res=max(res, sqrt(Born->nodeME()->factory()->scaleChoice()->renormalizationScale()));
return res;
}
double Merger::alphasUnlops( Energy next,Energy fixedScale) {
double betaZero = (11./6.)*SM().Nc() - (1./3.)*Nf(history[0].scale);
return (betaZero*log(sqr(fixedScale/next)))+
(theCMWScheme?(((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*Nf(history[0].scale)))):0.);
}
double Merger::pdfratio(NPtr Born,Energy nominator_scale, Energy denominator_scale,int side){
StdXCombPtr bXc = Born->xcomb();
assert (bXc->mePartonData()[side]->coloured());
double from,to;
if (side==0){
if (denominator_scale==nominator_scale) {
return 1.;
}
from=Born->nodeME()->pdf1(sqr(denominator_scale));
to =Born->nodeME()->pdf1(sqr( nominator_scale ));
if ( (to < 1e-8||from < 1e-8)&&(to/from>10000000.)){cout<<"\npdfratio to="<<to<<" from="<<from;return 0.;}
}
else{
if (denominator_scale==nominator_scale) {
return 1.;
}
from=Born->nodeME()->pdf2(sqr(denominator_scale));
to=Born->nodeME()->pdf2( sqr( nominator_scale ));
if ( (to < 1e-8||from < 1e-8)&&(to/from>10000000.)){cout<<"\npdfratio to="<<to<<" from="<<from;return 0.;}
}
return to/from;
}
bool Merger::dosudakov(NPtr Born,Energy running, Energy next, double& sudakov0_n) {
CKKW_PrepareSudakov(Born, running);
for ( list<DipoleChain>::iterator chain = DSH()->eventRecord().chains().begin() ;
chain != DSH()->eventRecord().chains().end() ; chain++ ) {
for ( list<Dipole>::iterator dip = (*chain).dipoles().begin() ; dip != (*chain).dipoles().end() ; ++dip ) {
sudakov0_n*=singlesudakov( dip, next,running,make_pair(true,false) );
sudakov0_n*=singlesudakov( dip, next,running,make_pair(false,true) );
if (sudakov0_n==0.0){
cleanup(Born);
return false;
}
}
}
cleanup(Born);
return true;
}
bool Merger::doUNLOPS(NPtr Born,Energy running, Energy next,Energy fixedScale, double& UNLOPS) {
CKKW_PrepareSudakov(Born, running);
for ( list<DipoleChain>::iterator chain = DSH()->eventRecord().chains().begin() ;
chain != DSH()->eventRecord().chains().end() ; chain++ ) {
for ( list<Dipole>::iterator dip = (*chain).dipoles().begin() ; dip != (*chain).dipoles().end() ; ++dip ) {
UNLOPS+=singleUNLOPS( dip, next,running,fixedScale,make_pair(true,false) );;
UNLOPS+=singleUNLOPS( dip, next,running,fixedScale,make_pair(false,true) );
}
}
cleanup(Born);
return true;
}
bool Merger::projectorStage(NPtr Born){
return (Born->deepHead()->nodeME()->projectorStage() ==
int((Born->deepHead()->legsize()
- Born->legsize())));
}
void Merger::cleanup(NPtr Born) {
DSH()->eventRecord().clear();
if(!Born->xcomb()->subProcess())return;
ParticleVector vecfirst = Born->xcomb()->subProcess()->incoming().first->children();
for ( ParticleVector::const_iterator it = vecfirst.begin() ; it != vecfirst.end() ; it++ )
Born->xcomb()->subProcess()->incoming().first->abandonChild(*it);
ParticleVector vecsecond = Born->xcomb()->subProcess()->incoming().second->children();
for ( ParticleVector::const_iterator it = vecsecond.begin() ; it != vecsecond.end() ; it++ )
Born->xcomb()->subProcess()->incoming().second->abandonChild(*it);
Born->xcomb()->subProcess(SubProPtr());
}
double Merger::singlesudakov(list<Dipole>::iterator dip ,Energy next,Energy running,pair<bool,bool> conf ){
double res=1.;
tPPtr emitter = dip->emitter(conf);
tPPtr spectator = dip->spectator(conf);
DipoleSplittingInfo candidate((*dip).index(conf),conf,(*dip).emitterX(conf),(*dip).spectatorX(conf),emitter,spectator);
if ( DSH()->generators().find(candidate.index()) == DSH()->generators().end() ) DSH()->getGenerators(candidate.index());
pair<GeneratorMap2::iterator,GeneratorMap2::iterator> gens = DSH()->generators().equal_range(candidate.index());
for ( GeneratorMap2::iterator gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
candidate.scale(dScale);
candidate.continuesEvolving();
Energy ptMax=(*gen).second->splittingKinematics()->ptMax(candidate.scale(),candidate.emitterX(), candidate.spectatorX(),
candidate.index(),*gen->second->splittingKernel());
candidate.hardPt(min(running,ptMax));
if (candidate.hardPt()>next){
res*=gen->second->sudakov(candidate,next);
}
}
return res;
}
double Merger::singleUNLOPS(list<Dipole>::iterator dip ,Energy next,Energy running,Energy fixedScale,pair<bool,bool> conf ){
double res=0.;
tPPtr emitter = dip->emitter(conf);
tPPtr spectator = dip->spectator(conf);
DipoleSplittingInfo candidate((*dip).index(conf),conf,(*dip).emitterX(conf),(*dip).spectatorX(conf),emitter,spectator);
if ( DSH()->generators().find(candidate.index()) == DSH()->generators().end() ) DSH()->getGenerators(candidate.index());
pair<GeneratorMap2::iterator,GeneratorMap2::iterator> gens = DSH()->generators().equal_range(candidate.index());
for ( GeneratorMap2::iterator gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
candidate.scale(dScale);
candidate.continuesEvolving();
Energy ptMax=gen->second->splittingKinematics()->ptMax(candidate.scale(),candidate.emitterX(), candidate.spectatorX(),
candidate.index(),*gen->second->splittingKernel());
candidate.hardPt(min(running,ptMax));
if (candidate.hardPt()>next){
res+=gen->second->unlops(candidate,next,fixedScale);
}
}
return res;
}
void Merger::firstNodeMap(Ptr<MatchboxMEBase>::ptr a,NPtr b){theFirstNodeMap.insert(make_pair(a,b));}
map<Ptr<MatchboxMEBase>::ptr,NPtr> Merger::firstNodeMap(){return theFirstNodeMap;}
void Merger::setXComb(Ptr<MatchboxMEBase>::ptr me,tStdXCombPtr xc,int st){
theFirstNodeMap[me]->setXComb(xc, st);
}
void Merger::setKinematics(Ptr<MatchboxMEBase>::ptr me){
theFirstNodeMap[me]->setKinematics();
}
void Merger::clearKinematics(Ptr<MatchboxMEBase>::ptr me){
theFirstNodeMap[me]->clearKinematics();
}
bool Merger::generateKinematics(Ptr<MatchboxMEBase>::ptr me,const double * r){
theFirstNodeMap[me]->firstgenerateKinematics(r, 0,me->lastXCombPtr()->lastSHat());
if (theFirstNodeMap[me]->cutStage()==0 ){
bool inAlphaPS=false;
NPtrVec children=theFirstNodeMap[me]->children();
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
treefactory()->setAlphaParameter(theGamma);
inAlphaPS|=theGamma!=1.?(*child)->dipol()->aboveAlpha():false;
treefactory()->setAlphaParameter(1.);
}
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for(SafeClusterMap::iterator
it=temp.begin();
it!=temp.end();++it){
if (!it->second.first&&!inAlphaPS)return false;
}
}
if (theFirstNodeMap[me]->cutStage()==1 ){
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for(SafeClusterMap::iterator
it=temp.begin();
it!=temp.end();++it){
if (!it->second.first && !it->second.second)return false;
}
}
return true;
}
bool Merger::calculateInNode() const{
return theCalculateInNode;
}
void Merger::fillProjectors(Ptr<MatchboxMEBase>::ptr me){
for (unsigned int i = 0; i < (theFirstNodeMap[me]->Projector()).size(); ++i) {
me->lastXCombPtr()->projectors().insert(
(theFirstNodeMap[me]->Projector())[i].first,
(theFirstNodeMap[me]->Projector())[i].second->xcomb());
}
}
pair<bool,bool> Merger::clusterSafe(Ptr<MatchboxMEBase>::ptr me,int emit,int emis,int spec){
return theFirstNodeMap[me]->clusterSafe().find(make_pair(make_pair(emit,emis),spec))->second;
}
bool Merger::matrixElementRegion(PVector particles,Energy winnerScale,Energy cutscale){
//cout<<"\nparticles s"<<particles.size()<<" "<<particles[0]<<" "<<particles[1]<<flush;
/*
if (defMERegionByJetAlg && !particles[0]->coloured()&& !particles[1]->coloured()) {
assert(false);
vector<fastjet::PseudoJet> input_particles;
for(size_t em=2; em < particles.size();em++){
input_particles.push_back(fastjet::PseudoJet(particles[em]->momentum().x()/GeV,
particles[em]->momentum().y()/GeV,
particles[em]->momentum().z()/GeV,
particles[em]->momentum().e()/GeV));
}
fastjet::JetDefinition jet_def(fastjet::ee_kt_algorithm);
fastjet::ClusterSequence clust_seq(input_particles, jet_def);
size_t n = particles.size()-2;
vector<fastjet::PseudoJet> exclusive_jets = clust_seq.exclusive_jets_ycut(ee_ycut);
return n==exclusive_jets.size();
}
if (defMERegionByJetAlg ) {
assert(false);
size_t noncol=0;
vector<fastjet::PseudoJet> input_particles;
for(size_t em=2; em < particles.size();em++){
if (particles[em]->coloured())
input_particles.push_back(fastjet::PseudoJet(particles[em]->momentum().x()/GeV,
particles[em]->momentum().y()/GeV,
particles[em]->momentum().z()/GeV,
particles[em]->momentum().e()/GeV));
else
noncol++;
}
double Rparam = 1.0;
fastjet::Strategy strategy = fastjet::Best;
fastjet::RecombinationScheme recomb_scheme = fastjet::E_scheme;
fastjet::JetDefinition jet_def(fastjet::kt_algorithm, Rparam, recomb_scheme, strategy);
fastjet::ClusterSequence clust_seq(input_particles, jet_def);
size_t n = particles.size()-2-noncol;
vector<fastjet::PseudoJet> exclusive_jets = clust_seq.exclusive_jets(pp_dcut);
// cout<<"\nn= "<<n<<" jets= "<<exclusive_jets.size();
//for (size_t i=0; i<exclusive_jets.size(); i++) {
//cout<<"\nn= "<<n<<" jetspt= "<<exclusive_jets[i].perp();
//}
return n==exclusive_jets.size();
}
*/
//assert(false);
Energy ptx=1000000.*GeV;
bool foundwinnerpt=false;
//FF
for(size_t em=2; em < particles.size();em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==emm) continue;
for(size_t spe=2; spe < particles.size();spe++){
if (!particles[spe]->coloured()) continue;
if (em==spe||emm==spe) continue;
if (!(particles[em]->id()==-particles[emm]->id()||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Energy pt=0*GeV;
- if (emittermom.m()==0.001*GeV&&emissionmom.m()==0.001*GeV&&spectatormom.m()==0.001*GeV) {
+ if (emittermom.m()<=0.001*GeV&&emissionmom.m()<=0.001*GeV&&spectatormom.m()<=0.001*GeV) {
pt=FFLTK->lastPt(emittermom,emissionmom,spectatormom);
}else{
pt=FFMTK->lastPt(emittermom,emissionmom,spectatormom);
}
//cout<<"\npt "<<pt/GeV<<" "<<winnerScale/GeV;
if (abs(pt-winnerScale)<0.001*GeV) {
foundwinnerpt=true;
}
// if(scale * sqrt(y*z*(1.-z))<optVeto&&winnerScale>optVeto)cout<<"\nFF "<<(scale * sqrt(y*z*(1.-z))/GeV);
ptx =min(ptx,pt);
}
}
}
//FI
for(size_t spe=0; spe < 2;spe++){
if (!particles[spe]->coloured()) continue;
for(size_t em=2; em < particles.size();em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==emm) continue;
if (!(particles[em]->id()==-particles[emm]->id()||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Energy pt=0*GeV;
- if (emittermom.m()==0.001*GeV&&emissionmom.m()==0.001*GeV&&spectatormom.m()==0.001*GeV) {
+ if (emittermom.m()<=0.001*GeV&&emissionmom.m()<=0.001*GeV&&spectatormom.m()<=0.001*GeV) {
pt=FILTK->lastPt(emittermom,emissionmom,spectatormom);
}else{
pt=FIMTK->lastPt(emittermom,emissionmom,spectatormom);
}
if (abs(pt-winnerScale)<0.001*GeV) {
foundwinnerpt=true;
}
if(pt>0.*GeV)
ptx =min(ptx,pt);
}
}
}
//IF
for(size_t em=0; em < 2;em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
for(size_t spe=2; spe < particles.size();spe++){
if (!particles[spe]->coloured()) continue;
if (emm==spe) continue;
if (!(particles[em]->id()>6|| particles[em]->id()==particles[emm]->id() ||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Energy pt=0*GeV;
if (emittermom.m()<=0.001*GeV&&emissionmom.m()<=0.001*GeV&&spectatormom.m()<=0.001*GeV) {
pt=IFLTK->lastPt(emittermom,emissionmom,spectatormom);
}else{
pt=IFMTK->lastPt(emittermom,emissionmom,spectatormom);
}
if (abs(pt-winnerScale)<0.01*GeV) {
foundwinnerpt=true;
}
ptx =min(ptx,pt);
}
}
}
//II
for(size_t em=0; em < 2;em++){
if (!particles[em]->coloured()) continue;
for(size_t spe=0; spe < 2;spe++){
if (!particles[spe]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==spe) continue;
if (!(particles[em]->id()>6||
particles[em]->id()==particles[emm]->id() ||
particles[emm]->id()>6))continue;
//assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Energy pt=IILTK->lastPt(emittermom,emissionmom,spectatormom);
if (abs(pt-winnerScale)<0.01*GeV) {
foundwinnerpt=true;
}
ptx =min(ptx, pt);
}
}
}
// cout<<"\n"<<cutscale/GeV<< " "<<foundwinnerpt<<" "<<ptx/GeV;
if(!foundwinnerpt){
cout<<"\nWarning: could not find winner with pt: "<<winnerScale/GeV;
for(size_t emm=0; emm < particles.size();emm++){
cout<<"\n"<<particles[emm]->id()<<" "<<particles[emm]->momentum()/GeV<<" "<<particles[emm]->momentum().m()/GeV<<flush;
}
}
return (ptx>cutscale) ;
}
int Merger::M()const{return theTreeFactory->M();}
int Merger::N()const{return theTreeFactory->N();}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void Merger::persistentOutput(PersistentOStream & os) const {
os << minusL<< Unlopsweights<<
theCMWScheme<< projected<<
isUnitarized<< isNLOUnitarized<<
defMERegionByJetAlg<<theOpenInitialStateZ<<
theChooseHistory<<
theN0<< theOnlyN<<
xiRenME<< xiFacME<<
xiRenSh<< xiFacSh<<
xiQSh<< weight<<weightCB<<theGamma<<ee_ycut<<pp_dcut<<theSmearing<<ounit(therenormscale, GeV)<<ounit(theIRSafePT, GeV)<<ounit(theMergePt, GeV)<<ounit(theCentralMergePt, GeV)<<theMergingJetFinder
<<theLargeNBasis
<< FFLTK
<< FILTK
<< IFLTK
<< IILTK
<< FFMTK
<< FIMTK
<< IFMTK
<<theDipoleShowerHandler<< theTreeFactory<<theFirstNodeMap ;
}
void Merger::persistentInput(PersistentIStream & is, int) {
is >> minusL>> Unlopsweights>>
theCMWScheme>> projected>>
isUnitarized>> isNLOUnitarized>>
defMERegionByJetAlg>>theOpenInitialStateZ>>
theChooseHistory>>
theN0>> theOnlyN>>
xiRenME>> xiFacME>>
xiRenSh>> xiFacSh>>
xiQSh>>
weight>>weightCB>>theGamma>>ee_ycut>>pp_dcut>>theSmearing>>iunit(therenormscale, GeV)>>iunit(theIRSafePT, GeV)>>iunit(theMergePt, GeV)>>iunit(theCentralMergePt, GeV)>>theMergingJetFinder>>theLargeNBasis>>
FFLTK
>> FILTK
>> IFLTK
>> IILTK
>> FFMTK
>> FIMTK
>> IFMTK>>
theDipoleShowerHandler>> theTreeFactory>>theFirstNodeMap ;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<Merger, Herwig::MergerBase>
describeHerwigMerger("Herwig::Merger", "HwDipoleShower.so");
//ClassDescription<Merger> Merger::initMerger;
// Definition of the static class description member.
void Merger::Init() {
static ClassDocumentation<Merger> documentation
("Merger.");
static Reference<Merger,DipoleShowerHandler> interfaceShowerHandler
("DipoleShowerHandler",
"",
&Merger::theDipoleShowerHandler, false, false, true, true, false);
static Switch<Merger,bool>
interfaceminusL ("minusL","",&Merger::minusL, false, false, false);
static SwitchOption interfaceminusLYes
(interfaceminusL,"Yes","",true);
static SwitchOption interfaceminusLNo
(interfaceminusL,"No","",false);
static Switch<Merger,bool>
interfaceUnlopsweights ("Unlopsweights","",&Merger::Unlopsweights, false, false, false);
static SwitchOption interfaceUnlopsweightsYes
(interfaceUnlopsweights,"Yes","",true);
static SwitchOption interfaceUnlopsweightsNo
(interfaceUnlopsweights,"No","",false);
static Switch<Merger,bool>
interfacetheCMWScheme ("CMWScheme","",&Merger::theCMWScheme, false, false, false);
static SwitchOption interfacetheCMWSchemeYes
(interfacetheCMWScheme,"Yes","",true);
static SwitchOption interfacetheCMWSchemeNo
(interfacetheCMWScheme,"No","",false);
static Parameter<Merger,Energy> interfaceMergerScale
("MergingScale",
"The MergingScale.",
&Merger::theCentralMergePt, GeV, 20.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Reference<Merger,MFactory> interfaceMergerHelper
("MFactory",
"",
&Merger::theTreeFactory, false, false, true, true, false);
static Parameter<Merger,double> interfacedcut
("pp_dcut",
"The MergingScale.",
&Merger::pp_dcut, 5.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<Merger,double> interfaceycut
("ee_ycut",
"The MergingScale.",
&Merger::ee_ycut, 0.2, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<Merger,double> interfacegamma
("gamma",
"gamma parameter.",
&Merger::theGamma, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Reference<Merger,JetFinder> interfaceMergingJetFinder
("MergingJetFinder",
"A reference to the jet finder used in Merging.",
&Merger::theMergingJetFinder, false, false, true, false, false);
static Reference<Merger,ColourBasis> interfaceLargeNBasis
("LargeNBasis",
"Set the large-N colour basis implementation.",
&Merger::theLargeNBasis, false, false, true, true, false);
static Switch<Merger,bool>
interfacedefMERegionByJetAlg
("MERegionByJetAlg","",&Merger::defMERegionByJetAlg, false, false, false);
static SwitchOption interfacedefMERegionByJetAlgYes
(interfacedefMERegionByJetAlg,"Yes","",true);
static SwitchOption interfacedefMERegionByJetAlgNo
(interfacedefMERegionByJetAlg,"No","",false);
static Switch<Merger,bool>
interfaceOpenInitialSateZ
("OpenInitialStateZ","",&Merger::theOpenInitialStateZ, false, false, false);
static SwitchOption interfaceOpenInitialSateZYes
(interfaceOpenInitialSateZ,"Yes","",true);
static SwitchOption interfaceOpenInitialSateZNo
(interfaceOpenInitialSateZ,"No","",false);
static Parameter<Merger, Energy>
interfaceIRSafePT
("IRSafePT", "Set the pt for which a matrixelement should be treated as IR-safe.",
&Merger::theIRSafePT,
GeV, 0.0 * GeV, ZERO, Constants::MaxEnergy, true, false, Interface::limited);
interfaceIRSafePT.setHasDefault(false);
static Parameter<Merger, double> interfacemergePtsmearing("MergingScaleSmearing", "Set the percentage the merging pt should be smeared.",
&Merger::theSmearing, 0., 0.,
0.0, 0.5, true, false, Interface::limited);
static Parameter<Merger, int> interfacechooseHistory("chooseHistory", "different ways of choosing the history", &Merger::theChooseHistory, 3, -1, 0,
false, false, Interface::lowerlim);
}
diff --git a/src/Matchbox/Merging.in b/src/Matchbox/Merging.in
--- a/src/Matchbox/Merging.in
+++ b/src/Matchbox/Merging.in
@@ -1,236 +1,236 @@
cd /Herwig/Merging
clear /Herwig/EventHandlers/EventHandler:SubProcessHandlers
insert /Herwig/EventHandlers/EventHandler:SubProcessHandlers[0] MFactory
set /Herwig/Cuts/MatchboxJetMatcher:Factory MFactory
set /Herwig/MatrixElements/Matchbox/MECorrectionHandler:Factory MFactory
set MFactory:Cuts /Herwig/Cuts/Cuts
set MFactory:ScaleChoice /Herwig/Merging/MScale
set Merger:MFactory MFactory
set Merger:DipoleShowerHandler /Herwig/DipoleShower/DipoleShowerHandler
set /Herwig/DipoleShower/DipoleShowerHandler:MergingHelper Merger
set MScale:MergingHelper Merger
set Merger:MergingJetFinder /Herwig/Cuts/JetFinder
set MFactory:MergingHelper Merger
set Merger:minusL No
cd /Herwig/Analysis
create ThePEG::RivetAnalysis Rivet RivetAnalysis.so
set /Herwig/Shower/ShowerHandler:HardScaleProfile NULL
set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleProfile NULL
set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF
set /Herwig/Couplings/NLOAlphaS:min_active_flavours 3
set /Herwig/Couplings/NLOAlphaS:max_active_flavours 5
set /Herwig/Couplings/NLOAlphaS:input_alpha_s 0.118
set /Herwig/Couplings/NLOAlphaS:input_scale 91.18*GeV
set /Herwig/Model:QCD/RunningAlphaS /Herwig/Couplings/NLOAlphaS
set /Herwig/DipoleShower/DipoleShowerHandler:GlobalAlphaS /Herwig/Couplings/NLOAlphaS
cd /Herwig/Generators
set EventGenerator:EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler
set /Herwig/DipoleShower/DipoleShowerHandler:MaxPtIsMuF Yes
create Herwig::MonacoSampler /Herwig/Samplers/Monaco
set /Herwig/Samplers/Sampler:BinSampler /Herwig/Samplers/Monaco
cd /Herwig/Particles
do c:UnsetHardProcessMass
do cbar:UnsetHardProcessMass
-do b:UnsetHardProcessMass
-do bbar:UnsetHardProcessMass
+#do b:UnsetHardProcessMass
+#do bbar:UnsetHardProcessMass
set c:NominalMass 0*GeV
set cbar:NominalMass 0*GeV
##################################################
## intrinsic pt
##################################################
set /Herwig/DipoleShower/IntrinsicPtGenerator:ValenceIntrinsicPtScale 2.0*GeV
set /Herwig/DipoleShower/IntrinsicPtGenerator:SeaIntrinsicPtScale 2.0*GeV
##################################################
## Dipole shower tune settings
##################################################
set /Herwig/Particles/g:ConstituentMass 0.7899848*GeV
##################################################
## cutoffs
##################################################
set /Herwig/DipoleShower/Kinematics/FFLightKinematics:IRCutoff 0.45*GeV
set /Herwig/DipoleShower/Kinematics/FFMassiveKinematics:IRCutoff 0.45*GeV
set /Herwig/DipoleShower/Kinematics/FILightKinematics:IRCutoff 0.45*GeV
set /Herwig/DipoleShower/Kinematics/FIMassiveKinematics:IRCutoff 0.45*GeV
set /Herwig/DipoleShower/Kinematics/IFLightKinematics:IRCutoff 0.45*GeV
set /Herwig/DipoleShower/Kinematics/IFMassiveKinematics:IRCutoff 0.45*GeV
set /Herwig/DipoleShower/Kinematics/IILightKinematics:IRCutoff 0.45*GeV
## needs to be synchronized with settings above
set /Herwig/MatrixElements/Matchbox/DipoleMatching:FFPtCut 0.45*GeV
set /Herwig/MatrixElements/Matchbox/DipoleMatching:FIPtCut 0.45*GeV
set /Herwig/MatrixElements/Matchbox/DipoleMatching:IIPtCut 0.45*GeV
##################################################
## hadronization parameters -- complete?
##################################################
cd /Herwig/Hadronization
set ClusterFissioner:ClMaxLight 3.30
set ClusterFissioner:ClPowLight 2.50
set ClusterFissioner:PSplitLight 1.29
set ClusterDecayer:ClDirLight 1
set ClusterDecayer:ClSmrLight 3.118342
set ClusterFissioner:ClMaxCharm 3.11*GeV
set ClusterFissioner:ClPowCharm 1.62
set ClusterFissioner:PSplitCharm 2.54
set ClusterDecayer:ClDirCharm 1
set ClusterDecayer:ClSmrCharm 0.
set HadronSelector:SingleHadronLimitCharm 0.0
set ClusterFissioner:ClMaxBottom 1.38*GeV
set ClusterFissioner:ClPowBottom 0.62
set ClusterFissioner:PSplitBottom 0.20
set ClusterDecayer:ClDirBottom 1
set ClusterDecayer:ClSmrBottom 0.
set HadronSelector:SingleHadronLimitBottom 0.12
set HadronSelector:PwtUquark 1.0
set HadronSelector:PwtDquark 1.0
set HadronSelector:PwtSquark 1.09
set HadronSelector:PwtCquark 1.0
set HadronSelector:PwtBquark 1.0
set HadronSelector:PwtDIquark 0.66
set HadronSelector:SngWt 1.0
set HadronSelector:DecWt 1.0
cd /
cd /Herwig/Samplers
set Monaco:EnhancementFactor 1.2
set Monaco:InitialPoints 5000
set Monaco:LuminosityMapperBins 8
set Monaco:NIterations 4
set Monaco:RemapChannelDimension Yes
set Monaco:RemapperMinSelection 0.0001
set Monaco:RemapperPoints 10000
set Monaco:UseAllIterations No
set Sampler:UpdateAfter 1000
set Sampler:AddUpSamplers Off
set Sampler:GlobalMaximumWeight Off
set Sampler:FlatSubprocesses Off
set Sampler:MinSelection 0.000001
set Sampler:AlmostUnweighted On
set Sampler:BinSampler:Kappa 1.
set Sampler:RunCombinationData Off
set Sampler:Verbose On
#######
# CMW scheme
#######
cd /Herwig/DipoleShower/Kernels
set FFgx2ggxDipoleKernel:CMWScheme On
set FFqx2qgxDipoleKernel:CMWScheme On
set FFgx2ddxDipoleKernel:CMWScheme On
set FFgx2uuxDipoleKernel:CMWScheme On
set FFgx2ccxDipoleKernel:CMWScheme On
set FFgx2ssxDipoleKernel:CMWScheme On
set FFgx2bbxDipoleKernel:CMWScheme On
set FFMgx2ggxDipoleKernel:CMWScheme On
set FFMdx2dgxDipoleKernel:CMWScheme On
set FFMux2ugxDipoleKernel:CMWScheme On
set FFMcx2cgxDipoleKernel:CMWScheme On
set FFMsx2sgxDipoleKernel:CMWScheme On
set FFMbx2bgxDipoleKernel:CMWScheme On
set FFMtx2tgxDipoleKernel:CMWScheme On
set FFMgx2ddxDipoleKernel:CMWScheme On
set FFMgx2uuxDipoleKernel:CMWScheme On
set FFMgx2ccxDipoleKernel:CMWScheme On
set FFMgx2ssxDipoleKernel:CMWScheme On
set FFMgx2bbxDipoleKernel:CMWScheme On
set FIgx2ggxDipoleKernel:CMWScheme On
set FIqx2qgxDipoleKernel:CMWScheme On
set FIgx2ddxDipoleKernel:CMWScheme On
set FIgx2uuxDipoleKernel:CMWScheme On
set FIgx2ccxDipoleKernel:CMWScheme On
set FIgx2ssxDipoleKernel:CMWScheme On
set FIgx2bbxDipoleKernel:CMWScheme On
set FIMdx2dgxDipoleKernel:CMWScheme On
set FIMux2ugxDipoleKernel:CMWScheme On
set FIMcx2cgxDipoleKernel:CMWScheme On
set FIMsx2sgxDipoleKernel:CMWScheme On
set FIMbx2bgxDipoleKernel:CMWScheme On
set FIMtx2tgxDipoleKernel:CMWScheme On
set FIMgx2ddxDipoleKernel:CMWScheme On
set FIMgx2uuxDipoleKernel:CMWScheme On
set FIMgx2ccxDipoleKernel:CMWScheme On
set FIMgx2ssxDipoleKernel:CMWScheme On
set FIMgx2bbxDipoleKernel:CMWScheme On
#set FIMgx2ttxDipoleKernel:CMWScheme On
set IFgx2ggxDipoleKernel:CMWScheme On
set IFqx2qgxDipoleKernel:CMWScheme On
set IFqx2gqxDipoleKernel:CMWScheme On
set IFgx2ddbarxDipoleKernel:CMWScheme On
set IFgx2dbardxDipoleKernel:CMWScheme On
set IFgx2uubarxDipoleKernel:CMWScheme On
set IFgx2ubaruxDipoleKernel:CMWScheme On
set IFgx2ccbarxDipoleKernel:CMWScheme On
set IFgx2cbarcxDipoleKernel:CMWScheme On
set IFgx2ssbarxDipoleKernel:CMWScheme On
set IFgx2sbarsxDipoleKernel:CMWScheme On
set IFMgx2ggxDipoleKernel:CMWScheme On
set IFMqx2qgxDipoleKernel:CMWScheme On
set IFMqx2gqxDipoleKernel:CMWScheme On
set IFMgx2ddbarxDipoleKernel:CMWScheme On
set IFMgx2dbardxDipoleKernel:CMWScheme On
set IFMgx2uubarxDipoleKernel:CMWScheme On
set IFMgx2ubaruxDipoleKernel:CMWScheme On
set IFMgx2ccbarxDipoleKernel:CMWScheme On
set IFMgx2cbarcxDipoleKernel:CMWScheme On
set IFMgx2ssbarxDipoleKernel:CMWScheme On
set IFMgx2sbarsxDipoleKernel:CMWScheme On
set IIgx2ggxDipoleKernel:CMWScheme On
set IIqx2qgxDipoleKernel:CMWScheme On
set IIqx2gqxDipoleKernel:CMWScheme On
set IIgx2ddbarxDipoleKernel:CMWScheme On
set IIgx2dbardxDipoleKernel:CMWScheme On
set IIgx2uubarxDipoleKernel:CMWScheme On
set IIgx2ubaruxDipoleKernel:CMWScheme On
set IIgx2ccbarxDipoleKernel:CMWScheme On
set IIgx2cbarcxDipoleKernel:CMWScheme On
set IIgx2ssbarxDipoleKernel:CMWScheme On
set IIgx2sbarsxDipoleKernel:CMWScheme On
diff --git a/src/Merging/3LO-2NLO-LEP.in b/src/Merging/3LO-2NLO-LEP.in
--- a/src/Merging/3LO-2NLO-LEP.in
+++ b/src/Merging/3LO-2NLO-LEP.in
@@ -1,138 +1,138 @@
# -*- ThePEG-repository -*-
##################################################
## Herwig/Merging example input file
##################################################
##################################################
## Collider type
##################################################
read Matchbox/EECollider.in
read Matchbox/Merging.in
##################################################
## Beam energy sqrt(s)
##################################################
cd /Herwig/EventHandlers
set EventHandler:LuminosityFunction:Energy 91.2*GeV
##################################################
## Process selection
##################################################
## Note that event generation may fail if no matching matrix element has
## been found. Coupling orders are with respect to the Born process,
## i.e. NLO QCD does not require an additional power of alphas.
## Model assumptions
read Matchbox/StandardModelLike.in
read Matchbox/DiagonalCKM.in
## Set the order of the couplings
cd /Herwig/Merging
set MFactory:OrderInAlphaS 0
set MFactory:OrderInAlphaEW 2
## Select the process
## You may use identifiers such as p, pbar, j, l, mu+, h0 etc.
do MFactory:Process e- e+ -> j j [ j j ]
set MFactory:NLOProcesses 2
set Merger:MergingScale 4.*GeV
set Merger:MergingScaleSmearing 0.1
set Merger:gamma 1.
cd /Herwig/MatrixElements/Matchbox/Utility
insert DiagramGenerator:ExcludeInternal 0 /Herwig/Particles/gamma
## Special settings required for on-shell production of unstable particles
## enable for on-shell top production
# read Matchbox/OnShellTopProduction.in
## enable for on-shell W, Z or h production
# read Matchbox/OnShellWProduction.in
# read Matchbox/OnShellZProduction.in
# read Matchbox/OnShellHProduction.in
##################################################
## Matrix element library selection
##################################################
## Select a generic tree/loop combination or a
## specialized NLO package
# read Matchbox/MadGraph-GoSam.in
# read Matchbox/MadGraph-MadGraph.in
# read Matchbox/MadGraph-NJet.in
-# read Matchbox/MadGraph-OpenLoops.in
-
+read Matchbox/MadGraph-OpenLoops.in
+set MadGraph:ProcessPath eejjjj
## Uncomment this to use ggh effective couplings
## currently only supported by MadGraph-GoSam
# read Matchbox/HiggsEffective.in
##################################################
## Cut selection
## See the documentation for more options
##################################################
## cuts on additional jets
# read Matchbox/DefaultEEJets.in
# set NJetsCut:NJetsMin 3
##################################################
## Scale choice
## See the documentation for more options
##################################################
cd /Herwig/MatrixElements/Matchbox/Scales/
set /Herwig/Merging/MScale:ScaleChoice SHatScale
##################################################
## Scale uncertainties
##################################################
# read Matchbox/MuDown.in
# read Matchbox/MuUp.in
##################################################
## Shower scale uncertainties
##################################################
# read Matchbox/MuQDown.in
# read Matchbox/MuQUp.in
##################################################
## Analyses
##################################################
cd /Herwig/Analysis
insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 Rivet
# insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 HepMC
read Matchbox/LEP91-Analysis.in
##################################################
## Do not apply profile scales for LEP as hard
## scale coincides with kinematic limit
##################################################
set /Herwig/Shower/ShowerHandler:HardScaleProfile NULL
set /Herwig/DipoleShower/DipoleShowerHandler:HardScaleProfile NULL
##################################################
## Save the generator
##################################################
do /Herwig/Merging/MFactory:ProductionMode
cd /Herwig/Generators
saverun 3LO-2NLO-LEP EventGenerator
diff --git a/src/Merging/Makefile.am b/src/Merging/Makefile.am
--- a/src/Merging/Makefile.am
+++ b/src/Merging/Makefile.am
@@ -1,59 +1,59 @@
BUILT_SOURCES = done-all-links
Matchboxdir = ${pkgdatadir}/Merging
INPUTFILES = \
1LO-0NLO-LEP.in \
1LO-0NLO-LEP_Qt.in \
1LO-0NLO-LEP_Matchbox.in \
1LO-0NLO-LEP_FixedOrder.in \
1LO-0NLO-LHC-W.in \
1LO-0NLO-LHC-Z.in \
2LO-0NLO-LEP.in \
2LO-0NLO-LHC-W.in \
2LO-0NLO-LHC-Z.in \
2LO-1NLO-LEP.in \
2LO-1NLO-LEP_Matchbox.in \
2LO-1NLO-LHC-Z.in \
2LO-1NLO-LHC-W.in \
3LO-0NLO-LEP.in \
3LO-0NLO-LEP_only0.in \
3LO-0NLO-LEP_only1.in \
3LO-0NLO-LEP_only2.in \
3LO-0NLO-LHC-W.in \
3LO-0NLO-LHC-Z.in \
3LO-1NLO-LEP.in \
3LO-1NLO-LHC-Z.in \
3LO-1NLO-LHC-W.in \
3LO-2NLO-LEP.in \
3LO-2NLO-LEP_light.in \
3LO-2NLO-LHC-Z.in \
3LO-2NLO-LHC-W.in \
4LO-0NLO-LHC-Z.in \
1LO-0NLO-LHC-Z_OISZ.in \
2LO-0NLO-LHC-Z_OISZ.in \
3LO-0NLO-LHC-Z_OISZ.in \
1LO-0NLO-LHC-J.in \
1LO-0NLO-LHC-T.in \
1LO-0NLO-LHC-T_QtDecay.in \
1LO-0NLO-LHC-J_OISZ.in \
2LO-0NLO-LHC-J.in \
2LO-0NLO-LHC-J_only0.in \
2LO-0NLO-LHC-J_only1.in \
2LO-0NLO-LHC-T.in \
2LO-0NLO-LHC-T_QtDecay.in \
3LO-0NLO-LHC-T.in \
-2LO-0NLO-LHC-J_OISZ.in \
-2LO-0NLO-LHC-J_OISZIC.in
+2LO-0NLO-LHC-J_OISZ.in
+#2LO-0NLO-LHC-J_OISZIC.in
dist_Matchbox_DATA = $(INPUTFILES)
CLEANFILES = done-all-links
done-all-links: $(INPUTFILES)
@echo "Linking input files"
@for i in $(INPUTFILES); do \
if test -f $(srcdir)/$$i -a ! -e $$i; then \
$(LN_S) -f $(srcdir)/$$i; fi; done
@touch done-all-links
diff --git a/src/Merging/test/Makefile b/src/Merging/test/Makefile
--- a/src/Merging/test/Makefile
+++ b/src/Merging/test/Makefile
@@ -1,67 +1,68 @@
#########################################################
# #
# usage: #
# #
# make NAME=LEP JOBS=4 LO=3 NLO=0 POINTS=10000 #
# #
# possible NAME: LEP, LHC-Z, LHC-W #
# #
#########################################################
NAME=default
JOBS=4
NLO=0
LO=2
POINTS=10000
all: $(NAME)
build$(NAME):
../../Herwig build Merging/$(LO)LO-$(NLO)NLO-$(NAME).in -y$(JOBS)
ls
integrate$(NAME): build$(NAME)
echo "$(shell for i in $$(seq 0 $$(( $(JOBS)-1 )) );do ../../Herwig integrate --jobid=$$i $(LO)LO-$(NLO)NLO-$(NAME).run & done)"
ls
run$(NAME): integrate$(NAME)
../../Herwig run $(LO)LO-$(NLO)NLO-$(NAME).run -N$(POINTS) -j$(JOBS)
yodamerge $(LO)LO-$(NLO)NLO-$(NAME)-*yoda > $(LO)LO-$(NLO)NLO-$(NAME).yoda
mergeyodas$(NAME): run$(NAME)
yodamerge $(LO)LO-$(NLO)NLO-$(NAME)-*yoda > $(LO)LO-$(NLO)NLO-$(NAME).yoda
mkdir -p $(LO)LO-$(NLO)NLO-$(NAME)
mkdir -p tars
tar$(NAME): mergeyodas$(NAME)
echo "$(shell mv $$(ls|grep "$(LO)LO-$(NLO)NLO-$(NAME)"|grep -v "$(LO)LO-$(NLO)NLO-$(NAME)_"| grep "\." |grep -v $(LO)LO-$(NLO)NLO-$(NAME).yoda) $(LO)LO-$(NLO)NLO-$(NAME))"
tar -czvf $(LO)LO-$(NLO)NLO-$(NAME).tar.gz $(LO)LO-$(NLO)NLO-$(NAME)
mv $(LO)LO-$(NLO)NLO-$(NAME).tar.gz tars
$(NAME): build$(NAME) integrate$(NAME) run$(NAME) mergeyodas$(NAME) tar$(NAME)
rm -fr $(LO)LO-$(NLO)NLO-$(NAME)
plots:
rivet-mkhtml *LEP.yoda *LEP-GAMMA-0.2.yoda -o LEP
rivet-mkhtml *LHC-Z.yoda *LHC-Z_OISZ.yoda -o LHC-Z
rivet-mkhtml *LHC-W.yoda -o LHC-W
rivet-mkhtml *LHC-J.yoda *LHC-J_OISZ.yoda -o LHC-J
help:
@echo "usage:"
@echo "make NAME=LEP JOBS=4 LO=3 NLO=0 POINTS=10000"
@echo ""
- @ls ../|grep NLO|cut -d"." -f1
+ @ls ../
+#|grep NLO|cut -d"." -f1

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 1:23 PM (19 h, 49 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022940
Default Alt Text
(110 KB)

Event Timeline