diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc b/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc --- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc +++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc @@ -1,959 +1,957 @@ // -*- C++ -*- // // MatchboxAmplitude.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 MatchboxAmplitude class. // #include "MatchboxAmplitude.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/MatrixElement/Matchbox/Utility/SpinorHelicity.h" #include "Herwig/MatrixElement/Matchbox/Utility/SU2Helper.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h" #include "ThePEG/Utilities/StringUtils.h" #include "MatchboxMEBase.h" #include #include #include using std::ostream_iterator; using namespace Herwig; MatchboxAmplitude::MatchboxAmplitude() : Amplitude(), theCleanupAfter(20), treeLevelHelicityPoints(0), oneLoopHelicityPoints(0), theTrivialColourLegs(false) { } -MatchboxAmplitude::~MatchboxAmplitude() {} - void MatchboxAmplitude::persistentOutput(PersistentOStream & os) const { os << theLastXComb << theColourBasis << theCleanupAfter << treeLevelHelicityPoints << oneLoopHelicityPoints << theTrivialColourLegs << theReshuffleMasses.size(); if ( !theReshuffleMasses.empty() ) { for (auto const & r : theReshuffleMasses ) os << r.first << ounit(r.second,GeV); } } void MatchboxAmplitude::persistentInput(PersistentIStream & is, int) { size_t reshuffleSize; is >> theLastXComb >> theColourBasis >> theCleanupAfter >> treeLevelHelicityPoints >> oneLoopHelicityPoints >> theTrivialColourLegs >> reshuffleSize; theReshuffleMasses.clear(); while ( reshuffleSize > 0 ) { long id; Energy m; is >> id >> iunit(m,GeV); theReshuffleMasses[id] = m; --reshuffleSize; } lastMatchboxXComb(theLastXComb); } Ptr::tptr MatchboxAmplitude::factory() const { return MatchboxFactory::currentFactory(); } void MatchboxAmplitude::doinit() { Amplitude::doinit(); if ( colourBasis() ) { colourBasis()->init(); } } void MatchboxAmplitude::doinitrun() { Amplitude::doinitrun(); if ( colourBasis()) colourBasis()->initrun(); } void MatchboxAmplitude::cloneDependencies(const std::string&,bool) {} MatchboxMEBasePtr MatchboxAmplitude::makeME(const PDVector&) const { return new_ptr(MatchboxMEBase()); } Selector MatchboxAmplitude::colourGeometries(tcDiagPtr d) const { if ( haveColourFlows() ) return colourBasis()->colourGeometries(d,lastLargeNAmplitudes()); return Selector(); } void MatchboxAmplitude::olpOrderFileHeader(ostream& os) const { os << "# OLP order file created by Herwig/Matchbox\n\n"; os << "InterfaceVersion BLHA2\n\n"; os << "Model SM\n" << "CorrectionType QCD\n" << "IRregularisation " << (isDR() ? "DRED" : "CDR") << "\n" << "Extra HelAvgInitial no\n" << "Extra ColAvgInitial no\n" << "Extra MCSymmetrizeFinal no\n"; os << "\n"; } void MatchboxAmplitude::olpOrderFileProcesses(ostream& os, const map,int>& proc) const { map > sorted; for (auto const & p : proc ) { sorted[p.second] = p.first; } unsigned int currentOrderInAlphaS = sorted.begin()->second.first.orderInAlphaS; unsigned int currentOrderInAlphaEW = sorted.begin()->second.first.orderInAlphaEW; int currentType = sorted.begin()->second.second; os << "AlphasPower " << currentOrderInAlphaS << "\n" << "AlphaPower " << currentOrderInAlphaEW << "\n" << "AmplitudeType "; if ( currentType == ProcessType::treeME2 ) { os << "tree\n"; } else if ( currentType == ProcessType::oneLoopInterference ) { os << "loop\n"; } else if ( currentType == ProcessType::colourCorrelatedME2 ) { os << "cctree\n"; } else if ( currentType == ProcessType::spinColourCorrelatedME2 ) { os << "sctree\n"; } else if ( currentType == ProcessType::loopInducedME2 ) { os << "loopinduced\n"; } else if ( currentType == ProcessType::spinCorrelatedME2 ) { os << "stree\n"; } else assert(false); for ( map >::const_iterator p = sorted.begin(); p != sorted.end(); ++p ) { if ( currentOrderInAlphaS != p->second.first.orderInAlphaS ) { currentOrderInAlphaS = p->second.first.orderInAlphaS; os << "AlphasPower " << currentOrderInAlphaS << "\n"; } if ( currentOrderInAlphaEW != p->second.first.orderInAlphaEW ) { currentOrderInAlphaEW = p->second.first.orderInAlphaEW; os << "AlphaPower " << currentOrderInAlphaEW << "\n"; } if ( currentType != p->second.second ) { currentType = p->second.second; os << "AmplitudeType "; if ( currentType == ProcessType::treeME2 ) { os << "tree\n"; } else if ( currentType == ProcessType::oneLoopInterference ) { os << "loop\n"; } else if ( currentType == ProcessType::colourCorrelatedME2 ) { os << "cctree\n"; } else if ( currentType == ProcessType::spinColourCorrelatedME2 ) { os << "sctree\n"; } else if ( currentType == ProcessType::spinCorrelatedME2 ) { os << "stree\n"; } else assert(false); } os << p->second.first.legs[0]->id() << " " << p->second.first.legs[1]->id() << " -> "; for ( PDVector::const_iterator o = p->second.first.legs.begin() + 2; o != p->second.first.legs.end(); ++o ) { os << (**o).id() << " "; } os << "\n"; } } bool MatchboxAmplitude::startOLP(const map,int>& procs) { string orderFileName = factory()->buildStorage() + name() + ".OLPOrder.lh"; ofstream orderFile(orderFileName.c_str()); olpOrderFileHeader(orderFile); olpOrderFileProcesses(orderFile,procs); string contractFileName = factory()->buildStorage() + name() + ".OLPContract.lh"; signOLP(orderFileName, contractFileName); // TODO check the contract file int status = 0; startOLP(contractFileName, status); if ( status != 1 ) return false; return true; } struct orderPartonData { bool operator()(const pair& a, const pair& b) const { if ( a.first == b.first ) return a.second < b.second; int acolour = a.first->iColour(); int bcolour = b.first->iColour(); if ( abs(acolour) != abs(bcolour) ) return abs(acolour) < abs(bcolour); if ( a.first->iSpin() != b.first->iSpin() ) return a.first->iSpin() < b.first->iSpin(); int acharge = a.first->iCharge(); int bcharge = b.first->iCharge(); if ( abs(acharge) != abs(bcharge) ) return abs(acharge) < abs(bcharge); if ( abs(a.first->id()) != abs(b.first->id()) ) return abs(a.first->id()) < abs(b.first->id()); return a.first->id() > b.first->id(); } }; void MatchboxAmplitude::setXComb(tStdXCombPtr xc) { theLastXComb = xc; lastMatchboxXComb(xc); fillCrossingMap(); if ( ( treeAmplitudes() || oneLoopAmplitudes() ) && hasAmplitudeMomenta() ) for ( size_t k = 0 ; k < meMomenta().size(); ++k ) amplitudeMomenta()[k] = amplitudeMomentum(k); } void MatchboxAmplitude::fillCrossingMap(size_t shift) { if ( !amplitudePartonData().empty() ) return; double csign = 1.; set,orderPartonData > processLegs; for ( unsigned int l = 0; l < mePartonData().size(); ++l ) { if ( l > 1 ) processLegs.insert(make_pair(mePartonData()[l],l)); else { if ( mePartonData()[l]->CC() ) { processLegs.insert(make_pair(mePartonData()[l]->CC(),l)); if ( mePartonData()[l]->iSpin() == PDT::Spin1Half ) csign *= -1.; } else { processLegs.insert(make_pair(mePartonData()[l],l)); } } } crossingSign(csign); set > amplitudeLegs; crossingMap().resize(mePartonData().size()); amplitudePartonData().resize(mePartonData().size()); if ( hasAmplitudeMomenta() ) amplitudeMomenta().resize(mePartonData().size()); int ampCount = 0; // process legs are already sorted, we only need to arrange for // adjacent particles and anti-particles while ( !processLegs.empty() ) { set,orderPartonData >::iterator next = processLegs.begin(); while ( next->first->id() < 0 ) { if ( ++next == processLegs.end() ) break; } //This happens for e.g. p p-> W- gamma & p p->W- W- j j //Still working for pp->W-H-W- e+ nue jj ??? if(next == processLegs.end()){ next = processLegs.begin(); for (;next!=processLegs.end();next++){ assert(next->first->id() < 0 ); crossingMap()[ampCount] = next->second - shift; amplitudeLegs.insert(make_pair(next->first,ampCount)); ++ampCount; processLegs.erase(next); } break; } crossingMap()[ampCount] = next->second - shift; amplitudeLegs.insert(make_pair(next->first,ampCount)); tcPDPtr check = next->first; processLegs.erase(next); ++ampCount; if ( check->CC() ) { set,orderPartonData>::iterator checkcc = processLegs.end(); for ( set,orderPartonData>::iterator c = processLegs.begin(); c != processLegs.end(); ++c ) { if ( c->first == check->CC() ) { checkcc = c; break; } } if ( checkcc == processLegs.end() ) for ( set,orderPartonData>::iterator c = processLegs.begin(); c != processLegs.end(); ++c ) { if ( !SU2Helper::SU2CC(check) ) continue; if ( c->first == SU2Helper::SU2CC(check)->CC() ) { checkcc = c; break; } } if ( checkcc == processLegs.end() ) { int f = SU2Helper::family(check); for ( int i = 1 - f; i < 5 - f; i++ ) { bool gotone = false; for ( set,orderPartonData>::iterator c = processLegs.begin(); c != processLegs.end(); ++c ) { if ( !SU2Helper::SU2CC(check,i) ) continue; if ( c->first == SU2Helper::SU2CC(check,i)->CC() ) { checkcc = c; gotone = true; break; } } if ( gotone ) break; } } // default to just pick the next available anti-particle if ( processLegs.empty() ) break; if ( checkcc == processLegs.end() ) { checkcc = processLegs.begin(); while ( checkcc->first->id() > 0 ) if ( ++checkcc == processLegs.end() ) break; } // if still not there, use whatever is available at the end if ( checkcc == processLegs.end() ) checkcc = processLegs.begin(); crossingMap()[ampCount] = checkcc->second - shift; amplitudeLegs.insert(make_pair(checkcc->first,ampCount)); processLegs.erase(checkcc); ++ampCount; } } for ( set >::const_iterator l = amplitudeLegs.begin(); l != amplitudeLegs.end(); ++l ) amplitudePartonData()[l->second] = l->first; if ( colourBasis() && !colourBasis()->indexMap().empty()) { assert(colourBasis()->indexMap().find(mePartonData()) != colourBasis()->indexMap().end()); const map colourCross = colourBasis()->indexMap().find(mePartonData())->second; for ( size_t k = 0; k < crossingMap().size(); ++k ) { if ( colourCross.find(crossingMap()[k]) != colourCross.end() ) { size_t ccross = colourCross.find(crossingMap()[k])->second; amplitudeToColourMap()[k] = ccross; colourToAmplitudeMap()[ccross] = k; } } } } const string& MatchboxAmplitude::colourOrderingString(size_t id) const { static string empty = ""; if ( !colourBasis() ) { return empty; } return colourBasis()->orderingString(mePartonData(),colourToAmplitudeMap(),id); } const set >& MatchboxAmplitude::colourOrdering(size_t id) const { static set > empty; if ( !colourBasis() ) { return empty; } return colourBasis()->ordering(mePartonData(),colourToAmplitudeMap(),id); } Lorentz5Momentum MatchboxAmplitude::amplitudeMomentum(int i) const { int iCrossed = crossingMap()[i]; Lorentz5Momentum res = meMomenta()[iCrossed]; if ( iCrossed < 2 ) res = -res; res.setMass(meMomenta()[iCrossed].mass()); Energy2 rho = res.t()*res.t() - res.mass2(); res.setRho(sqrt(abs(rho))); return res; } set > MatchboxAmplitude::generateHelicities() const { set > res; vector current(amplitudePartonData().size()); doGenerateHelicities(res,current,0); return res; } void MatchboxAmplitude::doGenerateHelicities(set >& res, vector& current, size_t pos) const { if ( pos == amplitudePartonData().size() ) { res.insert(current); return; } if ( amplitudePartonData()[pos]->iSpin() == PDT::Spin0 ) { current[pos] = 0; doGenerateHelicities(res,current,pos+1); } else if ( amplitudePartonData()[pos]->iSpin() == PDT::Spin1Half ) { current[pos] = 1; doGenerateHelicities(res,current,pos+1); current[pos] = -1; doGenerateHelicities(res,current,pos+1); }else if (amplitudePartonData()[pos]->iSpin() == PDT::Spin1 ) { if (amplitudePartonData()[pos]->hardProcessMass() != ZERO){ current[pos] = 0; doGenerateHelicities(res,current,pos+1); } current[pos] = 1; doGenerateHelicities(res,current,pos+1); current[pos] = -1; doGenerateHelicities(res,current,pos+1); } } vector MatchboxAmplitude::physicalHelicities(const vector&) const { throw Exception() << "MatchboxAmplitude::physicalHelicities(): The amplitude '" << name() << "' does not support the spin correlation algorithm" << Exception::runerror; static vector dummy; return dummy; } void MatchboxAmplitude::prepareAmplitudes(Ptr::tcptr) { if ( !calculateTreeAmplitudes() ) return; bool initialized = !lastAmplitudes().empty() && treeLevelHelicityPoints > theCleanupAfter; if ( !initialized ) { treeLevelHelicityPoints++; map,CVector> all; map,CVector> allLargeN; set > helicities = generateHelicities(); for ( set >::const_iterator h = helicities.begin(); h != helicities.end(); ++h ) { all.insert(make_pair(*h,CVector(max(colourBasisDim(),1)))); allLargeN.insert(make_pair(*h,CVector(max(colourBasisDim(),1)))); } AmplitudeIterator amp = all.begin(); AmplitudeIterator lamp = allLargeN.begin(); for ( ; amp != all.end(); ++amp, ++lamp ) { for ( size_t k = 0; k < max(colourBasisDim(),1); ++k ){ amp->second(k) = evaluate(k,amp->first,lamp->second(k)); if ( amp->second(k) != Complex(0.0) ) { if ( lastAmplitudes().find(amp->first)!=lastAmplitudes().end() ) { lastAmplitudes().find(amp->first)->second = amp->second; lastLargeNAmplitudes().find(lamp->first)->second = lamp->second; } else { lastAmplitudes().insert(*amp); lastLargeNAmplitudes().insert(*lamp); } } else if ( lastAmplitudes().find(amp->first)!=lastAmplitudes().end() ){ lastAmplitudes().find(amp->first)->second = amp->second; lastLargeNAmplitudes().find(lamp->first)->second = lamp->second; } } } } else { AmplitudeIterator amp = lastAmplitudes().begin(); AmplitudeIterator lamp = lastLargeNAmplitudes().begin(); for ( ;amp != lastAmplitudes().end(); ++amp, ++lamp ) { for ( size_t k = 0; k < max(colourBasisDim(),1); ++k ){ amp->second(k) = evaluate(k,amp->first,lamp->second(k)); } } } haveTreeAmplitudes(); } void MatchboxAmplitude::prepareOneLoopAmplitudes(Ptr::tcptr) { if ( !calculateOneLoopAmplitudes() ) return; bool initialized = !lastOneLoopAmplitudes().empty() && oneLoopHelicityPoints > theCleanupAfter; if ( !initialized ) { oneLoopHelicityPoints++; map,CVector> all; set > helicities = generateHelicities(); for ( set >::const_iterator h = helicities.begin(); h != helicities.end(); ++h ) { all.insert(make_pair(*h,CVector(max(colourBasisDim(),1)))); } AmplitudeIterator amp = all.begin(); for ( ; amp != all.end(); ++amp ) { for ( size_t k = 0; k < max(colourBasisDim(),1); ++k ){ amp->second(k) = evaluateOneLoop(k,amp->first); if ( amp->second(k) != Complex(0.0) ) { if ( lastOneLoopAmplitudes().find(amp->first)!=lastOneLoopAmplitudes().end() ) { lastOneLoopAmplitudes().find(amp->first)->second = amp->second; } else{ lastOneLoopAmplitudes().insert(*amp); } } else if ( lastOneLoopAmplitudes().find(amp->first)!=lastOneLoopAmplitudes().end() ){ lastOneLoopAmplitudes().find(amp->first)->second = amp->second; } } } } else { AmplitudeIterator amp = lastOneLoopAmplitudes().begin(); for ( ;amp != lastOneLoopAmplitudes().end(); ++amp ) { for ( size_t k = 0; k < max(colourBasisDim(),1); ++k ){ amp->second(k) = evaluateOneLoop(k,amp->first); } } } haveOneLoopAmplitudes(); } Complex MatchboxAmplitude::value(const tcPDVector&, const vector&, const vector&) { assert(false && "ThePEG::Amplitude interface is not sufficient at the moment."); throw Exception() << "MatchboxAmplitude::value(): ThePEG::Amplitude interface is not sufficient at the moment." << Exception::runerror; return 0.; } double MatchboxAmplitude::me2() const { if ( !calculateTreeME2() ) return lastTreeME2(); lastTreeME2(colourBasis()->me2(mePartonData(),lastAmplitudes())); return lastTreeME2(); } double MatchboxAmplitude::largeNME2(Ptr::tptr largeNBasis) const { if ( !calculateLargeNME2() ) return lastLargeNME2(); double res = 0.; if ( !trivialColourLegs() ) res = largeNBasis->me2(mePartonData(),lastLargeNAmplitudes()); else res = me2(); lastLargeNME2(res); return res; } double MatchboxAmplitude::oneLoopInterference() const { if ( !calculateOneLoopInterference() ) return lastOneLoopInterference(); lastOneLoopInterference(colourBasis()->interference(mePartonData(), lastOneLoopAmplitudes(),lastAmplitudes())); return lastOneLoopInterference(); } double MatchboxAmplitude::colourCharge(tcPDPtr data) const { double cfac = 1.; double Nc = generator()->standardModel()->Nc(); if ( data->iColour() == PDT::Colour8 ) { cfac = Nc; } else if ( data->iColour() == PDT::Colour3 || data->iColour() == PDT::Colour3bar ) { cfac = (sqr(Nc)-1.)/(2.*Nc); } else assert(false); return cfac; } double MatchboxAmplitude::largeNColourCharge(tcPDPtr data) const { double cfac = 1.; double Nc = generator()->standardModel()->Nc(); if ( data->iColour() == PDT::Colour8 ) { cfac = Nc; } else if ( data->iColour() == PDT::Colour3 || data->iColour() == PDT::Colour3bar ) { cfac = 0.5*Nc; } else assert(false); return cfac; } double MatchboxAmplitude::colourCorrelatedME2(pair ij) const { double cfac = colourCharge(mePartonData()[ij.first]); if ( !calculateColourCorrelator(ij) ) return lastColourCorrelator(ij)/cfac; double res = 0.; if ( !trivialColourLegs() ) res = colourBasis()->colourCorrelatedME2(ij,mePartonData(),lastAmplitudes()); else { // two partons TiTj = (-Ti^2-Tj^2)/2 // three partons TiTj = (-Ti^2-Tj^2+Tk^2)/2 int n = mePartonData().size(); for ( int i = 0; i < n; ++i ) { if ( !mePartonData()[i]->coloured() ) continue; if ( i == ij.first || i == ij.second ) res -= colourCharge(mePartonData()[i])*me2(); else res += colourCharge(mePartonData()[i])*me2(); } res *= 0.5; } lastColourCorrelator(ij,res); return res/cfac; } double MatchboxAmplitude::largeNColourCorrelatedME2(pair ij, Ptr::tptr largeNBasis) const { double cfac = largeNColourCharge(mePartonData()[ij.first]); if ( !calculateLargeNColourCorrelator(ij) ) return lastLargeNColourCorrelator(ij)/cfac; double res = 0.; if ( !trivialColourLegs() ) res = largeNBasis->colourCorrelatedME2(ij,mePartonData(),lastLargeNAmplitudes()); else { // two partons TiTj = (-Ti^2-Tj^2)/2 // three partons TiTj = (-Ti^2-Tj^2+Tk^2)/2 int n = mePartonData().size(); for ( int i = 0; i < n; ++i ) { if ( !mePartonData()[i]->coloured() ) continue; if ( i == ij.first || i == ij.second ) res -= largeNColourCharge(mePartonData()[i])*largeNME2(largeNBasis); else res += largeNColourCharge(mePartonData()[i])*largeNME2(largeNBasis); } res *= 0.5; } lastLargeNColourCorrelator(ij,res); return res/cfac; } // compare int vectors modulo certain element // which needs to differe between the two bool equalsModulo(unsigned int i, const vector& a, const vector& b) { assert(a.size()==b.size()); if ( a[i] == b[i] ) return false; for ( unsigned int k = 0; k < a.size(); ++k ) { if ( k == i ) continue; if ( a[k] != b[k] ) return false; } return true; } LorentzVector MatchboxAmplitude::plusPolarization(const Lorentz5Momentum& p, const Lorentz5Momentum& n, int) const { using namespace SpinorHelicity; LorentzVector > num = PlusSpinorCurrent(PlusConjugateSpinor(n),MinusSpinor(p)).eval(); complex den = sqrt(2.)*PlusSpinorProduct(PlusConjugateSpinor(n),PlusSpinor(p)).eval(); LorentzVector polarization(num.x()/den,num.y()/den,num.z()/den,num.t()/den); return polarization; } double MatchboxAmplitude::spinColourCorrelatedME2(pair ij, const SpinCorrelationTensor& c) const { Lorentz5Momentum p = meMomenta()[ij.first]; Lorentz5Momentum n = meMomenta()[ij.second]; LorentzVector polarization = plusPolarization(p,n,ij.first); Complex pFactor = (polarization*c.momentum())/sqrt(abs(c.scale())); double avg = colourCorrelatedME2(ij)*(-c.diagonal()+ (c.scale() > ZERO ? 1. : -1.)*norm(pFactor)); int iCrossed = -1; for ( unsigned int k = 0; k < crossingMap().size(); ++k ) if ( crossingMap()[k] == ij.first ) { iCrossed = k; break; } assert(iCrossed >= 0); Complex csCorr = 0.0; if ( calculateColourSpinCorrelator(ij) ) { set done; for ( AmplitudeConstIterator a = lastAmplitudes().begin(); a != lastAmplitudes().end(); ++a ) { if ( done.find(&(a->second)) != done.end() ) continue; AmplitudeConstIterator b = lastAmplitudes().begin(); while ( !equalsModulo(iCrossed,a->first,b->first) ) if ( ++b == lastAmplitudes().end() ) break; if ( b == lastAmplitudes().end() || done.find(&(b->second)) != done.end() ) continue; done.insert(&(a->second)); done.insert(&(b->second)); if ( a->first[iCrossed] == 1 ) swap(a,b); csCorr += colourBasis()->colourCorrelatedInterference(ij,mePartonData(),a->second,b->second); } lastColourSpinCorrelator(ij,csCorr); } else { csCorr = lastColourSpinCorrelator(ij); } double corr = 2.*real(csCorr*sqr(pFactor)); double Nc = generator()->standardModel()->Nc(); double cfac = 1.; if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) { cfac = Nc; } else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 || mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) { cfac = (sqr(Nc)-1.)/(2.*Nc); } else assert(false); return avg + (c.scale() > ZERO ? 1. : -1.)*corr/cfac; } double MatchboxAmplitude::spinCorrelatedME2(pair ij, const SpinCorrelationTensor& c) const { Lorentz5Momentum p = meMomenta()[ij.first]; Lorentz5Momentum n = meMomenta()[ij.second]; LorentzVector polarization = plusPolarization(p,n,ij.first); Complex pFactor = (polarization*c.momentum())/sqrt(abs(c.scale())); double avg = me2()*(-c.diagonal()+ (c.scale() > ZERO ? 1. : -1.)*norm(pFactor)); int iCrossed = -1; for ( unsigned int k = 0; k < crossingMap().size(); ++k ) if ( crossingMap()[k] == ij.first ) { iCrossed = k; break; } assert(iCrossed >= 0); Complex csCorr = 0.0; if ( calculateSpinCorrelator(ij) ) { set done; for ( AmplitudeConstIterator a = lastAmplitudes().begin(); a != lastAmplitudes().end(); ++a ) { if ( done.find(&(a->second)) != done.end() ) continue; AmplitudeConstIterator b = lastAmplitudes().begin(); while ( !equalsModulo(iCrossed,a->first,b->first) ) if ( ++b == lastAmplitudes().end() ) break; if ( b == lastAmplitudes().end() || done.find(&(b->second)) != done.end() ) continue; done.insert(&(a->second)); done.insert(&(b->second)); if ( a->first[iCrossed] == 1 ) swap(a,b); csCorr += colourBasis()->interference(mePartonData(),a->second,b->second); } lastSpinCorrelator(ij,csCorr); } else { csCorr = lastSpinCorrelator(ij); } double corr = 2.*real(csCorr*sqr(pFactor)); return avg + (c.scale() > ZERO ? 1. : -1.)*corr; } void MatchboxAmplitude::checkReshuffling(Ptr::tptr ps) { set noReshuffle; for ( map::const_iterator m = reshuffleMasses().begin(); m != reshuffleMasses().end(); ++m ) { tcPDPtr data = getParticleData(m->first); assert(data); bool needReshuffle = m->second != data->hardProcessMass(); needReshuffle |= (data->hardProcessWidth() != ZERO || data->massGenerator()) && ps->useMassGenerators(); if ( !needReshuffle ) noReshuffle.insert(m->first); } for ( set::const_iterator rm = noReshuffle.begin(); rm != noReshuffle.end(); ++rm ) theReshuffleMasses.erase(*rm); } string MatchboxAmplitude::doReshuffle(string in) { in = StringUtils::stripws(in); if ( in.empty() ) throw Exception() << "MatchboxAmplitude::doReshuffle(): Expecting PDG id and mass value" << Exception::runerror; istringstream ins(in); long id; ins >> id; if ( ins.eof() ) throw Exception() << "MatchboxAmplitude::doReshuffle(): expecting PDG id and mass value" << Exception::runerror; Energy m; ins >> iunit(m,GeV); theReshuffleMasses[id] = m; return ""; } string MatchboxAmplitude::doMassless(string in) { in = StringUtils::stripws(in); if ( in.empty() ) throw Exception() << "MatchboxAmplitude::doMassless(): Expecting PDG id" << Exception::runerror; istringstream ins(in); long id; ins >> id; theReshuffleMasses[id] = ZERO; return ""; } string MatchboxAmplitude::doOnShell(string in) { in = StringUtils::stripws(in); if ( in.empty() ) throw Exception() << "MatchboxAmplitude::doOnShell(): Expecting PDG id" << Exception::runerror; istringstream ins(in); long id; ins >> id; tcPDPtr data = getParticleData(id); assert(data); theReshuffleMasses[id] = data->hardProcessMass(); return ""; } string MatchboxAmplitude::doClearReshuffling(string) { theReshuffleMasses.clear(); return ""; } void MatchboxAmplitude::Init() { static ClassDocumentation documentation ("MatchboxAmplitude is the base class for amplitude " "implementations inside Matchbox."); static Reference interfaceColourBasis ("ColourBasis", "Set the colour basis implementation.", &MatchboxAmplitude::theColourBasis, false, false, true, true, false); static Parameter interfaceCleanupAfter ("CleanupAfter", "The number of points after which helicity combinations are cleaned up.", &MatchboxAmplitude::theCleanupAfter, 20, 1, 0, false, false, Interface::lowerlim); interfaceCleanupAfter.rank(-1); static Command interfaceReshuffle ("Reshuffle", "Reshuffle the mass for the given PDG id to a different mass shell for amplitude evaluation.", &MatchboxAmplitude::doReshuffle, false); interfaceReshuffle.rank(-1); static Command interfaceMassless ("Massless", "Reshuffle the mass for the given PDG id to be massless for amplitude evaluation.", &MatchboxAmplitude::doMassless, false); interfaceMassless.rank(-1); static Command interfaceOnShell ("OnShell", "Reshuffle the mass for the given PDG id to be the on-shell mass for amplitude evaluation.", &MatchboxAmplitude::doOnShell, false); interfaceOnShell.rank(-1); static Command interfaceClearReshuffling ("ClearReshuffling", "Do not perform any reshuffling.", &MatchboxAmplitude::doClearReshuffling, false); interfaceClearReshuffling.rank(-1); static Switch interfaceTrivialColourLegs ("TrivialColourLegs", "Assume the process considered has trivial colour correllations.", &MatchboxAmplitude::theTrivialColourLegs, false, false, false); static SwitchOption interfaceTrivialColourLegsYes (interfaceTrivialColourLegs, "Yes", "", true); static SwitchOption interfaceTrivialColourLegsNo (interfaceTrivialColourLegs, "No", "", false); interfaceTrivialColourLegs.rank(-1); } // *** 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). DescribeAbstractClass describeMatchboxAmplitude("Herwig::MatchboxAmplitude", "Herwig.so"); diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h --- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h +++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h @@ -1,714 +1,706 @@ // -*- C++ -*- // // MatchboxAmplitude.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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_MatchboxAmplitude_H #define HERWIG_MatchboxAmplitude_H // // This is the declaration of the MatchboxAmplitude class. // #include "ThePEG/MatrixElement/Amplitude.h" #include "ThePEG/Handlers/LastXCombInfo.h" #include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h" #include "Herwig/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h" #include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h" #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXComb.h" #include "Herwig/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief Process information with coupling order */ struct Process { PDVector legs; unsigned int orderInAlphaS; unsigned int orderInAlphaEW; Process() : orderInAlphaS(0), orderInAlphaEW(0) {} Process(const PDVector& p, unsigned int oas, unsigned int oae) : legs(p), orderInAlphaS(oas), orderInAlphaEW(oae) {} bool operator==(const Process& other) const { return legs == other.legs && orderInAlphaS == other.orderInAlphaS && orderInAlphaEW == other.orderInAlphaEW; } bool operator<(const Process& other) const { if ( orderInAlphaS != other.orderInAlphaS ) return orderInAlphaS < other.orderInAlphaS; if ( orderInAlphaEW != other.orderInAlphaEW ) return orderInAlphaEW < other.orderInAlphaEW; return legs < other.legs; } void persistentOutput(PersistentOStream & os) const { os << legs << orderInAlphaS << orderInAlphaEW; } void persistentInput(PersistentIStream & is) { is >> legs >> orderInAlphaS >> orderInAlphaEW; } }; /** * \ingroup Matchbox * \author Simon Platzer * * \brief Enumerate the type of calculation required */ namespace ProcessType { enum Types { treeME2 = 0, colourCorrelatedME2, spinColourCorrelatedME2, oneLoopInterference, loopInducedME2, spinCorrelatedME2 }; } /** * \ingroup Matchbox * \author Simon Platzer * * \brief MatchboxAmplitude is the base class for amplitude * implementations inside Matchbox. * * @see \ref MatchboxAmplitudeInterfaces "The interfaces" * defined for MatchboxAmplitude. */ class MatchboxAmplitude: public Amplitude, public LastXCombInfo, public LastMatchboxXCombInfo { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ MatchboxAmplitude(); - /** - * The destructor. - */ - virtual ~MatchboxAmplitude(); - //@} - public: /** * Return the amplitude. Needs to be implemented from * ThePEG::Amplitude but is actually ill-defined, as colours of the * external particles are not specified. To this extent, this * implementation just asserts. */ virtual Complex value(const tcPDVector & particles, const vector & momenta, const vector & helicities); /** * Return the factory which produced this matrix element */ Ptr::tptr factory() const; /** @name Subprocess information */ //@{ /** * Return true, if this amplitude can handle the given process. */ virtual bool canHandle(const PDVector& p, Ptr::tptr, bool) const { return canHandle(p); } /** * Return true, if this amplitude can handle the given process. */ virtual bool canHandle(const PDVector&) const { return false; } /** * Return the number of random numbers required to evaluate this * amplitude at a fixed phase space point. */ virtual int nDimAdditional() const { return 0; } /** * Return a ME instance appropriate for this amplitude and the given * subprocesses */ virtual Ptr::ptr makeME(const PDVector&) const; /** * Set the (tree-level) order in \f$g_S\f$ in which this matrix * element should be evaluated. */ virtual void orderInGs(unsigned int) {} /** * Return the (tree-level) order in \f$g_S\f$ in which this matrix * element is given. */ virtual unsigned int orderInGs() const = 0; /** * Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix * element should be evaluated. */ virtual void orderInGem(unsigned int) {} /** * Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix * element is given. */ virtual unsigned int orderInGem() const = 0; /** * Return the Herwig StandardModel object */ Ptr::tcptr standardModel() { if ( !hwStandardModel() ) hwStandardModel(dynamic_ptr_cast::tcptr>(HandlerBase::standardModel())); return hwStandardModel(); } /** * Return true, if this amplitude already includes averaging over * incoming parton's quantum numbers. */ virtual bool hasInitialAverage() const { return false; } /** * Return true, if this amplitude already includes symmetry factors * for identical outgoing particles. */ virtual bool hasFinalStateSymmetry() const { return false; } /** * Return true, if this amplitude is handled by a BLHA one-loop provider */ virtual bool isOLPTree() const { return false; } /** * Return true, if this amplitude is handled by a BLHA one-loop provider */ virtual bool isOLPLoop() const { return false; } /** * Return true, if colour and spin correlated matrix elements should * be ordered from the OLP */ virtual bool needsOLPCorrelators() const { return true; } /** * Write the order file header */ virtual void olpOrderFileHeader(ostream&) const; /** * Write the order file process list */ virtual void olpOrderFileProcesses(ostream&, const map,int>& procs) const; /** * Start the one loop provider, if appropriate, giving order and * contract files */ virtual void signOLP(const string&, const string&) { } /** * Start the one loop provider, if appropriate */ virtual void startOLP(const string&, int& status) { status = -1; } /** * Start the one loop provider, if appropriate. This default * implementation writes an BLHA 2.0 order file and starts the OLP */ virtual bool startOLP(const map,int>& procs); /** * Return true, if this amplitude needs to initialize an external * code. */ virtual bool isExternal() const { return false; } /** * Initialize this amplitude */ virtual bool initializeExternal() { return false; } /** * Return a generic process id for the given process */ virtual int externalId(const cPDVector&) { return 0; } /** * Return the map with masses to be used for amplitude evaluation */ const map& reshuffleMasses() const { return theReshuffleMasses; } /** * Check if reshuffling is needed at all */ void checkReshuffling(Ptr::tptr); /** * Return true, if this amplitude makes use of amplitudeMomenta */ virtual bool hasAmplitudeMomenta() const { return false; } //@} /** @name Colour basis. */ //@{ /** * Return the colour basis. */ virtual Ptr::tptr colourBasis() const { return theColourBasis; } /** * Return true, if the colour basis is capable of assigning colour * flows. */ virtual bool haveColourFlows() const { return colourBasis() ? colourBasis()->haveColourFlows() : false; } /** * Return a Selector with possible colour geometries for the selected * diagram weighted by their relative probabilities. */ virtual Selector colourGeometries(tcDiagPtr diag) const; /** * Return an ordering identifier for the current subprocess and * colour absis tensor index. */ const string& colourOrderingString(size_t id) const; /** * Return an ordering identifier for the current subprocess and * colour absis tensor index. */ const set >& colourOrdering(size_t id) const; //@} /** @name Phasespace point, crossing and helicities */ //@{ /** * Set the xcomb object. */ virtual void setXComb(tStdXCombPtr xc); /** * Return the momentum as crossed appropriate for this amplitude. */ Lorentz5Momentum amplitudeMomentum(int) const; /** * Perform a normal ordering of external legs and fill the * crossing information as. This default implementation sorts * lexicographically in (abs(colour)/spin/abs(charge)), putting pairs * of particles/anti-particles where possible. */ virtual void fillCrossingMap(size_t shift = 0); /** * Generate the helicity combinations. */ virtual set > generateHelicities() const; /** * Return the helicity combination of the physical process in the * conventions used by the spin correlation algorithm. */ virtual vector physicalHelicities(const vector&) const; //@} /** @name Tree-level amplitudes */ //@{ /** * Calculate the tree level amplitudes for the phasespace point * stored in lastXComb. */ virtual void prepareAmplitudes(Ptr::tcptr); /** * Return the matrix element squared. */ virtual double me2() const; /** * Return the colour charge of a given leg */ double colourCharge(tcPDPtr) const; /** * Return the large-N charge of a given leg */ double largeNColourCharge(tcPDPtr) const; /** * Return the largeN matrix element squared. */ virtual double largeNME2(Ptr::tptr largeNBasis) const; /** * Return the colour correlated matrix element. */ virtual double colourCorrelatedME2(pair ij) const; /** * Return the large-N colour correlated matrix element. */ virtual double largeNColourCorrelatedME2(pair ij, Ptr::tptr largeNBasis) const; /** * Return true if trivial colour configuration. */ bool trivialColourLegs() const { return theTrivialColourLegs; } /** * Return true, if this amplitude is capable of consistently filling * the rho matrices for the spin correllations */ virtual bool canFillRhoMatrix() const { return false; } /** * Return a positive helicity polarization vector for a gluon of * momentum p (with reference vector n) to be used when evaluating * spin correlations. */ virtual LorentzVector plusPolarization(const Lorentz5Momentum& p, const Lorentz5Momentum& n, int id = -1) const; /** * Return the colour and spin correlated matrix element. */ virtual double spinColourCorrelatedME2(pair emitterSpectator, const SpinCorrelationTensor& c) const; /** * Return the spin correlated matrix element. */ virtual double spinCorrelatedME2(pair emitterSpectator, const SpinCorrelationTensor& c) const; /** * Return true, if tree-level contributions will be evaluated at amplitude level. */ virtual bool treeAmplitudes() const { return true; } /** * Evaluate the amplitude for the given colour tensor id and * helicity assignment */ virtual Complex evaluate(size_t, const vector&, Complex&) { return 0.; } //@} /** @name One-loop amplitudes */ //@{ /** * Return the one-loop amplitude, if applicable. */ virtual Ptr::tptr oneLoopAmplitude() const { return Ptr::tptr(); } /** * Diasble one-loop functionality if not needed. */ virtual void disableOneLoop() {} /** * Return true, if this amplitude is capable of calculating one-loop * (QCD) corrections. */ virtual bool haveOneLoop() const { return false; } /** * Return true, if this amplitude only provides * one-loop (QCD) corrections. */ virtual bool onlyOneLoop() const { return false; } /** * Return true, if one-loop contributions will be evaluated at amplitude level. */ virtual bool oneLoopAmplitudes() const { return true; } /** * Return true, if one loop corrections have been calculated in * dimensional reduction. Otherwise conventional dimensional * regularization is assumed. Note that renormalization is always * assumed to be MSbar. */ virtual bool isDR() const { return false; } /** * Return true, if the amplitude is DRbar renormalized, otherwise * MSbar is assumed. */ virtual bool isDRbar() const { return true; } /** * Return true, if one loop corrections are given in the conventions * of the integrated dipoles. */ virtual bool isCS() const { return false; } /** * Return true, if one loop corrections are given in the conventions * of BDK. */ virtual bool isBDK() const { return false; } /** * Return true, if one loop corrections are given in the conventions * of everything expanded. */ virtual bool isExpanded() const { return false; } /** * Return the value of the dimensional regularization * parameter. Note that renormalization scale dependence is fully * restored in DipoleIOperator. */ virtual Energy2 mu2() const { return 0.*GeV2; } /** * Indicate that this amplitude is running alphas by itself. */ virtual bool hasRunningAlphaS() const { return false; } /** * Indicate that this amplitude is running alphaew by itself. */ virtual bool hasRunningAlphaEW() const { return false; } /** * If defined, return the coefficient of the pole in epsilon^2 */ virtual double oneLoopDoublePole() const { return 0.; } /** * If defined, return the coefficient of the pole in epsilon */ virtual double oneLoopSinglePole() const { return 0.; } /** * Calculate the one-loop amplitudes for the phasespace point * stored in lastXComb, if provided. */ virtual void prepareOneLoopAmplitudes(Ptr::tcptr); /** * Return the one-loop/tree interference. */ virtual double oneLoopInterference() const; /** * Evaluate the amplitude for the given colour tensor id and * helicity assignment */ virtual Complex evaluateOneLoop(size_t, const vector&) { return 0.; } //@} /** @name Caching and helpers to setup amplitude objects. */ //@{ /** * Flush all cashes. */ virtual void flushCaches() {} /** * Clone this amplitude. */ Ptr::ptr cloneMe() const { return dynamic_ptr_cast::ptr>(clone()); } /** * Clone the dependencies, using a given prefix. */ virtual void cloneDependencies(const std::string& prefix="" , bool slim=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(); // 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(); //@} private: /** * Recursively generate helicities */ void doGenerateHelicities(set >& res, vector& current, size_t pos) const; /** * The colour basis implementation to be used. */ Ptr::ptr theColourBasis; /** * The number of points after which helicity combinations wil be * cleaned up */ int theCleanupAfter; /** * The number of points that are calculated before a certain * helicity is excluded. Needed in pp->V */ int treeLevelHelicityPoints; /** * The number of points that are calculated before a certain * helicity is excluded. Needed in pp->V */ int oneLoopHelicityPoints; /** * The map with masses to be used for amplitude evaluation */ map theReshuffleMasses; /** * True if trivial colour configuration. */ bool theTrivialColourLegs; /** * A command to fill the reshuffle mass map */ string doReshuffle(string); /** * A command to fill the reshuffle mass map */ string doMassless(string); /** * A command to fill the reshuffle mass map */ string doOnShell(string); /** * Clear the reshuffling map */ string doClearReshuffling(string); /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MatchboxAmplitude & operator=(const MatchboxAmplitude &) = delete; }; inline PersistentOStream& operator<<(PersistentOStream& os, const Process& h) { h.persistentOutput(os); return os; } inline PersistentIStream& operator>>(PersistentIStream& is, Process& h) { h.persistentInput(is); return is; } } #endif /* HERWIG_MatchboxAmplitude_H */ diff --git a/MatrixElement/Matchbox/Base/MatchboxHybridAmplitude.cc b/MatrixElement/Matchbox/Base/MatchboxHybridAmplitude.cc --- a/MatrixElement/Matchbox/Base/MatchboxHybridAmplitude.cc +++ b/MatrixElement/Matchbox/Base/MatchboxHybridAmplitude.cc @@ -1,238 +1,236 @@ // -*- C++ -*- // // MatchboxHybridAmplitude.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 MatchboxHybridAmplitude class. // #include "MatchboxHybridAmplitude.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h" using namespace Herwig; MatchboxHybridAmplitude::MatchboxHybridAmplitude() : theUseOLPCorrelators(false) {} -MatchboxHybridAmplitude::~MatchboxHybridAmplitude() {} - IBPtr MatchboxHybridAmplitude::clone() const { return new_ptr(*this); } IBPtr MatchboxHybridAmplitude::fullclone() const { return new_ptr(*this); } bool MatchboxHybridAmplitude::isConsistent() const { assert(oneLoopAmplitude()); return !treeLevelAmplitude()->isOLPTree() && !treeLevelAmplitude()->isOLPLoop() && oneLoopAmplitude()->haveOneLoop() && treeLevelAmplitude()->orderInGs() == oneLoopAmplitude()->orderInGs() && treeLevelAmplitude()->orderInGem() == oneLoopAmplitude()->orderInGem() && treeLevelAmplitude()->hasRunningAlphaS() == oneLoopAmplitude()->hasRunningAlphaS() && treeLevelAmplitude()->hasRunningAlphaEW() == oneLoopAmplitude()->hasRunningAlphaEW() && !(treeLevelAmplitude()->nDimAdditional() != 0 && oneLoopAmplitude()->nDimAdditional() != 0); } bool MatchboxHybridAmplitude::canHandle(const PDVector& p, Ptr::tptr f, bool virt) const { if ( !virt ) return treeLevelAmplitude()->canHandle(p,f,false); if ( treeLevelAmplitude()->canHandle(p,f,false) && oneLoopAmplitude()->canHandle(p,f,true) ) { if ( !isConsistent() ) { generator()->log() << "Warning: Inconsistent settings encountered for MatchboxHybridAmplitude '" << name() << "'\n" << flush; return false; } return true; } return false; } void MatchboxHybridAmplitude::prepareAmplitudes(Ptr::tcptr me) { treeLevelAmplitude()->prepareAmplitudes(me); } void MatchboxHybridAmplitude::prepareOneLoopAmplitudes(Ptr::tcptr me) { assert(oneLoopAmplitude()); oneLoopAmplitude()->prepareOneLoopAmplitudes(me); } double MatchboxHybridAmplitude::symmetryRatio() const { assert(oneLoopAmplitude()); double ifact = 1.; if ( treeLevelAmplitude()->hasInitialAverage() && !oneLoopAmplitude()->hasInitialAverage() ) { ifact = 1./4.; if (lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour3 || lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour3bar ) ifact /= SM().Nc(); else if ( lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour8 ) ifact /= (SM().Nc()*SM().Nc()-1.); if ( lastMatchboxXComb()->matchboxME()->mePartonData()[1]->iColour() == PDT::Colour3 || lastMatchboxXComb()->matchboxME()->mePartonData()[1]->iColour() == PDT::Colour3bar ) ifact /= SM().Nc(); else if ( mePartonData()[1]->iColour() == PDT::Colour8 ) ifact /= (SM().Nc()*SM().Nc()-1.); } if ( !treeLevelAmplitude()->hasInitialAverage() && oneLoopAmplitude()->hasInitialAverage() ) { ifact = 4.; if ( lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour3 || lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour3bar ) ifact *= SM().Nc(); else if ( lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour8 ) ifact *= (SM().Nc()*SM().Nc()-1.); if ( lastMatchboxXComb()->matchboxME()->mePartonData()[1]->iColour() == PDT::Colour3 || lastMatchboxXComb()->matchboxME()->mePartonData()[1]->iColour() == PDT::Colour3bar ) ifact *= SM().Nc(); else if ( lastMatchboxXComb()->matchboxME()->mePartonData()[1]->iColour() == PDT::Colour8 ) ifact *= (SM().Nc()*SM().Nc()-1.); } if ( treeLevelAmplitude()->hasFinalStateSymmetry() && !oneLoopAmplitude()->hasFinalStateSymmetry() ) { assert(lastMatchboxXComb()->matchboxME()); ifact *= lastMatchboxXComb()->matchboxME()->finalStateSymmetry(); } if ( !treeLevelAmplitude()->hasFinalStateSymmetry() && oneLoopAmplitude()->hasFinalStateSymmetry() ) { assert(lastMatchboxXComb()->matchboxME()); ifact /= lastMatchboxXComb()->matchboxME()->finalStateSymmetry(); } return ifact; } void MatchboxHybridAmplitude::cloneDependencies(const std::string& prefix,bool slim) { if ( treeLevelAmplitude() ) { Ptr::ptr myTreeLevelAmplitude = treeLevelAmplitude()->cloneMe(); ostringstream pname; pname << (prefix == "" ? fullName() : prefix) << "/" << myTreeLevelAmplitude->name(); if ( ! (generator()->preinitRegister(myTreeLevelAmplitude,pname.str()) ) ) throw Exception() << "MatchboxHybridAmplitude::cloneDependencies(): Amplitude " << pname.str() << " already existing." << Exception::runerror; myTreeLevelAmplitude->cloneDependencies(pname.str()); treeLevelAmplitude(myTreeLevelAmplitude); } if ( oneLoopAmplitude() ) { Ptr::ptr myOneLoopAmplitude = oneLoopAmplitude()->cloneMe(); ostringstream pname; pname << (prefix == "" ? fullName() : prefix) << "/" << myOneLoopAmplitude->name(); if ( ! (generator()->preinitRegister(myOneLoopAmplitude,pname.str()) ) ) throw Exception() << "MatchboxHybridAmplitude::cloneDependencies(): Amplitude " << pname.str() << " already existing." << Exception::runerror; myOneLoopAmplitude->cloneDependencies(pname.str()); oneLoopAmplitude(myOneLoopAmplitude); } MatchboxAmplitude::cloneDependencies(prefix,slim); } void MatchboxHybridAmplitude::doinit() { MatchboxAmplitude::doinit(); if ( treeLevelAmplitude() ) treeLevelAmplitude()->init(); if ( oneLoopAmplitude() ) oneLoopAmplitude()->init(); } void MatchboxHybridAmplitude::doinitrun() { MatchboxAmplitude::doinitrun(); if ( treeLevelAmplitude() ) treeLevelAmplitude()->initrun(); if ( oneLoopAmplitude() ) oneLoopAmplitude()->initrun(); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void MatchboxHybridAmplitude::persistentOutput(PersistentOStream & os) const { os << theTreeLevelAmplitude << theOneLoopAmplitude << theUseOLPCorrelators; } void MatchboxHybridAmplitude::persistentInput(PersistentIStream & is, int) { is >> theTreeLevelAmplitude >> theOneLoopAmplitude >> theUseOLPCorrelators; } // *** 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 describeHerwigMatchboxHybridAmplitude("Herwig::MatchboxHybridAmplitude", "Herwig.so"); void MatchboxHybridAmplitude::Init() { static ClassDocumentation documentation ("MatchboxHybridAmplitude unifies two amplitude objects to " "provide tree and one-loop matrix elements."); static Reference interfaceTreeLevelAmplitude ("TreeLevelAmplitude", "Set the tree level amplitude to be used.", &MatchboxHybridAmplitude::theTreeLevelAmplitude, false, false, true, false, false); static Reference interfaceOneLoopAmplitude ("OneLoopAmplitude", "Set the one-loop amplitude to be used.", &MatchboxHybridAmplitude::theOneLoopAmplitude, false, false, true, true, false); static Switch interfaceUseOLPCorrelators ("UseOLPCorrelators", "Obtain correlated matrix elements from the OLP instead of " "the tree-level amplitude.", &MatchboxHybridAmplitude::theUseOLPCorrelators, false, false, false); static SwitchOption interfaceUseOLPCorrelatorsYes (interfaceUseOLPCorrelators, "Yes", "", true); static SwitchOption interfaceUseOLPCorrelatorsNo (interfaceUseOLPCorrelators, "No", "", false); interfaceUseOLPCorrelators.rank(-1); } diff --git a/MatrixElement/Matchbox/Base/MatchboxHybridAmplitude.h b/MatrixElement/Matchbox/Base/MatchboxHybridAmplitude.h --- a/MatrixElement/Matchbox/Base/MatchboxHybridAmplitude.h +++ b/MatrixElement/Matchbox/Base/MatchboxHybridAmplitude.h @@ -1,637 +1,629 @@ // -*- C++ -*- // // MatchboxHybridAmplitude.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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_MatchboxHybridAmplitude_H #define Herwig_MatchboxHybridAmplitude_H // // This is the declaration of the MatchboxHybridAmplitude class. // #include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief MatchboxHybridAmplitude unifies two amplitude objects to * provide tree and one-loop matrix elements. * * @see \ref MatchboxHybridAmplitudeInterfaces "The interfaces" * defined for MatchboxHybridAmplitude. */ class MatchboxHybridAmplitude: public Herwig::MatchboxAmplitude { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ MatchboxHybridAmplitude(); - /** - * The destructor. - */ - virtual ~MatchboxHybridAmplitude(); - //@} - public: /** * Return the amplitude object to provide tree-level amplitudes. */ Ptr::tptr treeLevelAmplitude() const { return theTreeLevelAmplitude; } /** * Set the amplitude object to provide tree-level amplitudes. */ void treeLevelAmplitude(Ptr::tptr amp) { theTreeLevelAmplitude = amp; } /** * Return the amplitude object to provide one-loop amplitudes. */ virtual Ptr::tptr oneLoopAmplitude() const { return theOneLoopAmplitude; } /** * Set the amplitude object to provide one-loop amplitudes. */ void oneLoopAmplitude(Ptr::tptr amp) { theOneLoopAmplitude = amp; } /** * Return true, if the two amplitude objects can be used in a * consistent way. */ bool isConsistent() const; public: /** @name Subprocess information */ //@{ /** * Return true, if this amplitude can handle the given process. */ virtual bool canHandle(const PDVector& p, Ptr::tptr f, bool virt) const; /** * Return the number of random numbers required to evaluate this * amplitude at a fixed phase space point. */ virtual int nDimAdditional() const { if ( !oneLoopAmplitude() ) return treeLevelAmplitude()->nDimAdditional(); return treeLevelAmplitude()->nDimAdditional() ? treeLevelAmplitude()->nDimAdditional() : oneLoopAmplitude()->nDimAdditional(); } /** * Return a ME instance appropriate for this amplitude and the given * subprocesses */ virtual Ptr::ptr makeME(const PDVector& p) const { return treeLevelAmplitude()->makeME(p); } /** * Set the (tree-level) order in \f$g_S\f$ in which this matrix * element should be evaluated. */ virtual void orderInGs(unsigned int n) { treeLevelAmplitude()->orderInGs(n); if ( oneLoopAmplitude() ) oneLoopAmplitude()->orderInGs(n); } /** * Return the (tree-level) order in \f$g_S\f$ in which this matrix * element is given. */ virtual unsigned int orderInGs() const { return treeLevelAmplitude()->orderInGs(); } /** * Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix * element should be evaluated. */ virtual void orderInGem(unsigned int n) { treeLevelAmplitude()->orderInGem(n); if ( oneLoopAmplitude() ) oneLoopAmplitude()->orderInGem(n); } /** * Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix * element is given. */ virtual unsigned int orderInGem() const { return treeLevelAmplitude()->orderInGem(); } /** * Return true, if this amplitude already includes averaging over * incoming parton's quantum numbers. */ virtual bool hasInitialAverage() const { return treeLevelAmplitude()->hasInitialAverage(); } /** * Return true, if this amplitude already includes symmetry factors * for identical outgoing particles. */ virtual bool hasFinalStateSymmetry() const { return treeLevelAmplitude()->hasFinalStateSymmetry(); } /** * Return true, if this amplitude is handled by a BLHA one-loop provider */ virtual bool isOLPTree() const { return false; } /** * Return true, if this amplitude is handled by a BLHA one-loop provider */ virtual bool isOLPLoop() const { if ( !oneLoopAmplitude() ) return false; return oneLoopAmplitude()->isOLPLoop(); } /** * Return true, if colour and spin correlated matrix elements should * be ordered from the OLP */ virtual bool needsOLPCorrelators() const { return theUseOLPCorrelators; } /** * Start the one loop provider, if appropriate. This default * implementation writes an BLHA 2.0 order file and starts the OLP */ virtual bool startOLP(const map,int>& procs) { assert(oneLoopAmplitude()); return oneLoopAmplitude()->startOLP(procs); } /** * Return true, if this amplitude needs to initialize an external * code. */ virtual bool isExternal() const { return treeLevelAmplitude()->isExternal(); } /** * Initialize this amplitude */ virtual bool initializeExternal() { return treeLevelAmplitude()->initializeExternal(); } /** * Return a generic process id for the given process */ virtual int externalId(const cPDVector& proc) { return treeLevelAmplitude()->externalId(proc); } //@} /** @name Colour basis. */ //@{ /** * Return the colour basis. */ virtual Ptr::tptr colourBasis() const { return treeLevelAmplitude()->colourBasis(); } /** * Return true, if the colour basis is capable of assigning colour * flows. */ virtual bool haveColourFlows() const { return treeLevelAmplitude()->haveColourFlows(); } /** * Return a Selector with possible colour geometries for the selected * diagram weighted by their relative probabilities. */ virtual Selector colourGeometries(tcDiagPtr diag) const { return treeLevelAmplitude()->colourGeometries(diag); } //@} /** @name Phasespace point, crossing and helicities */ //@{ /** * Set the xcomb object. */ virtual void setXComb(tStdXCombPtr xc) { treeLevelAmplitude()->setXComb(xc); if ( oneLoopAmplitude() ) oneLoopAmplitude()->setXComb(xc); lastMatchboxXComb(xc); } /** * Perform a normal ordering of external legs and fill the * crossing information as. This default implementation sorts * lexicographically in (abs(colour)/spin/abs(charge)), putting pairs * of particles/anti-particles where possible. */ virtual void fillCrossingMap(size_t shift = 0) { treeLevelAmplitude()->fillCrossingMap(shift); } //@} /** @name Tree-level amplitudes */ //@{ /** * Calculate the tree level amplitudes for the phasespace point * stored in lastXComb. */ virtual void prepareAmplitudes(Ptr::tcptr); /** * Return the matrix element squared. */ virtual double me2() const { return treeLevelAmplitude()->me2(); } /** * Return the largeN matrix element squared. */ virtual double largeNME2(Ptr::tptr largeNBasis) const { return treeLevelAmplitude()->largeNME2(largeNBasis); } /** * Return the colour correlated matrix element. */ virtual double colourCorrelatedME2(pair ij) const { return theUseOLPCorrelators ? oneLoopAmplitude()->colourCorrelatedME2(ij) : treeLevelAmplitude()->colourCorrelatedME2(ij); } /** * Return the large-N colour correlated matrix element. */ virtual double largeNColourCorrelatedME2(pair ij, Ptr::tptr largeNBasis) const { return treeLevelAmplitude()->largeNColourCorrelatedME2(ij,largeNBasis); } /** * Return a positive helicity polarization vector for a gluon of * momentum p (with reference vector n) to be used when evaluating * spin correlations. */ virtual LorentzVector plusPolarization(const Lorentz5Momentum& p, const Lorentz5Momentum& n, int id = -1) const { return theUseOLPCorrelators ? oneLoopAmplitude()->plusPolarization(p,n,id) : treeLevelAmplitude()->plusPolarization(p,n,id); } /** * Return the colour and spin correlated matrix element. */ virtual double spinColourCorrelatedME2(pair emitterSpectator, const SpinCorrelationTensor& c) const { return theUseOLPCorrelators ? oneLoopAmplitude()->spinColourCorrelatedME2(emitterSpectator,c) : treeLevelAmplitude()->spinColourCorrelatedME2(emitterSpectator,c); } /** * Return the spin correlated matrix element. */ virtual double spinCorrelatedME2(pair emitterSpectator, const SpinCorrelationTensor& c) const { return theUseOLPCorrelators ? oneLoopAmplitude()->spinCorrelatedME2(emitterSpectator,c) : treeLevelAmplitude()->spinCorrelatedME2(emitterSpectator,c); } /** * Return true, if this amplitude is capable of consistently filling * the rho matrices for the spin correllations */ virtual bool canFillRhoMatrix() const { return treeLevelAmplitude()->canFillRhoMatrix(); } /** * Return the helicity combination of the physical process in the * conventions used by the spin correlation algorithm. */ virtual vector physicalHelicities(const vector& hel) const { return treeLevelAmplitude()->physicalHelicities(hel); } /** * Return true, if tree-level contributions will be evaluated at amplitude level. */ virtual bool treeAmplitudes() const { return treeLevelAmplitude()->treeAmplitudes(); } /** * Evaluate the amplitude for the given colour tensor id and * helicity assignment */ virtual Complex evaluate(size_t a, const vector& hel, Complex& largeN) { return treeLevelAmplitude()->evaluate(a,hel,largeN); } //@} /** @name One-loop amplitudes */ //@{ /** * Diasble one-loop functionality if not needed. */ virtual void disableOneLoop() { oneLoopAmplitude(Ptr::ptr()); } /** * Return true, if this amplitude is capable of calculating one-loop * (QCD) corrections. */ virtual bool haveOneLoop() const { return oneLoopAmplitude(); } /** * Return true, if this amplitude only provides * one-loop (QCD) corrections. */ virtual bool onlyOneLoop() const { return false; } /** * Return true, if one-loop contributions will be evaluated at amplitude level. */ virtual bool oneLoopAmplitudes() const { assert(oneLoopAmplitude()); return oneLoopAmplitude()->oneLoopAmplitudes(); } /** * Return true, if the amplitude is DRbar renormalized, otherwise * MSbar is assumed. */ virtual bool isDRbar() const { assert(oneLoopAmplitude()); return oneLoopAmplitude()->isDRbar(); } /** * Return true, if one loop corrections have been calculated in * dimensional reduction. Otherwise conventional dimensional * regularization is assumed. Note that renormalization is always * assumed to be MSbar. */ virtual bool isDR() const { assert(oneLoopAmplitude()); return oneLoopAmplitude()->isDR(); } /** * Return true, if one loop corrections are given in the conventions * of the integrated dipoles. */ virtual bool isCS() const { assert(oneLoopAmplitude()); return oneLoopAmplitude()->isCS(); } /** * Return true, if one loop corrections are given in the conventions * of BDK. */ virtual bool isBDK() const { assert(oneLoopAmplitude()); return oneLoopAmplitude()->isBDK(); } /** * Return true, if one loop corrections are given in the conventions * of everything expanded. */ virtual bool isExpanded() const { assert(oneLoopAmplitude()); return oneLoopAmplitude()->isExpanded(); } /** * Return the value of the dimensional regularization * parameter. Note that renormalization scale dependence is fully * restored in DipoleIOperator. */ virtual Energy2 mu2() const { assert(oneLoopAmplitude()); return oneLoopAmplitude()->mu2(); } /** * Adjust the virtual symmetry factor conventions to the tree level * one */ double symmetryRatio() const; /** * Indicate that this amplitude is running alphas by itself. */ virtual bool hasRunningAlphaS() const { return treeLevelAmplitude()->hasRunningAlphaS(); } /** * Indicate that this amplitude is running alphaew by itself. */ virtual bool hasRunningAlphaEW() const { return treeLevelAmplitude()->hasRunningAlphaEW(); } /** * If defined, return the coefficient of the pole in epsilon^2 */ virtual double oneLoopDoublePole() const { assert(oneLoopAmplitude()); return symmetryRatio()*oneLoopAmplitude()->oneLoopDoublePole(); } /** * If defined, return the coefficient of the pole in epsilon */ virtual double oneLoopSinglePole() const { assert(oneLoopAmplitude()); return symmetryRatio()*oneLoopAmplitude()->oneLoopSinglePole(); } /** * Calculate the one-loop amplitudes for the phasespace point * stored in lastXComb, if provided. */ virtual void prepareOneLoopAmplitudes(Ptr::tcptr); /** * Return the one-loop/tree interference. */ virtual double oneLoopInterference() const { assert(oneLoopAmplitude()); return symmetryRatio()*oneLoopAmplitude()->oneLoopInterference(); } /** * Evaluate the amplitude for the given colour tensor id and * helicity assignment */ virtual Complex evaluateOneLoop(size_t a, const vector& hel) { assert(oneLoopAmplitude()); return oneLoopAmplitude()->evaluateOneLoop(a,hel); } //@} /** @name Caching and helpers to setup amplitude objects. */ //@{ /** * Flush all cashes. */ virtual void flushCaches() { treeLevelAmplitude()->flushCaches(); if ( oneLoopAmplitude() ) oneLoopAmplitude()->flushCaches(); } /** * Clone the dependencies, using a given prefix. */ virtual void cloneDependencies(const std::string& prefix= "" , bool slim=false); //@} 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(); //@} 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 amplitude object to provide tree-level amplitudes. */ Ptr::ptr theTreeLevelAmplitude; /** * The amplitude object to provide one-loop amplitudes. */ Ptr::ptr theOneLoopAmplitude; /** * True, if correlators should be used from the OLP amplitude */ bool theUseOLPCorrelators; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MatchboxHybridAmplitude & operator=(const MatchboxHybridAmplitude &) = delete; }; } #endif /* Herwig_MatchboxHybridAmplitude_H */ diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc --- a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc +++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc @@ -1,1680 +1,1678 @@ // -*- C++ -*- // // MatchboxMEBase.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 MatchboxMEBase class. // #include "MatchboxMEBase.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDF/PDF.h" #include "ThePEG/PDT/PDT.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/Handlers/StdXCombGroup.h" #include "ThePEG/EventRecord/SubProcess.h" #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h" #include "Herwig/MatrixElement/Matchbox/Utility/DiagramDrawer.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h" #include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h" #include "Herwig/API/RunDirectories.h" #include "Herwig/MatrixElement/ProductionMatrixElement.h" #include "Herwig/MatrixElement/HardVertex.h" #include #include using std::ostream_iterator; using namespace Herwig; MatchboxMEBase::MatchboxMEBase() : MEBase(), theOneLoop(false), theOneLoopNoBorn(false), theOneLoopNoLoops(false), theNoCorrelations(false), theHavePDFs(false,false), checkedPDFs(false) {} -MatchboxMEBase::~MatchboxMEBase() {} - Ptr::tptr MatchboxMEBase::factory() const { return MatchboxFactory::currentFactory(); } Ptr::tptr MatchboxMEBase::diagramGenerator() const { return factory()->diagramGenerator(); } Ptr::tptr MatchboxMEBase::processData() const { return factory()->processData(); } unsigned int MatchboxMEBase::getNLight() const { return factory()->nLight(); } vector MatchboxMEBase::getNLightJetVec() const { return factory()->nLightJetVec(); } vector MatchboxMEBase::getNHeavyJetVec() const { return factory()->nHeavyJetVec(); } vector MatchboxMEBase::getNLightProtonVec() const { return factory()->nLightProtonVec(); } double MatchboxMEBase::factorizationScaleFactor() const { return factory()->factorizationScaleFactor(); } double MatchboxMEBase::renormalizationScaleFactor() const { return factory()->renormalizationScaleFactor(); } bool MatchboxMEBase::fixedCouplings() const { return factory()->fixedCouplings(); } bool MatchboxMEBase::fixedQEDCouplings() const { return factory()->fixedQEDCouplings(); } bool MatchboxMEBase::checkPoles() const { return factory()->checkPoles(); } bool MatchboxMEBase::verbose() const { return factory()->verbose(); } bool MatchboxMEBase::initVerbose() const { return factory()->initVerbose(); } void MatchboxMEBase::getDiagrams() const { if ( diagramGenerator() && processData() ) { vector::ptr> diags; vector::ptr>& res = processData()->diagramMap()[subProcess().legs]; if ( res.empty() ) { res = diagramGenerator()->generate(subProcess().legs,orderInAlphaS(),orderInAlphaEW()); } copy(res.begin(),res.end(),back_inserter(diags)); processData()->fillMassGenerators(subProcess().legs); if ( diags.empty() ) return; for (auto const & d : diags ) add(d); return; } throw Exception() << "MatchboxMEBase::getDiagrams() expects a Tree2toNGenerator and ProcessData object.\n" << "Please check your setup." << Exception::runerror; } Selector MatchboxMEBase::diagrams(const DiagramVector & diags) const { if ( phasespace() ) { return phasespace()->selectDiagrams(diags); } throw Exception() << "MatchboxMEBase::diagrams() expects a MatchboxPhasespace object.\n" << "Please check your setup." << Exception::runerror; return Selector(); } Selector MatchboxMEBase::colourGeometries(tcDiagPtr diag) const { if ( matchboxAmplitude() ) { if ( matchboxAmplitude()->haveColourFlows() ) { if ( matchboxAmplitude()->treeAmplitudes() ) matchboxAmplitude()->prepareAmplitudes(this); return matchboxAmplitude()->colourGeometries(diag); } } Ptr::tcptr tdiag = dynamic_ptr_cast::tcptr>(diag); assert(diag && processData()); vector& flows = processData()->colourFlowMap()[tdiag]; if ( flows.empty() ) { list > > > cflows = ColourBasis::colourFlows(tdiag); for ( auto const & fit : cflows) flows.push_back(new ColourLines(ColourBasis::cfstring(fit))); } Selector res; for ( auto const & f : flows ) res.insert(1.0,f); return res; } void MatchboxMEBase::constructVertex(tSubProPtr sub, const ColourLines* cl) { if ( !canFillRhoMatrix() || !factory()->spinCorrelations() ) return; assert(matchboxAmplitude()); assert(matchboxAmplitude()->colourBasis()); // get the colour structure for the selected colour flow size_t cStructure = matchboxAmplitude()->colourBasis()->tensorIdFromFlow(lastXComb().lastDiagram(),cl); // hard process for processing the spin info tPVector hard; hard.push_back(sub->incoming().first); hard.push_back(sub->incoming().second); vector out; for ( auto const & p : sub->outgoing() ) { out.push_back(p->data().iSpin()); hard.push_back(p); } // calculate dummy wave functions to fill the spin info static vector dummyPolarizations; static vector dummySpinors; static vector dummyBarSpinors; for ( size_t k = 0; k < hard.size(); ++k ) { if ( hard[k]->data().iSpin() == PDT::Spin1Half ) { if ( hard[k]->id() > 0 && k > 1 ) { SpinorBarWaveFunction(dummyBarSpinors,hard[k], outgoing, true); } else if ( hard[k]->id() < 0 && k > 1 ) { SpinorWaveFunction(dummySpinors,hard[k], outgoing, true); } else if ( hard[k]->id() > 0 && k < 2 ) { SpinorWaveFunction(dummySpinors,hard[k], incoming, false); } else if ( hard[k]->id() < 0 && k < 2 ) { SpinorBarWaveFunction(dummyBarSpinors,hard[k], incoming, false); } } else if ( hard[k]->data().iSpin() == PDT::Spin1 ) { VectorWaveFunction(dummyPolarizations,hard[k], k > 1 ? outgoing : incoming, k > 1 ? true : false, hard[k]->data().hardProcessMass() == ZERO); } else if (hard[k]->data().iSpin() == PDT::Spin0 ) { ScalarWaveFunction(hard[k],k > 1 ? outgoing : incoming, k > 1 ? true : false); } else assert(false); } // fill the production matrix element ProductionMatrixElement pMe(mePartonData()[0]->iSpin(), mePartonData()[1]->iSpin(), out); for ( map,CVector>::const_iterator lamp = lastLargeNAmplitudes().begin(); lamp != lastLargeNAmplitudes().end(); ++lamp ) { vector pMeHelicities = matchboxAmplitude()->physicalHelicities(lamp->first); pMe(pMeHelicities) = lamp->second[cStructure]; } // set the spin information HardVertexPtr hardvertex = new_ptr(HardVertex()); hardvertex->ME(pMe); if ( sub->incoming().first->spinInfo() ) sub->incoming().first->spinInfo()->productionVertex(hardvertex); if ( sub->incoming().second->spinInfo() ) sub->incoming().second->spinInfo()->productionVertex(hardvertex); for ( auto const & p : sub->outgoing() ) if ( p->spinInfo() ) p->spinInfo()->productionVertex(hardvertex); } unsigned int MatchboxMEBase::orderInAlphaS() const { return subProcess().orderInAlphaS; } unsigned int MatchboxMEBase::orderInAlphaEW() const { return subProcess().orderInAlphaEW; } void MatchboxMEBase::setXComb(tStdXCombPtr xc) { MEBase::setXComb(xc); lastMatchboxXComb(xc); if ( phasespace() ) phasespace()->setXComb(xc); if ( scaleChoice() ) scaleChoice()->setXComb(xc); if ( matchboxAmplitude() ) matchboxAmplitude()->setXComb(xc); if (theMerger){ theMerger->setME(this); theMerger->setXComb( xc ); } } double MatchboxMEBase::generateIncomingPartons(const double* r1, const double* r2) { // shamelessly stolen from PartonExtractor.cc Energy2 shmax = lastCuts().sHatMax(); Energy2 shmin = lastCuts().sHatMin(); Energy2 sh = shmin*pow(shmax/shmin, *r1); double ymax = lastCuts().yHatMax(); double ymin = lastCuts().yHatMin(); double km = log(shmax/shmin); ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh))); ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh))); double y = ymin + (*r2)*(ymax - ymin); double x1 = exp(-0.5*log(lastS()/sh) + y); double x2 = exp(-0.5*log(lastS()/sh) - y); Lorentz5Momentum P1 = lastParticles().first->momentum(); LorentzMomentum p1 = lightCone((P1.rho() + P1.e())*x1, Energy()); p1.rotateY(P1.theta()); p1.rotateZ(P1.phi()); meMomenta()[0] = p1; Lorentz5Momentum P2 = lastParticles().second->momentum(); LorentzMomentum p2 = lightCone((P2.rho() + P2.e())*x2, Energy()); p2.rotateY(P2.theta()); p2.rotateZ(P2.phi()); meMomenta()[1] = p2; lastXCombPtr()->lastX1X2(make_pair(x1,x2)); lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2()); return km*(ymax - ymin); } bool MatchboxMEBase::generateKinematics(const double * r) { if ( phasespace() ) { jacobian(phasespace()->generateKinematics(r,meMomenta())); if ( jacobian() == 0.0 ) return false; setScale(); if (theMerger&&!theMerger->generateKinematics(r)){ return false; } logGenerateKinematics(r); assert(lastMatchboxXComb()); if ( nDimAmplitude() > 0 ) { amplitudeRandomNumbers().resize(nDimAmplitude()); copy(r + nDimPhasespace(), r + nDimPhasespace() + nDimAmplitude(), amplitudeRandomNumbers().begin()); } if ( nDimInsertions() > 0 ) { insertionRandomNumbers().resize(nDimInsertions()); copy(r + nDimPhasespace() + nDimAmplitude(), r + nDimPhasespace() + nDimAmplitude() + nDimInsertions(), insertionRandomNumbers().begin()); } return true; } throw Exception() << "MatchboxMEBase::generateKinematics() expects a MatchboxPhasespace object.\n" << "Please check your setup." << Exception::runerror; return false; } int MatchboxMEBase::nDim() const { if ( lastMatchboxXComb() ) return nDimPhasespace() + nDimAmplitude() + nDimInsertions(); int ampAdd = 0; if ( matchboxAmplitude() ) { ampAdd = matchboxAmplitude()->nDimAdditional(); } int insertionAdd = 0; for ( auto const & v : virtuals() ) { insertionAdd = max(insertionAdd,v->nDimAdditional()); } return nDimBorn() + ampAdd + insertionAdd; } int MatchboxMEBase::nDimBorn() const { if ( lastMatchboxXComb() ) return nDimPhasespace(); if ( phasespace() ) return phasespace()->nDim(diagrams().front()->partons()); throw Exception() << "MatchboxMEBase::nDim() expects a MatchboxPhasespace object.\n" << "Please check your setup." << Exception::runerror; return 0; } void MatchboxMEBase::setScale(Energy2 ren, Energy2 fac) const { if ( haveX1X2() ) { lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2()); } Energy2 fcscale = (fac == ZERO) ? factorizationScale() : fac; Energy2 fscale = fcscale*sqr(factorizationScaleFactor()); Energy2 rscale = (ren == ZERO ? renormalizationScale() : ren)*sqr(renormalizationScaleFactor()); Energy2 ewrscale = renormalizationScaleQED(); lastXCombPtr()->lastScale(fscale); lastXCombPtr()->lastCentralScale(fcscale); lastXCombPtr()->lastShowerScale(showerScale()); lastMatchboxXComb()->lastRenormalizationScale(rscale); if ( !fixedCouplings() ) { if ( rscale > lastCuts().scaleMin() ) lastXCombPtr()->lastAlphaS(SM().alphaS(rscale)); else lastXCombPtr()->lastAlphaS(SM().alphaS(lastCuts().scaleMin())); } else { lastXCombPtr()->lastAlphaS(SM().alphaS()); } if ( !fixedQEDCouplings() ) { lastXCombPtr()->lastAlphaEM(SM().alphaEMME(ewrscale)); } else { lastXCombPtr()->lastAlphaEM(SM().alphaEMMZ()); } logSetScale(); } Energy2 MatchboxMEBase::factorizationScale() const { if ( scaleChoice() ) { return scaleChoice()->factorizationScale(); } throw Exception() << "MatchboxMEBase::factorizationScale() expects a MatchboxScaleChoice object.\n" << "Please check your setup." << Exception::runerror; return ZERO; } Energy2 MatchboxMEBase::renormalizationScale() const { if ( scaleChoice() ) { return scaleChoice()->renormalizationScale(); } throw Exception() << "MatchboxMEBase::renormalizationScale() expects a MatchboxScaleChoice object.\n" << "Please check your setup." << Exception::runerror; return ZERO; } Energy2 MatchboxMEBase::renormalizationScaleQED() const { if ( scaleChoice() ) { return scaleChoice()->renormalizationScaleQED(); } return renormalizationScale(); } Energy2 MatchboxMEBase::showerScale() const { if ( scaleChoice() ) { return scaleChoice()->showerScale(); } throw Exception() << "MatchboxMEBase::showerScale() expects a MatchboxScaleChoice object.\n" << "Please check your setup." << Exception::runerror; return ZERO; } void MatchboxMEBase::setVetoScales(tSubProPtr) const {} bool MatchboxMEBase::havePDFWeight1() const { if ( checkedPDFs ) return theHavePDFs.first; theHavePDFs.first = factory()->isIncoming(mePartonData()[0]) && lastXCombPtr()->partonBins().first->pdf(); theHavePDFs.second = factory()->isIncoming(mePartonData()[1]) && lastXCombPtr()->partonBins().second->pdf(); checkedPDFs = true; return theHavePDFs.first; } bool MatchboxMEBase::havePDFWeight2() const { if ( checkedPDFs ) return theHavePDFs.second; theHavePDFs.first = factory()->isIncoming(mePartonData()[0]) && lastXCombPtr()->partonBins().first->pdf(); theHavePDFs.second = factory()->isIncoming(mePartonData()[1]) && lastXCombPtr()->partonBins().second->pdf(); checkedPDFs = true; return theHavePDFs.second; } void MatchboxMEBase::getPDFWeight(Energy2 factorizationScale) const { if ( !havePDFWeight1() && !havePDFWeight2() ) { lastMEPDFWeight(1.0); logPDFWeight(); return; } double w = 1.; if ( havePDFWeight1() ) w *= pdf1(factorizationScale); if ( havePDFWeight2() ) w *= pdf2(factorizationScale); lastMEPDFWeight(w); logPDFWeight(); } double MatchboxMEBase::pdf1(Energy2 fscale, double xEx, double xFactor) const { assert(lastXCombPtr()->partonBins().first->pdf()); if ( xEx < 1. && lastX1()*xFactor >= xEx ) { return ( ( 1. - lastX1()*xFactor ) / ( 1. - xEx ) ) * lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(), lastPartons().first->dataPtr(), fscale == ZERO ? lastScale() : fscale, xEx)/xEx; } return lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(), lastPartons().first->dataPtr(), fscale == ZERO ? lastScale() : fscale, lastX1()*xFactor)/lastX1()/xFactor; } double MatchboxMEBase::pdf2(Energy2 fscale, double xEx, double xFactor) const { assert(lastXCombPtr()->partonBins().second->pdf()); if ( xEx < 1. && lastX2()*xFactor >= xEx ) { return ( ( 1. - lastX2()*xFactor ) / ( 1. - xEx ) ) * lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(), lastPartons().second->dataPtr(), fscale == ZERO ? lastScale() : fscale, xEx)/xEx; } return lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(), lastPartons().second->dataPtr(), fscale == ZERO ? lastScale() : fscale, lastX2()*xFactor)/lastX2()/xFactor; } double MatchboxMEBase::me2() const { if ( matchboxAmplitude() ) { if ( matchboxAmplitude()->treeAmplitudes() ) matchboxAmplitude()->prepareAmplitudes(this); double res = matchboxAmplitude()->me2()* me2Norm(); return res; } throw Exception() << "MatchboxMEBase::me2() expects a MatchboxAmplitude object.\n" << "Please check your setup." << Exception::runerror; return 0.; } double MatchboxMEBase::largeNME2(Ptr::tptr largeNBasis) const { if ( matchboxAmplitude() ) { if ( matchboxAmplitude()->treeAmplitudes() ) { largeNBasis->prepare(mePartonData(),false); matchboxAmplitude()->prepareAmplitudes(this); } double res = matchboxAmplitude()->largeNME2(largeNBasis)* me2Norm(); return res; } throw Exception() << "MatchboxMEBase::largeNME2() expects a MatchboxAmplitude object.\n" << "Please check your setup." << Exception::runerror; return 0.; } double MatchboxMEBase::finalStateSymmetry() const { if ( symmetryFactor() > 0.0 ) return symmetryFactor(); double sFactor = 1.; map counts; cPDVector checkData; copy(mePartonData().begin()+2,mePartonData().end(),back_inserter(checkData)); cPDVector::iterator p = checkData.begin(); while ( !checkData.empty() ) { if ( counts.find((**p).id()) != counts.end() ) { counts[(**p).id()] += 1; } else { counts[(**p).id()] = 1; } checkData.erase(p); p = checkData.begin(); continue; } for ( auto const & c : counts) { if ( c.second == 1 ) continue; if ( c.second == 2 ) sFactor /= 2.; else if ( c.second == 3 ) sFactor /= 6.; else if ( c.second == 4 ) sFactor /= 24.; } symmetryFactor(sFactor); return symmetryFactor(); } double MatchboxMEBase::me2Norm(unsigned int addAlphaS) const { double fac = 1.; if ( !hasInitialAverage() ) { for ( size_t k = 0; k < 2; ++k ) { // Spin 0 does not involve any additional factors if ( mePartonData()[k]->iSpin() == PDT::Spin1Half ) fac *= 1./2.; else if ( mePartonData()[k]->iSpin() == PDT::Spin1 && mePartonData()[k]->hardProcessMass() == ZERO ) fac *= 1./2.; else if ( mePartonData()[k]->iSpin() == PDT::Spin1 && mePartonData()[k]->hardProcessMass() > ZERO ) fac *= 1./3.; else if ( mePartonData()[k]->iSpin() == PDT::Spin3Half ) fac *= 1./4.; else if ( mePartonData()[k]->iSpin() == PDT::Spin2 && mePartonData()[k]->hardProcessMass() == ZERO ) fac *= 1./2.; else if ( mePartonData()[k]->iSpin() == PDT::Spin2 && mePartonData()[k]->hardProcessMass() > ZERO ) fac *= 1./5.; } } double couplings = 1.0; if ( (orderInAlphaS() > 0 || addAlphaS != 0) && !hasRunningAlphaS() ) { fac *= pow(lastAlphaS()/SM().alphaS(),double(orderInAlphaS()+addAlphaS)); couplings *= pow(lastAlphaS(),double(orderInAlphaS()+addAlphaS)); } if ( orderInAlphaEW() > 0 && !hasRunningAlphaEW() ) { fac *= pow(lastAlphaEM()/SM().alphaEMMZ(),double(orderInAlphaEW())); couplings *= pow(lastAlphaEM(),double(orderInAlphaEW())); } lastMECouplings(couplings); if ( !hasInitialAverage() ) { if ( mePartonData()[0]->iColour() == PDT::Colour3 || mePartonData()[0]->iColour() == PDT::Colour3bar ) fac /= SM().Nc(); else if ( mePartonData()[0]->iColour() == PDT::Colour8 ) fac /= (SM().Nc()*SM().Nc()-1.); if ( mePartonData()[1]->iColour() == PDT::Colour3 || mePartonData()[1]->iColour() == PDT::Colour3bar ) fac /= SM().Nc(); else if ( mePartonData()[1]->iColour() == PDT::Colour8 ) fac /= (SM().Nc()*SM().Nc()-1.); } return !hasFinalStateSymmetry() ? finalStateSymmetry()*fac : fac; } CrossSection MatchboxMEBase::prefactor()const{ return (sqr(hbarc)/(2.*lastSHat())) *jacobian()* lastMEPDFWeight(); } CrossSection MatchboxMEBase::dSigHatDRB() const { getPDFWeight(); lastME2(me2()); return oneLoopNoBorn()?ZERO:prefactor() * lastME2(); } CrossSection MatchboxMEBase::dSigHatDRV() const { getPDFWeight(); lastME2(me2()); return ( oneLoop() && !oneLoopNoLoops() )?(prefactor() * oneLoopInterference()):ZERO; } CrossSection MatchboxMEBase::dSigHatDRI() const { getPDFWeight(); lastME2(me2()); CrossSection res=ZERO; if (oneLoop() &&!onlyOneLoop()) { for ( auto const & v : virtuals()) { v->setXComb(lastXCombPtr()); res += v->dSigHatDR(); } if ( checkPoles() && oneLoop() ) logPoles(); } return res; } CrossSection MatchboxMEBase::dSigHatDRAlphaDiff(double alpha) const { getPDFWeight(); lastME2(me2()); CrossSection res=ZERO; for ( auto const & v: virtuals() ) { v->setXComb(lastXCombPtr()); res+=v->dSigHatDRAlphaDiff( alpha); } return res; } CrossSection MatchboxMEBase::dSigHatDR() const { getPDFWeight(); if (theMerger){ lastMECrossSection(theMerger->MergingDSigDR()); return lastMECrossSection(); } else if (lastXCombPtr()->willPassCuts() ) { lastME2(me2()); CrossSection _dSigHatDRB, _dSigHatDRV, _dSigHatDRI, res = ZERO; // ----- dSigHatDRB ----- _dSigHatDRB = oneLoopNoBorn()?ZERO:prefactor() * lastME2(); // ----- dSigHatDRV ----- _dSigHatDRV = ( oneLoop() && !oneLoopNoLoops() )?(prefactor() * oneLoopInterference()):ZERO; // ----- dSigHatDRI ----- if (oneLoop() &&!onlyOneLoop()) { for ( auto const & v : virtuals()) { v->setXComb(lastXCombPtr()); res += v->dSigHatDR(); } if ( checkPoles() && oneLoop() ) logPoles(); } _dSigHatDRI = res; // ----- finalizing ----- lastMECrossSection(_dSigHatDRB + _dSigHatDRV + _dSigHatDRI); return lastMECrossSection(); } else { lastME2(ZERO); lastMECrossSection(ZERO); return lastMECrossSection(); } } double MatchboxMEBase::oneLoopInterference() const { if ( matchboxAmplitude() ) { if ( matchboxAmplitude()->oneLoopAmplitudes() ) matchboxAmplitude()->prepareOneLoopAmplitudes(this); double res = matchboxAmplitude()->oneLoopInterference()* me2Norm(1); return res; } throw Exception() << "MatchboxMEBase::oneLoopInterference() expects a MatchboxAmplitude object.\n" << "Please check your setup." << Exception::runerror; return 0.; } MatchboxMEBase::AccuracyHistogram::AccuracyHistogram(double low, double up, unsigned int nbins) : lower(low), upper(up), sameSign(0), oppositeSign(0), nans(0), overflow(0), underflow(0) { double step = (up-low)/nbins; for ( unsigned int k = 1; k <= nbins; ++k ) bins[lower + k*step] = 0.0; } void MatchboxMEBase::AccuracyHistogram::book(double a, double b) { if ( ! (isfinite(a) && isfinite(b)) ) { ++nans; return; } if ( a*b >= 0. ) ++sameSign; if ( a*b < 0. ) ++oppositeSign; double r = 1.; if ( abs(a) != 0.0 ) r = abs(1.-abs(b/a)); else if ( abs(b) != 0.0 ) r = abs(b); if ( log10(r) < lower || r == 0.0 ) { ++underflow; return; } if ( log10(r) > upper ) { ++overflow; return; } map::iterator bin = bins.upper_bound(log10(r)); if ( bin == bins.end() ) return; bin->second += 1.; } void MatchboxMEBase::AccuracyHistogram::dump(const std::string& folder, const std::string& prefix, const cPDVector& proc) const { ostringstream fname(""); for ( cPDVector::const_iterator p = proc.begin(); p != proc.end(); ++p ) fname << (**p).PDGName(); ofstream out((folder+"/"+prefix+fname.str()+".dat").c_str()); out << "# same sign : " << sameSign << " opposite sign : " << oppositeSign << " nans : " << nans << " overflow : " << overflow << " underflow : " << underflow << "\n"; for ( map::const_iterator b = bins.begin(); b != bins.end(); ++b ) { map::const_iterator bp = b; --bp; if ( b->second != 0. ) { if ( b != bins.begin() ) out << bp->first; else out << lower; out << " " << b->first << " " << b->second << "\n" << flush; } } ofstream gpout((folder+"/"+prefix+fname.str()+".gp").c_str()); gpout << "set terminal png\n" << "set xlabel 'accuracy of pole cancellation [decimal places]'\n" << "set ylabel 'counts\n" << "set xrange [-20:0]\n" << "set output '" << prefix << fname.str() << ".png'\n" << "plot '" << prefix << fname.str() << ".dat' using (0.5*($1+$2)):3 with linespoints pt 7 ps 1 not"; } void MatchboxMEBase::AccuracyHistogram::persistentOutput(PersistentOStream& os) const { os << lower << upper << bins << sameSign << oppositeSign << nans << overflow << underflow; } void MatchboxMEBase::AccuracyHistogram::persistentInput(PersistentIStream& is) { is >> lower >> upper >> bins >> sameSign >> oppositeSign >> nans >> overflow >> underflow; } void MatchboxMEBase::logPoles() const { double res2me = oneLoopDoublePole(); double res1me = oneLoopSinglePole(); double res2i = 0.; double res1i = 0.; for ( auto const & v : virtuals()) { res2i += v->oneLoopDoublePole(); res1i += v->oneLoopSinglePole(); } if (res2me != 0.0 || res2i != 0.0) epsilonSquarePoleHistograms[mePartonData()].book(res2me,res2i); if (res1me != 0.0 || res1i != 0.0) epsilonPoleHistograms[mePartonData()].book(res1me,res1i); } bool MatchboxMEBase::haveOneLoop() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->haveOneLoop(); return false; } bool MatchboxMEBase::onlyOneLoop() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->onlyOneLoop(); return false; } bool MatchboxMEBase::isDRbar() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->isDRbar(); return false; } bool MatchboxMEBase::isDR() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->isDR(); return false; } bool MatchboxMEBase::isCS() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->isCS(); return false; } bool MatchboxMEBase::isBDK() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->isBDK(); return false; } bool MatchboxMEBase::isExpanded() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->isExpanded(); return false; } Energy2 MatchboxMEBase::mu2() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->mu2(); return 0*GeV2; } double MatchboxMEBase::oneLoopDoublePole() const { if ( matchboxAmplitude() ) { return matchboxAmplitude()->oneLoopDoublePole()* me2Norm(1); } return 0.; } double MatchboxMEBase::oneLoopSinglePole() const { if ( matchboxAmplitude() ) { return matchboxAmplitude()->oneLoopSinglePole()* me2Norm(1); } return 0.; } vector MatchboxMEBase::getDipoles(const vector& dipoles, const vector & borns,bool slim) const { vector res; // keep track of the dipoles we already did set up set,int>,pair::tptr,Ptr::tptr> > > done; cPDVector rep = diagrams().front()->partons(); int nreal = rep.size(); // now loop over configs for ( int emitter = 0; emitter < nreal; ++emitter ) { list matchDipoles; for ( auto const & d : dipoles ) { if ( ! d->canHandleEmitter(rep,emitter) ) continue; matchDipoles.push_back(d); } if ( matchDipoles.empty() ) continue; for ( int emission = 2; emission < nreal; ++emission ) { if ( emission == emitter ) continue; list matchDipoles2; for ( auto const & d : matchDipoles ) { if ( !d->canHandleSplitting(rep,emitter,emission) ) continue; matchDipoles2.push_back(d); } if ( matchDipoles2.empty() ) continue; map::ptr,SubtractionDipole::MergeInfo> mergeInfo; for ( auto const & d : diagrams() ) { Ptr::ptr check = new_ptr(Tree2toNDiagram(*dynamic_ptr_cast::ptr>(d))); map theMergeLegs; for ( unsigned int i = 0; i < check->external().size(); ++i ) theMergeLegs[i] = -1; int theEmitter = check->mergeEmission(emitter,emission,theMergeLegs); // no underlying Born if ( theEmitter == -1 ) continue; SubtractionDipole::MergeInfo info; info.diagram = check; info.emitter = theEmitter; info.mergeLegs = theMergeLegs; mergeInfo[d] = info; } if ( mergeInfo.empty() ) continue; for ( int spectator = 0; spectator < nreal; ++spectator ) { if ( spectator == emitter || spectator == emission ) continue; list matchDipoles3; for ( auto const & d : matchDipoles2 ) { if ( ! d->canHandleSpectator(rep,spectator) ) continue; matchDipoles3.push_back(d); } if ( matchDipoles3.empty() ) continue; if ( noDipole(emitter,emission,spectator) ) continue; for ( auto const & d : matchDipoles3 ) { if ( !d->canHandle(rep,emitter,emission,spectator) ) continue; for ( auto const & b : borns ) { if ( b->onlyOneLoop() ) continue; if ( done.find(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(b,d))) != done.end() ) continue; // now get to work d->clearBookkeeping(); d->realEmitter(emitter); d->realEmission(emission); d->realSpectator(spectator); d->realEmissionME(const_cast(this)); d->underlyingBornME(b); d->setupBookkeeping(mergeInfo,slim); if ( ! d->empty() ) { res.push_back( d->cloneMe() ); Ptr::tptr nDipole = res.back(); done.insert(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(b,d))); if ( nDipole->isSymmetric() ) done.insert(make_pair(make_pair(make_pair(emission,emitter),spectator),make_pair(b,d))); ostringstream dname; if ( theMerger) { dname << fullName(); if (theOneLoopNoBorn) dname << ".virtual" << "." ; dname << b->name() << "." << d->name() << ".[(" << emitter << "," << emission << ")," << spectator << "]"; } else { dname << fullName() << "." << b->name() << "." << d->name() << ".[(" << emitter << "," << emission << ")," << spectator << "]"; } if ( ! (generator()->preinitRegister(nDipole,dname.str()) ) ) throw Exception() << "MatchboxMEBase::getDipoles(): Dipole " << dname.str() << " already existing." << Exception::runerror; if ( !factory()->reweighters().empty() ) { for ( auto const & rw : factory()->reweighters()) nDipole->addReweighter(rw); } if ( !factory()->preweighters().empty() ) { for ( auto const & rw : factory()->preweighters() ) nDipole->addPreweighter(rw); } nDipole->cloneDependencies(dname.str(),slim); } } } } } } vector::tptr> partners; copy(res.begin(),res.end(),back_inserter(partners)); for ( auto const & d : res ) d->partnerDipoles(partners); return res; } double MatchboxMEBase::colourCorrelatedME2(pair ij) const { if ( matchboxAmplitude() ) { if ( matchboxAmplitude()->treeAmplitudes() ) matchboxAmplitude()->prepareAmplitudes(this); double res = matchboxAmplitude()->colourCorrelatedME2(ij)* me2Norm(); return res; } throw Exception() << "MatchboxMEBase::colourCorrelatedME2() expects a MatchboxAmplitude object.\n" << "Please check your setup." << Exception::runerror; return 0.; } double MatchboxMEBase::largeNColourCorrelatedME2(pair ij, Ptr::tptr largeNBasis) const { if ( matchboxAmplitude() ) { if ( matchboxAmplitude()->treeAmplitudes() ) { largeNBasis->prepare(mePartonData(),false); matchboxAmplitude()->prepareAmplitudes(this); } double res = matchboxAmplitude()->largeNColourCorrelatedME2(ij,largeNBasis)* me2Norm(); return res; } throw Exception() << "MatchboxMEBase::largeNColourCorrelatedME2() expects a MatchboxAmplitude object.\n" << "Please check your setup." << Exception::runerror; return 0.; } double MatchboxMEBase::spinColourCorrelatedME2(pair ij, const SpinCorrelationTensor& c) const { if ( matchboxAmplitude() ) { if ( matchboxAmplitude()->treeAmplitudes() ) matchboxAmplitude()->prepareAmplitudes(this); double res = matchboxAmplitude()->spinColourCorrelatedME2(ij,c)* me2Norm(); return res; } throw Exception() << "MatchboxMEBase::spinColourCorrelatedME2() expects a MatchboxAmplitude object.\n" << "Please check your setup." << Exception::runerror; return 0.; } double MatchboxMEBase::spinCorrelatedME2(pair ij, const SpinCorrelationTensor& c) const { if ( matchboxAmplitude() ) { if ( matchboxAmplitude()->treeAmplitudes() ) matchboxAmplitude()->prepareAmplitudes(this); double res = matchboxAmplitude()->spinCorrelatedME2(ij,c)* me2Norm(); return res; } throw Exception() << "MatchboxMEBase::spinCorrelatedME2() expects a MatchboxAmplitude object.\n" << "Please check your setup." << Exception::runerror; return 0.; } void MatchboxMEBase::flushCaches() { if ( theMerger )theMerger->flushCaches(); MEBase::flushCaches(); if ( matchboxAmplitude() ) matchboxAmplitude()->flushCaches(); for ( auto const & r : reweights() ) r->flushCaches(); for ( auto const & v : virtuals()) v->flushCaches(); } void MatchboxMEBase::setKinematics() { MEBase::setKinematics(); if ( theMerger ) theMerger->setKinematics(); } void MatchboxMEBase::clearKinematics() { MEBase::clearKinematics(); if ( theMerger ) theMerger->clearKinematics(); } const MergerBasePtr MatchboxMEBase::merger() const { return theMerger; } MergerBasePtr MatchboxMEBase::merger() { return theMerger; } void MatchboxMEBase::merger(MergerBasePtr v) { theMerger = v; } void MatchboxMEBase::print(ostream& os) const { os << "--- MatchboxMEBase setup -------------------------------------------------------\n"; os << " '" << name() << "' for subprocess:\n"; os << " "; for ( PDVector::const_iterator pp = subProcess().legs.begin(); pp != subProcess().legs.end(); ++pp ) { os << (**pp).PDGName() << " "; if ( pp == subProcess().legs.begin() + 1 ) os << "-> "; } os << "\n"; os << " including " << (oneLoop() ? "" : "no ") << "virtual corrections"; if ( oneLoopNoBorn() ) os << " without Born contributions"; if ( oneLoopNoLoops() ) os << " without loop contributions"; os << "\n"; if ( oneLoop() && !onlyOneLoop() ) { os << " using insertion operators\n"; for ( vector::ptr>::const_iterator v = virtuals().begin(); v != virtuals().end(); ++v ) { os << " '" << (**v).name() << "' with " << ((**v).isDR() ? "" : "C") << "DR/"; if ( (**v).isCS() ) os << "CS"; if ( (**v).isBDK() ) os << "BDK"; if ( (**v).isExpanded() ) os << "expanded"; os << " conventions\n"; } } os << "--------------------------------------------------------------------------------\n"; os << flush; } void MatchboxMEBase::printLastEvent(ostream& os) const { os << "--- MatchboxMEBase last event information --------------------------------------\n"; os << " for matrix element '" << name() << "'\n"; os << " process considered:\n "; int in = 0; for ( cPDVector::const_iterator p = mePartonData().begin(); p != mePartonData().end(); ++p ) { os << (**p).PDGName() << " "; if ( ++in == 2 ) os << " -> "; } os << " kinematic environment as set by the XComb " << lastXCombPtr() << ":\n" << " sqrt(shat)/GeV = " << sqrt(lastSHat()/GeV2) << " x1 = " << lastX1() << " x2 = " << lastX2() << " alphaS = " << lastAlphaS() << "\n"; os << " momenta/GeV generated from random numbers\n "; copy(lastXComb().lastRandomNumbers().begin(), lastXComb().lastRandomNumbers().end(),ostream_iterator(os," ")); os << ":\n "; for ( vector::const_iterator p = meMomenta().begin(); p != meMomenta().end(); ++p ) { os << (*p/GeV) << "\n "; } os << "last cross section/nb calculated was:\n " << (lastMECrossSection()/nanobarn) << " (pdf weight " << lastMEPDFWeight() << ")\n"; os << "--------------------------------------------------------------------------------\n"; os << flush; } void MatchboxMEBase::logGenerateKinematics(const double * r) const { if ( !verbose() ) return; generator()->log() << "'" << name() << "' generated kinematics\nfrom " << nDim() << " random numbers:\n"; copy(r,r+nDim(),ostream_iterator(generator()->log()," ")); generator()->log() << "\n"; generator()->log() << "storing phase space information in XComb " << lastXCombPtr() << "\n"; generator()->log() << "generated phase space point (in GeV):\n"; vector::const_iterator pit = meMomenta().begin(); cPDVector::const_iterator dit = mePartonData().begin(); for ( ; pit != meMomenta().end() ; ++pit, ++dit ) generator()->log() << (**dit).PDGName() << " : " << (*pit/GeV) << "\n"; generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n" << "and Jacobian = " << jacobian() << " sHat/GeV2 = " << (lastSHat()/GeV2) << "\n" << flush; } void MatchboxMEBase::logSetScale() const { if ( !verbose() ) return; generator()->log() << "'" << name() << "' set scales using XComb " << lastXCombPtr() << ":\n" << "scale/GeV2 = " << (scale()/GeV2) << " xi_R = " << renormalizationScaleFactor() << " xi_F = " << factorizationScaleFactor() << "\n" << "alpha_s = " << lastAlphaS() << "\n" << flush; } void MatchboxMEBase::logPDFWeight() const { if ( !verbose() ) return; generator()->log() << "'" << name() << "' calculated pdf weight = " << lastMEPDFWeight() << " from XComb " << lastXCombPtr() << "\n" << "x1 = " << lastX1() << " (" << (mePartonData()[0]->coloured() ? "" : "not ") << "used) " << "x2 = " << lastX2() << " (" << (mePartonData()[1]->coloured() ? "" : "not ") << "used)\n" << flush; } void MatchboxMEBase::logME2() const { if ( !verbose() ) return; generator()->log() << "'" << name() << "' evaluated me2 using XComb " << lastXCombPtr() << "\n" << "and phase space point (in GeV):\n"; vector::const_iterator pit = meMomenta().begin(); cPDVector::const_iterator dit = mePartonData().begin(); for ( ; pit != meMomenta().end() ; ++pit, ++dit ) generator()->log() << (**dit).PDGName() << " : " << (*pit/GeV) << "\n"; generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n" << "sHat/GeV2 = " << (lastSHat()/GeV2) << "\n" << flush; } void MatchboxMEBase::logDSigHatDR() const { if ( !verbose() ) return; generator()->log() << "'" << name() << "' evaluated cross section using XComb " << lastXCombPtr() << "\n" << "Jacobian = " << jacobian() << " sHat/GeV2 = " << (lastSHat()/GeV2) << " dsig/nb = " << (lastMECrossSection()/nanobarn) << "\n" << flush; } void MatchboxMEBase::cloneDependencies(const std::string& prefix,bool slim) { if ( phasespace() && !slim ) { Ptr::ptr myPhasespace = phasespace()->cloneMe(); ostringstream pname; pname << (prefix == "" ? fullName() : prefix) << "/" << myPhasespace->name(); if ( ! (generator()->preinitRegister(myPhasespace,pname.str()) ) ) throw Exception() << "MatchboxMEBase::cloneDependencies(): Phasespace generator " << pname.str() << " already existing." << Exception::runerror; myPhasespace->cloneDependencies(pname.str()); phasespace(myPhasespace); } theAmplitude = dynamic_ptr_cast::ptr>(amplitude()); if ( matchboxAmplitude() ) { Ptr::ptr myAmplitude = matchboxAmplitude()->cloneMe(); ostringstream pname; pname << (prefix == "" ? fullName() : prefix) << "/" << myAmplitude->name(); if ( ! (generator()->preinitRegister(myAmplitude,pname.str()) ) ){ throw Exception() << "MatchboxMEBase::cloneDependencies(): Amplitude " << pname.str() << " already existing." << Exception::runerror; } myAmplitude->cloneDependencies(pname.str(),slim); matchboxAmplitude(myAmplitude); amplitude(myAmplitude); matchboxAmplitude()->orderInGs(orderInAlphaS()); matchboxAmplitude()->orderInGem(orderInAlphaEW()); } if ( scaleChoice() &&!slim ) { Ptr::ptr myScaleChoice = scaleChoice()->cloneMe(); ostringstream pname; pname << (prefix == "" ? fullName() : prefix) << "/" << myScaleChoice->name(); if ( ! (generator()->preinitRegister(myScaleChoice,pname.str()) ) ) throw Exception() << "MatchboxMEBase::cloneDependencies(): Scale choice " << pname.str() << " already existing." << Exception::runerror; scaleChoice(myScaleChoice); } for ( auto & rw : theReweights ) { Ptr::ptr myReweight = rw->cloneMe(); ostringstream pname; pname << (prefix == "" ? fullName() : prefix) << "/" << rw->name(); if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) ) throw Exception() << "MatchboxMEBase::cloneDependencies(): Reweight " << pname.str() << " already existing." << Exception::runerror; myReweight->cloneDependencies(pname.str()); rw = myReweight; } for ( auto & v : virtuals()) { Ptr::ptr myIOP = v->cloneMe(); ostringstream pname; pname << (prefix == "" ? fullName() : prefix) << "/" << v->name(); if ( ! (generator()->preinitRegister(myIOP,pname.str()) ) ) throw Exception() << "MatchboxMEBase::cloneDependencies(): Insertion operator " << pname.str() << " already existing." << Exception::runerror; v = myIOP; } } void MatchboxMEBase::prepareXComb(MatchboxXCombData& xc) const { // fixme We need to pass on the partons from the xcmob here, not // assuming one subprocess per matrix element if ( phasespace() ) xc.nDimPhasespace(phasespace()->nDim(diagrams().front()->partons())); if ( matchboxAmplitude() ) { xc.nDimAmplitude(matchboxAmplitude()->nDimAdditional()); if ( matchboxAmplitude()->colourBasis() ) { size_t cdim = matchboxAmplitude()->colourBasis()->prepare(diagrams(),noCorrelations()); xc.colourBasisDim(cdim); } if ( matchboxAmplitude()->isExternal() ) { xc.externalId(matchboxAmplitude()->externalId(diagrams().front()->partons())); } } int insertionAdd = 0; for ( auto const & v : virtuals() ) insertionAdd = max(insertionAdd,v->nDimAdditional()); xc.nDimInsertions(insertionAdd); xc.nLight(getNLight()); if(xc.nLightJetVec().empty()) for (auto const & id : getNLightJetVec()) xc.nLightJetVec( id ); if(xc.nHeavyJetVec().empty()) for (auto const & id :getNHeavyJetVec()) xc.nHeavyJetVec(id); if(xc.nLightProtonVec().empty()) for (auto const & id : getNLightProtonVec()) xc.nLightProtonVec(id); xc.olpId(olpProcess()); if ( initVerbose() ) { ostringstream fname_strm; // only allow alphanumeric, / and _ in filename for (const char c : name()) { switch (c) { case '+' : fname_strm << "+"; break; case '-' : fname_strm << "-"; break; case '~' : fname_strm << "_tilde"; break; case ']' : break; case ',' : fname_strm << "__"; break; default : fname_strm << (isalnum(c) ? c : '_'); break; } } fname_strm << ".diagrams"; const string fname = fname_strm.str(); ifstream test(fname.c_str()); if ( !test ) { test.close(); ofstream out(fname.c_str()); for ( vector::ptr>::const_iterator d = diagrams().begin(); d != diagrams().end(); ++d ) { DiagramDrawer::drawDiag(out,dynamic_cast(**d)); out << "\n"; } } } } StdXCombPtr MatchboxMEBase::makeXComb(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, const DiagramVector & newDiagrams, bool mir, const PartonPairVec&, tStdXCombPtr newHead, tMEPtr newME) { if ( !newME ) newME = this; Ptr::ptr xc = new_ptr(MatchboxXComb(newMaxEnergy, inc, newEventHandler, newSubProcessHandler, newExtractor, newCKKW, newPartonBins, newCuts, newME, newDiagrams, mir, newHead)); prepareXComb(*xc); return xc; } StdXCombPtr MatchboxMEBase::makeXComb(tStdXCombPtr newHead, const PBPair & newPartonBins, const DiagramVector & newDiagrams, tMEPtr newME) { if ( !newME ) newME = this; Ptr::ptr xc = new_ptr(MatchboxXComb(newHead, newPartonBins, newME, newDiagrams)); prepareXComb(*xc); return xc; } void MatchboxMEBase::persistentOutput(PersistentOStream & os) const { os << theLastXComb << thePhasespace << theAmplitude << theScaleChoice << theVirtuals << theReweights << theSubprocess << theOneLoop << theOneLoopNoBorn << theOneLoopNoLoops << epsilonSquarePoleHistograms << epsilonPoleHistograms << theMerger << theOLPProcess << theNoCorrelations << theHavePDFs << checkedPDFs; } void MatchboxMEBase::persistentInput(PersistentIStream & is, int) { is >> theLastXComb >> thePhasespace >> theAmplitude >> theScaleChoice >> theVirtuals >> theReweights >> theSubprocess >> theOneLoop >> theOneLoopNoBorn >> theOneLoopNoLoops >> epsilonSquarePoleHistograms >> epsilonPoleHistograms >> theMerger >> theOLPProcess >> theNoCorrelations >> theHavePDFs >> checkedPDFs; lastMatchboxXComb(theLastXComb); } void MatchboxMEBase::Init() { static ClassDocumentation documentation ("MatchboxMEBase is the base class for matrix elements " "in the context of the matchbox NLO interface."); } IBPtr MatchboxMEBase::clone() const { return new_ptr(*this); } IBPtr MatchboxMEBase::fullclone() const { return new_ptr(*this); } void MatchboxMEBase::doinit() { MEBase::doinit(); if ( !theAmplitude ) theAmplitude = dynamic_ptr_cast::ptr>(amplitude()); if ( matchboxAmplitude() ) matchboxAmplitude()->init(); if ( phasespace() ) { phasespace()->init(); matchboxAmplitude()->checkReshuffling(phasespace()); } if ( scaleChoice() ) { scaleChoice()->init(); } for (auto const & rw : theReweights) rw->init(); for (auto const & v : virtuals() ) v->init(); } void MatchboxMEBase::doinitrun() { MEBase::doinitrun(); if ( matchboxAmplitude() ) matchboxAmplitude()->initrun(); if ( phasespace() ) phasespace()->initrun(); if ( scaleChoice() ) scaleChoice()->initrun(); for (auto const & rw : theReweights) rw->initrun(); for (auto const & v : virtuals() ) v->initrun(); } void MatchboxMEBase::dofinish() { MEBase::dofinish(); for (auto const & b : epsilonSquarePoleHistograms ) { b.second.dump(factory()->poleData(),"epsilonSquarePoles-",b.first); } for (auto const & b : epsilonPoleHistograms ) { b.second.dump(factory()->poleData(),"epsilonPoles-",b.first); } } // *** 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 describeHerwigMatchboxMEBase("Herwig::MatchboxMEBase", "Herwig.so"); diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.h b/MatrixElement/Matchbox/Base/MatchboxMEBase.h --- a/MatrixElement/Matchbox/Base/MatchboxMEBase.h +++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.h @@ -1,1211 +1,1203 @@ // -*- C++ -*- // // MatchboxMEBase.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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_MatchboxMEBase_H #define HERWIG_MatchboxMEBase_H // // This is the declaration of the MatchboxMEBase class. // #include "ThePEG/MatrixElement/MEBase.h" #include "Herwig/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h" #include "Herwig/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h" #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h" #include "Herwig/MatrixElement/Matchbox/Utility/ProcessData.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxReweightBase.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh" #include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h" #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh" #include "Herwig/MatrixElement/Matchbox/InsertionOperators/MatchboxInsertionOperator.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh" #include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h" #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXComb.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief MatchboxMEBase is the base class for matrix elements * in the context of the matchbox NLO interface. * * @see \ref MatchboxMEBaseInterfaces "The interfaces" * defined for MatchboxMEBase. */ class MatchboxMEBase: public MEBase, public LastMatchboxXCombInfo { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ MatchboxMEBase(); - /** - * The destructor. - */ - virtual ~MatchboxMEBase(); - //@} - public: /** * Return the factory which produced this matrix element */ Ptr::tptr factory() const; /** @name Subprocess and diagram information. */ //@{ /** * Return the subprocess. */ const Process& subProcess() const { return theSubprocess; } /** * Access the subprocess. */ Process& subProcess() { return theSubprocess; } /** * Return the diagram generator. */ Ptr::tptr diagramGenerator() const; /** * Return the process data. */ Ptr::tptr processData() const; /** * Return true, if this matrix element does not want to * make use of mirroring processes; in this case all * possible partonic subprocesses with a fixed assignment * of incoming particles need to be provided through the diagrams * added with the add(...) method. */ virtual bool noMirror () const { return true; } /** * Add all possible diagrams with the add() function. */ virtual void getDiagrams() const; using MEBase::getDiagrams; /** * With the information previously supplied with the * setKinematics(...) method, a derived class may optionally * override this method to weight the given diagrams with their * (although certainly not physical) relative probabilities. */ virtual Selector diagrams(const DiagramVector &) const; using MEBase::diagrams; /** * Return a Selector with possible colour geometries for the selected * diagram weighted by their relative probabilities. */ virtual Selector colourGeometries(tcDiagPtr diag) const; /** * Return true, if this amplitude is capable of consistently filling * the rho matrices for the spin correllations */ virtual bool canFillRhoMatrix() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->canFillRhoMatrix(); return false; } /** * construct the spin information for the interaction */ virtual void constructVertex(tSubProPtr) {} /** * construct the spin information for the interaction */ virtual void constructVertex(tSubProPtr sub, const ColourLines* cl); /** * Return the order in \f$\alpha_S\f$ in which this matrix element * is given. */ virtual unsigned int orderInAlphaS() const; using MEBase::orderInAlphaS; /** * Return the order in \f$\alpha_{EM}\f$ in which this matrix * element is given. Returns 0. */ virtual unsigned int orderInAlphaEW() const; using MEBase::orderInAlphaEW; /** * Return true, if this amplitude already includes averaging over * incoming parton's quantum numbers. */ virtual bool hasInitialAverage() const { return matchboxAmplitude() ? matchboxAmplitude()->hasInitialAverage() : false; } /** * Return true, if this amplitude already includes symmetry factors * for identical outgoing particles. */ virtual bool hasFinalStateSymmetry() const { return matchboxAmplitude() ? matchboxAmplitude()->hasFinalStateSymmetry() : false; } /** * Return the number of light flavours, this matrix * element is calculated for. */ virtual unsigned int getNLight() const; /** * Return the vector that contains the PDG ids of * the light flavours, which are contained in the * jet particle group. */ virtual vector getNLightJetVec() const; /** * Return the vector that contains the PDG ids of * the heavy flavours, which are contained in the * jet particle group. */ virtual vector getNHeavyJetVec() const; /** * Return the vector that contains the PDG ids of * the light flavours, which are contained in the * proton particle group. */ virtual vector getNLightProtonVec() const; /** * Return true, if this matrix element is handled by a BLHA one-loop provider */ virtual bool isOLPTree() const { return matchboxAmplitude() ? matchboxAmplitude()->isOLPTree() : false; } /** * Return true, if this matrix element is handled by a BLHA one-loop provider */ virtual bool isOLPLoop() const { return matchboxAmplitude() ? matchboxAmplitude()->isOLPLoop() : false; } /** * Return true, if colour and spin correlated matrix elements should * be ordered from the OLP */ virtual bool needsOLPCorrelators() const { return matchboxAmplitude() ? matchboxAmplitude()->needsOLPCorrelators() : true; } /** * Return the process index, if this is an OLP handled matrix element */ const vector& olpProcess() const { return theOLPProcess; } /** * Set the process index, if this is an OLP handled matrix element */ void olpProcess(int pType, int id) { if ( theOLPProcess.empty() ) theOLPProcess.resize(5,0); theOLPProcess[pType] = id; } /** * Return true, if this is a real emission matrix element which does * not require colour correlators. */ bool noCorrelations() const { return theNoCorrelations; } /** * Indicate that this is a real emission matrix element which does * not require colour correlators. */ void needsNoCorrelations() { theNoCorrelations = true; } /** * Indicate that this is a virtual matrix element which does * require colour correlators. */ void needsCorrelations() { theNoCorrelations = false; } //@} /** @name Phasespace generation */ //@{ /** * Return the phase space generator to be used. */ Ptr::tptr phasespace() const { return thePhasespace; } /** * Set the phase space generator to be used. */ void phasespace(Ptr::ptr ps) { thePhasespace = ps; } /** * Set the XComb object to be used in the next call to * generateKinematics() and dSigHatDR(). */ virtual void setXComb(tStdXCombPtr xc); /** * Return true, if the XComb steering this matrix element * should keep track of the random numbers used to generate * the last phase space point */ virtual bool keepRandomNumbers() const { return true; } /** * Generate incoming parton momenta. This default * implementation performs the standard mapping * from x1,x2 -> tau,y making 1/tau flat; incoming * parton momenta are stored in meMomenta()[0,1], * only massless partons are supported so far; * return the Jacobian of the mapping */ double generateIncomingPartons(const double* r1, const double* r2); /** * Generate internal degrees of freedom given nDim() uniform random * numbers in the interval ]0,1[. To help the phase space generator, * the 'dSigHatDR' should be a smooth function of these numbers, * although this is not strictly necessary. The return value should * be true of the generation succeeded. If so the generated momenta * should be stored in the meMomenta() vector. Derived classes * must call this method once internal degrees of freedom are setup * and finally return the result of this method. */ virtual bool generateKinematics(const double * r); /** * Set the typed and momenta of the incoming and outgoing partons to * be used in subsequent calls to me() and colourGeometries() * according to the associated XComb object. If the function is * overridden in a sub class the new function must call the base * class one first. */ virtual void setKinematics(); /** * Clear the information previously provided by a call to * setKinematics(...). */ virtual void clearKinematics(); /** * The number of internal degreed of freedom used in the matrix * element. */ virtual int nDim() const; /** * The number of internal degrees of freedom used in the matrix * element for generating a Born phase space point */ virtual int nDimBorn() const; /** * Return true, if this matrix element will generate momenta for the * incoming partons itself. The matrix element is required to store * the incoming parton momenta in meMomenta()[0,1]. No mapping in * tau and y is performed by the PartonExtractor object, if a * derived class returns true here. The phase space jacobian is to * include a factor 1/(x1 x2). */ virtual bool haveX1X2() const { return (phasespace() ? phasespace()->haveX1X2() : false) || diagrams().front()->partons().size() == 3; } /** * Return true, if this matrix element expects * the incoming partons in their center-of-mass system */ virtual bool wantCMS() const { return (phasespace() ? phasespace()->wantCMS() : true) && diagrams().front()->partons().size() != 3; } /** * Return the meMomenta as generated at the last * phase space point. */ const vector& lastMEMomenta() const { return meMomenta(); } /** * Access the meMomenta. */ vector& lastMEMomenta() { return meMomenta(); } /** * leg size */ int legsize() const {return int(meMomenta().size());} //@} /** @name Scale choices, couplings and PDFs */ //@{ /** * Set the scale choice object */ void scaleChoice(Ptr::ptr sc) { theScaleChoice = sc; } /** * Return the scale choice object */ Ptr::tptr scaleChoice() const { return theScaleChoice; } /** * Set scales and alphaS */ void setScale(Energy2 ren=ZERO,Energy2 fac=ZERO) const; /** * Indicate that this matrix element is running alphas by itself. */ virtual bool hasRunningAlphaS() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->hasRunningAlphaS(); return false; } /** * Indicate that this matrix element is running alphaew by itself. */ virtual bool hasRunningAlphaEW() const { if ( matchboxAmplitude() ) return matchboxAmplitude()->hasRunningAlphaEW(); return false; } /** * Return the scale associated with the phase space point provided * by the last call to setKinematics(). */ virtual Energy2 scale() const { return lastScale(); } /** * Return the renormalization scale for the last generated phasespace point. */ virtual Energy2 factorizationScale() const; /** * Get the factorization scale factor */ virtual double factorizationScaleFactor() const; /** * Get the factorization scale factor */ virtual double facFac() const{return factorizationScaleFactor();} /** * Return the (QCD) renormalization scale for the last generated phasespace point. */ virtual Energy2 renormalizationScale() const; /** * Get the renormalization scale factor */ virtual double renormalizationScaleFactor() const; /** * Get the renormalization scale factor */ virtual double renFac() const{return renormalizationScaleFactor();} /** * Return the QED renormalization scale for the last generated phasespace point. */ virtual Energy2 renormalizationScaleQED() const; /** * Return the shower scale for the last generated phasespace point. */ virtual Energy2 showerScale() const; /** * Set veto scales on the particles at the given * SubProcess which has been generated using this * matrix element. */ virtual void setVetoScales(tSubProPtr) const; /** * Return true, if fixed couplings are used. */ bool fixedCouplings() const; /** * Return true, if fixed couplings are used. */ bool fixedQEDCouplings() const; /** * Return the value of \f$\alpha_S\f$ associated with the phase * space point provided by the last call to setKinematics(). This * versions returns SM().alphaS(scale()). */ virtual double alphaS() const { return lastAlphaS(); } /** * Return the value of \f$\alpha_EM\f$ associated with the phase * space point provided by the last call to setKinematics(). This * versions returns SM().alphaEM(scale()). */ virtual double alphaEM() const { return lastAlphaEM(); } /** * Return true, if this matrix element provides the PDF * weight for the first incoming parton itself. */ virtual bool havePDFWeight1() const; /** * Return true, if this matrix element provides the PDF * weight for the second incoming parton itself. */ virtual bool havePDFWeight2() const; /** * Set the PDF weight. */ void getPDFWeight(Energy2 factorizationScale = ZERO) const; /** * Supply the PDF weight for the first incoming parton. */ double pdf1(Energy2 factorizationScale = ZERO, double xEx = 1., double xFactor = 1.) const; /** * Supply the PDF weight for the second incoming parton. */ double pdf2(Energy2 factorizationScale = ZERO, double xEx = 1., double xFactor = 1.) const; //@} /** @name Amplitude information and matrix element evaluation */ //@{ /** * Return the amplitude. */ Ptr::tptr matchboxAmplitude() const { return theAmplitude; } /** * Set the amplitude. */ void matchboxAmplitude(Ptr::ptr amp) { theAmplitude = amp; } /** * Return the matrix element for the kinematical configuation * previously provided by the last call to setKinematics(), suitably * scaled by sHat() to give a dimension-less number. */ virtual double me2() const; /** * Return the matrix element for the kinematical configuation * previously provided by the last call to setKinematics(), suitably * scaled by sHat() to give a dimension-less number. */ virtual double largeNME2(Ptr::tptr largeNBasis) const; /** * Return the symmetry factor for identical final state particles. */ virtual double finalStateSymmetry() const; /** * Return the normalizing factor for the matrix element averaged * over quantum numbers and including running couplings. */ double me2Norm(unsigned int addAlphaS = 0) const; /** * Return the matrix element squared differential in the variables * given by the last call to generateKinematics(). */ virtual CrossSection dSigHatDR() const; /** * Same prefactor for all dSigHat **/ CrossSection prefactor()const; /** * Born part of the cross section **/ CrossSection dSigHatDRB() const ; /** * Virtual corrections of the cross section **/ CrossSection dSigHatDRV() const ; /** * Insertion operators of the cross section **/ CrossSection dSigHatDRI() const ; /** * If diffAlpha is not 1 and the matrix element has insertion operators * this routine adds the difference between the insertion operator calculated * with an alpha-Parameter to the insertion operator without alpha-parameter. */ CrossSection dSigHatDRAlphaDiff(double alpha) const ; //@} /** @name One-loop corrections */ //@{ /** * Return the one-loop/tree interference. */ virtual double oneLoopInterference() const; /** * Return true, if this matrix element is capable of calculating * one-loop (QCD) corrections. */ virtual bool haveOneLoop() const; /** * Return true, if this matrix element only provides * one-loop (QCD) corrections. */ virtual bool onlyOneLoop() const; /** * Return true, if the amplitude is DRbar renormalized, otherwise * MSbar is assumed. */ virtual bool isDRbar() const; /** * Return true, if one loop corrections have been calculated in * dimensional reduction. Otherwise conventional dimensional * regularization is assumed. Note that renormalization is always * assumed to be MSbar. */ virtual bool isDR() const; /** * Return true, if one loop corrections are given in the conventions * of the integrated dipoles. */ virtual bool isCS() const; /** * Return true, if one loop corrections are given in the conventions * of BDK. */ virtual bool isBDK() const; /** * Return true, if one loop corrections are given in the conventions * of everything expanded. */ virtual bool isExpanded() const; /** * Return the value of the dimensional regularization * parameter. Note that renormalization scale dependence is fully * restored in DipoleIOperator. */ virtual Energy2 mu2() const; /** * If defined, return the coefficient of the pole in epsilon^2 */ virtual double oneLoopDoublePole() const; /** * If defined, return the coefficient of the pole in epsilon */ virtual double oneLoopSinglePole() const; /** * Return true, if cancellationn of epsilon poles should be checked. */ bool checkPoles() const; /** * Simple histogram for accuracy checks */ struct AccuracyHistogram { /** * The lower bound */ double lower; /** * The upper bound */ double upper; /** * The bins, indexed by upper bound. */ map bins; /** * The number of points of same sign */ unsigned long sameSign; /** * The number of points of opposite sign */ unsigned long oppositeSign; /** * The number of points being nan or inf */ unsigned long nans; /** * The overflow */ unsigned long overflow; /** * The underflow */ unsigned long underflow; /** * Constructor */ AccuracyHistogram(double low = -40., double up = 0., unsigned int nbins = 80); /** * Book two values to be checked for numerical compatibility */ void book(double a, double b); /** * Write to file. */ void dump(const std::string& folder, const std::string& prefix, const cPDVector& proc) const; /** * Write to persistent ostream */ void persistentOutput(PersistentOStream&) const; /** * Read from persistent istream */ void persistentInput(PersistentIStream&); }; /** * Perform the check of epsilon pole cancellation. */ void logPoles() const; /** * Return the virtual corrections */ const vector::ptr>& virtuals() const { return theVirtuals; } /** * Return the virtual corrections */ vector::ptr>& virtuals() { return theVirtuals; } /** * Instruct this matrix element to include one-loop corrections */ void doOneLoop() { theOneLoop = true; } /** * Instruct this matrix element not to include one-loop corrections */ void noOneLoop() { theOneLoop = false; } /** * Return true, if this matrix element includes one-loop corrections */ bool oneLoop() const { return theOneLoop; } /** * Instruct this matrix element to include one-loop corrections but * no Born contributions */ void doOneLoopNoBorn() { theOneLoop = true; theOneLoopNoBorn = true; } void noOneLoopNoBorn() { theOneLoop = false; theOneLoopNoBorn = false; } /** * Return true, if this matrix element includes one-loop corrections * but no Born contributions */ bool oneLoopNoBorn() const { return theOneLoopNoBorn || onlyOneLoop(); } /** * Instruct this matrix element to include one-loop corrections but * no actual loop contributions */ void doOneLoopNoLoops() { theOneLoop = true; theOneLoopNoLoops = true; } /** * Return true, if this matrix element includes one-loop corrections * but no actual loop contributions */ bool oneLoopNoLoops() const { return theOneLoopNoLoops; } //@} /** @name Dipole subtraction */ //@{ /** * If this matrix element is considered a real * emission matrix element, return all subtraction * dipoles needed given a set of subtraction terms * and underlying Born matrix elements to choose * from. */ vector::ptr> getDipoles(const vector::ptr>&, const vector::ptr>&,bool slim=false) const; /** * If this matrix element is considered a real emission matrix * element, but actually neglecting a subclass of the contributing * diagrams, return true if the given emitter-emission-spectator * configuration should not be considered when setting up * subtraction dipoles. */ virtual bool noDipole(int,int,int) const { return false; } /** * If this matrix element is considered an underlying Born matrix * element in the context of a subtracted real emission, but * actually neglecting a subclass of the contributing diagrams, * return true if the given emitter-spectator configuration * should not be considered when setting up subtraction dipoles. */ virtual bool noDipole(int,int) const { return false; } /** * Return the colour correlated matrix element squared with * respect to the given two partons as appearing in mePartonData(), * suitably scaled by sHat() to give a dimension-less number. */ virtual double colourCorrelatedME2(pair) const; /** * Return the colour correlated matrix element squared in the * large-N approximation with respect to the given two partons as * appearing in mePartonData(), suitably scaled by sHat() to give a * dimension-less number. */ virtual double largeNColourCorrelatedME2(pair ij, Ptr::tptr largeNBasis) const; /** * Return the colour and spin correlated matrix element squared for * the gluon indexed by the first argument using the given * correlation tensor. */ virtual double spinColourCorrelatedME2(pair emitterSpectator, const SpinCorrelationTensor& c) const; /** * Return the spin correlated matrix element squared for * the vector boson indexed by the first argument using the given * correlation tensor. */ virtual double spinCorrelatedME2(pair emitterSpectator, const SpinCorrelationTensor& c) const; //@} /** @name Caching and diagnostic information */ //@{ /** * Inform this matrix element that a new phase space * point is about to be generated, so all caches should * be flushed. */ virtual void flushCaches(); /** * Return true, if verbose */ bool verbose() const; /** * Return true, if verbose */ bool initVerbose() const; /** * Dump the setup to an ostream */ void print(ostream&) const; /** * Print debug information on the last event */ virtual void printLastEvent(ostream&) const; /** * Write out diagnostic information for * generateKinematics */ void logGenerateKinematics(const double * r) const; /** * Write out diagnostic information for * setting scales */ void logSetScale() const; /** * Write out diagnostic information for * pdf evaluation */ void logPDFWeight() const; /** * Write out diagnostic information for * me2 evaluation */ void logME2() const; /** * Write out diagnostic information * for dsigdr evaluation */ void logDSigHatDR() const; //@} /** @name Reweight objects */ //@{ /** * Insert a reweight object */ void addReweight(Ptr::ptr rw) { theReweights.push_back(rw); } /** * Return the reweights */ const vector::ptr>& reweights() const { return theReweights; } /** * Access the reweights */ vector::ptr>& reweights() { return theReweights; } /** * Return the theMerger. */ const MergerBasePtr merger() const; /** * Return the theMerger. */ MergerBasePtr merger() ; /** * Set the theMerger. */ void merger(MergerBasePtr v); //@} /** @name Methods used to setup MatchboxMEBase objects */ //@{ /** * Return true if this object needs to be initialized before all * other objects (except those for which this function also returns * true). This default version always returns false, but subclasses * may override it to return true. */ virtual bool preInitialize() const { return true; } /** * Clone this matrix element. */ Ptr::ptr cloneMe() const { return dynamic_ptr_cast::ptr>(clone()); } /** * Clone the dependencies, using a given prefix. */ void cloneDependencies(const std::string& prefix = "",bool slim = false ); /** * Prepare an xcomb */ void prepareXComb(MatchboxXCombData&) const; /** * For the given event generation setup return a xcomb object * appropriate to this matrix element. */ virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, const DiagramVector & newDiagrams, bool mir, const PartonPairVec& allPBins, tStdXCombPtr newHead = tStdXCombPtr(), tMEPtr newME = tMEPtr()); /** * For the given event generation setup return a dependent xcomb object * appropriate to this matrix element. */ virtual StdXCombPtr makeXComb(tStdXCombPtr newHead, const PBPair & newPartonBins, const DiagramVector & newDiagrams, tMEPtr newME = tMEPtr()); //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * 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 phase space generator to be used. */ Ptr::ptr thePhasespace; /** * The amplitude to be used */ Ptr::ptr theAmplitude; /** * The scale choice object */ Ptr::ptr theScaleChoice; /** * The virtual corrections. */ vector::ptr> theVirtuals; /** * A vector of reweight objects the sum of which * should be applied to reweight this matrix element */ vector::ptr> theReweights; private: /** * The subprocess to be considered. */ Process theSubprocess; /** * True, if this matrix element includes one-loop corrections */ bool theOneLoop; /** * True, if this matrix element includes one-loop corrections * but no Born contributions */ bool theOneLoopNoBorn; /** * True, if this matrix element includes one-loop corrections * but no actual loop contributions (e.g. finite collinear terms) */ bool theOneLoopNoLoops; /** * The process index, if this is an OLP handled matrix element */ vector theOLPProcess; /** * Histograms of epsilon^2 pole cancellation */ mutable map epsilonSquarePoleHistograms; /** * Histograms of epsilon pole cancellation */ mutable map epsilonPoleHistograms; /** * True, if this is a real emission matrix element which does * not require colour correlators. */ bool theNoCorrelations; /** * Flag which pdfs should be included. */ mutable pair theHavePDFs; /** * True, if already checked for which PDFs to include. */ mutable bool checkedPDFs; /** * The merging helper to be used. * Only the head ME has a pointer to this helper. */ MergerBasePtr theMerger; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MatchboxMEBase & operator=(const MatchboxMEBase &) = delete; }; inline PersistentOStream& operator<<(PersistentOStream& os, const MatchboxMEBase::AccuracyHistogram& h) { h.persistentOutput(os); return os; } inline PersistentIStream& operator>>(PersistentIStream& is, MatchboxMEBase::AccuracyHistogram& h) { h.persistentInput(is); return is; } } #endif /* HERWIG_MatchboxMEBase_H */ diff --git a/MatrixElement/Matchbox/Base/MatchboxOLPME.cc b/MatrixElement/Matchbox/Base/MatchboxOLPME.cc --- a/MatrixElement/Matchbox/Base/MatchboxOLPME.cc +++ b/MatrixElement/Matchbox/Base/MatchboxOLPME.cc @@ -1,355 +1,353 @@ // -*- C++ -*- // // MatchboxOLPME.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 MatchboxOLPME class. // #include "MatchboxOLPME.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "MatchboxMEBase.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; MatchboxOLPME::MatchboxOLPME() : theOrderInGs(0), theOrderInGem(0), theSetMuToMuR(false), theUseRunningAlphaS(false), theUseRunningAlphaEW(false) {} -MatchboxOLPME::~MatchboxOLPME() {} - bool MatchboxOLPME::canHandle(const PDVector& p, Ptr::tptr factory, bool) const { if ( factory->processData()->diagramMap().find(p) != factory->processData()->diagramMap().end() ) return true; vector::ptr> diags = factory->diagramGenerator()->generate(p,orderInGs(),orderInGem()); if ( diags.empty() ) return false; factory->processData()->diagramMap()[p] = diags; return true; } void MatchboxOLPME::setXComb(tStdXCombPtr xc) { theLastXComb = xc; lastMatchboxXComb(xc); } double MatchboxOLPME::me2() const { if ( !calculateTreeME2() ) return lastTreeME2(); evalSubProcess(); return lastTreeME2(); } double MatchboxOLPME::colourCorrelatedME2(pair ij) const { double cfac = 1.; double Nc = generator()->standardModel()->Nc(); if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) { cfac = Nc; } else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 || mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) { cfac = (sqr(Nc)-1.)/(2.*Nc); } else assert(false); if ( !calculateColourCorrelator(ij) ) return lastColourCorrelator(ij)/cfac; evalColourCorrelator(ij); return lastColourCorrelator(ij)/cfac; } double MatchboxOLPME::spinColourCorrelatedME2(pair ij, const SpinCorrelationTensor& c) const { Lorentz5Momentum p = meMomenta()[ij.first]; Lorentz5Momentum n = meMomenta()[ij.second]; LorentzVector polarization = plusPolarization(p,n,ij.first); Complex pFactor = (polarization*c.momentum())/sqrt(abs(c.scale())); double avg = colourCorrelatedME2(ij)*(-c.diagonal()+ (c.scale() > ZERO ? 1. : -1.)*norm(pFactor)); Complex csCorr = 0.0; if ( calculateColourSpinCorrelator(ij) ) evalSpinColourCorrelator(ij); csCorr = lastColourSpinCorrelator(ij); double corr = 2.*real(csCorr*sqr(pFactor)); double Nc = generator()->standardModel()->Nc(); double cfac = 1.; if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) { cfac = Nc; } else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 || mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) { cfac = (sqr(Nc)-1.)/(2.*Nc); } else assert(false); return avg + (c.scale() > ZERO ? 1. : -1.)*corr/cfac; } double MatchboxOLPME::spinCorrelatedME2(pair ij, const SpinCorrelationTensor& c) const { Lorentz5Momentum p = meMomenta()[ij.first]; Lorentz5Momentum n = meMomenta()[ij.second]; LorentzVector polarization = plusPolarization(p,n,ij.first); Complex pFactor = (polarization*c.momentum())/sqrt(abs(c.scale())); double avg = me2()*(-c.diagonal()+ (c.scale() > ZERO ? 1. : -1.)*norm(pFactor)); Complex csCorr = 0.0; if ( calculateSpinCorrelator(ij) ) evalSpinCorrelator(ij); csCorr = lastSpinCorrelator(ij); double corr = 2.*real(csCorr*sqr(pFactor)); return avg + (c.scale() > ZERO ? 1. : -1.)*corr; } void MatchboxOLPME::evalSpinCorrelator(pair) const { throw Exception() << "MatchboxOLPME::spinCorrelatedME2() is not implemented.\n" << "Please check your setup." << Exception::runerror; } double MatchboxOLPME::oneLoopDoublePole() const { if ( !calculateOneLoopPoles() ) return lastOneLoopPoles().first; evalSubProcess(); return lastOneLoopPoles().first; } double MatchboxOLPME::oneLoopSinglePole() const { if ( !calculateOneLoopPoles() ) return lastOneLoopPoles().second; evalSubProcess(); return lastOneLoopPoles().second; } double MatchboxOLPME::oneLoopInterference() const { if ( !calculateOneLoopInterference() ) return lastOneLoopInterference(); evalSubProcess(); return lastOneLoopInterference(); } double MatchboxOLPME::largeNColourCorrelatedME2(pair ij, Ptr::tptr basis) const { if ( trivialColourLegs() ) return MatchboxAmplitude::largeNColourCorrelatedME2(ij,basis); throw Exception() << "MatchboxOLPME::largeNColourCorrelatedME2(): not supported" << Exception::runerror; return 0.; } double MatchboxOLPME::largeNME2(Ptr::tptr basis) const { if ( trivialColourLegs() ) return MatchboxAmplitude::largeNME2(basis); throw Exception() << "MatchboxOLPME::largeNME2(): not supported" << Exception::runerror; return 0.; } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void MatchboxOLPME::doinit() { if ( theUseRunningAlphaS && !theSetMuToMuR ) { throw Exception() << "MatchboxOLPME::doinit(): " << "Amplitude '" << name() << "' " << "uses a running alpha_s but fixed renormalization scale!\n" << Exception::runerror; } if ( !theUseRunningAlphaS && theSetMuToMuR ) { throw Exception() << "MatchboxOLPME::doinit(): " << "Amplitude '" << name() << "' " << "uses a fixed alpha_s but running renormalization scale!\n" << Exception::runerror; } if ( !didStartOLP() ) { string contractFileName = optionalContractFile().empty() ? factory()->buildStorage() + name() + ".OLPContract.lh" : optionalContractFile(); int status = -1; startOLP(contractFileName,status); didStartOLP()=true; if ( status != 1 ) { throw Exception() << "MatchboxOLPME::doinit(): " << "Failed to restart one loop provider for amplitude '" << name() << "'\n" << Exception::runerror; } } MatchboxAmplitude::doinit(); } void MatchboxOLPME::doinitrun() { if ( theUseRunningAlphaS && !theSetMuToMuR ) { throw Exception() << "MatchboxOLPME::doinitrun(): " << "Amplitude '" << name() << "' " << "uses a running alpha_s but fixed renormalization scale!\n" << Exception::runerror; } if ( !theUseRunningAlphaS && theSetMuToMuR ) { throw Exception() << "MatchboxOLPME::doinitrun(): " << "Amplitude '" << name() << "' " << "uses a fixed alpha_s but running renormalization scale!\n" << Exception::runerror; } if ( !didStartOLP() ) { string contractFileName = optionalContractFile().empty() ? factory()->buildStorage() + name() + ".OLPContract.lh" : optionalContractFile(); int status = -1; startOLP(contractFileName,status); didStartOLP()=true; if ( status != 1 ) { throw Exception() << "MatchboxOLPME::doinitrun(): " << "Failed to restart one loop provider for amplitude '" << name() << "'\n" << Exception::runerror; } } MatchboxAmplitude::doinitrun(); } Energy2 MatchboxOLPME::mu2() const { if (theSetMuToMuR) { return lastMatchboxXComb()->lastRenormalizationScale(); } return lastSHat(); } bool MatchboxOLPME::hasRunningAlphaS() const { if (theUseRunningAlphaS) { return true; } return false; } bool MatchboxOLPME::hasRunningAlphaEW() const { if (theUseRunningAlphaEW) { return true; } return false; } void MatchboxOLPME::persistentOutput(PersistentOStream & os) const { os << theOrderInGs << theOrderInGem << theSetMuToMuR << theUseRunningAlphaS << theUseRunningAlphaEW; } void MatchboxOLPME::persistentInput(PersistentIStream & is, int) { is >> theOrderInGs >> theOrderInGem >> theSetMuToMuR >> theUseRunningAlphaS >> theUseRunningAlphaEW; } // *** 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). DescribeAbstractClass describeHerwigMatchboxOLPME("Herwig::MatchboxOLPME", "Herwig.so"); void MatchboxOLPME::Init() { static ClassDocumentation documentation ("MatchboxOLPME implements OLP interfaces."); static Switch interfaceSetMuToMuR ("SetMuToMuR", "Switch On to set the value of the dimensional regularization parameter mu2 for this OLP" "to the value of the renormalization scale muR2. Default is Off. The restoration for the " "full renormalization scale dependence in the DipoleIOperator isn't needed in this case.", &MatchboxOLPME::theSetMuToMuR, false, false, false); static SwitchOption interfaceSetMuToMuRYes (interfaceSetMuToMuR, "Yes", "Yes", true); static SwitchOption interfaceSetMuToMuRNo (interfaceSetMuToMuR, "No", "No", false); interfaceSetMuToMuR.rank(-1); static Switch interfaceUseRunningAlphaS ("UseRunningAlphaS", "Switch On to set the value of alpha_s for this OLP to the value of the running alpha_s " "instead of to the value of the reference alpha_s. Default is Off. This also sets the value " "for hasRunningAlphaS() to true.", &MatchboxOLPME::theUseRunningAlphaS, false, false, false); static SwitchOption interfaceUseRunningAlphaSYes (interfaceUseRunningAlphaS, "Yes", "Yes", true); static SwitchOption interfaceUseRunningAlphaSNo (interfaceUseRunningAlphaS, "No", "No", false); interfaceUseRunningAlphaS.rank(-1); static Switch interfaceUseRunningAlphaEW ("UseRunningAlphaEW", "Switch On to set the value of alpha_ew for this OLP to the value of the running alpha_ew " "instead of to the value of the reference alpha_ew. Default is Off. This also sets the value " "for hasRunningAlphaEW() to true.", &MatchboxOLPME::theUseRunningAlphaEW, false, false, false); static SwitchOption interfaceUseRunningAlphaEWYes (interfaceUseRunningAlphaEW, "Yes", "Yes", true); static SwitchOption interfaceUseRunningAlphaEWNo (interfaceUseRunningAlphaEW, "No", "No", false); interfaceUseRunningAlphaEW.rank(-1); } diff --git a/MatrixElement/Matchbox/Base/MatchboxOLPME.h b/MatrixElement/Matchbox/Base/MatchboxOLPME.h --- a/MatrixElement/Matchbox/Base/MatchboxOLPME.h +++ b/MatrixElement/Matchbox/Base/MatchboxOLPME.h @@ -1,332 +1,324 @@ // -*- C++ -*- // // MatchboxOLPME.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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_MatchboxOLPME_H #define Herwig_MatchboxOLPME_H // // This is the declaration of the MatchboxOLPME class. // #include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief MatchboxOLPME implements OLP interfaces. */ class MatchboxOLPME: public MatchboxAmplitude { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ MatchboxOLPME(); - /** - * The destructor. - */ - virtual ~MatchboxOLPME(); - //@} - public: /** * Return true, if this amplitude can handle the given process. */ virtual bool canHandle(const PDVector& p, Ptr::tptr, bool) const; /** * Set the (tree-level) order in \f$g_S\f$ in which this matrix * element should be evaluated. */ virtual void orderInGs(unsigned int ogs) { theOrderInGs = ogs; } /** * Return the (tree-level) order in \f$g_S\f$ in which this matrix * element is given. */ virtual unsigned int orderInGs() const { return theOrderInGs; } /** * Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix * element should be evaluated. */ virtual void orderInGem(unsigned int oge) { theOrderInGem = oge; } /** * Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix * element is given. */ virtual unsigned int orderInGem() const { return theOrderInGem; } /** * Return true, if this amplitude is handled by a BLHA one-loop provider */ virtual bool isOLPTree() const { return true; } /** * Return true, if this amplitude is handled by a BLHA one-loop provider */ virtual bool isOLPLoop() const { return true; } /** * Return true, if the colour basis is capable of assigning colour * flows. */ virtual bool haveColourFlows() const { return false; } /** * Set the xcomb object. */ virtual void setXComb(tStdXCombPtr xc); /** * Calculate the tree level amplitudes for the phasespace point * stored in lastXComb. */ virtual void prepareAmplitudes(Ptr::tcptr) {} /** * Return the matrix element squared. */ virtual double me2() const; /** * Return the colour correlated matrix element. */ virtual double colourCorrelatedME2(pair ij) const; /** * Return the large-N colour correlated matrix element. */ virtual double largeNColourCorrelatedME2(pair, Ptr::tptr) const; /** * Return the largeN matrix element squared. */ virtual double largeNME2(Ptr::tptr largeNBasis) const; /** * Return the colour and spin correlated matrix element. */ virtual double spinColourCorrelatedME2(pair ij, const SpinCorrelationTensor& c) const; /** * Return the spin correlated matrix element. */ virtual double spinCorrelatedME2(pair ij, const SpinCorrelationTensor& c) const; /** * Return true, if tree-level contributions will be evaluated at amplitude level. */ virtual bool treeAmplitudes() const { return false; } /** * Return true, if this amplitude is capable of calculating one-loop * (QCD) corrections. */ virtual bool haveOneLoop() const { return true; } /** * Return true, if this amplitude only provides * one-loop (QCD) corrections. */ virtual bool onlyOneLoop() const { return false; } /** * Return true, if one-loop contributions will be evaluated at amplitude level. */ virtual bool oneLoopAmplitudes() const { return false; } /** * Return true, if one loop corrections are given in the conventions * of everything expanded. */ virtual bool isExpanded() const { return true; } /** * Return the value of the dimensional regularization * parameter. Note that renormalization scale dependence is fully * restored in DipoleIOperator. */ virtual Energy2 mu2() const; /** * Indicate that this amplitude is running alphas by itself. */ virtual bool hasRunningAlphaS() const; /** * Indicate that this amplitude is running alphaew by itself. */ virtual bool hasRunningAlphaEW() const; /** * If defined, return the coefficient of the pole in epsilon^2 */ virtual double oneLoopDoublePole() const; /** * If defined, return the coefficient of the pole in epsilon */ virtual double oneLoopSinglePole() const; /** * Calculate the one-loop amplitudes for the phasespace point * stored in lastXComb, if provided. */ virtual void prepareOneLoopAmplitudes(Ptr::tcptr) {} /** * Return the one-loop/tree interference. */ virtual double oneLoopInterference() const; public: /** * Call OLP_EvalSubProcess and fill in the results */ virtual void evalSubProcess() const = 0; /** * Fill in results for the given colour correlator */ virtual void evalColourCorrelator(pair ij) const = 0; /** * Fill in results for the given colour/spin correlator */ virtual void evalSpinColourCorrelator(pair ij) const = 0; /** * Fill in results for the given spin correlator; may not be supported */ virtual void evalSpinCorrelator(pair ij) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); // 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(); //@} /** * Set an optional contract file name to be used */ static string& optionalContractFile() { static string s = ""; return s; } /** * Indicate that the OLP has been started */ static bool& didStartOLP() { static bool f = false; return f; } private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MatchboxOLPME & operator=(const MatchboxOLPME &) = delete; /** * The (tree-level) order in \f$g_S\f$ in which this matrix * element is given. */ unsigned int theOrderInGs; /** * The (tree-level) order in \f$g_{EM}\f$ in which this matrix * element is given. */ unsigned int theOrderInGem; /** * Set the value of the dimensional regularization parameter * to the value of the renormalization scale */ bool theSetMuToMuR; /** * Use the running alpha_s instead of the reference alpha_s. * This also sets hasRunningAlphaS() to true. */ bool theUseRunningAlphaS; /** * Use the running alpha_ew instead of the reference alpha_ew. * This also sets hasRunningAlphaEW() to true. */ bool theUseRunningAlphaEW; }; } #endif /* Herwig_MatchboxOLPME_H */ diff --git a/MatrixElement/Matchbox/Base/MatchboxReweightBase.cc b/MatrixElement/Matchbox/Base/MatchboxReweightBase.cc --- a/MatrixElement/Matchbox/Base/MatchboxReweightBase.cc +++ b/MatrixElement/Matchbox/Base/MatchboxReweightBase.cc @@ -1,45 +1,38 @@ // -*- C++ -*- // // MatchboxReweightBase.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 MatchboxReweightBase class. // #include "MatchboxReweightBase.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; -MatchboxReweightBase::MatchboxReweightBase() {} - -MatchboxReweightBase::~MatchboxReweightBase() {} - void MatchboxReweightBase::cloneDependencies(const std::string&) {} void MatchboxReweightBase::persistentOutput(PersistentOStream &) const {} void MatchboxReweightBase::persistentInput(PersistentIStream &, int) {} void MatchboxReweightBase::Init() { static ClassDocumentation documentation ("MatchboxReweightBase"); } -// *** 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). +// The following static variable is needed for the type +// description system in ThePEG. DescribeAbstractClass describeMatchboxReweightBase("Herwig::MatchboxReweightBase", "Herwig.so"); diff --git a/MatrixElement/Matchbox/Base/MatchboxReweightBase.h b/MatrixElement/Matchbox/Base/MatchboxReweightBase.h --- a/MatrixElement/Matchbox/Base/MatchboxReweightBase.h +++ b/MatrixElement/Matchbox/Base/MatchboxReweightBase.h @@ -1,146 +1,131 @@ // -*- C++ -*- // // MatchboxReweightBase.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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_MatchboxReweightBase_H #define HERWIG_MatchboxReweightBase_H // // This is the declaration of the MatchboxReweightBase class. // #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/Handlers/StandardXComb.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief MatchboxReweightBase is the base class * for reweighting MatchboxMEBase matrix elements * as |M|^2 ( w_1 + ... + w_n ) * */ class MatchboxReweightBase: public HandlerBase { public: - /** @name Standard constructors and destructors. */ - //@{ - /** - * The default constructor. - */ - MatchboxReweightBase(); - - /** - * The destructor. - */ - virtual ~MatchboxReweightBase(); - //@} - -public: - /** * Clone this reweight. */ Ptr::ptr cloneMe() const { return dynamic_ptr_cast::ptr>(clone()); } /** * Clone the dependencies, using a given prefix. */ virtual void cloneDependencies(const std::string& prefix = ""); /** * Set the XComb object. */ virtual void setXComb(tStdXCombPtr) = 0; /** * Return true, if applies to the process in the xcomb. */ virtual bool apply() const = 0; /** * Inform this matrix element that a new phase space * point is about to be generated, so all caches should * be flushed. */ virtual void flushCaches() = 0; /** * Evaluate the reweight. */ virtual double evaluate() const = 0; /** * Set veto scales on the particles at the given * SubProcess which has been generated using this * matrix element. */ virtual void setVetoScales(tSubProPtr) const {} public: /** * Dump the setup to an ostream */ virtual void print(ostream&) const {} /** * Print debug information on the last event */ virtual void printLastEvent(ostream&) const {} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); // 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. */ MatchboxReweightBase & operator=(const MatchboxReweightBase &) = delete; }; } #endif /* HERWIG_MatchboxReweightBase_H */ diff --git a/MatrixElement/Matchbox/Base/MergerBase.cc b/MatrixElement/Matchbox/Base/MergerBase.cc --- a/MatrixElement/Matchbox/Base/MergerBase.cc +++ b/MatrixElement/Matchbox/Base/MergerBase.cc @@ -1,44 +1,39 @@ // -*- C++ -*- // // MergerBase.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 MergerBase class. // #include "MergerBase.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; MergerBase::MergerBase() : HandlerBase() {} -MergerBase::~MergerBase() {} - void MergerBase::Init() { static ClassDocumentation documentation ("MergerBase is the base class for merging helpers."); } -// *** 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). +// The following static variable is needed for the type +// description system in ThePEG. DescribeAbstractNoPIOClass describeHerwigMergerBase("Herwig::MergerBase", "Herwig.so"); 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,100 +1,94 @@ // -*- C++ -*- // // MergerBase.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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; //@} 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/MatrixElement/Matchbox/Base/SubtractedME.cc b/MatrixElement/Matchbox/Base/SubtractedME.cc --- a/MatrixElement/Matchbox/Base/SubtractedME.cc +++ b/MatrixElement/Matchbox/Base/SubtractedME.cc @@ -1,790 +1,788 @@ // -*- C++ -*- // // SubtractedME.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 SubtractedME class. // #include "SubtractedME.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Repository/Repository.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/Handlers/StdXCombGroup.h" #include "ThePEG/Utilities/Rebinder.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h" #include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h" #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h" using namespace Herwig; SubtractedME::SubtractedME() : MEGroup(), theRealShowerSubtraction(false), theVirtualShowerSubtraction(false), theLoopSimSubtraction(false) {} -SubtractedME::~SubtractedME() {} - Ptr::tcptr SubtractedME::factory() const { return MatchboxFactory::currentFactory(); } bool SubtractedME::subProcessGroups() const { return (factory()->subProcessGroups() && !showerApproximation()) || factory()->subtractionData() != ""; } Ptr::tptr SubtractedME::showerApproximation() const { return factory()->showerApproximation(); } const vector::ptr>& SubtractedME::borns() const { return theBorns.empty() ? factory()->bornMEs() : theBorns; } bool SubtractedME::verbose() const { return factory()->verbose(); } bool SubtractedME::initVerbose() const { return factory()->initVerbose(); } IBPtr SubtractedME::clone() const { return new_ptr(*this); } IBPtr SubtractedME::fullclone() const { return new_ptr(*this); } StdXCombPtr SubtractedME::makeXComb(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, const DiagramVector & newDiagrams, bool mir, const PartonPairVec& allPBins, tStdXCombPtr newHead, tMEPtr newME) { tMEGroupPtr newMEGroup = dynamic_ptr_cast(newME); if ( !newMEGroup ) newMEGroup = this; Ptr::ptr res = new_ptr(MatchboxXCombGroup(newMaxEnergy, inc, newEventHandler, newSubProcessHandler, newExtractor, newCKKW, newPartonBins, newCuts, newMEGroup, newDiagrams, mir, newHead)); res->build(allPBins); theReal->prepareXComb(*res); if ( factory()->subtractionData() != "" ) { set procs; for ( DiagramVector::const_iterator d = head()->diagrams().begin(); d != head()->diagrams().end(); ++d ) { if ( procs.find((**d).partons()) == procs.end() ) procs.insert((**d).partons()); } for ( set::const_iterator p = procs.begin(); p != procs.end(); ++p ) { for ( size_t i = 0; i < (*p).size(); ++i ) { if ( !(*p)[i]->coloured() ) continue; if ( i > 1 && (*p)[i]->id() == ParticleID::g ) { softHistograms[SoftSubtractionIndex(*p,i)] = SubtractionHistogram(0.00001,1000.); ostringstream fname(""); fname << factory()->subtractionData(); const cPDVector& myproc = SoftSubtractionIndex(*p,i).first; for (cPDVector::const_iterator pp = myproc.begin(); pp != myproc.end(); ++pp) fname << (**pp).PDGName(); fname << "-" << i << "-" << i << "-scatter.dat"; fnamesSoftSubtraction[SoftSubtractionIndex(*p,i)] = fname.str(); if ( theReal->phasespace() ) res->singularLimits().insert(make_pair(i,i)); } for ( size_t j = i+1; j < (*p).size(); ++j ) { if ( !(*p)[j]->coloured() ) continue; long iid = (*p)[i]->id(); long jid = (*p)[j]->id(); if ( i < 2 && j < 2 ) continue; if ( i < 2 && j > 1 ) { if ( abs(iid) < 7 && abs(jid) < 7 && iid != jid ) continue; } if ( i > 1 && j > 1 ) { if ( abs(iid) < 7 && abs(jid) < 7 && iid + jid != 0 ) continue; } bool haveDipole = false; for ( MEVector::const_iterator k = dependent().begin(); k != dependent().end(); ++k ) { const SubtractionDipole& dip = dynamic_cast(**k); if ( ( (size_t)(dip.realEmitter()) == i && (size_t)(dip.realEmission()) == j ) || ( (size_t)(dip.realEmitter()) == j && (size_t)(dip.realEmission()) == i ) ) { haveDipole = true; break; } } if ( !haveDipole ) continue; collinearHistograms[CollinearSubtractionIndex(*p,make_pair(i,j))] = SubtractionHistogram(0.00001,1000.); ostringstream fname(""); fname << factory()->subtractionData(); const cPDVector& myproc = CollinearSubtractionIndex(*p,make_pair(i,j)).first; for (cPDVector::const_iterator pp = myproc.begin(); pp != myproc.end(); ++pp) fname << (**pp).PDGName(); fname << "-" << i << "-" << j << "-scatter.dat"; fnamesCollinearSubtraction[CollinearSubtractionIndex(*p,make_pair(i,j))] = fname.str(); if ( theReal->phasespace() ) res->singularLimits().insert(make_pair(i,j)); } } } } return res; } void SubtractedME::setXComb(tStdXCombPtr xc) { MEGroup::setXComb(xc); lastMatchboxXComb(xc); } MEBase::DiagramVector SubtractedME::dependentDiagrams(const cPDVector& proc, tMEPtr depME) const { Ptr::tptr dipole = dynamic_ptr_cast::tptr>(depME); if ( !dipole ) { throw Exception() << "SubtractedME: A dependent matrix element of SubtractedME " << "has not been derived from SubtractionDipole. " << "Please check the corresponding input file." << Exception::runerror; } return dipole->underlyingBornDiagrams(proc); } vector::ptr> SubtractedME::dipoles() { if ( dependent().empty() ) getDipoles(); vector::ptr> res; for ( MEVector::const_iterator k = dependent().begin(); k != dependent().end(); ++k ) res.push_back(dynamic_ptr_cast::ptr>(*k)); return res; } void SubtractedME::getDipoles() { if ( !dependent().empty() ) return; Ptr::tptr real = dynamic_ptr_cast::tptr>(head()); if ( borns().empty() || !real ) throw Exception() << "SubtractedME: The SubtractedME '" << name() << "' could not generate " << "subtraction terms for the real emission " << "matrix element '" << real->name() << "'. " << "Please check the corresponding input file." << Exception::runerror; Ptr::ptr myRealEmissionME = real->cloneMe(); ostringstream pname; pname << fullName() << "/" << myRealEmissionME->name(); if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) ) throw Exception() << "SubtractedME: Matrix element " << pname.str() << " already existing." << Exception::runerror; myRealEmissionME->cloneDependencies(pname.str()); head(myRealEmissionME); real = myRealEmissionME; dependent().clear(); vector::ptr> genDipoles = real->getDipoles(DipoleRepository::dipoles(factory()->dipoleSet()),borns()); if ( factory()->subtractionData() != "" ) { for ( vector::ptr>::const_iterator d = genDipoles.begin(); d != genDipoles.end(); ++d ) (**d).doTestSubtraction(); } if ( genDipoles.empty() && factory()->initVerbose() ) { // probably finite real contribution, but warn generator()->log() << "\nWarning: No subtraction dipoles could be found for the process:\n"; generator()->log() << real->subProcess().legs[0]->PDGName() << " " << real->subProcess().legs[1]->PDGName() << " -> "; for ( PDVector::const_iterator p = real->subProcess().legs.begin() + 2; p != real->subProcess().legs.end(); ++p ) generator()->log() << (**p).PDGName() << " "; generator()->log() << "\n" << flush; generator()->log() << "Assuming finite tree-level O(alphaS) correction.\n"; } dependent().resize(genDipoles.size()); copy(genDipoles.begin(),genDipoles.end(),dependent().begin()); if ( !factory()->reweighters().empty() ) { for ( MEVector::const_iterator d = dependent().begin(); d != dependent().end(); ++d ) { for ( vector::const_iterator rw = factory()->reweighters().begin(); rw != factory()->reweighters().end(); ++rw ) (**d).addReweighter(*rw); } } if ( !factory()->preweighters().empty() ) { for ( MEVector::const_iterator d = dependent().begin(); d != dependent().end(); ++d ) { for ( vector::const_iterator rw = factory()->preweighters().begin(); rw != factory()->preweighters().end(); ++rw ) (**d).addPreweighter(*rw); } } } void SubtractedME::cloneRealME(const string& prefix) { theReal = dynamic_ptr_cast::tptr>(head()); if ( theReal ) { Ptr::ptr myRealEmissionME = theReal->cloneMe(); ostringstream pname; pname << (prefix == "" ? fullName() : prefix) << "/" << myRealEmissionME->name(); if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) ) throw Exception() << "SubtractedME: Matrix element " << pname.str() << " already existing." << Exception::runerror; myRealEmissionME->cloneDependencies(pname.str()); theReal = myRealEmissionME; } head(theReal); } void SubtractedME::cloneDipoles(const string& prefix) { MEVector dipMEs; for ( MEVector::const_iterator m = dependent().begin(); m != dependent().end(); ++m ) { Ptr::tptr dip = dynamic_ptr_cast::tptr>(*m); assert(dip); Ptr::ptr cloned = dip->cloneMe(); ostringstream pname; pname << (prefix == "" ? fullName() : prefix) << "/" << cloned->name(); if ( ! (generator()->preinitRegister(cloned,pname.str()) ) ) throw Exception() << "SubtractedME: Subtraction dipole " << pname.str() << " already existing." << Exception::runerror; cloned->cloneDependencies(pname.str()); dipMEs.push_back(cloned); } dependent() = dipMEs; } vector::ptr> SubtractedME::splitDipoles(const cPDVector& born) { vector::ptr> dips = dipoles(); vector::ptr> res; for ( vector::ptr>::iterator d = dips.begin(); d != dips.end(); ++d ) { for ( DiagramVector::const_iterator p = (**d).underlyingBornME()->diagrams().begin(); p != (**d).underlyingBornME()->diagrams().end(); ++p ) if ( born == (**p).partons() ) { res.push_back(*d); break; } } return res; } void SubtractedME::doRealEmissionScales() { for ( MEVector::const_iterator m = dependent().begin(); m != dependent().end(); ++m ) { Ptr::tptr dip = dynamic_ptr_cast::tptr>(*m); assert(dip); dip->doRealEmissionScales(); } } void SubtractedME::doRealShowerSubtraction() { theRealShowerSubtraction = true; for ( MEVector::const_iterator m = dependent().begin(); m != dependent().end(); ++m ) { Ptr::tptr dip = dynamic_ptr_cast::tptr>(*m); assert(dip); dip->showerApproximation(showerApproximation()); dip->doRealShowerSubtraction(); } } void SubtractedME::doVirtualShowerSubtraction() { theVirtualShowerSubtraction = true; for ( MEVector::const_iterator m = dependent().begin(); m != dependent().end(); ++m ) { Ptr::tptr dip = dynamic_ptr_cast::tptr>(*m); assert(dip); dip->showerApproximation(showerApproximation()); dip->doVirtualShowerSubtraction(); } } void SubtractedME::doLoopSimSubtraction() { theLoopSimSubtraction = true; for ( MEVector::const_iterator m = dependent().begin(); m != dependent().end(); ++m ) { Ptr::tptr dip = dynamic_ptr_cast::tptr>(*m); assert(dip); dip->showerApproximation(showerApproximation()); dip->doLoopSimSubtraction(); } } void SubtractedME::setVetoScales(tSubProPtr) const {} void SubtractedME::fillProjectors() { if ( !virtualShowerSubtraction() && !loopSimSubtraction() ) return; Ptr::tptr group = dynamic_ptr_cast::tptr>(lastXCombPtr()); for ( vector::const_iterator d = group->dependent().begin(); d != group->dependent().end(); ++d ) { if ( !(**d).matrixElement()->apply() || !(**d).kinematicsGenerated() ) continue; if ( (**d).willPassCuts() && (**d).lastMECrossSection()/picobarn != 0.0 ) { lastXCombPtr()->projectors().insert(abs((**d).cutWeight()*(**d).lastMECrossSection()/picobarn),*d);// } } } double SubtractedME::reweightHead(const vector&) { if ( showerApproximation() ) { if ( realShowerSubtraction() ) return 1.; if ( virtualShowerSubtraction() || loopSimSubtraction() ) return 0.; } return 1.; } double SubtractedME::reweightDependent(tStdXCombPtr xc, const vector& dep) { if ( showerApproximation() ) { if ( realShowerSubtraction() ) return 1.0; if ( virtualShowerSubtraction() || loopSimSubtraction() ) { if ( !lastXComb().lastProjector() ) return 0.0; if ( xc != lastXComb().lastProjector() ) return 0.0; double invPAlpha = 0.; for ( vector::const_iterator d = dep.begin(); d != dep.end(); ++d ) { if ( !(**d).matrixElement()->apply() || !(**d).kinematicsGenerated() ) continue; if ( (**d).willPassCuts() && (**d).lastMECrossSection()/picobarn != 0.0 ) { invPAlpha += abs((**d).cutWeight()*(**d).lastMECrossSection()/picobarn); } } assert(invPAlpha != 0.0 && xc->cutWeight() != 0.0 && xc->lastMECrossSection()/picobarn != 0.0); double palpha = abs((xc->cutWeight())*(xc->lastMECrossSection()/picobarn))/invPAlpha; return 1./palpha; } } return 1.; } void SubtractedME::doinit() { // has been deactivated by the factory if ( !head() ) { MEBase::doinit(); return; } theReal = dynamic_ptr_cast::tptr>(head()); if ( theReal ) { getDipoles(); } for ( vector::ptr>::iterator b = theBorns.begin(); b != theBorns.end(); ++b ) (**b).init(); if ( initVerbose() ) print(Repository::clog()); MEGroup::doinit(); } void SubtractedME::doinitrun() { // has been deactivated by the factory if ( !head() ) { MEBase::doinitrun(); return; } theReal = dynamic_ptr_cast::tptr>(head()); for ( vector::ptr>::iterator b = theBorns.begin(); b != theBorns.end(); ++b ) (**b).initrun(); MEGroup::doinitrun(); } void SubtractedME::dofinish() { // has been deactivated by the factory if ( !head() ) { MEBase::dofinish(); return; } MEGroup::dofinish(); for ( map:: const_iterator b = collinearHistograms.begin(); b != collinearHistograms.end(); ++b ) { b->second.dump(factory()->subtractionData(), factory()->subtractionPlotType(), factory()->subtractionScatterPlot(), b->first.first, b->first.second.first, b->first.second.second); } for ( map:: const_iterator b = softHistograms.begin(); b != softHistograms.end(); ++b ) { b->second.dump(factory()->subtractionData(), factory()->subtractionPlotType(), factory()->subtractionScatterPlot(), b->first.first, b->first.second, b->first.second); } } void SubtractedME::print(ostream& os) const { os << "--- SubtractedME setup ---------------------------------------------------------\n"; os << " '" << name() << "' subtracting real emission\n '" << head()->name() << "' using the dipoles:\n"; for ( MEVector::const_iterator d = dependent().begin(); d != dependent().end(); ++d ) dynamic_ptr_cast::tptr>(*d)->print(os); os << "--------------------------------------------------------------------------------\n"; os << flush; } void SubtractedME::printLastEvent(ostream& os) const { os << "--- SubtractedME last event information ----------------------------------------\n"; os << " for subtracted matrix element '" << name() << "'\n"; os << " real emission event information:\n"; dynamic_ptr_cast::tptr>(head())->printLastEvent(os); os << " dipoles event information:\n"; for ( MEVector::const_iterator d = dependent().begin(); d != dependent().end(); ++d ) dynamic_ptr_cast::tptr>(*d)->printLastEvent(os); os << "--- end SubtractedME last event information ------------------------------------\n\n\n"; os << flush; } void SubtractedME::lastEventStatistics() { MEGroup::lastEventStatistics(); if ( !generator() ) return; /* if ( verbose() ) printLastEvent(generator()->log()); */ if ( !collinearHistograms.empty() ) lastEventSubtraction(); } SubtractedME::SubtractionHistogram:: SubtractionHistogram(double low, double up, unsigned int nbins) : lower(low) { nbins = nbins + 1; double c = log10(up/low) / (nbins-1.); for ( unsigned int k = 1; k < nbins; ++k ) { bins[low*pow(10.0,k*c)] = make_pair(Constants::MaxDouble,0.); } } void SubtractedME::SubtractionHistogram::persistentOutput(PersistentOStream& os) const { os << lower << bins; } void SubtractedME::SubtractionHistogram::persistentInput(PersistentIStream& is) { is >> lower >> bins; } void SubtractedME::SubtractionHistogram:: dump(const std::string& prefix, const int& plottype, const bool& scatterplot, const cPDVector& proc, int i, int j) const { bool bbmin = true; double bmin = bins.begin()->first; double bmax = bins.begin()->first; ostringstream fname(""); for ( cPDVector::const_iterator p = proc.begin(); p != proc.end(); ++p ) fname << (**p).PDGName(); fname << "-" << i << "-" << j; ofstream out((prefix+fname.str()+".dat").c_str()); for ( map >::const_iterator b = bins.begin(); b != bins.end(); ++b ) { map >::const_iterator bp = b; if (bp== bins.begin())continue; --bp; if ( b->second.first != Constants::MaxDouble || b->second.second != 0.0 ) { if ( b != bins.begin() ){ out << bp->first; if (bbmin){ bmin = bp->first; bbmin = false; } } else { out << lower; if (bbmin){ bmin = lower; bbmin = false; } } bmax = b->first; out << " " << b->first << " " << b->second.first << " " << b->second.second << "\n" << flush; } } double xmin = pow(10.0, floor(log10(bmin))); double xmax = pow(10.0, ceil(log10(bmax))); ofstream gpout((prefix+fname.str()+".gp").c_str()); gpout << "set terminal epslatex color solid\n" << "set output '" << fname.str() << "-plot.tex'\n" << "set format x '$10^{%T}$'\n" << "set logscale x\n" << "set xrange [" << xmin << ":" << xmax << "]\n"; if ( i != j ) { gpout << "set xlabel '$\\sqrt{s_{" << i << j << "}}/{\\rm GeV}$'\n"; } else { gpout << "set xlabel '$E_{" << i << "}/{\\rm GeV}$'\n"; } if (plottype == 1){ gpout << "set size 0.5,0.6\n" << "set yrange [0:2]\n"; gpout << "plot 1 w lines lc rgbcolor \"#DDDDDD\" notitle, '" << fname.str() << ".dat' u (($1+$2)/2.):3:($4 < 4. ? $4 : 4.) w filledcurves lc rgbcolor \"#00AACC\" t '$"; } else if (plottype == 2){ gpout << "set key left top Left reverse\n" << "set logscale y\n" << "set format y '$10^{%T}$'\n" << "set size 0.7,0.8\n" << "set yrange [1e-6:1e1]\n" << "set ylabel '$\\max\\left\\{\\left|\\mathcal{D}-\\mathcal{M}\\right|/\\left|\\mathcal{M}\\right|\\right\\}$'\n" << "unset bars\n"; gpout << "plot '"; if (scatterplot) gpout << fname.str() << "-scatter.dat' w points pt 7 ps 0.5 lc rgbcolor \"#00AACC\" not, \\\n'"; gpout << fname.str() << ".dat' u (($1+$2)/2.):4 w lines lw 4 lc rgbcolor \"#00AACC\" t '$"; } for ( size_t k = 0; k < proc.size(); k++ ) { if ( k == 2 ) gpout << "\\to "; gpout << (proc[k]->id() < 0 ? "\\bar{" : "") << (proc[k]->id() < 0 ? proc[k]->CC()->PDGName() : proc[k]->PDGName()) << (proc[k]->id() < 0 ? "}" : "") << " "; } gpout << "$'\n"; gpout << "reset\n"; } void SubtractedME::lastEventSubtraction() { tStdXCombGroupPtr xc = dynamic_ptr_cast(lastXCombPtr()); CrossSection xcme2 = xc->lastHeadCrossSection(); CrossSection xcdip = ZERO; for ( vector::const_iterator d = xc->dependent().begin(); d != xc->dependent().end(); ++d ) { if ( !(*d) ) continue; if ( !(**d).matrixElement()->apply() ) continue; if ( !(**d).willPassCuts() ) continue; xcdip += (**d).lastCrossSection(); } // want a real emission safely above the cut if ( xc->cutWeight() < 1.0 ) return; double delta; if (factory()->subtractionPlotType() == 2) delta = abs(xcdip+xcme2)/abs(xcme2); else delta = abs(xcdip)/abs(xcme2); if ( theReal->phasespace() ) { size_t i = lastSingularLimit()->first; size_t j = lastSingularLimit()->second; if ( i == j && softHistograms.find(SoftSubtractionIndex(head()->mePartonData(),i)) != softHistograms.end() ) { softHistograms[SoftSubtractionIndex(head()->mePartonData(),i)]. book(meMomenta()[i].t()/GeV,delta); if ( factory()->subtractionScatterPlot() ){ ofstream outstream((fnamesSoftSubtraction[SoftSubtractionIndex(head()->mePartonData(),i)]).c_str(),ofstream::app); outstream << meMomenta()[i].t()/GeV << " " << delta << "\n"; } } if ( i != j && collinearHistograms.find(CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))) != collinearHistograms.end() ) { double s = sqrt(2.*meMomenta()[i]*meMomenta()[j])/GeV; collinearHistograms[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))]. book(s,delta); if ( factory()->subtractionScatterPlot() ){ ofstream outstream((fnamesCollinearSubtraction[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))]).c_str(),ofstream::app); outstream << s << " " << delta << "\n"; } } return; } for ( size_t i = 0; i < meMomenta().size(); ++i ) { if ( i > 1 ) { if ( softHistograms.find(SoftSubtractionIndex(head()->mePartonData(),i)) != softHistograms.end() ) { softHistograms[SoftSubtractionIndex(head()->mePartonData(),i)]. book(meMomenta()[i].t()/GeV,delta); if ( factory()->subtractionScatterPlot() ){ ofstream outstream((fnamesSoftSubtraction[SoftSubtractionIndex(head()->mePartonData(),i)]).c_str(),ofstream::app); outstream << meMomenta()[i].t()/GeV << " " << delta << "\n"; } } } for ( size_t j = i+1; j < meMomenta().size(); ++j ) { if ( collinearHistograms.find(CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))) == collinearHistograms.end() ) continue; double s = sqrt(2.*meMomenta()[i]*meMomenta()[j])/GeV; collinearHistograms[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))]. book(s,delta); if ( factory()->subtractionScatterPlot() ){ ofstream outstream((fnamesCollinearSubtraction[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))]).c_str(),ofstream::app); outstream << s << " " << delta << "\n"; } } } } void SubtractedME::persistentOutput(PersistentOStream & os) const { os << theLastXComb << theBorns << theReal << collinearHistograms << softHistograms << fnamesSoftSubtraction << theRealShowerSubtraction << theVirtualShowerSubtraction << theLoopSimSubtraction; } void SubtractedME::persistentInput(PersistentIStream & is, int) { is >> theLastXComb >> theBorns >> theReal >> collinearHistograms >> softHistograms >> fnamesSoftSubtraction >> theRealShowerSubtraction >> theVirtualShowerSubtraction >> theLoopSimSubtraction; lastMatchboxXComb(theLastXComb); } void SubtractedME::Init() { static ClassDocumentation documentation ("SubtractedME represents a subtracted real emission matrix element."); } // *** 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 describeHerwigSubtractedME("Herwig::SubtractedME", "Herwig.so"); diff --git a/MatrixElement/Matchbox/Base/SubtractedME.h b/MatrixElement/Matchbox/Base/SubtractedME.h --- a/MatrixElement/Matchbox/Base/SubtractedME.h +++ b/MatrixElement/Matchbox/Base/SubtractedME.h @@ -1,536 +1,528 @@ // -*- C++ -*- // // SubtractedME.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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_SubtractedME_H #define HERWIG_SubtractedME_H // // This is the declaration of the SubtractedME class. // #include "ThePEG/MatrixElement/MEGroup.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h" #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h" #include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief SubtractedME represents a subtracted real emission matrix element. * * @see \ref SubtractedMEInterfaces "The interfaces" * defined for SubtractedME. */ class SubtractedME: public MEGroup, public LastMatchboxXCombInfo { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ SubtractedME(); - /** - * The destructor. - */ - virtual ~SubtractedME(); - //@} - public: /** * Return the factory which produced this matrix element */ Ptr::tcptr factory() const; /** @name Phasespace and subprocess information */ //@{ /** * For the given event generation setup return an xcomb object * appropriate to this matrix element. */ virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, const DiagramVector & newDiagrams, bool mir, const PartonPairVec& allPBins, tStdXCombPtr newHead = tStdXCombPtr(), tMEPtr newME = tMEPtr()); /** * Set the XComb object to be used in the next call to * generateKinematics() and dSigHatDR(). */ virtual void setXComb(tStdXCombPtr); /** * Return true, if the same additional random numbers * should be presented to any of the dependent * matrix elements. */ virtual bool uniformAdditional() const { return true; } /** * Return true, if the XComb steering this matrix element * should keep track of the random numbers used to generate * the last phase space point */ virtual bool keepRandomNumbers() const { return true; } /** * Given a process from the head matrix element, * return a list of diagrams which should be considered for * the given dependent matrix element. */ virtual MEBase::DiagramVector dependentDiagrams(const cPDVector& proc, tMEPtr depME) const; /** * Return true, if SubProcessGroups should be * setup from this MEGroup. If not, a single SubProcess * is constructed from the data provided by the * head matrix element. */ virtual bool subProcessGroups() const; /** * Return true, if one of the dependent subprocesses should be * constructed in place of the one driven by the head matrix element * or a full subprocess group. */ virtual bool selectDependentSubProcess() const { return false; } /** * Fill the projectors object of xcombs to choose subprocesses * different than the one currently integrated. */ virtual void fillProjectors(); /** * Return true, if projectors will be used */ virtual bool willProject() const { return virtualShowerSubtraction() || loopSimSubtraction(); } /** * Return true, if this MEGroup will reweight the contributing cross * sections. */ virtual bool groupReweighted() const { return showerApproximation(); } /** * Reweight the head cross section */ virtual double reweightHead(const vector&); /** * Reweight the dependent cross section */ virtual double reweightDependent(tStdXCombPtr, const vector&); /** * Switch on or off that scales should be calculated from real emission kinematics */ void doRealEmissionScales(); //@} /** @name Methods relevant to matching */ //@{ /** * Inform this matrix element that a new phase space * point is about to be generated, so all caches should * be flushed. */ virtual void flushCaches() { MEGroup::flushCaches(); if ( showerApproximation() ) showerApproximation()->resetBelowCutoff(); } /** * Return the shower approximation. */ Ptr::tptr showerApproximation() const; /** * Indicate that the shower real emission contribution should be subtracted. */ void doRealShowerSubtraction(); /** * Return true, if the shower real emission contribution should be subtracted. */ bool realShowerSubtraction() const { return theRealShowerSubtraction; } /** * Indicate that the shower virtual contribution should be subtracted. */ void doVirtualShowerSubtraction(); /** * Return true, if the shower virtual contribution should be subtracted. */ bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; } /** * Indicate that the loopsim matched virtual contribution should be subtracted. */ void doLoopSimSubtraction(); /** * Return true, if the loopsim matched virtual contribution should be subtracted. */ bool loopSimSubtraction() const { return theLoopSimSubtraction; } /** * Return true, if this configuration of cross sections should not * be included due to their relative magnitude. Arguments are head * cross section and dependent cross section, including all * reweights. */ virtual bool discard(const CrossSection&, const CrossSection&) const { return false; } //@} /** @name Matrix element and dipole information */ //@{ /** * Return the subtraction dipoles. */ vector::ptr> dipoles(); /** * Return the underlying born matrix elements. */ const vector::ptr>& borns() const; /** * Access the underlying born matrix elements, * overriding the ones contained in the factory object. */ void setBorns(const vector::ptr>& newBorns) { theBorns = newBorns; } /** * Build up dipoles needed. */ void getDipoles(); /** * Clone all dipoles. */ void cloneDipoles(const string& prefix = ""); /** * Clone the real emission matrix element. */ void cloneRealME(const string& prefix = ""); /** * Clone all dependencies. */ void cloneDependencies(const string& prefix = "") { cloneDipoles(prefix); cloneRealME(prefix); } /** * Return all dipoles matching the given Born process */ vector::ptr> splitDipoles(const cPDVector&); //@} /** @name Veto scale settings */ //@{ /** * Set veto scales on the particles at the given * SubProcess which has been generated using this * matrix element. */ virtual void setVetoScales(tSubProPtr) const; //@} /** @name Diagnostic information */ //@{ /** * Dump the setup to an ostream */ void print(ostream&) const; /** * Collect information on the last evaluated phasespace * point for verification or debugging purposes. This * only called, if the StdXCombGroup did accumulate * a non-zero cross section from this ME group. */ virtual void lastEventStatistics(); /** * Print debug information on the last event */ void printLastEvent(ostream&) const; /** * Check the subtraction for the last event */ void lastEventSubtraction(); /** * Return true, if verbose */ bool verbose() const; /** * Return true, if verbose */ bool initVerbose() const; //@} /** @name Setup of Subtracted ME objects */ //@{ /** * Return true if this object needs to be initialized before all * other objects (except those for which this function also returns * true). This default version always returns false, but subclasses * may override it to return true. */ virtual bool preInitialize() const { return true; } /** * Simple envelope histogram to keep track of subtraction */ struct SubtractionHistogram { /** * The lower bound */ double lower; /** * The bins, indexed by upper bound. */ map > bins; /** * Constructor */ SubtractionHistogram(double low = 0.001, double up = 10., unsigned int nbins = 100); /** * Book an event. */ void book(double inv, double diff) { map >::iterator b = bins.upper_bound(inv); if ( b == bins.end() ) return; b->second.first = min(b->second.first,diff); b->second.second = max(b->second.second,diff); } /** * Write to file given name and invariant. */ void dump(const std::string& prefix, const int& plottype, const bool& scatterplot, const cPDVector& proc, int i, int j) const; /** * Write to persistent ostream */ void persistentOutput(PersistentOStream&) const; /** * Read from persistent istream */ void persistentInput(PersistentIStream&); }; //@} 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 underlying born matrix elements, overriding the ones * contained in the factory object. */ vector::ptr> theBorns; /** * Pointer to the head real emission ME casted to a MatchboxMEBase * object. */ Ptr::ptr theReal; /** * Define the key for the collinear subtraction data. */ typedef pair > CollinearSubtractionIndex; /** * subtraction data for collinear limits. */ map collinearHistograms; /** * names of files to which subtraction data is written for all phase space points individually */ map fnamesCollinearSubtraction; /** * Define the key for the soft subtraction data. */ typedef pair SoftSubtractionIndex; /** * subtraction data for soft limits. */ map softHistograms; /** * names of files to which subtraction data is written for all phase space points individually */ map fnamesSoftSubtraction; /** * True, if the shower real emission contribution should be subtracted. */ bool theRealShowerSubtraction; /** * True, if the shower virtual contribution should be subtracted. */ bool theVirtualShowerSubtraction; /** * True, if the loopsim matched virtual contribution should be subtracted. */ bool theLoopSimSubtraction; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ SubtractedME & operator=(const SubtractedME &) = delete; }; inline PersistentOStream& operator<<(PersistentOStream& os, const SubtractedME::SubtractionHistogram& h) { h.persistentOutput(os); return os; } inline PersistentIStream& operator>>(PersistentIStream& is, SubtractedME::SubtractionHistogram& h) { h.persistentInput(is); return is; } } #endif /* HERWIG_SubtractedME_H */