diff --git a/MatrixElement/Matchbox/Base/MergerBase.h b/MatrixElement/Matchbox/Base/MergerBase.h --- a/MatrixElement/Matchbox/Base/MergerBase.h +++ b/MatrixElement/Matchbox/Base/MergerBase.h @@ -1,104 +1,100 @@ // -*- C++ -*- // // MergerBase.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_MergerBase_H #define HERWIG_MergerBase_H // // This is the declaration of the MergerBase class. // #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/MatrixElement/MEBase.h" #include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h" #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXComb.h" #include "MatchboxMEBase.fh" namespace Herwig { using namespace ThePEG; class MergerBase; ThePEG_DECLARE_POINTERS(MergerBase,MergerBasePtr); /** * \ingroup Matchbox * \author Johannes Bellm * * \brief MergerBase is the base class MergingHelpers. * * @see \ref MergerBaseInterfaces "The interfaces" * defined for MergerBase. */ class MergerBase: public HandlerBase { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ MergerBase(); /** * The destructor. */ virtual ~MergerBase(); /// define the ME region for a particle vector. virtual bool matrixElementRegion(PVector incoming, PVector outcoming, Energy winnerScale, Energy cutscale) const = 0; /// cross section of as given by the merging virtual CrossSection MergingDSigDR() = 0; /// set the current xcomb, called from ME virtual void setXComb( tStdXCombPtr) = 0; /// set the current ME virtual void setME(Ptr::ptr) = 0; /// set kinematics, called from ME virtual void setKinematics() = 0; /// clear kinematics, called from ME virtual void clearKinematics() = 0; /// generate kinematics, called from ME virtual bool generateKinematics( const double * ) = 0; /** * flush all chaches of the subleading nodes. * Note: this is called not from the ME but before the first * kinematics is generated. **/ virtual void flushCaches() = 0; /// return the current maximum legs, the shower should veto virtual size_t maxLegs() const = 0; /// return the current merging scale, virtual Energy mergingScale() const = 0; - /// Potential additional emission probability for unitarising LO contributions. - virtual double emissionProbability() const = 0; - /// set the potential additional emission probability. - virtual void setEmissionProbability( double ) = 0; //@} public: /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); }; } #endif /* HERWIG_MergerBase_H */ diff --git a/Shower/Dipole/DipoleShowerHandler.cc b/Shower/Dipole/DipoleShowerHandler.cc --- a/Shower/Dipole/DipoleShowerHandler.cc +++ b/Shower/Dipole/DipoleShowerHandler.cc @@ -1,1918 +1,1927 @@ // -*- C++ -*- // // DipoleShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the DipoleShowerHandler class. // #include #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/ParVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" // include theses to have complete types #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.h" #include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h" #include using namespace Herwig; bool DipoleShowerHandler::firstWarn = true; DipoleShowerHandler::DipoleShowerHandler() : ShowerHandler(), chainOrderVetoScales(true), nEmissions(0), discardNoEmissions(false), firstMCatNLOEmission(false), thePowhegDecayEmission(true), //theAnalyseSpinCorrelations(false), realignmentScheme(0), doSubleadingNc(false),subleadingNcEmissionsLimit(0), densityOperatorEvolution(0),densityOperatorCutoff(1.0*GeV2), doPartialUnweightingAtEmission(false), doPartialUnweighting(false),referenceWeight(0.1), cmecReweightFactor(1.0),negCMECScaling(1.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), theInputColouredOffShellInShower(), theZBoundaries(1) {} DipoleShowerHandler::~DipoleShowerHandler() {} IBPtr DipoleShowerHandler::clone() const { return new_ptr(*this); } IBPtr DipoleShowerHandler::fullclone() const { return new_ptr(*this); } void DipoleShowerHandler::cascade(tPVector ) { throw Exception() << "DipoleShowerHandler: Dipoleshower not implemented as second shower." << "Check your setup or contact Herwig authors." << Exception::runerror; } tPPair DipoleShowerHandler::cascade(tSubProPtr sub, XCombPtr, Energy optHardPt, Energy optCutoff) { useMe(); prepareCascade(sub); resetWeights(); if ( !doFSR() && ! doISR() ) return sub->incoming(); eventRecord().setSubleadingNc(doSubleadingNc, subleadingNcEmissionsLimit); eventRecord().clear(); eventRecord().prepare(sub,dynamic_ptr_cast(lastXCombPtr()),newStep(),pdfs(), ShowerHandler::currentHandler()->generator()->currentEvent()->incoming(), firstInteraction(), offShellPartons(), !doSubleadingNc); if ( doSubleadingNc ) { if ( !theSplittingReweight ) { throw Exception() << "No splitting reweight was found. " << "A ColourMatrixElementCorrection " << "splitting reweight is required " << "for the subleading colour shower." << Exception::runerror; } //Set the evolution scheme for the density operator eventRecord().setDensityOperatorEvolution( densityOperatorEvolution, densityOperatorCutoff ); //Set the CMEC reweight factor theSplittingReweight->reweightFactor(cmecReweightFactor); theSplittingReweight->negativeScaling(negCMECScaling); theSplittingReweight->updateCurrentHandler(); } // SW: Removed simple test on doFSR and doISR and moved // here to account for the case of a hard event involving // no coloured particles but with unstable outgoing particles if ( !doFSR() && ! doISR() && eventRecord().decays().empty() ) return sub->incoming(); if ( !doISR() && eventRecord().outgoing().empty() && eventRecord().decays().empty() ) return sub->incoming(); if ( !doFSR() && !eventRecord().incoming().first->coloured() && !eventRecord().incoming().second->coloured() && eventRecord().decays().empty() ) return sub->incoming(); nTries = 0; // Clear the vertex record for spin correlations if ( spinCorrelations() ) //|| theAnalyseSpinCorrelations ) vertexRecord().clear(); 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() ) { LorentzRotation rot = intrinsicPtGenerator->kick(eventRecord().incoming(), eventRecord().intermediates()); eventRecord().transform(rot); } } didRealign = realign(); constituentReshuffle(); // backup subleading switch if decays fail bool doneSubleadingNc = doSubleadingNc; // subleading N can't handle decays doSubleadingNc = false; try { // Decay and shower any particles that require decaying while ( !eventRecord().decays().empty() ) { map::const_iterator decayIt = eventRecord().decays().begin(); if ( eventRecord().nextDecay() ) { decayIt = eventRecord().decays().find(eventRecord().nextDecay() ); } else { // 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 = !( nEmissions && nEmitted==nEmissions) ? thePowhegDecayEmission : false; // 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(), offShellPartons()); // 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 ( auto & ch : eventRecord().chains()) { for ( auto & dip : ch.dipoles()) { assert ( showerScale > ZERO ); dip.leftScale( showerScale ); dip.rightScale( showerScale ); } } // Prepare vertex record for spin correlations in decay shower if ( spinCorrelations() ) vertexRecord().prepareParticleDecay(incoming); // Perform the cascade doCascade(nEmitted,optHardPt,optCutoff,true); if ( spinCorrelations() ) vertexRecord().updateParticleDecay(); // Do the constituent mass shell reshuffling decayConstituentReshuffle(eventRecord().currentDecay()); } // Update the decays, adding any decays and updating momenta eventRecord().updateDecays(eventRecord().currentDecay()); eventRecord().decays().erase(decayIt); } } catch(...) { // reset flag doSubleadingNc = doneSubleadingNc; throw; } doSubleadingNc = doneSubleadingNc; break; } catch (RedoShower&) { resetWeights(); if ( ++nTries > maxtry() ) throw ShowerTriesVeto(maxtry()); eventRecord().clear(); eventRecord().prepare(sub, dynamic_ptr_cast(lastXCombPtr()), newStep(), pdfs(), ShowerHandler::currentHandler()->generator()->currentEvent()->incoming(), firstInteraction(), offShellPartons(), !doSubleadingNc); if ( doSubleadingNc ) { theSplittingReweight->updateCurrentHandler(); } continue; } catch (...) { throw; } } tPPair incoming=eventRecord().fillEventRecord(newStep(),firstInteraction(),didRealign); setDidRunCascade(true); return incoming; } // Reshuffle the outgoing partons from the hard process onto their constituent mass shells void DipoleShowerHandler::constituentReshuffle() { if ( constituentReshuffler && ShowerHandler::currentHandler()->retConstituentMasses() ) { if ( eventRecord().decays().empty() ) { constituentReshuffler->reshuffle(eventRecord().outgoing(), eventRecord().incoming(), eventRecord().intermediates()); return; } else { PList decaying; for(auto const & dec : eventRecord().decays()) decaying.push_back(dec.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 > decays; for(auto const & dec : eventRecord().decays() ) decays.push_back({dec.first,dec.second}); for(auto const & dec : decays) { PPtr unstable = dec.first; PList::iterator pos = find(eventRecord().intermediates().begin(), eventRecord().intermediates().end(), dec.first); // Update the PPtr in theDecays if(pos!=eventRecord().intermediates().end()) { unstable = *pos; while(!unstable->children().empty()) { unstable = unstable->children()[0]; } eventRecord().decays().erase(dec.first); eventRecord().decays()[unstable] = dec.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) { if ( Debug::level > 2 ){ // 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 ( auto const & testDecayOutItBefore : decayProc->outgoing() ) { testOutMomBefore += testDecayOutItBefore.first->momentum(); } testInvMassBefore = testOutMomBefore.m(); // decayReshuffle updates both the event record and the decay perturbative process if ( constituentReshuffler && ShowerHandler::currentHandler()->retConstituentMasses()) { constituentReshuffler->decayReshuffle(decayProc, eventRecord().outgoing(), eventRecord().hard(), eventRecord().intermediates()); } Lorentz5Momentum testOutMomAfter (ZERO,ZERO,ZERO,ZERO); Energy testInvMassAfter = ZERO; for ( auto const & testDecayOutItAfter : decayProc->outgoing() ) { 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); }else{ // decayReshuffle updates both the event record and the decay perturbative process if ( constituentReshuffler && ShowerHandler::currentHandler()->retConstituentMasses() ) { constituentReshuffler->decayReshuffle(decayProc, eventRecord().outgoing(), eventRecord().hard(), eventRecord().intermediates()); } return; } } // 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 ( auto const & p : eventRecord().outgoing() ) 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 ( auto const & p : eventRecord().hard()) phard += p->momentum(); Energy mhard = phard.m(); maxPt = mhard; } maxPt *= hardScaleFactor(); } else { maxPt = hardScaleFactor()*sqrt(muf); } muPt = maxPt; } else { muPt = hardScaleFactor()*sqrt(muf); } if ( doSubleadingNc ) { return; } for ( auto & ch : eventRecord().chains()) { // 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 ( auto & dip : ch.dipoles()) { // 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 ( auto const & k : kernels) { pair conf = {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 = {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({true,false},maxFirst); maxSecond = min(maxPt,maxSecond); dip.emitterScale({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 ( auto & dip : ch.dipoles() ) { dip.leftScale(min(dip.leftScale(),minVetoScale)); dip.rightScale(min(dip.rightScale(),minVetoScale)); } } } } void DipoleShowerHandler::hardScalesSubleading(list candidates, Energy hardPt) { maxPt = hardPt;//generator()->maximumCMEnergy(); // Note that minVetoScale is a value for each competing dipole (i.e. all dipoles // for the subleading shower. // It will contain the minimum veto scale from all of the dipoles Energy minVetoScale = -1.*GeV; for ( list::iterator cand = candidates.begin(); cand != candidates.end(); ++cand ) { // max scale Energy maxScale = ZERO; // Loop over kernels for ( vector::ptr>::iterator k = kernels.begin(); k != kernels.end(); ++k ) { if ( (**k).canHandle(cand->index()) ) { Energy scale = evolutionOrdering()->hardScale(cand->emitter(),cand->spectator(), cand->emitterX(),cand->spectatorX(), **k,cand->index()); maxScale = max(maxScale,scale); } } if ( cand->emitter()->vetoScale() >= ZERO ) { maxScale = min(maxScale,sqrt(cand->emitter()->vetoScale())); if ( minVetoScale >= ZERO ) minVetoScale = min(minVetoScale,sqrt(cand->emitter()->vetoScale())); else minVetoScale = sqrt(cand->emitter()->vetoScale()); } maxScale = min(maxPt,maxScale); cand->scale(maxScale); } if ( !evolutionOrdering()->independentDipoles() && chainOrderVetoScales && minVetoScale >= ZERO ) { for ( list::iterator cand = candidates.begin(); cand != candidates.end(); ++cand ) { cand->scale(min(cand->scale(),minVetoScale)); } } } void DipoleShowerHandler::addCandidates(PPair particles, list& clist) const { DipoleSplittingInfo candidate; Energy2 scale = ZERO; pair is(particles.first == eventRecord().incoming().first, particles.second == eventRecord().incoming().second); if ( (is.first && !is.second) || (!is.first && is.second) ) { scale = -(particles.first->momentum() - particles.second->momentum()).m2(); } else { scale = (particles.first->momentum() + particles.second->momentum()).m2(); } DipoleIndex index(particles.first->dataPtr(),particles.second->dataPtr(), is.first ? eventRecord().pdfs().first : PDF(), is.second ? eventRecord().pdfs().second : PDF()); candidate.scale(sqrt(scale)); candidate.index(index); candidate.configuration(make_pair(true,false)); candidate.emitter(particles.first); candidate.emitterX(is.first ? eventRecord().fractions().first : 1.0); candidate.spectator(particles.second); candidate.spectatorX(is.second ? eventRecord().fractions().second : 1.0); clist.push_back(candidate); index.swap(); candidate.index(index); candidate.configuration(make_pair(false,true)); candidate.emitter(particles.second); candidate.emitterX(is.second ? eventRecord().fractions().second : 1.0); candidate.spectator(particles.first); candidate.spectatorX(is.first ? eventRecord().fractions().first : 1.0); clist.push_back(candidate); } void DipoleShowerHandler::getCandidates(list& clist) const { clist.clear(); for ( PList::const_iterator i = eventRecord().outgoing().begin(); i != eventRecord().outgoing().end(); ++i ) { PList::const_iterator j = i; ++j; for ( ; j != eventRecord().outgoing().end(); ++j ) { addCandidates(make_pair(*i,*j),clist); } // Changed order of *i and inc().first if ( eventRecord().incoming().first->coloured() ) addCandidates(make_pair(eventRecord().incoming().first,*i),clist); if ( eventRecord().incoming().second->coloured() ) addCandidates(make_pair(*i,eventRecord().incoming().second),clist); } if ( eventRecord().incoming().first->coloured() && eventRecord().incoming().second->coloured() ) { addCandidates(eventRecord().incoming(),clist); } } void DipoleShowerHandler::performSplitting(DipoleSplittingInfo& split) const { Ptr::tptr kinematics = split.splittingKinematics(); kinematics->generateKinematics(split.emitter()->momentum(), split.spectator()->momentum(), split); split.splitEmitter(split.emitterData()->produceParticle(kinematics->lastEmitterMomentum())); split.splitSpectator(split.spectatorData()->produceParticle(kinematics->lastSpectatorMomentum())); split.emission(split.emissionData()->produceParticle(kinematics->lastEmissionMomentum())); // Setting resolution scales for the particles split.emission()->scale(sqr(split.lastPt())); split.splitEmitter()->scale(sqr(split.lastPt())); split.splitSpectator()->scale(split.spectator()->scale()); PVector neighbours; if ( DipolePartonSplitter::colourConnected(split.emitter(), eventRecord().incoming().first) && split.emitter() != eventRecord().incoming().first ) neighbours.push_back(eventRecord().incoming().first); if ( DipolePartonSplitter::colourConnected(split.emitter(), eventRecord().incoming().second) && split.emitter() != eventRecord().incoming().second ) neighbours.push_back(eventRecord().incoming().second); for ( PList::const_iterator p = eventRecord().outgoing().begin(); p != eventRecord().outgoing().end(); ++p ) { if ( *p == split.emitter() ) continue; if ( DipolePartonSplitter::colourConnected(split.emitter(),*p) ) neighbours.push_back(*p); } assert(neighbours.size() == 1 || neighbours.size() == 2 ); if ( neighbours.size() == 2 ) { if ( UseRandom::rnd() < 0.5 ) swap(neighbours[0],neighbours[1]); } DipolePartonSplitter::split(split.emitter(),split.splitEmitter(),split.emission(), neighbours.front(),split.index().initialStateEmitter(),false); DipolePartonSplitter::change(split.spectator(),split.splitSpectator(), split.index().initialStateSpectator(),false); } Energy DipoleShowerHandler::nextSubleadingSplitting(Energy hardPt, Energy optHardPt, Energy optCutoff, const bool decay) { list candidates; getCandidates(candidates); hardScalesSubleading(candidates,hardPt); for ( list::iterator cand = candidates.begin(); cand != candidates.end(); cand++ ) { cand->scale(hardPt); } list::iterator split = candidates.end(); // Winner of all dipoles DipoleSplittingInfo winner; // Winner for the current iteration of the for loop DipoleSplittingInfo candWinner; Energy winnerScale = 0.0*GeV; Energy nextScale = 0.0*GeV; for ( list::iterator cand = candidates.begin(); cand != candidates.end(); cand++ ) { nextScale = getWinner(candWinner, cand->index(), cand->emitterX(),cand->spectatorX(), make_pair(true,false), cand->emitter(),cand->spectator(), hardPt, optHardPt, optCutoff); if ( nextScale > winnerScale ) { winnerScale = nextScale; winner = candWinner; split = cand; winnerIndex = winningKernelIndex;//check } } if ( split == candidates.end() ) return ZERO; if ( decay ) winner.isDecayProc( true ); split->fill(winner); performSplitting(*split); eventRecord().update(*split); for ( list::iterator dip = candidates.begin(); dip != candidates.end(); ++dip ) { if ( dip == split ) continue; dip->emission(split->emission()); if ( dip->emitter() == split->emitter() ) { dip->splitEmitter(split->splitEmitter()); } else { dip->splitEmitter(dip->emitter()); } if ( dip->spectator() == split->spectator() ) { dip->splitSpectator(split->splitSpectator()); } else { dip->splitSpectator(dip->spectator()); } } // Update the ShowerHandler of the splitting reweight. if ( doSubleadingNc ) { theSplittingReweight->updateCurrentHandler(); } return split->lastPt(); } Energy DipoleShowerHandler::getWinner(DipoleSplittingInfo& winner, const Dipole& dip, pair 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 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; } if ( index.incomingDecaySpectator() && !doFSR() ) { 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); candidate.emitter(emitter); candidate.spectator(spectator); 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 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 ) { if ( doPartialUnweighting ) gen->second->doPartialUnweighting(referenceWeight); // (*) 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().incomingDecaySpectator() || candidate.index().incomingDecayEmitter() ) { Energy recoilMass = gen->second->splittingKinematics()->recoilMassKin(emitter->momentum(), spectator->momentum()); candidate.recoilMass(recoilMass); } // Store emitter and spectator masses, needed in kinematics if ( candidate.index().emitterData()->mass() != ZERO ) { if ( !candidate.index().offShellEmitter() ) candidate.emitterMass( emitter->nominalMass() ); else candidate.emitterMass( emitter->mass() ); } if ( candidate.index().spectatorData()->mass() != ZERO ) { if ( !candidate.index().offShellSpectator() ) candidate.spectatorMass( spectator->nominalMass() ); else candidate.spectatorMass( spectator->mass() ); } 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; if ( continueSubleadingNc() ) winningKernelIndex = kernelIndex+1;//check } if ( continueSubleadingNc() ) { kernelIndex++;//check scales.push_back(nextScale);//check theWeightsVector.push_back(gen->second->splittingWeightVector()); } 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; if ( doSubleadingNc ) { unsigned int subEmDone = 0; // Set the starting scale Energy hardPt = muPt; double wref = referenceWeight; while ( subEmDone < subleadingNcEmissionsLimit && hardPt != ZERO && continueSubleadingNc() ) { // Clear out the weights from the earlier step theWeightsVector.clear(); kernelIndex = 0;//check scales.clear();//check hardPt = nextSubleadingSplitting( hardPt, optHardPt, optCutoff, decay ); // Partial unweighting if ( doPartialUnweightingAtEmission ) { const double w = reweight(); if ( abs(w) < wref ) { if ( abs(w)/wref < UseRandom::rnd() ) { // Set weight to zero and end this event reweight(0.0); return; } else reweight( wref*w/abs(w) ); } // Update the reference weight after emission wref *= referenceWeight; } // When the winning scale is larger than the cutoff // remove the added weights that are under the winning scale if ( hardPt != ZERO ) { Energy maxq = 0.0*GeV; size_t iwinner = theWeightsVector.size();//check for ( size_t i = 0; i < theWeightsVector.size(); i++ ) { if ( theWeightsVector[i].size() > 0 ) { // get<2> is true for an accept step. if ( std::get<2>(theWeightsVector[i].back()) && std::get<0>(theWeightsVector[i].back()) > maxq) { maxq = std::get<0>(theWeightsVector[i].back()); iwinner = i;//check } } } assert(winnerIndex-1 == iwinner);//check double correctionWeight = 1.0; for ( size_t i = 0; i < theWeightsVector.size(); i++ ) { for ( size_t j = 0; j < theWeightsVector[i].size(); j++ ) { if ( std::get<0>(theWeightsVector[i][j]) < maxq ) correctionWeight *= std::get<1>(theWeightsVector[i][j]); } } reweight(reweight()/correctionWeight); } // Increment the number of subleading Nc emissions done subEmDone++; // Stop if the limit of emissions is reached if ( nEmissions ) if ( ++emDone == nEmissions ) return; } // Subleading shower done, prepare chains for the standard // dipole shower eventRecord().prepareChainsSubleading( decay ); // Set scales for ( list::iterator ch = eventRecord().chains().begin(); ch != eventRecord().chains().end(); ch++ ) { for ( list::iterator dp = ch->dipoles().begin(); dp != ch->dipoles().end(); dp++ ) { dp->emitterScale(make_pair(true,false),hardPt); dp->emitterScale(make_pair(false,true),hardPt); } } } DipoleSplittingInfo winner; DipoleSplittingInfo dipoleWinner; while ( eventRecord().haveChain() ) { // allow the dipole chain to be rearranged according to arXiv:1801.06113 if( _rearrange && ( _rearrangeNEmissions < 0 || _rearrangeNEmissions >= int(emDone) ) ){ eventRecord().currentChain().rearrange(_dipmax,_diplong); } if ( verbosity > 2 ) { generator()->log() << "DipoleShowerHandler selecting splittings for the chain:\n" << eventRecord().currentChain() << flush; } list::iterator winnerDip = eventRecord().currentChain().dipoles().end(); Energy winnerScale = 0.0*GeV; Energy nextLeftScale = 0.0*GeV; Energy nextRightScale = 0.0*GeV; for ( list::iterator dip = eventRecord().currentChain().dipoles().begin(); dip != eventRecord().currentChain().dipoles().end(); ++dip ) { nextLeftScale = getWinner(dipoleWinner,*dip,{true,false},optHardPt,optCutoff); if ( nextLeftScale > winnerScale ) { winnerScale = nextLeftScale; winner = dipoleWinner; winnerDip = dip; } nextRightScale = getWinner(dipoleWinner,*dip,{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 - // but first see if the emission would produce a configuration in the ME region. - if ( theMergingHelper - && eventHandler()->currentCollision() - && !decay - && firstInteraction() ) { - if (theMergingHelper->maxLegs()>eventRecord().outgoing().size()+ - eventRecord().hard().size() - +2){//incoming - - if (theMergingHelper->mergingScale()emissionProbability() < UseRandom::rnd()) { - - theMergingHelper->setEmissionProbability(0.); - const bool transparent=true; - if (transparent) { - pair::iterator,list::iterator> tmpchildren; - DipoleSplittingInfo tmpwinner=winner; - DipoleChain* tmpfirstChain = nullptr; - DipoleChain* tmpsecondChain = nullptr; - - auto New=eventRecord().tmpsplit(winnerDip,tmpwinner, - tmpchildren,tmpfirstChain, - tmpsecondChain); - - - if (theMergingHelper->matrixElementRegion(New.first, - New.second, - winnerScale, - theMergingHelper->mergingScale())) { - optHardPt=winnerScale; - continue; - } - }else{ - optHardPt=winnerScale; - continue; - - } - } + // Check if the produced splitting would be part + // of a matrix element region in a multi jet merging. + // If so, the optHardPt is set to the current winner scale and the slitting + // is vetoed. Note: It is possible that the current scale is larger than the + // merging scale but another shower configuration shows a scale that is + // below the required scale. + if(firstInteraction() && !decay){ + if ( isMERegion( winnerScale , winner , winnerDip) ){ + optHardPt=winnerScale; // Vetoed shower. + continue; + }else if( theMergingHelper ){ + optHardPt=ZERO; } } - if(theMergingHelper&&firstInteraction()) - optHardPt=ZERO; - - + // otherwise perform the splitting didRadiate = true; eventRecord().isMCatNLOSEvent(false); eventRecord().isMCatNLOHEvent(false); pair::iterator,list::iterator> children; DipoleChain* firstChain = nullptr; DipoleChain* secondChain = nullptr; // Generate the azimuthal angle if ( spinCorrelations() ) vertexRecord().generatePhi(winner,*winnerDip); if ( decay ) winner.isDecayProc( true ); // 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); // Update the vertex record following the splitting if ( spinCorrelations() ) vertexRecord().update(winner); 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::isMERegion(const Energy winnerScale, + const DipoleSplittingInfo & winner, + list::iterator winnerDip ){ + // First check if we have a merging setup adn additional conditions. + if(!theMergingHelper)return false; + if(!eventHandler()->currentCollision())return false; + if(theMergingHelper->maxLegs()<=eventRecord().outgoing().size()+ + eventRecord().hard().size() + +2)//incoming + return false; + + if( theMergingHelper->mergingScale() > winnerScale ) return false; + + // Now produce a temporary splitting to check if the final configuration is + // part of the matrix element region + pair::iterator,list::iterator> tmpchildren; + DipoleSplittingInfo tmpwinner=winner; + DipoleChain* tmpfirstChain = nullptr; + DipoleChain* tmpsecondChain = nullptr; + + auto New=eventRecord().tmpsplit(winnerDip,tmpwinner, + tmpchildren,tmpfirstChain, + tmpsecondChain); + // This is the same function that is used in the ME calculations of the + // multi jet merging. + if (theMergingHelper->matrixElementRegion(New.first, + New.second, + winnerScale, + theMergingHelper->mergingScale())) { + return true; + }else return false; +} + + + + + + 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 inMomenta (eventRecord().incoming().first->momentum(), eventRecord().incoming().second->momentum()); LorentzRotation 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 masses(eventRecord().incoming().first->mass(), eventRecord().incoming().second->mass()); pair 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::tptr as) { for ( auto & k : kernels) { 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::tptr rw) { for ( auto & g : generators() ) g.second->splittingReweight(rw); } void DipoleShowerHandler::getGenerators(const DipoleIndex& ind, Ptr::tptr rw) { bool gotone = false; for ( auto & k : kernels ) { 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::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({ind,nGenerator}); } } if ( !gotone ) { throw 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::runerror; } } // 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); // copy off-shell particle ids before showering from input vector to the // set used in the simulation if ( theColouredOffShellInShower.empty() ) { for(unsigned int ix=0;ixsplittingKinematics()->openZBoundaries() == 0 ) zChoice0 = true; else if ( k->splittingKinematics()->openZBoundaries() == 1 ) zChoice1 = true; else zChoiceOther = true; // either inconsistent or other option which cannot be handled by the matching if ( zChoice0 && zChoice1 ) { zChoiceOther = true; break; } } if ( zChoiceOther ) theZBoundaries = 2; else if ( zChoice1 ) theZBoundaries = 1; else if ( zChoice0 ) theZBoundaries = 0; } 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 //<< theAnalyseSpinCorrelations << doSubleadingNc << subleadingNcEmissionsLimit << densityOperatorEvolution << ounit(densityOperatorCutoff,GeV2) << doPartialUnweightingAtEmission << doPartialUnweighting << referenceWeight << cmecReweightFactor << negCMECScaling << realignmentScheme << verbosity << printEvent << ounit(theRenormalizationScaleFreeze,GeV) << ounit(theFactorizationScaleFreeze,GeV) << theShowerApproximation << theDoCompensate << theFreezeGrid << theDetuning << theEventReweight << theSplittingReweight << ounit(maxPt,GeV) << ounit(muPt,GeV)<< theMergingHelper << theColouredOffShellInShower << theInputColouredOffShellInShower << _rearrange << _dipmax << _diplong << _rearrangeNEmissions << theZBoundaries; } void DipoleShowerHandler::persistentInput(PersistentIStream & is, int) { is >> kernels >> theEvolutionOrdering >> constituentReshuffler >> intrinsicPtGenerator >> theGlobalAlphaS >> chainOrderVetoScales >> nEmissions >> discardNoEmissions >> firstMCatNLOEmission >> thePowhegDecayEmission //>> theAnalyseSpinCorrelations >> doSubleadingNc >> subleadingNcEmissionsLimit >> densityOperatorEvolution >> iunit(densityOperatorCutoff,GeV2) >> doPartialUnweightingAtEmission >> doPartialUnweighting >> referenceWeight >> cmecReweightFactor >> negCMECScaling >> realignmentScheme >> verbosity >> printEvent >> iunit(theRenormalizationScaleFreeze,GeV) >> iunit(theFactorizationScaleFreeze,GeV) >> theShowerApproximation >> theDoCompensate >> theFreezeGrid >> theDetuning >> theEventReweight >> theSplittingReweight >> iunit(maxPt,GeV) >> iunit(muPt,GeV)>>theMergingHelper >> theColouredOffShellInShower >> theInputColouredOffShellInShower >> _rearrange >> _dipmax >> _diplong >> _rearrangeNEmissions >> theZBoundaries; } ClassDescription DipoleShowerHandler::initDipoleShowerHandler; // Definition of the static class description member. void DipoleShowerHandler::Init() { static ClassDocumentation 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 interfaceKernels ("Kernels", "Set the splitting kernels to be used by the dipole shower.", &DipoleShowerHandler::kernels, -1, false, false, true, false, false); static Reference interfaceEvolutionOrdering ("EvolutionOrdering", "Set the evolution ordering to be used.", &DipoleShowerHandler::theEvolutionOrdering, false, false, true, false, false); static Reference interfaceConstituentReshuffler ("ConstituentReshuffler", "The object to be used to reshuffle partons to their constitutent mass shells.", &DipoleShowerHandler::constituentReshuffler, false, false, true, true, false); static Reference interfaceIntrinsicPtGenerator ("IntrinsicPtGenerator", "Set the object in charge to generate intrinsic pt for incoming partons.", &DipoleShowerHandler::intrinsicPtGenerator, false, false, true, true, false); static Reference interfaceGlobalAlphaS ("GlobalAlphaS", "Set a global strong coupling for all splitting kernels.", &DipoleShowerHandler::theGlobalAlphaS, false, false, true, true, false); // Start: Trying to add interface for the subleading Nc static Switch interfaceDoSubleadingNc ("DoSubleadingNc", "Switch on or off subleading Nc corrections.", &DipoleShowerHandler::doSubleadingNc, true, false, false); static SwitchOption interfaceDoSubleadingNcOn (interfaceDoSubleadingNc, "On", "Switch on subleading Nc corrections.", true); static SwitchOption interfaceDoSubleadingNcOff (interfaceDoSubleadingNc, "Off", "Switch off subleading Nc corrections.", false); // Limit for how many subleading Nc emissions should be calculated static Parameter interfaceSubleadingNcEmissionsLimit ("SubleadingNcEmissionsLimit", "Number of emissions to calculate subleading Nc corrections for.", &DipoleShowerHandler::subleadingNcEmissionsLimit,0,0,0, false, false, Interface::lowerlim); static Parameter interfaceDensityOperatorEvolution ("DensityOperatorEvolution", "Scheme for evolving the density operator.", &DipoleShowerHandler::densityOperatorEvolution,0,0,0, false, false, Interface::lowerlim); static Parameter interfaceDensityOperatorCutoff ("DensityOperatorCutoff", "Cutoff for momentum invariants for the density operator evolution.", &DipoleShowerHandler::densityOperatorCutoff,GeV2,1.0*GeV2,0.0*GeV2,0*GeV2, false, false, Interface::lowerlim); static Switch interfaceDoPartialUnweightingAtEmission ("DoPartialUnweightingAtEmission", "Switch on or off partial unweighting at the emission level.", &DipoleShowerHandler::doPartialUnweightingAtEmission,true,false,false); static SwitchOption interfaceDoPartialUnweightingAtEmissionOn (interfaceDoPartialUnweightingAtEmission, "On", "Switch on partial unweighting.", true); static SwitchOption interfaceDoPartialUnweightingAtEmissionOff (interfaceDoPartialUnweightingAtEmission, "Off", "Switch off partial unweighting.", false); static Switch interfaceDoPartialUnweighting ("DoPartialUnweighting", "Switch on or off partial unweighting at the dipole splitting level.", &DipoleShowerHandler::doPartialUnweighting,true,false,false); static SwitchOption interfaceDoPartialUnweightingOn (interfaceDoPartialUnweighting, "On", "Switch on partial unweighting.", true); static SwitchOption interfaceDoPartialUnweightingOff (interfaceDoPartialUnweighting, "Off", "Switch off partial unweighting.", false); static Parameter interfaceReferenceWeight ("ReferenceWeight", "Reference weight for the partial unweighting.", &DipoleShowerHandler::referenceWeight,0.1,0.0,0, false, false, Interface::lowerlim); static Parameter interfaceCMECReweightFactor ("CMECReweightFactor", "Factor used in the reweighting algorithm.", &DipoleShowerHandler::cmecReweightFactor,1.0,0.0,0, false, false, Interface::lowerlim); static Parameter interfaceNegCMECScaling ("NegCMECScaling", "Scaling factor for the negative colour matrix element corrections (CMECs).", &DipoleShowerHandler::negCMECScaling,0.0,0.0,0, false, false, Interface::lowerlim); static Switch 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 interfaceChainOrderVetoScales ("ChainOrderVetoScales", "[experimental] Switch the chain ordering for veto scales on or off.", &DipoleShowerHandler::chainOrderVetoScales, true, false, false); static SwitchOption interfaceChainOrderVetoScalesYes (interfaceChainOrderVetoScales, "Yes", "Switch on chain ordering for veto scales.", true); static SwitchOption interfaceChainOrderVetoScalesNo (interfaceChainOrderVetoScales, "No", "Switch off chain ordering for veto scales.", false); interfaceChainOrderVetoScales.rank(-1); static Parameter 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 interfaceDiscardNoEmissions ("DiscardNoEmissions", "[debug option] Discard events without radiation.", &DipoleShowerHandler::discardNoEmissions, false, false, false); static SwitchOption interfaceDiscardNoEmissionsYes (interfaceDiscardNoEmissions, "Yes", "Discard events without radiation.", true); static SwitchOption interfaceDiscardNoEmissionsNo (interfaceDiscardNoEmissions, "No", "Do not discard events without radiation.", false); interfaceDiscardNoEmissions.rank(-1); static Switch interfaceFirstMCatNLOEmission ("FirstMCatNLOEmission", "[debug option] Only perform the first MC@NLO emission.", &DipoleShowerHandler::firstMCatNLOEmission, false, false, false); static SwitchOption interfaceFirstMCatNLOEmissionYes (interfaceFirstMCatNLOEmission, "Yes", "Perform only the first MC@NLO emission.", true); static SwitchOption interfaceFirstMCatNLOEmissionNo (interfaceFirstMCatNLOEmission, "No", "Produce all emissions.", false); interfaceFirstMCatNLOEmission.rank(-1); static Parameter 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 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 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 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 interfaceDoCompensate ("DoCompensate", "", &DipoleShowerHandler::theDoCompensate, false, false, false); static SwitchOption interfaceDoCompensateYes (interfaceDoCompensate, "Yes", "", true); static SwitchOption interfaceDoCompensateNo (interfaceDoCompensate, "No", "", false); static Parameter interfaceFreezeGrid ("FreezeGrid", "", &DipoleShowerHandler::theFreezeGrid, 500000, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceDetuning ("Detuning", "A value to detune the overestimate kernel.", &DipoleShowerHandler::theDetuning, 1.0, 1.0, 0, false, false, Interface::lowerlim); static Reference interfaceEventReweight ("EventReweight", "", &DipoleShowerHandler::theEventReweight, false, false, true, true, false); static Reference interfaceSplittingReweight ("SplittingReweight", "Set the splitting reweight.", &DipoleShowerHandler::theSplittingReweight, false, false, true, true, false); static Switch 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 ParVector interfaceOffShellInShower ("OffShellInShower", "PDG codes of the coloured particles that can be off-shell in the process.", &DipoleShowerHandler::theInputColouredOffShellInShower, -1, 0l, -10000000l, 10000000l, false, false, Interface::limited); /* static Switch interfaceAnalyseSpinCorrelations ("AnalyseSpinCorrelations", "Record the information required for the spin correlation analyis.", &DipoleShowerHandler::theAnalyseSpinCorrelations, false, false, false); static SwitchOption interfaceAnalyseSpinCorrelationsYes (interfaceAnalyseSpinCorrelations,"Yes","Record the information for analysing the spin correlations.", true); static SwitchOption interfaceAnalyseSpinCorrelationsNo (interfaceAnalyseSpinCorrelations,"No","Do not record extra information.", false); */ static Switch interfacerearrange ("Rearrange", "Allow rearranging of dipole chains according to arXiv:1801.06113", &DipoleShowerHandler::_rearrange, false, false, false); static SwitchOption interfacerearrangeYes (interfacerearrange,"Yes","_rearrange on", true); static SwitchOption interfacerearrangeNo (interfacerearrange,"No","_rearrange off", false); static Parameter interfacedipmax ("DipMax", "Allow rearrangment of color chains with ME including dipmax dipoles.", &DipoleShowerHandler::_dipmax, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfacediplong ("DipLong", "Dipole chains with more than dipmax dipoles are treated as long. \ diplong=3 rearranges these chains with eeuugg MEs, \ diplong=4 rearranges these chains with eeuuggg MEs (slower), \ diplong=5 rearranges these chains with eeuugggg MEs (slow).\ Note: Numerically there is no difference between the options. ", &DipoleShowerHandler::_diplong, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfacedcorrectNemissions ("RearrangeNEmissions", "Allow rearrangment of color chains up to the nth emission.", &DipoleShowerHandler::_rearrangeNEmissions, 0, 0, 0, false, false, Interface::lowerlim); } diff --git a/Shower/Dipole/DipoleShowerHandler.h b/Shower/Dipole/DipoleShowerHandler.h --- a/Shower/Dipole/DipoleShowerHandler.h +++ b/Shower/Dipole/DipoleShowerHandler.h @@ -1,777 +1,785 @@ // -*- C++ -*- // // DipoleShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_DipoleShowerHandler_H #define HERWIG_DipoleShowerHandler_H // // This is the declaration of the DipoleShowerHandler class. // #include "Herwig/Shower/ShowerHandler.h" #include "Herwig/Shower/Dipole/DipoleShowerHandler.fh" #include "Herwig/Shower/Dipole/Base/DipoleSplittingInfo.h" #include "Herwig/Shower/Dipole/Base/DipoleSplittingReweight.h" #include "Herwig/Shower/Dipole/Kernels/DipoleSplittingKernel.h" #include "Herwig/Shower/Dipole/Base/DipoleSplittingGenerator.h" #include "Herwig/Shower/Dipole/Base/DipoleEventRecord.h" #include "Herwig/Shower/Dipole/Base/DipoleEvolutionOrdering.h" #include "Herwig/Shower/Dipole/Base/DipoleEventReweight.h" #include "Herwig/Shower/Dipole/Utility/ConstituentReshuffler.h" #include "Herwig/Shower/Dipole/Utility/IntrinsicPtGenerator.h" #include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h" #include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h" #include "Herwig/Shower/Dipole/SpinCorrelations/DipoleVertexRecord.h" #include "Herwig/MatrixElement/Matchbox/Utility/DensityOperator.h" #include namespace Herwig { using namespace ThePEG; /** * \ingroup DipoleShower * \author Simon Platzer, Stephen Webster * * \brief The DipoleShowerHandler class manages the showering using * the dipole shower algorithm. * * @see \ref DipoleShowerHandlerInterfaces "The interfaces" * defined for DipoleShowerHandler. */ class DipoleShowerHandler: public ShowerHandler { friend class Merger; public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ DipoleShowerHandler(); /** * The destructor. */ virtual ~DipoleShowerHandler(); //@} public: inline void colourPrint(); /** * Indicate a problem in the shower. */ struct RedoShower {}; /** * Insert an additional splitting kernel. */ void addSplitting(Ptr::ptr sp) { kernels.push_back(sp); } /** * Reset the alpha_s for all splitting kernels. */ void resetAlphaS(Ptr::tptr); virtual void cascade(tPVector); /** * Reset the splitting reweight for all splitting kernels. */ void resetReweight(Ptr::tptr); /** * Return true, if the shower handler can generate a truncated * shower for POWHEG style events generated using Matchbox */ virtual bool canHandleMatchboxTrunc() const { return false; } /** * Return true, if this cascade handler will perform reshuffling from hard * process masses. */ virtual bool isReshuffling() const { return false; } /** * Return the relevant hard scale to be used in the profile scales */ virtual Energy hardScale() const { return muPt; } /** * Calculate the alpha_s value the shower uses for Q. */ double as(Energy Q)const{return theGlobalAlphaS->value(sqr(Q));} /** * Return the number of scale dependent active flavours from * the alpha_s object. */ double Nf(Energy Q)const{return theGlobalAlphaS->Nf(sqr(Q));} /* * Access the vertex record */ DipoleVertexRecord& vertexRecord() { return theVertexRecord; } /** * Return the event record */ const DipoleVertexRecord& vertexRecord() const { return theVertexRecord; } /** * Set the pointer to the Merging Helper. * Used by the merging factory. */ void setMerger(Ptr::ptr mh){theMergingHelper=mh;} public: /** * Return the dictionary for the incoming particles and outgoing partons * in the event, will be empty if subleading Nc is turned off. */ const map& particleIndices() const { return eventRecord().particleIndices(); } /** * Return the particle data vector of the particles in the event, will * be empty if subleading Nc is turned off. */ const cPDVector& particlesAfter() const { return eventRecord().particlesAfter(); } /** * Return the density operator from the event record. */ DensityOperator& densityOperator() { return eventRecord().densityOperator(); } /** * Return the colour matrix element correction map, will be empty if * subleading Nc is turned off. NOT USED, REMOVE */ const map,pair >,double>& correlatorMap() const { return eventRecord().densityOperator().correlatorMap(); } /** * Return the colour basis used in the density operator. NOT USED, REMOVE */ Ptr::tptr colourBasis() { return eventRecord().densityOperator().colourBasis(); } /** * Return the continue subleading Nc flag from the event record. */ bool continueSubleadingNc() const { return eventRecord().getContinueSubleadingNc(); } protected: /** * Add the possible splitting candidates given a pair of emitting particles */ void addCandidates(PPair particles, list& clist) const; /** * Get all possible radiating dipoles */ void getCandidates(list& clist) const; /** * Perform a splitting independent of any chains */ void performSplitting(DipoleSplittingInfo&) const; /** * Generate the next subleading Nc improved splitting and return the scale */ Energy nextSubleadingSplitting(Energy hardPt, Energy optHardPt, Energy optCutoff, const bool decay); protected: typedef multimap::ptr> GeneratorMap; /** * The main method which manages the showering of a subprocess. */ virtual tPPair cascade(tSubProPtr sub, XCombPtr xcomb) { return cascade(sub,xcomb,ZERO,ZERO); } /** * The main method which manages the showering of a subprocess. */ tPPair cascade(tSubProPtr sub, XCombPtr xcomb, Energy optHardPt, Energy optCutoff); /** * Build splitting generators for the given * dipole index. */ void getGenerators(const DipoleIndex&, Ptr::tptr rw = Ptr::tptr()); /** * Setup the hard scales. */ void hardScales(Energy2 scale); /** * Setup the hard scales of the dipoles after subleading emissions.//debug */ void hardScalesSubleading(list candidates,Energy hardPt); /** * Return the evolution ordering */ Ptr::tptr evolutionOrdering() const { return theEvolutionOrdering; } /** * Reshuffle to constituent mass shells */ void constituentReshuffle(); /** * Reshuffle to constituent mass shells */ void decayConstituentReshuffle( PerturbativeProcessPtr decayProc); /** * Access the generator map */ GeneratorMap& generators() { return theGenerators; } /** * Access the event record */ DipoleEventRecord& eventRecord() { return theEventRecord; } /** * Return the event record */ const DipoleEventRecord& eventRecord() const { return theEventRecord; } /** * Return the splitting kernels. */ const vector::ptr>& splittingKernels() const { return kernels; } /** * Return the set of offshell parton ids. **/ const set& offShellPartons() { return theColouredOffShellInShower; } /** + * In a merging setup this function checks if the next shower + * configuration is part of the matrix element region. + */ + bool isMERegion(const Energy winnerScale, + const DipoleSplittingInfo & winner, + const list::iterator winnerDip); + + /** * Realign the event such as to have the incoming partons along thre * beam axes. */ bool realign(); /** * The choice of z boundaries; 0 = restricted, 1 = open, 2 = mixed/other */ virtual int showerPhaseSpaceOption() const { return theZBoundaries; } protected: /** * Perform the cascade. */ void doCascade(unsigned int& emDone, Energy optHardPt = ZERO, Energy optCutoff = ZERO, const bool decay = false); /** * Set the number of emissions **/ void setNEmissions(unsigned int n){nEmissions=n;} /** * Get the winning splitting for the * given dipole and configuration. */ Energy getWinner(DipoleSplittingInfo& winner, const Dipole& dip, pair conf, Energy optHardPt = ZERO, Energy optCutoff = ZERO); /** * Get the winning splitting for the * given dipole and configuration. */ Energy getWinner(DipoleSplittingInfo& winner, Energy optHardPt = ZERO, Energy optCutoff = ZERO); /** * Get the winning splitting for the * given dipole and configuration. */ Energy getWinner(SubleadingSplittingInfo& winner, Energy optHardPt = ZERO, Energy optCutoff = ZERO); /** * Get the winning splitting for the * given dipole and configuration. */ Energy getWinner(DipoleSplittingInfo& winner, const DipoleIndex& index, double emitterX, double spectatorX, pair conf, tPPtr emitter, tPPtr spectator, Energy startScale, Energy optHardPt = ZERO, Energy optCutoff = ZERO); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); //@} private: /** * The splitting kernels to be used. */ vector::ptr> kernels; /** * The evolution ordering considered */ Ptr::ptr theEvolutionOrdering; /** * The ConstituentReshuffler to be used */ Ptr::ptr constituentReshuffler; /** * The intrinsic pt generator to be used. */ Ptr::ptr intrinsicPtGenerator; /** * A global alpha_s to be used for all splitting kernels. */ Ptr::ptr theGlobalAlphaS; /** * Apply chain ordering to events from matrix * element corrections. */ bool chainOrderVetoScales; /** * Limit the number of emissions. * Limit applied if > 0. */ unsigned int nEmissions; /** * Discard events which did not radiate. */ bool discardNoEmissions; /** * Perform the first MC@NLO emission only. */ bool firstMCatNLOEmission; /** * True if powheg style emissions are to be used in the decays */ bool thePowhegDecayEmission; /** * Switch to record information required for the * nearest neighbour analysis. */ //bool theAnalyseSpinCorrelations; /** * The realignment scheme */ int realignmentScheme; /** * Switch on or off subleading Nc corrections */ bool doSubleadingNc; /** * Number of emissions to do subleading Nc corrections to. */ size_t subleadingNcEmissionsLimit; /** * Current reference weight used for partial unweighting of the * subleading colour shower (updated each emission). */ int currentReferenceWeight; /** * Integer used to set which method of evolving the density operator * to use: * 0 - Vijk is Eikonal but there is a cutoff. * 1 - Vijk is Eikonal. * 2 - Vijk=1 for all i,j,k. * 3 - Semi-leading Nc, Vijk=0 for all */ int densityOperatorEvolution; /** * Cutoff scale for the invariants (e.g. pEmitter*pEmission) in the * Eikonal dipole kernel, Vijk. */ Energy2 densityOperatorCutoff; /** * Switch on or off partial unweighting in after each subleading emission. */ bool doPartialUnweightingAtEmission; /** * Switch on or off partial unweighting in the splitting generator. */ bool doPartialUnweighting; /** * Reference weight for the partial unweighting. */ double referenceWeight; /** * Factor changing the acceptance probability for the veto algorithm. */ double cmecReweightFactor; /** * Scaling factor for the negative colour matrix element corrections. */ double negCMECScaling; private: /** * The verbosity level. * 0 - print no info * 1 - print diagnostic information on setting up * splitting generators etc. * 2 - print detailed event information for up to * printEvent events. * 3 - print dipole chains after each splitting. */ int verbosity; /** * See verbosity. */ int printEvent; private: /** * The splitting generators indexed by the dipole * indices they can work on. */ GeneratorMap theGenerators; /** * The evnt record used. */ DipoleEventRecord theEventRecord; /** * The vertex record. **/ DipoleVertexRecord theVertexRecord; /** * The number of shoer tries so far. */ unsigned int nTries; /** * Whether or not we did radiate anything */ bool didRadiate; /** * Whether or not we did realign the event */ bool didRealign; /** * Vector of candidate splittings containing a vector of the * weights, scale and a bool for every step in the reweighted * veto algorithm. The bool is true for an accept step. */ vector > > theWeightsVector; /** * Winning candidate index. */ size_t winnerIndex;//debug size_t kernelIndex;//debug size_t winningKernelIndex;//debug vector scales;//debug private: /** * A freezing value for the renormalization scale */ Energy theRenormalizationScaleFreeze; /** * A freezing value for the factorization scale */ Energy theFactorizationScaleFreeze; /** * The matching subtraction, if appropriate */ Ptr::tptr theShowerApproximation; /** * True, if sampler should apply compensation */ bool theDoCompensate; /** * Return the number of accepted points after which the grid should * be frozen */ unsigned long theFreezeGrid; /** * The detuning factor applied to the sampling overestimate kernel */ double theDetuning; /** * A pointer to the dipole event reweight object */ Ptr::ptr theEventReweight; /** * A pointer to a global dipole splitting reweight */ Ptr::ptr theSplittingReweight; /** * True if no warnings have been issued yet */ static bool firstWarn; /** * The shower starting scale for the last event encountered */ Energy maxPt; /** * The shower hard scale for the last event encountered */ Energy muPt; /** * The merging helper takes care of merging multiple LO and NLO * cross sections. Here we need to check if an emission would * radiate in the matrix element region of an other multipicity. * If so, the emission is vetoed. */ Ptr::ptr theMergingHelper; /** * PDG codes of the partons which can have an off-shell mass, * this is fast storage for use during running */ set theColouredOffShellInShower; /** * PDG codes of the partons which can have an off-shell mass, * this is a vector that is interfaced so they can be changed */ vector theInputColouredOffShellInShower; /** * Allow the dipole chains to be rearranged */ bool _rearrange=false; /** * number of maximal ME dipoles in the rearrangement. */ unsigned int _dipmax=3; /** * If a chain is considered long (more than dipmax dipoles) * ME with diplong dipoles are used to test for rearrangement. */ unsigned int _diplong=3; /** * Number of emissions to be rearranged. */ int _rearrangeNEmissions=-1; /** * The choice of z boundaries; 0 = restricted, 1 = open, 2 = mixed/other */ int theZBoundaries; private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initDipoleShowerHandler; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ DipoleShowerHandler & operator=(const DipoleShowerHandler &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of DipoleShowerHandler. */ template <> struct BaseClassTrait { /** Typedef of the first base class of DipoleShowerHandler. */ typedef Herwig::ShowerHandler NthBase; }; /** This template specialization informs ThePEG about the name of * the DipoleShowerHandler class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::DipoleShowerHandler"; } /** * The name of a file containing the dynamic library where the class * DipoleShowerHandler is implemented. It may also include several, space-separated, * libraries if the class DipoleShowerHandler depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_DipoleShowerHandler_H */ diff --git a/Shower/Dipole/Kernels/DipoleSplittingKernel.h b/Shower/Dipole/Kernels/DipoleSplittingKernel.h --- a/Shower/Dipole/Kernels/DipoleSplittingKernel.h +++ b/Shower/Dipole/Kernels/DipoleSplittingKernel.h @@ -1,535 +1,540 @@ // -*- C++ -*- // // DipoleSplittingKernel.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_DipoleSplittingKernel_H #define HERWIG_DipoleSplittingKernel_H // // This is the declaration of the DipoleSplittingKernel class. // #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/StandardModel/AlphaSBase.h" #include "ThePEG/PDF/PDF.h" #include "Herwig/Shower/Dipole/Utility/PDFRatio.h" #include "Herwig/Shower/Dipole/Base/DipoleSplittingInfo.h" #include "Herwig/Shower/Dipole/Kinematics/DipoleSplittingKinematics.h" #include "ThePEG/EventRecord/RhoDMatrix.h" #include "Herwig/Decay/DecayMatrixElement.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" namespace Herwig { using namespace ThePEG; /** * \ingroup DipoleShower * \author Simon Platzer * * \brief DipoleSplittingKernel is the base class for all kernels * used within the dipole shower. * * @see \ref DipoleSplittingKernelInterfaces "The interfaces" * defined for DipoleSplittingKernel. */ class DipoleSplittingKernel: public HandlerBase { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ DipoleSplittingKernel(); /** * The destructor. */ virtual ~DipoleSplittingKernel(); //@} public: /** * Return the alpha_s to be used */ Ptr::tptr alphaS() const { return theAlphaS; } /** * Set the alpha_s to be used */ void alphaS(Ptr::tptr ap) { theAlphaS = ap; } /** * Return the splitting kinematics object */ Ptr::tptr splittingKinematics() const { return theSplittingKinematics; } /** * Return the mc check object */ Ptr::ptr mcCheck() const { return theMCCheck; } /** * Set the splitting kinematics object */ void splittingKinematics(Ptr::tptr sp) { theSplittingKinematics = sp; } /** * Return the PDFRatio object */ Ptr::tptr pdfRatio() const { return thePDFRatio; } /** * Set the PDFRatio object */ void pdfRatio(Ptr::tptr sp) { thePDFRatio = sp; } /** * Return the number of additional parameter * random variables needed to evaluate this kernel * except the momentum fractions of incoming partons. * These will be accessible through the * lastSplittingParameters() container of the splitting * info object. */ virtual int nDimAdditional() const { return 0; } /** * Set the freezing value for the renormalization scale */ void renormalizationScaleFreeze(Energy s) { theRenormalizationScaleFreeze = s; } /** * Set the freezing value for the factorization scale */ void factorizationScaleFreeze(Energy s) { theFactorizationScaleFreeze = s; } /** * Get the freezing value for the renormalization scale */ Energy renormalizationScaleFreeze() const { return theRenormalizationScaleFreeze; } /** * Get the freezing value for the factorization scale */ Energy factorizationScaleFreeze() const { return theFactorizationScaleFreeze; } public: /** * Return true, if this splitting kernel * applies to the given dipole index. */ virtual bool canHandle(const DipoleIndex&) const = 0; /** * Return true, if this splitting kernel is * the same for the given index a, as the given * splitting kernel for index b. */ virtual bool canHandleEquivalent(const DipoleIndex& a, const DipoleSplittingKernel& sk, const DipoleIndex& b) const = 0; /** * Return the emitter data after splitting, given * a dipole index. */ virtual tcPDPtr emitter(const DipoleIndex&) const = 0; /** * Return the emission data after splitting, given * a dipole index. */ virtual tcPDPtr emission(const DipoleIndex&) const = 0; /** * Return the spectator data after splitting, given * a dipole index. */ virtual tcPDPtr spectator(const DipoleIndex&) const = 0; /** * Return the flavour produced, if this cannot * be determined from the dipole. */ PDPtr flavour() const { return theFlavour; } /** * Return true, if this splitting kernel is supposed to work in a * strict large-N limit, i.e. replacing C_F by C_A/2 */ bool strictLargeN() const { return theStrictLargeN; } public: /** * Inform this splitting kernel, that it is being * presampled until a call to stopPresampling */ virtual void startPresampling(const DipoleIndex&) { presampling = true; } /** * Inform this splitting kernel, that it is not being * presampled until a call to startPresampling */ virtual void stopPresampling(const DipoleIndex&) { presampling = false; } /** * Return the number of points to presample this * splitting generator. */ unsigned long presamplingPoints() const { return thePresamplingPoints; } /** * Return the maximum number of trials * to generate a splitting. */ unsigned long maxtry() const { return theMaxtry; } /** * Return the number of accepted points after which the grid should * be frozen */ unsigned long freezeGrid() const { return theFreezeGrid; } /** * Set the number of accepted points after which the grid should * be frozen */ void freezeGrid(unsigned long n) { theFreezeGrid = n; } /** * Set a detuning factor to be applied to the sampling overestimate kernel */ void detuning(double d) { theDetuning = d; } /** * Return the detuning factor applied to the sampling overestimate kernel */ double detuning() const { return theDetuning; } /** * Evaluate this splitting kernel for the given * dipole splitting. */ virtual double evaluate(const DipoleSplittingInfo&) const = 0; /** * Evaluate rho_ii' V_ijk V*_i'jk / equivalent for initial-state splitting, * required for generating spin-correlated azimuthal angles. **/ virtual vector< pair > generatePhi( const DipoleSplittingInfo& dInfo, const RhoDMatrix& rho) const = 0; /** * Return the completely spin-unaveraged (i.e. spin-indexed) splitting kernel. **/ virtual DecayMEPtr matrixElement(const DipoleSplittingInfo& dInfo) const = 0; /** * Clear the alphaPDF cache */ void clearAlphaPDFCache() const { theAlphaSCache.clear(); thePDFCache.clear(); } /** * Update the variations vector at the given splitting using the indicated * kernel and overestimate values. */ virtual void accept(const DipoleSplittingInfo&, double, double, map&) const; /** * Update the variations vector at the given splitting using the indicated * kernel and overestimate values. */ virtual void veto(const DipoleSplittingInfo&, double, double, map&) const; /** * Return true, if this kernel is capable of * delivering an overestimate to the kernel, and * of inverting the integral over the overestimate * w.r.t. the phasepsace provided by the given * DipoleSplittingInfo object. */ virtual bool haveOverestimate(const DipoleSplittingInfo&) const { return false; } /** * Return the overestimate to this splitting kernel * for the given dipole splitting. */ virtual double overestimate(const DipoleSplittingInfo&) const { return -1.; } /** * Invert the integral over the overestimate * w.r.t. the phasepsace provided by the given * DipoleSplittingInfo object to equal * the given value. */ virtual double invertOverestimateIntegral(const DipoleSplittingInfo&, double) const { return -1.; } /** * . */ bool useThisKernel() const { return theUseThisKernel; } public: /** * Get the factorization scale factor */ double factorizationScaleFactor() const { return theFactorizationScaleFactor; } /** * Set the factorization scale factor */ void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; } /** * Get the renormalization scale factor */ double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; } /** * Set the renormalization scale factor */ void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; } + /** + * Return the CMW sheme used by the kernel + */ + unsigned int cmwScheme() const {return theCMWScheme;} + protected: /** * Return the common factor of (alphas/2pi)*(pdf ratio) */ double alphaPDF(const DipoleSplittingInfo&, Energy optScale = ZERO, double rScaleFactor = 1.0, double fScaleFactor = 1.0) const; /** * Return true, if the virtuality of the splitting should be used as the * argument of alphas rather than the pt */ bool virtualitySplittingScale() const { return theVirtualitySplittingScale; } public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The alpha_s to be used. */ Ptr::ptr theAlphaS; /** * An optional 'colour screening' scale * for alternative intrinsic pt generation. */ Energy theScreeningScale; /** * The splitting kinematics to be used. */ Ptr::ptr theSplittingKinematics; /** * An optional PDF ratio object to be used * when evaluating this kernel. */ Ptr::ptr thePDFRatio; /** * The number of points to presample this * splitting generator. */ unsigned long thePresamplingPoints; /** * The maximum number of trials * to generate a splitting. */ unsigned long theMaxtry; /** * Return the number of accepted points after which the grid should * be frozen */ unsigned long theFreezeGrid; /** * The detuning factor applied to the sampling overestimate kernel */ double theDetuning; /** * The flavour produced, if this cannot * be determined from the dipole. */ PDPtr theFlavour; /** * Pointer to a check histogram object */ Ptr::ptr theMCCheck; /** * True, if this splitting kernel is supposed to work in a * strict large-N limit, i.e. replacing C_F by C_A/2 */ bool theStrictLargeN; /** * The factorization scale factor. */ double theFactorizationScaleFactor; /** * The renormalization scale factor. */ double theRenormalizationScaleFactor; /** * A freezing value for the renormalization scale */ Energy theRenormalizationScaleFreeze; /** * A freezing value for the factorization scale */ Energy theFactorizationScaleFreeze; /** * True, if the virtuality of the splitting should be used as the * argument of alphas rather than the pt */ bool theVirtualitySplittingScale; /** * Implementing CMW in the kernels. **/ unsigned int theCMWScheme=0; /** * Cache for alphas evaluations */ mutable map theAlphaSCache; /** * Cache for PDF evaluations */ mutable map thePDFCache; /** * True, if we are presampling */ bool presampling; /** * True, if the kernel should be used */ bool theUseThisKernel = true; private: /** * The static object used to initialize the description of this class. * Indicates that this is an abstract class with persistent data. */ static AbstractClassDescription initDipoleSplittingKernel; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ DipoleSplittingKernel & operator=(const DipoleSplittingKernel &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of DipoleSplittingKernel. */ template <> struct BaseClassTrait { /** Typedef of the first base class of DipoleSplittingKernel. */ typedef HandlerBase NthBase; }; /** This template specialization informs ThePEG about the name of * the DipoleSplittingKernel class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::DipoleSplittingKernel"; } /** * The name of a file containing the dynamic library where the class * DipoleSplittingKernel is implemented. It may also include several, space-separated, * libraries if the class DipoleSplittingKernel depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_DipoleSplittingKernel_H */ diff --git a/Shower/Dipole/Merging/Merger.cc b/Shower/Dipole/Merging/Merger.cc --- a/Shower/Dipole/Merging/Merger.cc +++ b/Shower/Dipole/Merging/Merger.cc @@ -1,1614 +1,1479 @@ // -*- C++ -*- // // Merger.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL , see COPYING for details. // Please respect the MCnet academic guidelines , see GUIDELINES for details. // // // This is the implementation of the non-inlined , non-templated member // functions of the Merger class. // #include "Merger.h" #include "Node.h" #include "MergingFactory.h" // other includes when needed below. using namespace Herwig; IBPtr Merger::clone() const { return new_ptr( *this ); } IBPtr Merger::fullclone() const { return new_ptr( *this ); } - - - - namespace { double decideClustering(const NodePtr sub,const NodePtr head,bool& pro){ if( sub != head ){// at least one history step -> unitarisation if ( UseRandom::rndbool() ){ pro = true; return -2.; } else{ pro = false; return 2.; } } // no ordered history -> no projection else{ pro = false; return 1.; } } } CrossSection Merger::MergingDSigDRBornStandard( ){ // get the history for the process const NodePtr productionNode = currentNode()-> getHistory( true, DSH()->hardScaleFactor() ); // decide if to cluster - weight = decideClustering(productionNode, currentNode(), projected); + double weight = decideClustering(productionNode, currentNode(), projected); // check if we only want to calculate the current multiplicity. if(notOnlyMulti()) return ZERO; // Check for cuts on the production proces. if ( !productionNode->xcomb()->willPassCuts() ) return ZERO; // calculate the staring scale for the production node Energy startscale = CKKW_StartScale( productionNode ); // fill history with caluclation of sudakov supression fillHistory( startscale , productionNode , currentNode() ); // fill the projector -> return the scale of the last splitting currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) ); // the weight has three components to get shower history weight weight *= history.back().weight* // Sudakov suppression alphaReweight()* // alpha_s reweight pdfReweight(); // pdf reweight // If weight is zero return. if( weight == ZERO ) return ZERO; //calculate the cross section - return weight*TreedSigDR( startscale , 1. ); + return weight*TreedSigDR( startscale ); } CrossSection Merger::MergingDSigDRVirtualStandard( ){ // get the history for the process const NodePtr productionNode = currentNode()-> getHistory( true , DSH()->hardScaleFactor() ); // decide if to cluster - weight = decideClustering(productionNode,currentNode(),projected); + double weight = decideClustering(productionNode,currentNode(),projected); // Check for cuts on the production proces. if ( !productionNode->xcomb()->willPassCuts() )return ZERO; // calculate the staring scale Energy startscale = CKKW_StartScale( productionNode ); // fill history with caluclation of sudakov supression fillHistory( startscale , productionNode , currentNode() ); // fill the projector -> return the scale of the last splitting currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) ); // the weight has three components to get shower history weight double ww1 = history.back().weight; double ww2 = alphaReweight(true); double ww3 = pdfReweight(); weight *= ww1*ww2*ww3; // If weight is zero return. if( weight == 0. )return ZERO; // calculate the cross section for virtual contribution // and various insertion operators. CrossSection matrixElement = LoopdSigDR( startscale ); // Now calculate the expansion of the shower history. // first: the born contibution: CrossSection bornWeight = currentME()->dSigHatDRB(); // second: expansion of pdf ,alpha_s-ratio and sudakov suppression. double w1 = -sumPdfReweightExpansion(); double w2 = -sumAlphaSReweightExpansion(); double w3 = -sumFillHistoryExpansion(); // put together the expansion weights. CrossSection expansionweight = bornWeight*SM().alphaS()/( 2.*ThePEG::Constants::pi ); if (theShowerExpansionWeights == 0){ expansionweight *=0.; }else if ( theShowerExpansionWeights == 1 ){ expansionweight *=w1+w2+w3; }else if ( theShowerExpansionWeights == 2 ){ expansionweight *=w1+w2+w3*pow(as( startscale*DSH()->renFac() )/SM().alphaS(),currentME()->orderInAlphaS())/ww2; }else if ( theShowerExpansionWeights == 3 ){ expansionweight *=(w1+w2+w3)*pow(as( startscale*DSH()->renFac() )/SM().alphaS(),currentME()->orderInAlphaS())/ww2; }else if ( theShowerExpansionWeights == 4 ){ expansionweight *= w1+w3+w2*pow(as( startscale*DSH()->renFac() )/SM().alphaS(),currentME()->orderInAlphaS())/ww2; }else assert(false && theShowerExpansionWeights); // [ DEBUG ] if ( currentNode()->legsize() == 5 && Debug::level > 2 ) debugVirt(weight,w1,w2,w3,matrixElement,ww1,ww2,ww3,productionNode,bornWeight); // return with correction that ME was calculated with fixed alpha_s return weight* as( startscale*DSH()->renFac() )/ SM().alphaS()* ( matrixElement+expansionweight ); } CrossSection Merger::MergingDSigDRRealStandard(){ if ( currentNode()->children().empty() ) { throw Exception() << "Real emission contribution without underlying born." << "These are finite contibutions already handled in LO merging." << Exception::abortnow; } // check for IR Safe Cutoff if( !currentNode()->allAbove( theIRSafePT ) )return ZERO; auto inOutPair = currentNode()->getInOut(); NodePtr randomChild = currentNode()->randomChild(); bool meRegion =matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ); if ( meRegion )return MergingDSigDRRealAllAbove( ); if ( UseRandom::rndbool() ) return 2.*MergingDSigDRRealBelowSubReal( ); return 2.*MergingDSigDRRealBelowSubInt( ); } CrossSection Merger::MergingDSigDRRealAllAbove( ){ //If all dipoles pts are above , we cluster to this dipole. NodePtr CLNode = currentNode()->randomChild(); // Check if phase space poing is in ME region--> else rm PSP if ( !CLNode->children().empty() ) { auto inOutPair = CLNode->getInOut(); NodePtr randomChild = CLNode->randomChild(); if( !matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ) )return ZERO; } // first find the history for the acctual Node NodePtr productionNode = currentNode()-> getHistory( true , DSH()->hardScaleFactor() ); // If CLNode is not part of the history , dont calculate the Real contribution // else multiply the real contribution with N ( number of children ). // this returns the sudakov suppression according to // the clustering of the born parts. bool inhist = CLNode->isInHistoryOf( productionNode ); if(productionNode== currentNode())assert(!inhist); // get the history for the clustered Node. productionNode = CLNode-> getHistory( false , DSH()->hardScaleFactor() ); // decide if to cluster - weight = decideClustering(productionNode,CLNode,projected); + double weight = decideClustering(productionNode,CLNode,projected); // Check for cuts on the production process. if ( !productionNode->xcomb()->willPassCuts() )return ZERO; // calculate the staring scale Energy startscale = CKKW_StartScale( productionNode ); // fill history with caluclation of sudakov supression fillHistory( startscale , productionNode , CLNode ); // fill the projector -> return the scale of the last splitting currentNode()->runningPt( fillProjector( projected ? 2 : 1 ) ); // the weight has three components to get shower history weight weight *= history.back().weight*alphaReweight(true)*pdfReweight(); if( weight == 0. )return ZERO; // The inhist flag produces the correct cluster density. CrossSection me = ( inhist?TreedSigDR( startscale ):ZERO ); // calculate the dipole CrossSection dip = CLNode->calcDip( startscale* currentME()->renFac() ); CrossSection res = weight*as( startscale*DSH()->renFac() )/SM().alphaS()* currentNode()->children().size()* ( me - dip ); // [ DEBUG ] if ( currentNode()->legsize() == 6&&Debug::level > 2 ) debugReal("RAA",weight,me,dip); return res; } CrossSection Merger::MergingDSigDRRealBelowSubReal( ){ // Choose a random child to produce history from. NodePtr HistNode = currentNode()->randomChild(); // Check that this subleading point is in ME region. if ( !HistNode->children().empty() ) { auto inOutPair = HistNode->getInOut(); NodePtr randomChild = HistNode->randomChild(); // Here we make sure that clustering and splitting are invertible if( !matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ) )return ZERO; } // As the HistNode is now in ME region, we can cluster according to LO merging. // In this clustering we do not require a orering for the last step to the currentNode. const NodePtr productionNode = HistNode-> getHistory( false , DSH()->hardScaleFactor() ); // If the real emission contribution should be unitarised, we decide here if we cluster. // This applies to NLO corrections for the first emission w.r.t. the production process the first time. - weight = decideClustering(productionNode,HistNode,projected); + double weight = decideClustering(productionNode,HistNode,projected); // The production node needs to fulfill the cut criterion of the production process. if ( !productionNode->xcomb()->willPassCuts() && !currentNode()->xcomb()->willPassCuts() )return ZERO; // Calculate the starting scale w.r.t. the production process. Energy startscale = CKKW_StartScale( productionNode ); // If the production process does not fullfill the cut criterion use the real emission point (DEBUG) if (!productionNode->xcomb()->willPassCuts()) startscale = CKKW_StartScale( currentNode() ); // DEBUG trial //currentNode()->xcomb()->lastProjector( productionNode->xcomb()); // Calculate the sudakov weights starting from the production node to the histNode fillHistory( startscale , productionNode , HistNode ); // Set the running Pt of the process as it is used in the vetoed parton shower. currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) ); currentNode()->runningPt(max(HistNode->pT(),theMergePt)); // Calculate the alpha_S ratios and pdf ratios. weight *= history.back().weight*alphaReweight(true)*pdfReweight(); if( weight == 0. )return ZERO; // Start calculation of subtraction contribution. CrossSection sumPS = ZERO; // Iterate over all subtraction contributions. for( auto const & child : currentNode()->children() ){ if ( child->allAbove( mergePt() ) && child->xcomb()->willPassCuts() ){ Energy relevantScale=child->children().empty()?CKKW_StartScale( child ):child->maxChildPt(); if( ( child )->pT()>mergePt()/theRealSubtractionRatio ){ if( child ->pT()inShowerPS(relevantScale))){ //DEBUG: CKKW_StartScale( child);???? sumPS += child->calcPs( startscale* currentME()->renFac() ); } }else{ if( child ->pT()calcDip( startscale* currentME()->renFac() ); } } } } CrossSection me = ZERO; if(currentNode()->xcomb()->willPassCuts()){ me = TreedSigDR( startscale ); } // [ DEBUG ] if ( currentNode()->legsize() == 6&&Debug::level > 2 ) debugReal("RBSR",weight,me,sumPS); //Here we subtract the PS ( and below the dynamical cutoff the Dip ) return weight*as( startscale*DSH()->renFac() )/SM().alphaS()* ( me-sumPS ); } CrossSection Merger::MergingDSigDRRealBelowSubInt( ){ if( currentNode()->children().empty() )return ZERO; NodePtr CLNode = currentNode()->randomChild(); if( CLNode->pT()children().empty() ) { auto inOutPair = CLNode->getInOut( ); NodePtr randomChild = CLNode->randomChild(); // Here we make sure that clustering and splitting are invertible if( !matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ) )return ZERO; } const NodePtr productionNode = CLNode-> getHistory( false , DSH()->hardScaleFactor() ); - weight = decideClustering(productionNode,CLNode,projected); + double weight = decideClustering(productionNode,CLNode,projected); if ( !CLNode->allAbove( mergePt() ) )return ZERO; if ( !productionNode->xcomb()->willPassCuts() )return ZERO; Energy startscale = CKKW_StartScale( productionNode ); fillHistory( startscale , productionNode , CLNode ); currentNode()->runningPt( fillProjector( projected ? 2 : 1 ) ); weight *= history.back().weight*alphaReweight(true)*pdfReweight(); if( weight == 0. )return ZERO; pair DipAndPs = CLNode->calcDipandPS( startscale* currentME()->renFac() ); // [ DEBUG ] if ( currentNode()->legsize() == 6&&Debug::level > 2 ) debugReal("RBSI",weight,DipAndPs.second,DipAndPs.first); Energy relevantScale=CLNode->children().empty()?CKKW_StartScale( CLNode ):CLNode->maxChildPt(); bool calcPScontribution=CLNode->pT()inShowerPS(relevantScale)); //Here we add the PS and subtrac the Dip ( only above the dynamical cutoff ) return weight*as( startscale*DSH()->renFac() )/SM().alphaS()* currentNode()->children().size()*( (calcPScontribution?DipAndPs.second:ZERO)-DipAndPs.first ); } -CrossSection Merger::MergingDSigDRBornGamma( ){ - - double weightCL = 0.; - weight = 1.; - - if ( !currentNode()->children().empty() ) { - auto const inOutPair = currentNode()->getInOut(); - // Here we make sure that clustering and splitting are invertible. - NodePtr randomChild = currentNode()->randomChild(); - // Check if point is part of the ME region. - if( !matrixElementRegion( inOutPair.first , - inOutPair.second , - randomChild->pT() , - theMergePt ) )weight *= 0.; - } - - const NodePtr productionNode = currentNode()->getHistory( true , DSH()->hardScaleFactor() ); - NodePtr CLNode; - NodePtr BornCL; - - - if( !currentNode()->children().empty() ){ - if ( UseRandom::rndbool() ){ - CLNode = currentNode()->randomChild(); - bool inhist = CLNode->isInHistoryOf( productionNode ); - weight *= inhist?( -2.*int( currentNode()->children().size() ) ):0.; - projected = true; - weightCL = 2.*int( currentNode()->children().size() ); - BornCL = CLNode-> getHistory( false , DSH()->hardScaleFactor() ); - }else{ - weight = 2.; - projected = false; - } - }else{ - weight = 1.; - projected = false; - } - - - if ( treefactory()->onlymulti() != -1&& - treefactory()->onlymulti() != - int( currentNode()->legsize()-(projected ? 1 : 0) ) ) - return ZERO; - - - if( !currentNode()->allAbove( mergePt()*(1.-1e-6) ) )weight = 0.; - if( CLNode&&!CLNode->allAbove( mergePt()*(1.-1e-6) ) )weightCL = 0.; - if ( !productionNode->xcomb()->willPassCuts() ){ - return ZERO; - } - - CrossSection res = ZERO; - bool maxMulti = currentNode()->legsize() == int( maxLegsLO() ); - - - if( weight != 0. ){ - Energy startscale = CKKW_StartScale( productionNode ); - fillHistory( startscale , productionNode , currentNode() ); - currentNode()->runningPt( fillProjector( (projected ? 1 : 0) ) ); - weight *= history.back().weight*alphaReweight()*pdfReweight(); - if( weight == 0.&&weightCL == 0. )return ZERO; - - res += weight*TreedSigDR( startscale , ( !maxMulti&&!projected )?theGamma:1. ); - } - - if( CLNode&&theGamma != 1. ){ - Energy startscale = CKKW_StartScale( BornCL ); - fillHistory( startscale , BornCL , CLNode ); - currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) ); - weightCL *= history.back().weight*alphaReweight()*pdfReweight(); - CrossSection diff = ZERO; - currentME()->factory()->setAlphaParameter( 1. ); - diff -= weightCL*CLNode->dipole()->dSigHatDR( sqr( startscale* currentME()->renFac() ) ); - currentME()->factory()->setAlphaParameter( theGamma ); - - string alp = ( CLNode->dipole()->aboveAlpha()?"Above":"Below" ); - - diff += weightCL*CLNode->dipole()->dSigHatDR( sqr( startscale* currentME()->renFac() ) ); - currentME()->factory()->setAlphaParameter( 1. ); - - res += diff; - } - return res; -} - - - - -CrossSection Merger::TreedSigDR( Energy startscale , double diffAlpha ){ +CrossSection Merger::TreedSigDR( Energy startscale ){ currentME()->setScale( sqr( startscale ) , sqr( startscale ) ); CrossSection res = currentME()->dSigHatDRB(); - /*bool useDipolesForME=false; - if (useDipolesForME && !currentNode()->children().empty()){ - res=ZERO; - for (auto const & child : currentNode()->children() ) - res-=child->dipole()->dSigHatDR(sqr( startscale )); - } - */ - if ( projected && emitDipoleMEDiff ) { - CrossSection resDip=ZERO; - for (auto const & child : currentNode()->children() ) - resDip-=child->dipole()->dSigHatDR(sqr( startscale )); - setEmissionProbability(1.-min(resDip/res,res/resDip)); - }else{ - setEmissionProbability(0.); - } - - if ( diffAlpha != 1. ) { - res += currentME()->dSigHatDRAlphaDiff( diffAlpha ); - } if( std::isnan( double( res/nanobarn ) ) ){ generator()->logWarning(Exception() << "Merger: TreedSigDR is nan" << Exception::warning); res = ZERO;}; return res; } CrossSection Merger::LoopdSigDR( Energy startscale ){ currentME()->setScale( sqr( startscale ) , sqr( startscale ) ); currentME()->doOneLoopNoBorn(); CrossSection res = currentME()->dSigHatDRV()+ currentME()->dSigHatDRI(); currentME()->noOneLoopNoBorn(); return res; } Energy Merger::fillProjector( int pjs ){ // in the shower handler the scale is multiplied // by DSH()->hardScaleFactor() so here we need // to devide by the factor. double xiQSh = history.begin()->node->legsize() == N0()?DSH()->hardScaleFactor():1.; if( pjs == 0 ){ return ( history.size() == 1?1.:( 1./xiQSh ) )*history.back().scale; } for( auto const & hs : history ) if ( isProjectorStage( hs.node , pjs )&&pjs != 0 ){ currentNode()->xcomb()->lastProjector( hs.node->xcomb() ); return ( hs.node == history[0].node?1.:( 1./xiQSh ) )*hs.scale; } throw Exception() << "Could not fill projector." << Exception::abortnow; return ZERO; } double Merger::pdfReweight(){ // TODO factorization scale inside double res = 1.; // consider both sides. for( int side : {0 , 1} ){ // The side scale defines the scale that the leg is changig thou the history. // We start reweighting at the seed process. // We only need to calculate the pdf if the emission whould change the leg, // otherwise the leg remains the same. // If no emission is prduced from this leg the pdf ratios from this leg are always 1. Energy sidescale=history[0].scale*( history[0].node->legsize() == N0() ? currentME()->facFac(): DSH()->facFac()); bool sidechanged=false; // only if the incoming parton is coloured. if( history[0].node->xcomb()->mePartonData()[side]->coloured() ){ // go though the history. for ( auto const & hs : history ){ //pdf-ratio only to the last step if ( !hs.node->parent() ) continue; if ( hs.node == history.back().node ) continue; if ( !hs.node->dipole() ){ throw Exception() << "\nMerger: pdfReweight: history step has no dipol. " << Exception::abortnow; return 0.; } // if the emitter is the side the emission changes the momentum fraction. if(!(hs.node->dipole()->bornEmitter()==side|| // II dipoles dont change the momentum fraction of the spectator only FI ( hs.node->dipole()->bornSpectator()==side && hs.node->dipole()->bornEmitter()>1))) continue; const bool fromIsME = false; const bool toIsME = history[0].node == hs.node; res *= pdfratio( hs.node, //numerator sidescale, // denominator DSH()->facFac()*hs.node->pT(), side, fromIsME, toIsME); sidescale= DSH()->facFac()*hs.node->pT(); sidechanged=true; } const bool fromIsME = true; const bool toIsME = !sidechanged && history[0].node->legsize() == N0(); res *= pdfratio( history.back().node, sidescale, history[0].scale * currentME()->facFac(), side, fromIsME, toIsME ); } } if ( std::isnan( res ) ) generator()->logWarning(Exception() << "Merger: pdfReweight is nan." << Exception::warning); return res; } double Merger::cmwAlphaS(Energy q)const{ using Constants::pi; // No cmw-scheme if (theCMWScheme==0) return as( q ); // Linear cmw-scheme else if(theCMWScheme==1){ double als=as( q ); return als * (1.+(3.*(67./18.-1./6.*sqr(pi)) -5./9.*Nf(q))* als/2./pi); } // cmw-scheme as factor in argument. else if(theCMWScheme==2){ double kg=exp(-(67.-3.*sqr(pi)-10/3*Nf(q)) /( 2. *(33.-2.*Nf(q)))); //Note factor 2 since we here dealing with Energy return as(max(kg*q,1_GeV)); }else{ throw Exception() << "This CMW-Scheme is not implemented." << Exception::abortnow; } return -1; } double Merger::alphaReweight(bool nocmw){ double res = 1.; Energy Q_R = ( history[0].node->legsize() == N0()? currentME()->renFac(): DSH()->renFac() )* history[0].scale; using Constants::pi; const auto Q_qed=history[0].node->nodeME()->factory()->scaleChoice()->renormalizationScaleQED(); const auto Oew=history[0].node->nodeME()->orderInAlphaEW(); const auto Oqcd=history[0].node->nodeME()->orderInAlphaS(); if (!history[0].node->children().empty()) { assert(Oqcd!=0); } res *= pow( SM().alphaEMME( Q_qed )/ SM().alphaEMMZ() , Oew ); res *= pow( (nocmw?as(Q_R):cmwAlphaS(Q_R)) / SM().alphaS() , Oqcd ); for ( auto const & hs : history ) if ( hs.node!= history.back().node ){ Energy q_i = DSH()->renFac()* hs.node->pT(); res *= cmwAlphaS(q_i) / SM().alphaS(); } if ( std::isnan( res ) ) generator()->logWarning(Exception() << "Merger: alphaReweight is nan. "<< Exception::warning); return res; } void Merger::fillHistory( Energy scale , NodePtr begin , NodePtr endNode ){ history.clear(); double sudakov0_n = 1.; history.push_back( {begin , sudakov0_n , scale} ); double xiQSh = history.begin()->node->legsize() == N0()? DSH()->hardScaleFactor():1.; scale *= xiQSh; if ( begin->parent()||!isUnitarized ) { while ( begin->parent() && ( begin != endNode ) ) { if ( !dosudakov( begin , scale , begin->pT() , sudakov0_n ) ){ history.push_back( { begin->parent() , 0. , scale } ); } if ( std::isnan( sudakov0_n ) ) generator()->logWarning(Exception() << "Merger: sudakov"<pT()/GeV<<"0_n is nan. " << Exception::warning); scale = begin->pT(); history.push_back( { begin->parent() , sudakov0_n , begin->pT() } ); begin = begin->parent(); } Energy notunirunning = scale; if ( !isUnitarized&&N()+N0() > int( currentNode()->legsize() ) ) { if ( !dosudakov( begin , notunirunning , mergePt() , sudakov0_n ) ){ history.back().weight = 0.; }else{ history.back().weight = sudakov0_n; } } } if( history.size() == 1 )scale /= DSH()->hardScaleFactor(); } double Merger::sumPdfReweightExpansion()const{ double res = 0.; Energy beam1Scale = history[0].scale* ( history[0].node->legsize() == N0()? currentME()->facFac(): DSH()->facFac() ); Energy beam2Scale = history[0].scale* ( history[0].node->legsize() == N0()? currentME()->facFac(): DSH()->facFac() ); for ( auto const & hs : history ){ //pdf expansion only to the last step if( !hs.node->parent() )continue; if( hs.node->xcomb()->mePartonData()[0]->coloured()&& hs.node->nodeME()->lastX1()>0.99 )return 0.; if( hs.node->xcomb()->mePartonData()[1]->coloured()&& hs.node->nodeME()->lastX2()>0.99 )return 0.; if( hs.node->nodeME()->lastX1()<0.00001 )return 0.; if( hs.node->nodeME()->lastX2()<0.00001 )return 0.; if ( hs.node->dipole()->bornEmitter() == 0 ){ res += pdfExpansion( hs.node , 0 , beam1Scale , ( hs.node->pT() ) , hs.node->nodeME()->lastX1() , Nf( history[0].scale ) , history[0].scale ); beam1Scale = ( hs.node->pT() )*DSH()->facFac(); } else if ( hs.node->dipole()->bornEmitter() == 1 ){ res += pdfExpansion( hs.node , 1 , beam2Scale , ( hs.node->pT() ) , hs.node->nodeME()->lastX2() , Nf( history[0].scale ) , history[0].scale ); beam2Scale = ( hs.node->pT() )*DSH()->facFac(); } // if we're here we know hs.node->dipole()->bornEmitter() > 1 // works only in collinear scheme else if ( hs.node->dipole()->bornSpectator() == 0 ){ res += pdfExpansion( hs.node , 0 , beam1Scale , ( hs.node->pT() ) , hs.node->nodeME()->lastX1() , Nf( history[0].scale ) , history[0].scale ); beam1Scale = ( hs.node->pT() )*DSH()->facFac(); } else if ( hs.node->dipole()->bornSpectator() == 1 ){ res += pdfExpansion( hs.node , 1 , beam2Scale , ( hs.node->pT() ) , hs.node->nodeME()->lastX2() , Nf( history[0].scale ) , history[0].scale ); beam2Scale = ( hs.node->pT() )*DSH()->facFac(); } } if ( currentNode()->xcomb()->mePartonData()[0]->coloured() ){ res += pdfExpansion( history.back().node , 0 , beam1Scale , history[0].scale* currentME()->facFac() , ( history.back() ).node->nodeME()->lastX1() , Nf( history[0].scale ) , history[0].scale ); } if ( currentNode()->xcomb()->mePartonData()[1]->coloured() ) { res += pdfExpansion( history.back().node , 1 , beam2Scale , history[0].scale* currentME()->facFac() , ( history.back() ).node->nodeME()->lastX2() , Nf( history[0].scale ) , history[0].scale ); } return res; } #include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h" double Merger::pdfExpansion( NodePtr node , int side , Energy running , Energy next , double x , int nlp , Energy fixedScale ) const { tcPDPtr particle , parton; tcPDFPtr pdf; if ( side == 0 ) { particle = node->nodeME()->lastParticles().first->dataPtr(); parton = node->nodeME()->lastPartons().first->dataPtr(); pdf = node->xcomb()->partonBins().first->pdf(); }else{ assert( side == 1 ); particle = node->nodeME()->lastParticles().second->dataPtr(); parton = node->nodeME()->lastPartons().second->dataPtr(); pdf = node->xcomb()->partonBins().second->pdf(); } //copied from PKOperator double res = 0.; int number = 10; for ( int nr = 0;nr 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; static double CA = SM().Nc(); static 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; // Pqg 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; } } // Pgg 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.; else res += 1*restmp*log( sqr( running/next ) )/PDFxparton*mapz; } return res/number; } double Merger::sumAlphaSReweightExpansion()const{ double res = 0.; const auto Oqcd=history[0].node->nodeME()->orderInAlphaS(); res += alphasExpansion( history[0].scale* DSH()->renFac() , history[0].scale* currentME()->renFac() )* Oqcd; // dsig is calculated with fixed alpha_s for ( auto const & hs : history ){ //expansion only to the last step if( !hs.node->parent() )continue; res += alphasExpansion( hs.node->pT()*DSH()->renFac() ,currentME()->renFac()*history[0].scale ); } return res; } double Merger::sumFillHistoryExpansion(){ double res = 0.; double xiQSh = history[0].node->legsize() == N0()?DSH()->hardScaleFactor():1.; for ( auto const & hs : history ){ if( !hs.node->parent() )continue; doHistExpansion( hs.node , ( hs.node == history[0].node?xiQSh:1. )*hs.scale , hs.node->pT() , history[0].scale , res ); } return res; } MergingFactoryPtr Merger::treefactory()const{return theTreeFactory;} void Merger::doinit(){ if ( !DSH()->hardScaleIsMuF() ) { throw Exception() << "Merger: Merging is currently only sensible " << "if we are using the hardScale as MuF." << Exception::abortnow; } DSH()->init(); if ( !( DSH()->showerPhaseSpaceOption() == 0 || DSH()->showerPhaseSpaceOption() == 1) ){ throw InitException() << "Merger::doinit(): Choice of shower phase space cannot be handled by the merging"; } + + /** + * Make sure all kernels use the same CMW scheme and set the scheme used in + * merging. Should be solved by allowing the kernels to have a pointer to + * the Shower Handler. + */ + bool setCMW=false; + assert(DSH()->splittingKernels().size() > 0); + for ( const auto & sk : DSH()->splittingKernels() ){ + if ( not setCMW ){ + theCMWScheme=sk->cmwScheme(); + setCMW=true; + } + else if( theCMWScheme != sk->cmwScheme() ){ + throw InitException() << "Merger::doinit(): The merging currently only " + <<"supports if all kernels emit with the same" + <<"CMW scheme. The kernel "<name()<<" has a different convention."; + } + } + } namespace{ void setXCombScales(StdXCombPtr xc,Energy2 scale){ xc->lastShowerScale ( scale ); xc->partonBinInstances().first->scale ( scale ); xc->partonBinInstances().second->scale ( scale ); } } CrossSection Merger::MergingDSigDR() { history.clear(); assert(currentNode()==theFirstNodeMap[ currentME()]); if(DSH()->doesSplitHardProcess()){ throw Exception() << "Merger: The splithardprocess option is currently not supported." << Exception::abortnow; } //get the PDF's (from ShowerHandler.cc) if (!DSH()->getPDFA()||!DSH()->firstPDF().particle()){ tSubProPtr sub = currentNode()->xcomb()->construct(); const auto pb=currentNode()->xcomb()->partonBins(); tcPDFPtr first = DSH()->getPDFA() ? tcPDFPtr(DSH()->getPDFA()) : DSH()->firstPDF().pdf(); tcPDFPtr second = DSH()->getPDFB() ? tcPDFPtr(DSH()->getPDFB()) : DSH()->secondPDF().pdf(); DSH()->resetPDFs( {first,second },pb ); } DSH()->eventHandler( generator()->eventHandler() ); CrossSection res = ZERO; if( currentNode()->subtractedReal() ){ res = MergingDSigDRRealStandard(); theCurrentMaxLegs = maxLegsNLO(); }else if( currentNode()->virtualContribution() ){ res = MergingDSigDRVirtualStandard(); theCurrentMaxLegs = maxLegsNLO(); - }else if( theGamma!= 1. ){ - res = MergingDSigDRBornGamma(); - theCurrentMaxLegs = maxLegsLO(); }else{ res = MergingDSigDRBornStandard(); theCurrentMaxLegs = maxLegsLO(); } auto lxc= currentME()->lastXCombPtr(); setXCombScales(lxc,sqr( currentNode()->runningPt())); auto lp= currentME()->lastXCombPtr()->lastProjector(); if( lp ) setXCombScales( lp, sqr( currentNode()->runningPt())); if ( res == ZERO ){ history.clear(); return ZERO; } cleanup( currentNode() ); lxc->subProcess( SubProPtr() ); history.clear(); if( !std::isfinite( double( res/nanobarn ) ) ){ generator()->logWarning(Exception() << "Merger weight is " << res/nanobarn<< " -> setting to 0" << Exception::warning); return ZERO; } return res; } #include "Herwig/PDF/HwRemDecayer.h" void Merger::CKKW_PrepareSudakov( NodePtr node , Energy running ){ tSubProPtr sub = node->xcomb()->construct(); const auto pb=node->xcomb()->partonBins(); DSH()->setCurrentHandler(); DSH()->currentHandler()->generator()->currentEventHandler(currentNode()->xcomb()->eventHandlerPtr() ); DSH()->currentHandler()->remnantDecayer()->setHadronContent( currentNode()->xcomb()->lastParticles() ); DSH()->eventRecord().clear(); DSH()->eventRecord().slimprepare( sub , dynamic_ptr_cast( node->xcomb() ) , DSH()->pdfs() , currentNode()->xcomb()->lastParticles(), DSH()->offShellPartons() ); DSH()->hardScales( sqr( running ) ); } Energy Merger::CKKW_StartScale( NodePtr node ) const { Energy res = generator()->maximumCMEnergy(); const int N=node->legsize(); const auto & data=node->nodeME()->mePartonData(); const auto & momenta=node->nodeME()->lastMEMomenta(); if( !node->children().empty() ){ for ( int i = 0; icoloured() )continue; for ( int j = 2; jcoloured() )continue; for ( int k = 0; kcoloured() || i == k || j == k )continue; if(i<2){ if(k<2) { res = min( res , IILTK->lastPt( momenta[i] , momenta[j] , momenta[k] )); } else { res = min( res , (data[k]->mass()+data[j]->mass()+data[i]->mass()>ZERO)? IFMTK->lastPt( momenta[i] , momenta[j] , momenta[k] ): IFLTK->lastPt( momenta[i] , momenta[j] , momenta[k] )); } }else{ if(k<2) { res = min( res , (data[k]->mass()+data[j]->mass()+data[i]->mass()>ZERO)? FIMTK->lastPt( momenta[i] , momenta[j] , momenta[k] ): FILTK->lastPt( momenta[i] , momenta[j] , momenta[k] )); } else { res = min( res , (data[k]->mass()+data[j]->mass()+data[i]->mass()>ZERO)? FFMTK->lastPt( momenta[i] , momenta[j] , momenta[k] ): FFLTK->lastPt( momenta[i] , momenta[j] , momenta[k] )); } } } } } }else{ node->nodeME()->factory()->scaleChoice()->setXComb( node->xcomb() ); res = sqrt( node->nodeME()->factory()->scaleChoice()->renormalizationScale() ); } node->nodeME()->factory()->scaleChoice()->setXComb( node->xcomb() ); res = max( res , sqrt( node->nodeME()->factory()->scaleChoice()->renormalizationScale() ) ); return res; } double Merger::alphasExpansion( Energy next , Energy fixedScale ) const { double betaZero = ( 11./6. )*SM().Nc() - ( 1./3. )*Nf( history[0].scale ); double K=3.*( 67./18.-1./6.*sqr(Constants::pi) ) -5./9.*Nf( history[0].scale); return ( betaZero*log( sqr( fixedScale/next ) ) )+( theCMWScheme>0?K:0. ); } double Merger::pdfratio( NodePtr node , Energy numerator_scale , Energy denominator_scale , int side , bool fromIsME, bool toIsME ){ StdXCombPtr bXc = node->xcomb(); if( !bXc->mePartonData()[side]->coloured() ) throw Exception() << "Merger: pdf-ratio required for non-coloured particle." << Exception::abortnow; double from = 1.; double to = 1.; if ( side == 0 ){ if ( denominator_scale == numerator_scale && fromIsME==toIsME ) { return 1.; } if (fromIsME) { from = node->nodeME()->pdf1( sqr( denominator_scale ) ); }else{ from = DSH()->firstPDF().xfx(node->xcomb()->lastPartons().first->dataPtr(), sqr( denominator_scale ), node->xcomb()->lastX1())/node->xcomb()->lastX1(); } if (toIsME) { to = node->nodeME()->pdf1( sqr( numerator_scale ) ); }else{ to = DSH()->firstPDF().xfx(node->xcomb()->lastPartons().first->dataPtr(), sqr( numerator_scale ), node->xcomb()->lastX1())/node->xcomb()->lastX1(); } if ( ( to < 1e-8||from < 1e-8 )&&( to/from>10000000. ) ){ generator()->logWarning(Exception() << "Merger: pdfratio to = " << to << " from = " << from << Exception::warning); return 0.; } } else{ if ( denominator_scale == numerator_scale && fromIsME==toIsME ) { return 1.; } if (fromIsME) { from = node->nodeME()->pdf2( sqr( denominator_scale ) ); }else{ from =DSH()->secondPDF().xfx(node->xcomb()->lastPartons().second->dataPtr(), sqr( denominator_scale ), node->xcomb()->lastX2())/node->xcomb()->lastX2(); } if (toIsME) { to = node->nodeME()->pdf2( sqr( numerator_scale ) ); }else{ to = DSH()->secondPDF().xfx(node->xcomb()->lastPartons().second->dataPtr(), sqr( numerator_scale ), node->xcomb()->lastX2())/node->xcomb()->lastX2(); } if ( ( to < 1e-8||from < 1e-8 )&&( to/from>10000000. ) ){ generator()->logWarning(Exception() << "Merger: pdfratio to = " << to << " from = " << from << Exception::warning); return 0.;} } return to/from; } bool Merger::dosudakov( NodePtr node , Energy running , Energy next , double& sudakov0_n ) { CKKW_PrepareSudakov( node , running ); for( DipoleChain const & chain : DSH()->eventRecord().chains() ){ for( Dipole const & dip : chain.dipoles() ){ sudakov0_n *= singlesudakov( dip , next , running , { true , false } ); sudakov0_n *= singlesudakov( dip , next , running , { false , true } ); if ( sudakov0_n == 0.0 ){ cleanup( node ); return false; } } } cleanup( node ); return true; } bool Merger::doHistExpansion( NodePtr node , Energy running , Energy next , Energy fixedScale , double& histExpansion ) { CKKW_PrepareSudakov( node , running ); for( DipoleChain const & chain : DSH()->eventRecord().chains() ){ for( Dipole const & dip : chain.dipoles() ){ histExpansion += singleHistExpansion( dip , next , running , fixedScale , { true , false } );; histExpansion += singleHistExpansion( dip , next , running , fixedScale , { false , true } ); } } cleanup( node ); return true; } bool Merger::isProjectorStage( NodePtr node , int pjs )const{ return ( pjs == int( ( currentNode()->legsize() - node->legsize() ) ) ); } void Merger::cleanup( NodePtr node ) { DSH()->eventRecord().clear(); if( !node->xcomb()->subProcess() )return; ParticleVector vecfirst = node->xcomb()->subProcess()->incoming().first->children(); for( auto const & particle : vecfirst ) node->xcomb()->subProcess()->incoming().first->abandonChild( particle ); ParticleVector vecsecond = node->xcomb()->subProcess()->incoming().second->children(); for( auto const & particle : vecsecond ) node->xcomb()->subProcess()->incoming().second->abandonChild( particle ); node->xcomb()->subProcess( SubProPtr() ); } double Merger::singlesudakov( Dipole dip , Energy next , Energy running , pair 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() ); auto const & gens = DSH()->generators().equal_range( candidate.index() ); for ( auto 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 , *gen->second->splittingKernel() ); candidate.hardPt( min( running , ptMax ) ); if ( candidate.hardPt()>next ){ res *= gen->second->sudakov( candidate , next ); } } return res; } double Merger::singleHistExpansion( Dipole dip , Energy next , Energy running , Energy fixedScale , pair 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() ); auto const & gens = DSH()->generators().equal_range( candidate.index() ); for ( auto 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 , *gen->second->splittingKernel() ); candidate.hardPt( min( running , ptMax ) ); if ( candidate.hardPt()>next ){ res += gen->second->sudakovExpansion( candidate , next , fixedScale ); } } return res; } void Merger::firstNodeMap( MatchboxMEBasePtr a , NodePtr b ){ theFirstNodeMap.insert( { a , b } ); } map Merger::firstNodeMap()const{return theFirstNodeMap;} void Merger::setXComb( tStdXCombPtr xc ){ currentNode()->setXComb( xc ); } void Merger::setKinematics( ){ currentNode()->setKinematics(); } void Merger::clearKinematics( ){ currentNode()->clearKinematics(); } void Merger::flushCaches(){ if (currentNode()&¤tNode()->xcomb()->lastParticles().first) { currentNode()->flushCaches(); } } bool Merger::generateKinematics( const double * r ){ return currentNode()->firstgenerateKinematics( r , ! currentNode()->subtractedReal() ); } bool Merger::matrixElementRegion( PVector incoming , PVector outgoing , Energy winnerScale , Energy cutscale )const{ Energy ptx = Constants::MaxEnergy; bool foundwinnerpt = false; using namespace boost; //FF for( auto const & em : outgoing ){ if ( ! em->coloured() ) continue; for( auto const & emm : outgoing ){ if ( !emm->coloured() ) continue; if ( em == emm ) continue; for( auto const & spe : outgoing ){ if ( !spe->coloured() ) continue; if ( em == spe||emm == spe ) continue; if ( !( em->id() == -emm->id()||emm->id()>6 ) )continue; Energy pt = ZERO; if ( em->momentum().m()<= 10_MeV && emm->momentum().m()<= 10_MeV && spe->momentum().m()<= 10_MeV ) { pt = FFLTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() ); }else{ pt = FFMTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() ); } if ( abs( pt-winnerScale ) < 10_MeV ) { foundwinnerpt = true; } ptx = min( ptx , pt ); } } } //FI for( auto const & spe : incoming ){ if ( ! spe->coloured() ) continue; for( auto const & emm : outgoing ){ if ( ! emm->coloured() ) continue; for( auto const & em : outgoing ){ if ( ! em->coloured() ) continue; if ( em == emm ) continue; if ( !( em->id() == -emm->id() || emm->id()>6 ) )continue; Energy pt = ZERO; if ( em->momentum().m()<= 10_MeV && emm->momentum().m()<= 10_MeV && spe->momentum().m()<= 10_MeV ) { pt = FILTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() ); }else{ pt = FIMTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() ); } if ( abs( pt-winnerScale )<10_MeV ) { foundwinnerpt = true; } if( pt > ZERO ) ptx = min( ptx , pt ); } } } //IF for( auto const & em : incoming ){ if ( ! em->coloured() ) continue; for( auto const & emm : outgoing ){ if ( ! emm->coloured() ) continue; for( auto const & spe : outgoing ){ if ( ! spe->coloured() ) continue; if ( emm == spe ) continue; if ( !( em->id()>6|| em->id() == emm->id() ||emm->id()>6 ) )continue; Energy pt = ZERO; if ( em->momentum().m()<= 10_MeV && emm->momentum().m()<= 10_MeV && spe->momentum().m()<= 10_MeV ) { //massless pt = IFLTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() ); }else{ //massiv pt = IFMTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() ); } if ( abs( pt-winnerScale )< 10_MeV ) { foundwinnerpt = true; } ptx = min( ptx , pt ); } } } //II for( auto const & em : incoming ){ if ( ! em->coloured() ) continue; for( auto const & spe : incoming ){ if ( ! spe->coloured() ) continue; if ( em == spe ) continue; for( auto const & emm : outgoing ){ if ( ! emm->coloured() ) continue; if ( !( em->id()>6||em->id() == emm->id() ||emm->id()>6 ) )continue; Energy pt = IILTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() ); if ( abs( pt-winnerScale )< 10_MeV ) { foundwinnerpt = true; } ptx = min( ptx , pt ); } } } if( !foundwinnerpt ){ generator()->logWarning( Exception() << "Merger: Could not find winner with pt." << "Run with -d3 to get phase space points. " << Exception::warning ); if( Debug::level > 2 ) { generator()->log() << "\nWarning: could not find winner with pt: " << winnerScale/GeV; for( auto const & emm : incoming ){ generator()->log() << "\n" << emm->id() << " " << emm->momentum()/GeV << " " << emm->momentum().m()/GeV << flush; } for( auto const & emm : outgoing ){ generator()->log() <<"\n" << emm->id() << " " << emm->momentum()/GeV << " " << emm->momentum().m()/GeV << flush; } } } return ( ptx>cutscale ) ; } bool Merger::notOnlyMulti() const { return ( treefactory()->onlymulti() != -1&& treefactory()->onlymulti() != int( currentNode()->legsize()-( projected ? 1 : 0 ) ) ); } int Merger::M()const{return theTreeFactory->M();} int Merger::N()const{return theTreeFactory->N();} void Merger::debugVirt(double weight,double w1,double w2,double w3, CrossSection matrixElement,double ww1,double ww2, double ww3, NodePtr node,CrossSection bornWeight)const{ Energy minPT = Constants::MaxEnergy; for( auto const & no :currentNode()->children() )minPT = min( no->pT() , minPT ); generator()->log() << "\nVIRT " << minPT/GeV << " " << weight << " " << w1; generator()->log() << " " << w2; generator()->log() << " " << w3; generator()->log() << " " << ( matrixElement/nanobarn ) << " " << ww1 << " " << ww2 << " " << ww3 << " " << node->pT()/GeV << " " << node->nodeME()->mePartonData()[3]->mass()/GeV << " " << ( bornWeight*SM().alphaS()/( 2.*ThePEG::Constants::pi )/nanobarn ); } void Merger::debugReal(string realstr, double weight, CrossSection me, CrossSection dip)const { Energy minPT = Constants::MaxEnergy; for( auto const & no :currentNode()->children() )minPT = min( no->pT() , minPT ); generator()->log() << "\n"<> theShowerExpansionWeights >> theCMWScheme >> projected >> isUnitarized >> isNLOUnitarized >> theChooseHistory >> theN0 >> theOnlyN >> - weight >> - theGamma >> theSmearing >> iunit( theIRSafePT , GeV ) >> + theSmearing >> iunit( theIRSafePT , GeV ) >> iunit( theMergePt , GeV ) >> iunit( theCentralMergePt , GeV ) >> theLargeNBasis >> FFLTK >> FILTK >> IFLTK >> IILTK >> FFMTK >> FIMTK >> IFMTK >> theDipoleShowerHandler >> theTreeFactory >> - theFirstNodeMap >> theRealSubtractionRatio >> emitDipoleMEDiff; + theFirstNodeMap >> theRealSubtractionRatio ; } // *** 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 ). #include "ThePEG/Utilities/DescribeClass.h" DescribeClass describeHerwigMerger( "Herwig::Merger" , "HwDipoleShower.so" ); //ClassDescription Merger::initMerger; // Definition of the static class description member. #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" void Merger::Init() { static ClassDocumentation documentation ( "The Merger class takes care of merging multiple LO & NLO cross sections." ); ////////////////////////////////////////////////// static Switch interfaceShowerExpansionWeights ( "ShowerExpansionWeights" , "Calculate the expansions of the shower history to be NLO accurate, in different schemes." , &Merger::theShowerExpansionWeights , 3 , false , false ); static SwitchOption interfaceShowerExpansionWeightsScheme0 ( interfaceShowerExpansionWeights , "NoWeights" , "Switch off the expansion." , 0 ); static SwitchOption interfaceShowerExpansionWeightsScheme1 ( interfaceShowerExpansionWeights , "FlatAndHistoryReweighted" , "Sum alphaS expansion and weight with same Historyweight as LO." , 1 ); static SwitchOption interfaceShowerExpansionWeightsScheme2 ( interfaceShowerExpansionWeights , "NoAlphaSReweightForSudakovExpansion" , "Switch off the expansion." , 2 ); static SwitchOption interfaceShowerExpansionWeightsScheme3 ( interfaceShowerExpansionWeights , "NoAlphaSReweightForAllExpansions" , "Switch off the expansion." , 3 ); static SwitchOption interfaceShowerExpansionWeightsScheme4 ( interfaceShowerExpansionWeights , "NoAlphaSReweightForAlphaSExpansion" , "Switch off the expansion." , 4 ); ///////////////////////////////////////////////////////////////////// - static Switch - interfacetheCMWScheme - ( "CMWScheme" , - "Use CMW-Scheme to calculate the alpha_s for the shower expressions." , - &Merger::theCMWScheme , - 0 , - false , false ); - static SwitchOption interfacetheCMWSchemeNo - (interfacetheCMWScheme, - "No", - "No CMW-Scheme", - 0); - static SwitchOption interfacetheCMWSchemeLinear - (interfacetheCMWScheme, - "Linear", - "Linear CMW multiplication: alpha_s(q) -> alpha_s(q)(1+K_g*alpha_s(q)/2pi )", - 1); - static SwitchOption interfacetheCMWSchemeFactor - (interfacetheCMWScheme, - "Factor", - "Use factor in alpha_s argument: alpha_s(q) -> alpha_s(fac*q) with fac=exp(-(67-3pi^2-10/3*Nf)/(33-2Nf)) ", - 2); static Parameter interfaceMergerScale ( "MergingScale" , "The MergingScale." , &Merger::theCentralMergePt , GeV , 20.0*GeV , 0.0*GeV , 0*GeV , false , false , Interface::lowerlim ); - - static Parameter interfacegamma - ( "gamma" , - "gamma parameter." , - &Merger::theGamma , 1.0 , 0.0 , 0 , - false , false , Interface::lowerlim ); - interfacegamma.rank(-1); - - static Switch - interfaceemitDipoleMEDiff - ( "emitDipoleMEDiff" , - "Allow emissions of the unitarisation contribution with prob. 1-min(B_n/sum Dip_n,sum Dip_n/B_n)" , - &Merger::emitDipoleMEDiff , false , false , false ); - static SwitchOption interfaceemitDipoleMEDiffYes - ( interfaceemitDipoleMEDiff , "Yes" , "" , true ); - static SwitchOption interfaceemitDipoleMEDiffNo - ( interfaceemitDipoleMEDiff , "No" , "" , false ); - interfaceemitDipoleMEDiff.rank(-1); static Parameter 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 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 interfacechooseHistory( "chooseHistory" , "Various ways of choosing the history weights: 0(default):dipole xs, 1:dipole/born, 2:flat, 3: 1/pt_dip" , &Merger::theChooseHistory , 3 , -1 , 0 , false , false , Interface::lowerlim ); static Parameter interfacetheRealSubtractionRatio( "RealSubtractionRatio" , "If the pt of the Dipole divided by the merging scale is lower than the ratio, the dipole is used to subtract the real emission. Otherwise the shower approximation is used, " , &Merger::theRealSubtractionRatio , 0. , 0. , 1. , 10.0 , true , false , Interface::limited ); } diff --git a/Shower/Dipole/Merging/Merger.h b/Shower/Dipole/Merging/Merger.h --- a/Shower/Dipole/Merging/Merger.h +++ b/Shower/Dipole/Merging/Merger.h @@ -1,408 +1,404 @@ /// -*- C++ -*- // /// Merger.h is a part of Herwig - A multi-purpose Monte Carlo event generator /// Copyright (C) 2002-2017 The Herwig Collaboration // /// Herwig is licenced under version 3 of the GPL, see COPYING for details. /// Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_Merger_H #define HERWIG_Merger_H // /// This is the declaration of the Merger class. // #include "MergingFactory.fh" #include "Node.fh" #include "ThePEG/Handlers/HandlerBase.h" #include "Herwig/Shower/Dipole/DipoleShowerHandler.h" //#include "Herwig/Shower/Dipole/Base/DipoleSplittingGenerator.h" #include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.h" #include "ThePEG/Cuts/JetFinder.h" #include "ThePEG/Cuts/Cuts.h" namespace Herwig { using namespace ThePEG; class Merger; ThePEG_DECLARE_POINTERS(Merger , MergerPtr ); typedef vector NodePtrVec; //definition of a history step struct HistoryStep { /// containing the full information NodePtr node; /// current sudakov weight of the history double weight; /// current scale of the history Energy scale; }; typedef vector< HistoryStep > History; typedef multimap::ptr> GeneratorMap2; /** * \ingroup DipoleShower * \author Johannes Bellm * - * \brief Merger handles the Merger ....... //TODO . + * \brief This class is responsible for the handling of the merging process + * after the setup stage, performed by the merging factory. + * The class is inherited from the MergerBase class that is visible + * (and the interface) to the shower and matrix elements. + * + * The main functions of the Merger are: + * - matrixElementRegion: + * Given vectors of incoming and outgoing particles, the + * function defines the matrix element region. + * This function is used by the shower after an emission + * that produces a multiplicity still handeld by the merging + * and the Matrix element to determine the region to calculate the + * full ME contribution. + * - mergingScale: + * This is the possibly smeared version of the merging pt + * (centralMergePt) given as input to the Merger. This scaleis used + * in the definitions of ME regions. + * - N and M: + * are the number of additional legs of the highest LO or NLO (virtual) + * contribution. + * - MergingDSigDR: + * The main function to call the calculation of the underlying + * process and the history reweights. + * Here it is also decided if the process should be + * clustered/projected to a process with less legs. + * * * @see \ref MergerInterfaces "The interfaces" * defined for Merger. */ class Merger: public MergerBase { friend class MergingFactory; friend class Node; public: // define the ME region for a particle vector. bool matrixElementRegion(PVector incoming, PVector outgoing, Energy winnerScale = ZERO, Energy cutscale = ZERO)const; /// return the current merging scale, /// gets smeared around the central merging scale in generate kinematics. Energy mergingScale()const{return theMergePt;} /// return the current merging pt (should be unified with mergingScale) Energy mergePt()const {return theMergePt;} /// legsize of highest process with NLO corrections int M()const; /// legsize of the highest LO merged process int N()const; /// legsize of the production process int N0()const{return theN0;} /// cross section of as given by the merging CrossSection MergingDSigDR(); /// ***** virtual functions of the base class ****/// /// set the current xcomb, called from ME void setXComb( tStdXCombPtr ); /// set kinematics, called from ME void setKinematics(); /// clear kinematics, called from ME void clearKinematics(); /// generate kinematics, called from ME bool generateKinematics( const double * ); /// generate kinematics, called from ME void flushCaches(); /// return the current maximum legs, the shower should veto size_t maxLegs() const {return theCurrentMaxLegs;} /// set the current ME void setME(MatchboxMEBasePtr me){ theCurrentME=me; assert(theFirstNodeMap.count(theCurrentME)); theCurrentNode=theFirstNodeMap[theCurrentME]; } - /// allow emissions with a given probability.in the ME region - double emissionProbability() const{ return theEmissionProbability; } - /// set a probability of allowed emission into ME region. - void setEmissionProbability(double x){theEmissionProbability=x;} - protected: /// the merging factory needs to set the legsize of the production process void N0(int n){ theN0=n;} /// return the large-N basis (TODO: implement check if born ME works with the choice) Ptr::ptr largeNBasis()const{return theLargeNBasis;} /// smear the merging pt void smearMergePt(){ const double factor = 1. + (-1. + 2.*UseRandom::rnd() ) * smear(); theMergePt = factor * centralMergePt(); } /// true if the phase space for initial emissions should not be restricted in z. int openZBoundaries()const{return DSH()->showerPhaseSpaceOption();} /// return the current ME MatchboxMEBasePtr currentME() const { return theCurrentME; } /// return the current Node NodePtr currentNode() const { return theCurrentNode; } - /// the gamma parameter to subtract dipoles above a alpha parameter - /// and subtract the corresponding IPK operator - double gamma()const{return theGamma;} private: /// calculate a single sudakov step for a given dipole double singlesudakov(Dipole, Energy, Energy, pair); /// calculate the sudakov supression for a clusternode between /// the current running scale and next scale bool dosudakov(NodePtr Born, Energy running, Energy next, double& sudakov0_n); /// cleanup void cleanup(NodePtr); /// return true if the cluster node has the matching number of /// legs to the current projector stage bool isProjectorStage( NodePtr , int )const; /** * Calculate the staring scale: * if Node is part of the production process, calculate according to the * scale choice object in the merging scale objekt, else * return max(scale as scalechoice , min(Mass(i, j))) */ Energy CKKW_StartScale(NodePtr) const; /// prepare the sudakov calculation void CKKW_PrepareSudakov(NodePtr, Energy); /// number of active flavours as given by the shower double Nf(Energy scale)const{return DSH()->Nf(scale);} /// pointer to the factory MergingFactoryPtr treefactory() const; /// map from ME to first clusternode map firstNodeMap() const ; /// set the current merging pt, smeared in generate kinematics void mergePt(Energy x) {theMergePt = x;} /// return the central merging pt Energy centralMergePt() const {return theCentralMergePt;} private: /// calculate the history weighted born cross section CrossSection MergingDSigDRBornStandard(); - /** - * calculate the history weighted born cross section - * add the difference of IPK with and without alpha parameter - * subtract the dipoles above the alpha parameter - */ - CrossSection MergingDSigDRBornGamma(); /// calculate the history weighted virtual contribution CrossSection MergingDSigDRVirtualStandard(); /** * calculate the history weighted real contribution * splitted into 3 differnt contibutions */ CrossSection MergingDSigDRRealStandard(); /// calculate the history weighted real contribution /// all dipoles above: /// N*(R rnd(i)-Dip_i) history_i U(\phi^n_i) CrossSection MergingDSigDRRealAllAbove(); /// calculate the history weighted real contribution /// not all dipoles above: /// (R - sum PS_i) history_rnd U(\phi^n+1) CrossSection MergingDSigDRRealBelowSubReal(); /// calculate the history weighted real contribution /// not all dipoles above: /// rnd(i)-> N*(PS_i - Dip_i) history_i U(\phi^n_i) CrossSection MergingDSigDRRealBelowSubInt(); /// max legssize the shower should veto for LO size_t maxLegsLO() const {return N0()+N();} /// Calculate the LO partonic cross section. - /// if diffalpha != 1, add the difference of IPK(1)-IPK(diffalpha) - CrossSection TreedSigDR(Energy startscale, double diffalpha=1.); + CrossSection TreedSigDR(Energy startscale); /// fill the projecting xcomb Energy fillProjector(int); /// fill the history, including calculation of sudakov supression void fillHistory(Energy, NodePtr, NodePtr ); /// calculate the pdf ratio for the given clusternode double pdfratio(NodePtr, Energy, Energy, int, bool fromIsME, bool toIsME); /// return the pdf-ratio reweight for the history double pdfReweight(); /// return the alpha_s reweight for the history double alphaReweight(bool nocmw=false); /// max legssize the shower should veto for NLO size_t maxLegsNLO()const {return N0()+M();} /// calculate the virtual contribution. CrossSection LoopdSigDR(Energy startscale ); /// calculate alpha_s expansion of the pdf-ratios double sumPdfReweightExpansion()const; /// calculate alpha_s expansion of the alpha_s-ratios, including K_g double sumAlphaSReweightExpansion()const; /// calculate alpha_s expansion of the sudakov exponents double sumFillHistoryExpansion(); /// calculate alpha_s expansion of the single step alpha_s-ratio, including K_g double alphasExpansion( Energy next, Energy fixedScale)const; /// calculate alpha_s expansion of the single step pdf-ratio double pdfExpansion(NodePtr, int, Energy, Energy, double, int, Energy)const; /// calculate alpha_s expansion of the single step sudakov exponent bool doHistExpansion(NodePtr Born, Energy running, Energy next, Energy fixedScale, double& HistExpansion); /// calculate alpha_s expansion of the single dipole sudakov exponent double singleHistExpansion(Dipole, Energy, Energy, Energy, pair); //alpha_s as given in the shower double as(Energy q)const{return DSH()->as(q);} // set the pointer to the Mergingfactory. void setFactory(MergingFactoryPtr f){theTreeFactory=f;} // set the pointer to the DipoleShower. void setDipoleShower(DipoleShowerHandlerPtr dsh){theDipoleShowerHandler=dsh;} //return the dipole shower handler DipoleShowerHandlerPtr DSH(){return theDipoleShowerHandler;} //return the const dipole shower handler cDipoleShowerHandlerPtr DSH()const{return theDipoleShowerHandler;} /// insert map from ME to first clusternode void firstNodeMap(MatchboxMEBasePtr, NodePtr); /// history choice: weighted history choice int chooseHistory()const {return theChooseHistory;} /// the smearing factor for the merging scale double smear()const{return theSmearing;} /// return the large-N colour basis void largeNBasis(Ptr::ptr x){theLargeNBasis=x;} /// helper function to check the only multi condition. bool notOnlyMulti()const; /// Calculate the CMW AlphaS double cmwAlphaS(Energy q)const; /// debug output for virtual void debugVirt(double, double, double, double, CrossSection, double, double, double, NodePtr,CrossSection) const; /// debug output for reals void debugReal( string, double, CrossSection, CrossSection) const; - private: /// calculate the history expansion unsigned int theShowerExpansionWeights = 2; /// use CMW scheme unsigned int theCMWScheme = 0; /// true if current point should be projected bool projected = true; /// true if LO cross sections should be unitarised bool isUnitarized = true; /// true if NLO contributions should be unitarised bool isNLOUnitarized = true; /// history weight choice int theChooseHistory = 0; /// legsize of production process int theN0 = 0; /// calculate only the N particle contribution int theOnlyN = -1; /// the current maxlegs (either LO or NLO maxlegs) int theCurrentMaxLegs = -1; - /// current weight and weight of clustered born - double weight = 1.0; - double weightCB = 1.0; - /// subtract the dipole contribution above a given gamma - double theGamma = 1.0; /// smearing factor for merging scale double theSmearing = 0.; - - /** - * Allow emissions for the unitarising LO contributions - * according to the prob 1-min(B_n/sum Dip_n,Dip_n/B_n) - */ - bool emitDipoleMEDiff = false; - /// The conditional emission probability if emitDipoleMEDiff is true. - double theEmissionProbability = 0.; /// cutoff for real emission contribution Energy theIRSafePT = 1_GeV; /// current merging scale Energy theMergePt = 4_GeV; /// central merging scale Energy theCentralMergePt = 4_GeV; /// below mergingscale/theRealSubtractionRatio the dipoles are used to subtract. /// above the shower approximation is in use. double theRealSubtractionRatio=3.; /// current cluster histoy including sudakov weights History history; /// pointer to the large-N basis Ptr::ptr theLargeNBasis; /// current Node NodePtr theCurrentNode; /// current ME MatchboxMEBasePtr theCurrentME; /// Tilde kinematics pointers, only to use lastPt(emitter, emission, spectator) Ptr::ptr FFLTK = new_ptr( FFLightTildeKinematics() ); Ptr::ptr FILTK = new_ptr( FILightTildeKinematics() ); Ptr::ptr IFLTK = new_ptr( IFLightTildeKinematics() ); Ptr::ptr IILTK = new_ptr( IILightTildeKinematics() ); Ptr::ptr FFMTK = new_ptr( FFMassiveTildeKinematics() ); Ptr::ptr FIMTK = new_ptr( FIMassiveTildeKinematics() ); Ptr::ptr IFMTK = new_ptr( IFMassiveTildeKinematics() ); //pointer to the shower handler DipoleShowerHandlerPtr theDipoleShowerHandler; /// pointer to the MergingFactory MergingFactoryPtr theTreeFactory; /// map from ME to first Node map theFirstNodeMap; /// map from ME to highest ME weight so far map theHighMeWeightMap; protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ Merger & operator=(const Merger &) = delete; }; } #endif /* HERWIG_Merger_H */ diff --git a/Shower/Dipole/Merging/MergingFactory.cc b/Shower/Dipole/Merging/MergingFactory.cc --- a/Shower/Dipole/Merging/MergingFactory.cc +++ b/Shower/Dipole/Merging/MergingFactory.cc @@ -1,683 +1,680 @@ // -*- C++ -*- // // MergeboxFactory.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the MergeboxFactory class. // #include "MergingFactory.h" #include "Node.h" #include "ThePEG/Repository/Repository.h" #include "ThePEG/Utilities/ColourOutput.h" using namespace Herwig; using std::ostream_iterator; IBPtr MergingFactory::clone() const { return new_ptr(*this); } IBPtr MergingFactory::fullclone() const { return new_ptr(*this); } void MergingFactory::doinit(){ MatchboxFactory::doinit(); if (subProcessGroups()) { throw InitException() << "There are no subprocess groups in merging!"; } } void MergingFactory::productionMode() { if(M()<0) for ( vector::ptr>::iterator amp = amplitudes().begin(); amp != amplitudes().end(); ++amp ) { Repository::clog() << "One-loop contributions from '" << (**amp).name() << "' are not required and will be disabled.\n" << flush; (**amp).disableOneLoop(); } MatchboxFactory::productionMode(); } void MergingFactory::fillMEsMap() { olpProcesses().clear(); assert( getProcesses().size() == 1 ); processMap[0] = getProcesses()[0]; if ( MH()->M() >= 0 ) setHighestVirt(processMap[0].size()+MH()->M()); MH()->N0(processMap[0].size()); for ( int i = 1 ; i <= MH()->N() ; ++i ) { processMap[i] = processMap[i - 1]; processMap[i].push_back("j"); } for ( int i = 0 ; i <= MH()->N() ; ++i ) { const bool below_maxNLO = i < MH()->M() + 1; vector ames = makeMEs(processMap[i], orderInAlphaS() + i, below_maxNLO ); copy(ames.begin(), ames.end(), back_inserter(pureMEsMap()[i])); } } #include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h" void MergingFactory::prepare_BV(int i) { // check if we have virtual contributions bool haveVirtuals = true; for ( auto born : pureMEsMap()[i]) { prepareME(born); if ( born->isOLPTree() ) { int id = orderOLPProcess(born->subProcess(), born->matchboxAmplitude(), ProcessType::treeME2); born->olpProcess(ProcessType::treeME2,id); id = orderOLPProcess(born->subProcess(), born->matchboxAmplitude(), ProcessType::colourCorrelatedME2); born->olpProcess(ProcessType::colourCorrelatedME2,id); bool haveGluon = false; for ( const auto & p : born->subProcess().legs ) if ( p->id() == 21 ) { haveGluon = true; break; } if ( haveGluon ) { id = orderOLPProcess(born->subProcess(), born->matchboxAmplitude(), ProcessType::spinColourCorrelatedME2); born->olpProcess(ProcessType::spinColourCorrelatedME2,id); } } if ( born->isOLPLoop() && i <= MH()->M() ) { int id = orderOLPProcess(born->subProcess(), born->matchboxAmplitude(), ProcessType::oneLoopInterference); born->olpProcess(ProcessType::oneLoopInterference,id); if ( !born->onlyOneLoop() && born->needsOLPCorrelators() ) { id = orderOLPProcess(born->subProcess(), born->matchboxAmplitude(), ProcessType::colourCorrelatedME2); born->olpProcess(ProcessType::colourCorrelatedME2,id); } } haveVirtuals &= born->haveOneLoop(); } // check for consistent conventions on virtuals, if we are to include MH()->M() if (!(i > MH()->M()||haveVirtuals)) throw InitException() << MH()->M() << " NLO corrections requested,\n" << "but no virtual contributions are found."; } void MergingFactory::prepare_R(int i) { for ( auto real : pureMEsMap()[i]) prepareME(real); } #include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h" void MergingFactory::getVirtuals(MatchboxMEBasePtr nlo, bool clone){ const auto & partons = nlo->diagrams().front()->partons(); for ( auto I : DipoleRepository::insertionIOperators(dipoleSet()) ) if ( I->apply(partons) ){ auto myI = I; if ( clone ) myI = I->cloneMe(); nlo->virtuals().push_back(myI); } for ( auto PK : DipoleRepository::insertionPKOperators(dipoleSet()) ) if ( PK->apply(partons) ){ auto myPK = PK; if ( clone ) myPK = PK->cloneMe(); nlo->virtuals().push_back(myPK); } } void MergingFactory::pushB(MatchboxMEBasePtr born, int i) { MatchboxMEBasePtr bornme = born->cloneMe(); bornme->maxMultCKKW(1); bornme->minMultCKKW(0); string pname = fullName() + "/" + bornme->name() + ".Born"; if ( !(generator()->preinitRegister(bornme, pname)) ) throw InitException() << "Born ME "<< pname << " already existing."; - if (MH()->gamma()!=1.) - getVirtuals(bornme,false); - NodePtr clusternode = new_ptr(Node(bornme, 0, MH())); clusternode->deepHead(clusternode); MH()->firstNodeMap(bornme,clusternode); bornme->merger(MH()); vector current = {{clusternode}}; vector children; unsigned int k = 1; while ( ! thePureMEsMap[i - k].empty() ) { for ( auto tmp : current ){//j tmp->birth(thePureMEsMap[i - k]); for ( auto tmpchild : tmp->children() ) {//m children.push_back(tmpchild); } } current = children; children.clear(); ++k; } if ( MH()->N() > i ) bornme->needsCorrelations(); else bornme->needsNoCorrelations(); bornme->cloneDependencies(); MEs().push_back(bornme); } void MergingFactory::pushV(MatchboxMEBasePtr born, int i) { MatchboxMEBasePtr nlo = born->cloneMe(); string pname = fullName() + "/" + nlo->name() + ".Virtual"; if ( !(generator()->preinitRegister(nlo, pname)) ) throw InitException() << "Virtual ME "<< pname << " already existing."; ////////////////////////////////////NLO/////////////////////////// nlo->virtuals().clear(); getVirtuals(nlo , false); if ( nlo->virtuals().empty() ) throw InitException() << "No insertion operators have been found for " << born->name() << ".\n"; nlo->doOneLoopNoBorn(); ////////////////////////////////////NLO/////////////////////////// NodePtr clusternode = new_ptr(Node(nlo, 0,MH())); clusternode->deepHead(clusternode); clusternode->virtualContribution(true); MH()->firstNodeMap(nlo,clusternode); nlo->merger(MH()); vector current = {{clusternode}}; vector children; unsigned int k = 1; while ( ! thePureMEsMap[i - k].empty() ) { for ( auto tmp : current ){ tmp->birth(thePureMEsMap[i - k]); for ( auto tmpchild : tmp->children()) children.push_back(tmpchild); } current = children; children.clear(); ++k; } if ( nlo->isOLPLoop() ) { int id = orderOLPProcess(nlo->subProcess(), born->matchboxAmplitude(), ProcessType::oneLoopInterference); nlo->olpProcess(ProcessType::oneLoopInterference,id); if ( !nlo->onlyOneLoop() && nlo->needsOLPCorrelators() ) { id = orderOLPProcess(nlo->subProcess(), born->matchboxAmplitude(), ProcessType::colourCorrelatedME2); nlo->olpProcess(ProcessType::colourCorrelatedME2,id); } } nlo->needsCorrelations(); nlo->cloneDependencies(); MEs().push_back(nlo); } void MergingFactory::pushR(MatchboxMEBasePtr born, int i) { MatchboxMEBasePtr bornme = born->cloneMe(); string pname = fullName() + "/" + bornme->name() + ".Real"; if ( !(generator()->preinitRegister(bornme, pname)) ) throw InitException() << "Subtracted ME " << pname << " already existing."; NodePtr clusternode = new_ptr(Node(bornme, 1, MH())); clusternode->deepHead(clusternode); clusternode->subtractedReal(true); MH()->firstNodeMap(bornme,clusternode); bornme->merger(MH()); vector current = {{clusternode}}; vector children; unsigned int k = 1; while ( ! thePureMEsMap[i - k].empty() ) { for ( auto tmp : current ){ tmp->birth(thePureMEsMap[i - k]); for ( auto tmpchild : tmp->children()) children.push_back(tmpchild); } current = children; children.clear(); ++k; } if(clusternode->children().empty()){ // This is a finite real contribution. // This process is included in the LO merging. return; } if ( MH()->N() > i ) bornme->needsCorrelations(); else bornme->needsNoCorrelations(); bornme->cloneDependencies(pname); MEs().push_back(bornme); } // MergingFactory should never order OLPs here, // they're done elsewhere. void MergingFactory::orderOLPs() {} #include "ThePEG/Utilities/StringUtils.h" vector MergingFactory::parseProcess(string in) { vector process = StringUtils::split(in); if ( process.size() < 3 ) throw Exception() << "MatchboxFactory: Invalid process."<< Exception::runerror; for ( string & p : process) { p = StringUtils::stripws(p); } theN = 0; bool prodprocess = true; vector result; for ( const string & p : process ) { if ( p == "->" ) continue; if (p=="[") { prodprocess = false; } else if (p=="]") { prodprocess = false; // TODO what if there's stuff after the bracket? assert( p == process.back() ); break; } else if (p=="[j") { prodprocess = false; ++theN; } else if (p=="j" && !prodprocess) { ++theN; prodprocess = false; } else if (p=="j]") { ++theN; prodprocess = false; // TODO what if there's stuff after the bracket? assert( p == process.back() ); break; } else if ( prodprocess ) { result.push_back(p); } else { throw InitException() << "Unknown particle class \"" << p << '"' << " in the process definition merging bracket.\n" << "Only jets (\"j\") are supported at the moment."; } } return result; } #include "Herwig/Utilities/Progress.h" void MergingFactory::setup() { useMe(); DipoleShowerHandlerPtr dsh=dynamic_ptr_cast(this->CKKWHandler()); if(! dsh )throw InitException() << "The showerhandler for the MergingFactory must be the DipoleShower. "; dsh->setMerger(MH()); MH()->setFactory(this); MH()->setDipoleShower(dsh); if(!ransetup){ generator()->log() <<"\nStarting merging setup.\n\n"; olpProcesses().clear(); externalAmplitudes().clear(); // We set the couplings in the ME to be fixed // and reweight in the history weight for this. setFixedCouplings(true); setFixedQEDCouplings(true); // rebind the particle data objects, can't use rebind() function for ( auto & g : particleGroups()) { for ( auto & p : g.second) { p = getParticleData(p->id()); } } const PDVector& partons = particleGroups()["j"]; unsigned int nl = 0; for ( const auto p : partons ) { const Energy mass = p->hardProcessMass(); const long pid = p->id(); if ( abs(pid) < 7 && mass == ZERO ) ++nl; if ( pid > 0 && pid < 7 && mass == ZERO ) nLightJetVec( pid ); if ( pid > 0 && pid < 7 && mass != ZERO ) nHeavyJetVec( pid ); } nLight(nl/2); const PDVector& partonsInP = particleGroups()["p"]; for ( const auto pip : partonsInP ) if ( pip->id() > 0 && pip->id() < 7 && pip->hardProcessMass() == ZERO ) nLightProtonVec( pip->id() ); // fill the amplitudes if ( !amplitudes().empty() ) fillMEsMap(); // Use the colour basis of the first element of amplitudes // to set the large N colour basis for the MergingHelper assert(!amplitudes().empty() ); if ( !amplitudes()[0]->colourBasis() ) throw Exception() << "MergingFactory::setup(): Expecting a colour basis object." << Exception::runerror; auto largeNBasis = amplitudes()[0]->colourBasis()->cloneMe(); largeNBasis->clear(); largeNBasis->doLargeN(); MH()->largeNBasis(largeNBasis); // prepare the Born and virtual matrix elements for ( int i = 0 ; i <= max(0, MH()->N()) ; ++i ) prepare_BV(i); // prepare the real emission matrix elements for ( int i = 0 ; i <= MH()->N() ; ++i ) prepare_R(i); if (MH()->N()<=MH()->M()) { throw InitException() << "Merging: The number of NLOs need to be" << "\nsmaller than the number of LO processes.\n"; } // Order the external Amplitudes. orderOLPs(); // start creating matrix elements MEs().clear(); // count the subprocesses size_t numb = 0; size_t numv = 0; size_t numr = 0; for (int i = 0; i <= max(0, MH()->N()) ; ++i ) { for ( auto born : thePureMEsMap[i] ) { if (bornContributions() ) { numb++; } } } for (int i = 0 ; i <=max(0, MH()->N()); ++i ) for ( auto virt : thePureMEsMap[i] ) if ( virtualContributions() && i <= MH()->M()) { numv++; } for (int i = 1; i <= max(0, MH()->N()) ; ++i ) for ( auto real : thePureMEsMap[i] ) if (realContributions() && i <= MH()->M() + 1 ){ numr++; } if(int(numb+numv+numr) < theChunk){ throw InitException() << "You try to chunk (Chunk="< 0 && theChunkPart == 0 ) throw InitException() <<" Set the ChunkPart ( = "<log() << ANSI::yellow << "\n\nWarning: \nYou split up the runs into theChunks. This is no standard feature." << "\nYou are now responsible to make sure to run all theChunk parts." << "\nThis setup run is for: " << theChunkPart << "/" << theChunk<<"\n\n "; generator()->log() << ANSI::red << "Preparing Merging: "; generator()->log() << ANSI::green << numb << " x Born " << ANSI::red; if (MH()->M()>-1) { generator()->log() << ANSI::yellow << numv << " x Virtual "; generator()->log() << ANSI::blue << numr << " x Real " << ANSI::red << flush; } int countchunk=0; progress_display progressBar{ numb+numv+numr, generator()->log() }; // insert the born contributions to ME vector for (int i = 0; i <= max(0, MH()->N()) ; ++i ) for ( auto born : thePureMEsMap[i]) if (bornContributions() ){ countchunk++; // theChunkPart is in [1,...,theChunk] if ( theChunk == 0 || theChunkPart == countchunk ) pushB( born , i ); if ( countchunk == theChunk ) countchunk=0; generator()->log() << ANSI::green; ++progressBar; generator()->log() << ANSI::reset; } // insert the virtual contributions to ME vector for (int i = 0 ; i <=max(0, MH()->N()); ++i ) for ( auto virt : thePureMEsMap[i]) if ( virtualContributions() && i <= MH()->M()){ countchunk++;// theChunkPart is in [1,...,theChunk] if ( theChunk == 0 || theChunkPart == countchunk ) pushV(virt, i); if ( countchunk == theChunk ) countchunk=0; generator()->log() << ANSI::yellow; ++progressBar; } // insert the real contributions to ME vector for (int i = 1; i <= max(0, MH()->N()) ; ++i ) for ( auto real : thePureMEsMap[i] ) if (realContributions()&& i <= MH()->M() + 1 ){ countchunk++;// theChunkPart is in [1,...,theChunk] if ( theChunk == 0 || theChunkPart == countchunk ) pushR(real, i); if ( countchunk == theChunk ) countchunk=0; generator()->log() << ANSI::blue; ++progressBar; } generator()->log() << ANSI::reset; if ( !externalAmplitudes().empty() ) { generator()->log() << "Initializing external amplitudes." << endl; progress_display progressBar{ externalAmplitudes().size(), generator()->log() }; for ( const auto ext : externalAmplitudes() ) { if ( ! ext->initializeExternal() ) { throw InitException() << "error: failed to initialize amplitude '" << ext->name() << "'\n"; } ++progressBar; } generator()->log() << "---------------------------------------------------" << endl; } if ( !olpProcesses().empty() ) { generator()->log() << "Initializing one-loop provider(s)." << endl; map::tptr, map, int> > olps; for (const auto oit : olpProcesses()) { olps[oit.first] = oit.second; } progress_display progressBar{ olps.size(), generator()->log() }; for ( const auto olpit : olps ) { if ( !olpit.first->startOLP(olpit.second) ) { throw InitException() << "error: failed to start OLP for amplitude '" << olpit.first->name() << "'\n"; } ++progressBar; } generator()->log() << "---------------------------------------------------\n" << flush; } generator()->log() <<"\nGenerated "<log() << "---------------------------------------------------\n" << flush; generator()->log() <<"\n\n" << ANSI::red <<"Note: Due to the unitarization of the higher " <<"\nmultiplicities, the individual cross sections " <<"\ngiven in the integration and run step are not" <<"\nmeaningful without merging." << ANSI::reset << endl; ransetup=true; } } #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" void MergingFactory::persistentOutput(PersistentOStream & os) const { os << theonlymulti << ransetup << processMap << theMergingHelper <> theonlymulti >> ransetup >> processMap >> theMergingHelper >>theM>>theN>>theNonQCDCuts>>theChunk>>theChunkPart; } #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Command.h" void MergingFactory::Init() { static Parameter interfaceonlymulti("onlymulti", "Calculate only the ME with k additional partons.", &MergingFactory::theonlymulti, -1, -1, 0, false, false, Interface::lowerlim); static Switch interface_Unitarized("Unitarized", "Unitarize the cross section (default is unitarised. NLO merging must be unitarised).", &MergingFactory::unitarized, true, false, false); static SwitchOption interface_UnitarizedYes(interface_Unitarized, "Yes", "Switch on the unitarized cross section.", true); static SwitchOption interface_UnitarizedNo(interface_Unitarized, "No", "Switch off the unitarized cross section.", false); static Reference interfaceMergingHelper("MergingHelper", "Pointer to the Merging Helper.", &MergingFactory::theMergingHelper, false, false, true, true, false); static Parameter interfaceaddNLOLegs("NLOProcesses", "Set the number of virtual corrections to consider. 0 is default for no virtual correction.", &MergingFactory::theM, 0, 0, 0, false, false, Interface::lowerlim); static Reference interfaceNonQcdCuts("NonQCDCuts", "Cut on non-QCD modified observables. Be carefull!", &MergingFactory::theNonQCDCuts, false, false, true, true, false); static Parameter interfacetheChunk("Chunk", "Cut the number of subprocesses into n theChunks.", &MergingFactory::theChunk, -1, -1, 0, false, false, Interface::lowerlim); static Parameter interfacetheChunkPart("ChunkPart", "If theChunk is larger then 0, set this parameter to the n'th part. Make sure to add the ChunksParts afterwards.", &MergingFactory::theChunkPart, -1, -1, 0, false, false, Interface::lowerlim); } // *** 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 describeHerwigMergingFactory("Herwig::MergingFactory", "HwDipoleShower.so"); diff --git a/Shower/Dipole/Merging/Node.cc b/Shower/Dipole/Merging/Node.cc --- a/Shower/Dipole/Merging/Node.cc +++ b/Shower/Dipole/Merging/Node.cc @@ -1,489 +1,457 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the Node class. // #include "Node.h" #include "MergingFactory.h" #include "Merger.h" using namespace Herwig; Node::Node(MatchboxMEBasePtr nodeME, int cutstage, MergerPtr mh) :Interfaced(), thenodeMEPtr(nodeME), thedipol(), theparent(), theCutStage(cutstage), - isOrdered(true), theSubtractedReal(false), theVirtualContribution(false), theMergingHelper(mh) { nodeME->maxMultCKKW(1); nodeME->minMultCKKW(0); } Node::Node(NodePtr deephead, NodePtr head, SubtractionDipolePtr dipol, MatchboxMEBasePtr nodeME, int cutstage) :Interfaced(), thenodeMEPtr(nodeME), thedipol(dipol), theparent(head), theDeepHead(deephead), theCutStage(cutstage), -isOrdered(true), theSubtractedReal(false), theVirtualContribution(false), theMergingHelper() //The subnodes have no merging helper { } Node::~Node() { } SubtractionDipolePtr Node::dipole() const { return thedipol; } /** returns the matrix element pointer */ const MatchboxMEBasePtr Node::nodeME() const { return thenodeMEPtr; } /** access the matrix element pointer */ MatchboxMEBasePtr Node::nodeME() { return thenodeMEPtr; } pair Node::getInOut( ){ PVector in; const auto me= nodeME(); const auto pd=me->mePartonData(); for( auto i : {0 , 1} ) in.push_back(pd[i]->produceParticle( me->lastMEMomenta()[i] ) ); PVector out; for ( size_t i = 2;i< pd.size();i++ ){ PPtr p = pd[i]->produceParticle( me->lastMEMomenta()[i] ); out.push_back( p ); } return { in , out }; } int Node::legsize() const {return nodeME()->legsize();} NodePtr Node::randomChild() { return thechildren[UseRandom::irnd(thechildren.size())]; } bool Node::allAbove(Energy pt) { for (NodePtr child : thechildren) if ( child->pT() < pt ) return false; return true; } Energy Node::maxChildPt(){ Energy maxi=-1*GeV; for (NodePtr child : thechildren)maxi=max(child->pT(),maxi); return maxi; } bool Node::isInHistoryOf(NodePtr other) { while (other->parent()) { if (other == this) return true; other = other->parent(); } return false; } void Node::flushCaches() { if (didflush) return; didflush=true; for ( auto const & ch: thechildren) { ch->xcomb()->clean(); ch->nodeME()->flushCaches(); ch->flushCaches(); } } void Node::setKinematics() { for (auto const & ch: thechildren) { ch->dipole()->setXComb(ch->xcomb()); ch->dipole()->setKinematics(); ch->nodeME()->setKinematics(); ch->setKinematics(); } } void Node::clearKinematics() { for (auto const & ch: thechildren) { ch->dipole()->setXComb(ch->xcomb()); ch->nodeME()->clearKinematics(); ch->dipole()->clearKinematics(); ch->clearKinematics(); } } bool Node::generateKinematics(const double *r, bool directCut) { didflush=false; // If there are no children to the child process we are done. if(children().empty()) return true; assert(parent()); if ( ! directCut && pT() < deepHead()->MH()->mergePt()) { // Real emission: // If there are children to the child process, // we now require that all subsequent children // with pt < merging scale are in their ME region. // Since the possible children of the real emission // contribution are now in their ME region, // it is clear that a second clustering is possible // -- modulo phase space restrictions. // Therefore the real emission contribution // are unitarised and the cross section is // hardly modified. auto inOutPair = getInOut(); NodePtr rc = randomChild(); rc->dipole()->setXComb(rc->xcomb()); if(!rc->dipole()->generateKinematics(r))assert(false); // If not in ME -> return false if(!deepHead()->MH()->matrixElementRegion( inOutPair.first , inOutPair.second , rc->pT() , deepHead()->MH()->mergePt() ) )return false; } for (auto & ch : children() ) { ch->dipole()->setXComb(ch->xcomb()); if ( !ch->dipole()->generateKinematics(r) ) { assert(false); } ch->generateKinematics( r, true); } return true; } bool Node::firstgenerateKinematics(const double *r, bool directCut) { didflush=false; // This is called form the merging helper for the first node. So: assert(!parent()); assert(xcomb()); if(MH()->treefactory()->nonQCDCuts()){ tcPDVector outdata(xcomb()->mePartonData().begin()+2, xcomb()->mePartonData().end()); vector outmomenta(xcomb()->meMomenta().begin()+2, xcomb()->meMomenta().end()); if ( !MH()->treefactory()->nonQCDCuts()->passCuts(outdata,outmomenta, xcomb()->mePartonData()[0], xcomb()->mePartonData()[1]) ) return false; } ///// This should not be needed!!! ///// ( Warning inMerger::matrixElementRegion gets triggered.) flushCaches(); //Set here the new merge Pt for the next phase space point.( Smearing!!!) MH()->smearMergePt(); // If there are no children to this node, we are done here: if (children().empty()) return true; // directCut is for born and for virtual contributions. // if directCut is true, then cut on the first ME region. // call recursiv generate kinematics for subsequent nodes. if ( directCut ){ auto inOutPair = getInOut(); NodePtr rc = randomChild(); rc->dipole()->setXComb(rc->xcomb()); if ( !rc->dipole()->generateKinematics(r) ) { return false; } rc->nodeME()->setXComb(rc->xcomb()); - if(MH()->gamma() == 1.){ + if(!MH()->matrixElementRegion( inOutPair.first , inOutPair.second , rc->pT() , MH()->mergePt() ) ){ return false; - } - - }else{ - - // Different treatment if gamma is not 1. - // Since the dipoles need to be calculated always - // their alpha region is touched. - // AlphaRegion != MERegion !!! - bool inAlphaPS = false; - for (auto const & ch: thechildren) { - ch->dipole()->setXComb(ch->xcomb()); - if ( !ch->dipole()->generateKinematics(r) ) - return false; - MH()->treefactory()->setAlphaParameter( MH()->gamma() ); - inAlphaPS |= ch->dipole()->aboveAlpha(); - MH()->treefactory()->setAlphaParameter( 1. ); - } - NodePtr rc = randomChild(); - if(!inAlphaPS&& - !MH()->matrixElementRegion( inOutPair.first , - inOutPair.second , - rc->pT() , - MH()->mergePt() ) ) - return false; - - - - } } for (auto const & ch: thechildren) { ch->dipole()->setXComb(ch->xcomb()); if ( !ch->dipole()->generateKinematics(r) ) { cout<<"\nCould not generate dipole kinematics";;return false; } if( ! ch->generateKinematics(r,directCut) )return false; } return true; } StdXCombPtr Node::xcomb() const { assert(thexcomb); return thexcomb; } StdXCombPtr Node::xcomb(){ if(thexcomb)return thexcomb; assert(parent()); thexcomb=dipole()->makeBornXComb(parent()->xcomb()); xcomb()->head(parent()->xcomb()); dipole()->setXComb(thexcomb); return thexcomb; } void Node::setXComb(tStdXCombPtr xc) { assert ( !parent() ); thexcomb=xc; assert(thexcomb->lastParticles().first); } #include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h" void Node::birth(const vector & vec) { // produce the children vector dipoles = nodeME()->getDipoles(DipoleRepository::dipoles( nodeME()->factory()->dipoleSet()), vec, true); for ( auto const & dip : dipoles ) { dip->doSubtraction(); NodePtr node = new_ptr(Node(theDeepHead, this, dip, dip->underlyingBornME(), theDeepHead->cutStage())); thechildren.push_back(node); } } vector Node::getNextOrderedNodes(bool normal, double hardScaleFactor) const { vector temp = children(); vector res; for (NodePtr const & child : children()) { if(deepHead()->MH()->mergePt()>child->pT()) { res.clear(); return res; } } for (NodePtr const & child: children()) { if (parent()&& normal) { if ( child->pT() < pT() ) { continue; } } if ( child->children().size() != 0 ) { for (NodePtr itChild: child->children()) { if( itChild->pT() > child->pT()&&child->inShowerPS(itChild->pT()) ) { res.push_back(child); break; } } } else { const auto sc=child->nodeME()->factory()->scaleChoice(); sc->setXComb(child->xcomb()); if ( sqr(hardScaleFactor)* sc->renormalizationScale() >= sqr(child->pT()) && child->inShowerPS(hardScaleFactor*sqrt(sc->renormalizationScale()))) { res.push_back(child); } } } return res; } bool Node::inShowerPS(Energy hardpT)const { // Here we decide if the current phase space // point can be reached from the underlying Node. // Full phase space available -> Tilde Kinematic is always fine. if(deepHead()->MH()->openZBoundaries()==1) return true; double z_ = dipole()->lastZ(); // restrict according to hard scale if(deepHead()->MH()->openZBoundaries()==0){ pair zbounds = dipole()->tildeKinematics()->zBounds(pT(), hardpT); return (zbounds.first temp = getNextOrderedNodes(normal, hardScaleFactor); Energy minpt = Constants::MaxEnergy; Selector subprosel; while (temp.size() != 0) { minpt = Constants::MaxEnergy; subprosel.clear(); for (NodePtr const & child : temp) { assert(deepHead()->MH()->largeNBasis()); if( child->dipole()->underlyingBornME()->largeNColourCorrelatedME2( {child->dipole()->bornEmitter(), child->dipole()->bornSpectator()}, deepHead()->MH()->largeNBasis()) != 0. ) { double weight = 1.; if ( deepHead()->MH()->chooseHistory() == 0 ) weight = abs(child->dipole()->dSigHatDR()/nanobarn); else if ( deepHead()->MH()->chooseHistory() == 1 ) weight = abs(child->dipole()->dSigHatDR()/child->nodeME()->dSigHatDRB()); else if ( deepHead()->MH()->chooseHistory() == 2 ) weight = 1.; else if ( deepHead()->MH()->chooseHistory() == 3 ) weight = 1_GeV/child->pT(); else assert(false); if(weight != 0.) { subprosel.insert(weight , child); minpt = min(minpt, child->pT()); } } } if (subprosel.empty()) return res; res = subprosel.select(UseRandom::rnd()); temp = res->getNextOrderedNodes(true, hardScaleFactor); } return res; } pair Node::calcDipandPS(Energy scale)const { return dipole()->dipandPs(sqr(scale), deepHead()->MH()->largeNBasis()); } CrossSection Node::calcPs(Energy scale)const { return dipole()->ps(sqr(scale), deepHead()->MH()->largeNBasis()); } CrossSection Node::calcDip(Energy scale)const { return dipole()->dip(sqr(scale)); } IBPtr Node::clone() const { return new_ptr(*this); } IBPtr Node::fullclone() const { return new_ptr(*this); } #include "ThePEG/Persistency/PersistentOStream.h" void Node::persistentOutput(PersistentOStream & os) const { os << thexcomb<< thenodeMEPtr<< thedipol<< thechildren<< theparent<< theProjector<< theDeepHead<< theCutStage<< ounit(theRunningPt, GeV)<< theSubtractedReal<< theVirtualContribution<< theMergingHelper; } #include "ThePEG/Persistency/PersistentIStream.h" void Node::persistentInput(PersistentIStream & is, int) { is >> thexcomb>> thenodeMEPtr>> thedipol>> thechildren>> theparent>> theProjector>> theDeepHead>> theCutStage>> iunit(theRunningPt, GeV)>> theSubtractedReal>> theVirtualContribution>> theMergingHelper; } // *** 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). #include "ThePEG/Utilities/DescribeClass.h" DescribeClass describeHerwigNode("Herwig::Node", "HwDipoleShower.so"); void Node::Init() { static ClassDocumentation documentation("There is no documentation for the Node class"); } diff --git a/Shower/Dipole/Merging/Node.h b/Shower/Dipole/Merging/Node.h --- a/Shower/Dipole/Merging/Node.h +++ b/Shower/Dipole/Merging/Node.h @@ -1,237 +1,242 @@ // -*- C++ -*- #ifndef Herwig_Node_H #define Herwig_Node_H // // This is the declaration of the Node class. // #include "Node.fh" #include "MergingFactory.fh" #include "Merger.h" #include "ThePEG/Config/ThePEG.h" #include "ThePEG/Config/std.h" #include "ThePEG/Interface/Interfaced.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh" #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh" #include "Herwig/Shower/Dipole/Base/DipoleEventRecord.h" #include "ThePEG/MatrixElement/MEBase.h" #include namespace Herwig { using namespace ThePEG; /** - * Here is the documentation of the Node class. + * The Node class represents a part of the shower history + * in the merging procedure. + * Each node can search for its "childen" in the "birth" step. + * The children processes that can be created by clustering two legs + * with a given spectator -- just like in the shower. + * To perform the "birth" a vector of processes is given to the node + * and the possible dipoles are given by the factory object. + * * * @see \ref NodeInterfaces "The interfaces" * defined for Node. */ class Node : public Interfaced { public: /** @name Standard constructors and destructors. */ //@{ Node (){}; // constructor for first nodes Node(MatchboxMEBasePtr nodeME , int cutstage , MergerPtr mh ); // another constructor for underlying nodes Node(NodePtr deephead, NodePtr head, SubtractionDipolePtr dipol, MatchboxMEBasePtr nodeME, int cutstage); /// The destructor. virtual ~Node(); //@} public: // get children from vector void birth(const vector & vec); /// recursive setXComb. proStage is the number of clusterings /// before the projectors get filled. void setXComb(tStdXCombPtr xc); /// calculate the dipole and ps approximation pair calcDipandPS(Energy scale)const; /// calculate the ps approximation CrossSection calcPs(Energy scale)const; /// calculate the dipole CrossSection calcDip(Energy scale)const; /// recursive flush caches and clean up XCombs. void flushCaches(); /// recursive clearKinematics void clearKinematics(); /// recursive setKinematics void setKinematics(); /// recursive generateKinematics using tilde kinematics of the dipoles bool generateKinematics(const double *r, bool directCut); /// generate the kinamatics of the first node bool firstgenerateKinematics(const double *r, bool directCut); //return the ME const MatchboxMEBasePtr nodeME() const; //return the node ME MatchboxMEBasePtr nodeME(); //return the parent Node NodePtr parent() const {return theparent;} /// vector of children nodes created in birth vector< NodePtr > children() const {return thechildren;} //pick a random child (flat) NodePtr randomChild(); /// true if all children show scales above pt bool allAbove(Energy pt); /// return maximum of all child pts. Energy maxChildPt(); /// true if the node is in the history of other. bool isInHistoryOf(NodePtr other); /// legsize of the node ME int legsize() const; /// set the first node (first men). only use in factory void deepHead(NodePtr deephead) {theDeepHead = deephead;} /// return the first node NodePtr deepHead() const {return theDeepHead;} /// returns the dipol of the node. SubtractionDipolePtr dipole() const; /// return the xcomb StdXCombPtr xcomb() const; /// return the xcomb (if not created, create one from head) StdXCombPtr xcomb() ; /// return the current running pt Energy runningPt() const { return theRunningPt; } /// set the current running pt void runningPt(Energy x) { theRunningPt=x; } /// return the cut stage to cut on merging pt in generate kinematics int cutStage() const { return theCutStage; } /// get a vector of the next nodes, ordered in pt (and in parton shower phace space) vector getNextOrderedNodes(bool normal=true, double hardscalefactor=1.) const; //true if the node is in shower history for a given pt bool inShowerPS(Energy hardpt)const; //get the history NodePtr getHistory(bool normal=true, double hardscalefactor=1.); //true if node correspond to a subtracted real. bool subtractedReal() const {return theSubtractedReal;} /// set if node correspont to a subtracted real. void subtractedReal(bool x) { theSubtractedReal = x;} //true if node correspond to a virtual contribution. bool virtualContribution() const { return theVirtualContribution ;} /// set if node correspont to a virtual contribution. void virtualContribution(bool x) {theVirtualContribution = x;} //pointer to the merging helper MergerPtr MH()const{return theMergingHelper;} /// set the merging helper void MH(MergerPtr a){theMergingHelper=a;} /// pT of the dipole Energy pT()const{return dipole()->lastPt();} /// get incoming and outgoing particles (TODO: expensive) pair getInOut(); private: /// the Matrixelement representing this node. MatchboxMEBasePtr thenodeMEPtr; /// the dipol used to substract /// and generate kinematics using tilde kinematics SubtractionDipolePtr thedipol; /// the parent node NodePtr theparent; /// The godfather node of whole tree.(Firstnode) NodePtr theDeepHead; /** * The CutStage is number of clusterings which are possible without * introducing a merging scale to cut away singularities. * -> subtracted MEs have the CutStage 1. * -> virtual and normal tree level ME get 0. */ int theCutStage; - /// tell if node belongs to an ordered history - bool isOrdered; /// flag to tell if node is subtracted real bool theSubtractedReal; /// flag to tell if node is virtual contribution bool theVirtualContribution; - /// the merging helper + /// the merging helper (should be static) MergerPtr theMergingHelper; //the xcomb of the node StdXCombPtr thexcomb; /// vector of the children node vector< NodePtr > thechildren; /// the current running pt Energy theRunningPt; /// The nodes of the projection stage. NodePtr theProjector; /// flag not to enter infinite loop. (There should be a better solution...) bool didflush=false; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ Node & operator=(const Node &) = delete; }; } #endif /* Herwig_Node_H */ diff --git a/src/Merging/FactorCMWScheme.in b/src/Merging/FactorCMWScheme.in --- a/src/Merging/FactorCMWScheme.in +++ b/src/Merging/FactorCMWScheme.in @@ -1,90 +1,97 @@ ####### # CMW scheme ####### -cd /Herwig/Merging -set Merger:CMWScheme Factor - cd /Herwig/DipoleShower/Kernels set FFgx2ggxDipoleKernel:CMWScheme Factor set FFqx2qgxDipoleKernel:CMWScheme Factor set FFgx2ddxDipoleKernel:CMWScheme Factor set FFgx2uuxDipoleKernel:CMWScheme Factor set FFgx2ccxDipoleKernel:CMWScheme Factor set FFgx2ssxDipoleKernel:CMWScheme Factor set FFgx2bbxDipoleKernel:CMWScheme Factor set FFMgx2ggxDipoleKernel:CMWScheme Factor set FFMdx2dgxDipoleKernel:CMWScheme Factor set FFMux2ugxDipoleKernel:CMWScheme Factor set FFMcx2cgxDipoleKernel:CMWScheme Factor set FFMsx2sgxDipoleKernel:CMWScheme Factor set FFMbx2bgxDipoleKernel:CMWScheme Factor set FFMtx2tgxDipoleKernel:CMWScheme Factor set FFMgx2ddxDipoleKernel:CMWScheme Factor set FFMgx2uuxDipoleKernel:CMWScheme Factor set FFMgx2ccxDipoleKernel:CMWScheme Factor set FFMgx2ssxDipoleKernel:CMWScheme Factor set FFMgx2bbxDipoleKernel:CMWScheme Factor set FIgx2ggxDipoleKernel:CMWScheme Factor set FIqx2qgxDipoleKernel:CMWScheme Factor set FIgx2ddxDipoleKernel:CMWScheme Factor set FIgx2uuxDipoleKernel:CMWScheme Factor set FIgx2ccxDipoleKernel:CMWScheme Factor set FIgx2ssxDipoleKernel:CMWScheme Factor set FIgx2bbxDipoleKernel:CMWScheme Factor set FIMdx2dgxDipoleKernel:CMWScheme Factor set FIMux2ugxDipoleKernel:CMWScheme Factor set FIMcx2cgxDipoleKernel:CMWScheme Factor set FIMsx2sgxDipoleKernel:CMWScheme Factor set FIMbx2bgxDipoleKernel:CMWScheme Factor set FIMtx2tgxDipoleKernel:CMWScheme Factor set FIMgx2ddxDipoleKernel:CMWScheme Factor set FIMgx2uuxDipoleKernel:CMWScheme Factor set FIMgx2ccxDipoleKernel:CMWScheme Factor set FIMgx2ssxDipoleKernel:CMWScheme Factor set FIMgx2bbxDipoleKernel:CMWScheme Factor #set FIMgx2ttxDipoleKernel:CMWScheme Factor set IFgx2ggxDipoleKernel:CMWScheme Factor set IFqx2qgxDipoleKernel:CMWScheme Factor set IFqx2gqxDipoleKernel:CMWScheme Factor set IFgx2ddbarxDipoleKernel:CMWScheme Factor set IFgx2dbardxDipoleKernel:CMWScheme Factor set IFgx2uubarxDipoleKernel:CMWScheme Factor set IFgx2ubaruxDipoleKernel:CMWScheme Factor set IFgx2ccbarxDipoleKernel:CMWScheme Factor set IFgx2cbarcxDipoleKernel:CMWScheme Factor set IFgx2ssbarxDipoleKernel:CMWScheme Factor set IFgx2sbarsxDipoleKernel:CMWScheme Factor set IFMgx2ggxDipoleKernel:CMWScheme Factor set IFMqx2qgxDipoleKernel:CMWScheme Factor set IFMqx2gqxDipoleKernel:CMWScheme Factor set IFMgx2ddbarxDipoleKernel:CMWScheme Factor set IFMgx2dbardxDipoleKernel:CMWScheme Factor set IFMgx2uubarxDipoleKernel:CMWScheme Factor set IFMgx2ubaruxDipoleKernel:CMWScheme Factor set IFMgx2ccbarxDipoleKernel:CMWScheme Factor set IFMgx2cbarcxDipoleKernel:CMWScheme Factor set IFMgx2ssbarxDipoleKernel:CMWScheme Factor set IFMgx2sbarsxDipoleKernel:CMWScheme Factor set IIgx2ggxDipoleKernel:CMWScheme Factor set IIqx2qgxDipoleKernel:CMWScheme Factor set IIqx2gqxDipoleKernel:CMWScheme Factor set IIgx2ddbarxDipoleKernel:CMWScheme Factor set IIgx2dbardxDipoleKernel:CMWScheme Factor set IIgx2uubarxDipoleKernel:CMWScheme Factor set IIgx2ubaruxDipoleKernel:CMWScheme Factor set IIgx2ccbarxDipoleKernel:CMWScheme Factor set IIgx2cbarcxDipoleKernel:CMWScheme Factor set IIgx2ssbarxDipoleKernel:CMWScheme Factor set IIgx2sbarsxDipoleKernel:CMWScheme Factor set IFgx2bbbarxDipoleKernel:CMWScheme Factor set IFgx2bbarbxDipoleKernel:CMWScheme Factor set IFMgx2bbbarxDipoleKernel:CMWScheme Factor set IFMgx2bbarbxDipoleKernel:CMWScheme Factor set IIgx2bbbarxDipoleKernel:CMWScheme Factor set IIgx2bbarbxDipoleKernel:CMWScheme Factor +set FIMDecaygx2bbxDipoleKernel:CMWScheme Factor +set FIMDecaygx2ccxDipoleKernel:CMWScheme Factor +set FIMDecaygx2ssxDipoleKernel:CMWScheme Factor +set FIMDecaygx2uuxDipoleKernel:CMWScheme Factor +set FIMDecaygx2ddxDipoleKernel:CMWScheme Factor +set FIMDecaygx2ggxDipoleKernel:CMWScheme Factor +set FIMDecaycx2cgxDipoleKernel:CMWScheme Factor +set FIMDecaysx2sgxDipoleKernel:CMWScheme Factor +set FIMDecayux2ugxDipoleKernel:CMWScheme Factor +set FIMDecaydx2dgxDipoleKernel:CMWScheme Factor +set FIMDecaybx2bgxDipoleKernel:CMWScheme Factor - diff --git a/src/Merging/LEP-Merging.in b/src/Merging/LEP-Merging.in --- a/src/Merging/LEP-Merging.in +++ b/src/Merging/LEP-Merging.in @@ -1,155 +1,155 @@ # -*- ThePEG-repository -*- ################################################## ## Herwig/Merging example input file ################################################## ################################################## ## Collider type ################################################## read snippets/DipoleMerging.in read snippets/EECollider.in read snippets/MonacoSampler.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 MergingFactory:OrderInAlphaS 0 set MergingFactory:OrderInAlphaEW 2 ## Select the process ## You may use identifiers such as p, pbar, j, l, mu+, h0 etc. do MergingFactory:Process e- e+ -> j j [ j j ] set MergingFactory:NLOProcesses 2 set Merger:MergingScale 4.*GeV set Merger:MergingScaleSmearing 0.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 ## As massive b-quarks are currently not supported by the ## build in ME an external ME-Provider is needed (e.g. MG+OL). # 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 ################################################## ## Cut selection ## See the documentation for more options ################################################## ## cuts on additional jets cd /Herwig/Cuts # read Matchbox/DefaultEEJets.in # set NJetsCut:NJetsMin 3 ################################################## ## Scale choice ## See the documentation for more options ################################################## cd /Herwig/MatrixElements/Matchbox/Scales/ set /Herwig/Merging/MergingFactory:ScaleChoice SHatScale ################################################## ## Scale uncertainties ################################################## # read Matchbox/MuDown.in # read Matchbox/MuUp.in ################################################## ## Shower scale uncertainties ################################################## # read Matchbox/MuQDown.in # read Matchbox/MuQUp.in ################################################## ## CMW - Scheme ################################################## read snippets/Dipole_AutoTune_prel.in ### Use factor in alpha_s argument: alpha_s(q) -> alpha_s(fac*q) ### with fac=exp(-(67-3pi^2-10/3*Nf)/(33-2Nf)) read Merging/FactorCMWScheme.in ### Linear CMW multiplication: ### alpha_s(q) -> alpha_s(q)(1+K_g*alpha_s(q)/2pi ) # read Merging/LinearCMWScheme.in ################################################## ## Analyses ################################################## cd /Herwig/Analysis ## Write HepMC events. Modify the PrintEvent interface for your needs. # insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 HepMCFile ## Setup the Rivet analysis: #read snippets/Rivet.in #insert Rivet:Analyses 0 XXX_2017_ABC123 ## Here we collected a various Rivet analysis for LEP ## at the Z Mass. (The collection might not be complete.) # read Merging/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/MergingFactory:ProductionMode cd /Herwig/Generators saverun LEP-Merging EventGenerator diff --git a/src/Merging/LinearCMWScheme.in b/src/Merging/LinearCMWScheme.in --- a/src/Merging/LinearCMWScheme.in +++ b/src/Merging/LinearCMWScheme.in @@ -1,85 +1,96 @@ ####### # CMW scheme ####### -cd /Herwig/Merging -set Merger:CMWScheme Linear - cd /Herwig/DipoleShower/Kernels set FFgx2ggxDipoleKernel:CMWScheme Linear set FFqx2qgxDipoleKernel:CMWScheme Linear set FFgx2ddxDipoleKernel:CMWScheme Linear set FFgx2uuxDipoleKernel:CMWScheme Linear set FFgx2ccxDipoleKernel:CMWScheme Linear set FFgx2ssxDipoleKernel:CMWScheme Linear set FFgx2bbxDipoleKernel:CMWScheme Linear set FFMgx2ggxDipoleKernel:CMWScheme Linear set FFMdx2dgxDipoleKernel:CMWScheme Linear set FFMux2ugxDipoleKernel:CMWScheme Linear set FFMcx2cgxDipoleKernel:CMWScheme Linear set FFMsx2sgxDipoleKernel:CMWScheme Linear set FFMbx2bgxDipoleKernel:CMWScheme Linear set FFMtx2tgxDipoleKernel:CMWScheme Linear set FFMgx2ddxDipoleKernel:CMWScheme Linear set FFMgx2uuxDipoleKernel:CMWScheme Linear set FFMgx2ccxDipoleKernel:CMWScheme Linear set FFMgx2ssxDipoleKernel:CMWScheme Linear set FFMgx2bbxDipoleKernel:CMWScheme Linear set FIgx2ggxDipoleKernel:CMWScheme Linear set FIqx2qgxDipoleKernel:CMWScheme Linear set FIgx2ddxDipoleKernel:CMWScheme Linear set FIgx2uuxDipoleKernel:CMWScheme Linear set FIgx2ccxDipoleKernel:CMWScheme Linear set FIgx2ssxDipoleKernel:CMWScheme Linear set FIgx2bbxDipoleKernel:CMWScheme Linear set FIMdx2dgxDipoleKernel:CMWScheme Linear set FIMux2ugxDipoleKernel:CMWScheme Linear set FIMcx2cgxDipoleKernel:CMWScheme Linear set FIMsx2sgxDipoleKernel:CMWScheme Linear set FIMbx2bgxDipoleKernel:CMWScheme Linear set FIMtx2tgxDipoleKernel:CMWScheme Linear set FIMgx2ddxDipoleKernel:CMWScheme Linear set FIMgx2uuxDipoleKernel:CMWScheme Linear set FIMgx2ccxDipoleKernel:CMWScheme Linear set FIMgx2ssxDipoleKernel:CMWScheme Linear set FIMgx2bbxDipoleKernel:CMWScheme Linear #set FIMgx2ttxDipoleKernel:CMWScheme Linear set IFgx2ggxDipoleKernel:CMWScheme Linear set IFqx2qgxDipoleKernel:CMWScheme Linear set IFqx2gqxDipoleKernel:CMWScheme Linear set IFgx2ddbarxDipoleKernel:CMWScheme Linear set IFgx2dbardxDipoleKernel:CMWScheme Linear set IFgx2uubarxDipoleKernel:CMWScheme Linear set IFgx2ubaruxDipoleKernel:CMWScheme Linear set IFgx2ccbarxDipoleKernel:CMWScheme Linear set IFgx2cbarcxDipoleKernel:CMWScheme Linear set IFgx2ssbarxDipoleKernel:CMWScheme Linear set IFgx2sbarsxDipoleKernel:CMWScheme Linear set IFMgx2ggxDipoleKernel:CMWScheme Linear set IFMqx2qgxDipoleKernel:CMWScheme Linear set IFMqx2gqxDipoleKernel:CMWScheme Linear set IFMgx2ddbarxDipoleKernel:CMWScheme Linear set IFMgx2dbardxDipoleKernel:CMWScheme Linear set IFMgx2uubarxDipoleKernel:CMWScheme Linear set IFMgx2ubaruxDipoleKernel:CMWScheme Linear set IFMgx2ccbarxDipoleKernel:CMWScheme Linear set IFMgx2cbarcxDipoleKernel:CMWScheme Linear set IFMgx2ssbarxDipoleKernel:CMWScheme Linear set IFMgx2sbarsxDipoleKernel:CMWScheme Linear set IIgx2ggxDipoleKernel:CMWScheme Linear set IIqx2qgxDipoleKernel:CMWScheme Linear set IIqx2gqxDipoleKernel:CMWScheme Linear set IIgx2ddbarxDipoleKernel:CMWScheme Linear set IIgx2dbardxDipoleKernel:CMWScheme Linear set IIgx2uubarxDipoleKernel:CMWScheme Linear set IIgx2ubaruxDipoleKernel:CMWScheme Linear set IIgx2ccbarxDipoleKernel:CMWScheme Linear set IIgx2cbarcxDipoleKernel:CMWScheme Linear set IIgx2ssbarxDipoleKernel:CMWScheme Linear set IIgx2sbarsxDipoleKernel:CMWScheme Linear set IFgx2bbbarxDipoleKernel:CMWScheme Linear set IFgx2bbarbxDipoleKernel:CMWScheme Linear set IFMgx2bbbarxDipoleKernel:CMWScheme Linear set IFMgx2bbarbxDipoleKernel:CMWScheme Linear set IIgx2bbbarxDipoleKernel:CMWScheme Linear set IIgx2bbarbxDipoleKernel:CMWScheme Linear + + +set FIMDecaygx2bbxDipoleKernel:CMWScheme Linear +set FIMDecaygx2ccxDipoleKernel:CMWScheme Linear +set FIMDecaygx2ssxDipoleKernel:CMWScheme Linear +set FIMDecaygx2uuxDipoleKernel:CMWScheme Linear +set FIMDecaygx2ddxDipoleKernel:CMWScheme Linear +set FIMDecaygx2ggxDipoleKernel:CMWScheme Linear +set FIMDecaycx2cgxDipoleKernel:CMWScheme Linear +set FIMDecaysx2sgxDipoleKernel:CMWScheme Linear +set FIMDecayux2ugxDipoleKernel:CMWScheme Linear +set FIMDecaydx2dgxDipoleKernel:CMWScheme Linear +set FIMDecaybx2bgxDipoleKernel:CMWScheme Linear +