diff --git a/MatrixElement/Matchbox/Utility/DensityOperator.cc b/MatrixElement/Matchbox/Utility/DensityOperator.cc --- a/MatrixElement/Matchbox/Utility/DensityOperator.cc +++ b/MatrixElement/Matchbox/Utility/DensityOperator.cc @@ -1,479 +1,477 @@ // -*- C++ -*- // // DensityOperator.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 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the DensityOperator class. // #include "DensityOperator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXCombData.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; DensityOperator::DensityOperator() : Nc(3.0), TR(0.5) { } -DensityOperator::~DensityOperator() {} - void DensityOperator::clear() { theCorrelatorMap.clear(); } // Prepare density operator if not done before void DensityOperator::prepare(const cPDVector& mePartonData) { // Take parton data and create key of type 33bar888 to use as key for basis vector mePartonDataColoured = theColourBasis->normalOrderMap(mePartonData); if( theDensityOperatorMap.count( mePartonDataColoured ) == 0 ){ // Get basis dimension size_t dim = theColourBasis->prepare( mePartonData, false ); // Allocate space for density matrix theDensityOperatorMap.insert(make_pair(mePartonDataColoured,matrix (dim,dim))); } } // Fill the density matrix for the first time void DensityOperator::fill(const Ptr::ptr MBXCombPtr, const cPDVector& partons, const vector& momenta){ // Map from helicity structure to amplitude vector in the color basis const map,CVector>& amplitudeMap = MBXCombPtr->lastAmplitudes(); const cPDVector& mePartonData = MBXCombPtr->mePartonData(); // Get the dimension of the basis for the indexChange method size_t dim = theColourBasis->prepare(mePartonData,false); // Normal order partons according to ColourBasis const vector mePartonDataColoured = theColourBasis->normalOrderMap(mePartonData); // Get the colour basis to colour basis map for this hard subprocess map >::iterator cb2cbit = theColourBasisToColourBasisMap.find(mePartonData); // Fill the map with this hard subprocess if it doesn't have it yet if ( cb2cbit == theColourBasisToColourBasisMap.end() ) { // Get the index map for the partons as they are ordered in the MatchboxXComb // object const map& mbIndexMap = theColourBasis->indexMap().at(mePartonData); // Get the index map for the partons as they are ordered // in the shower. // Use prepare, as it might not have been done for this order of the partons theColourBasis->prepare(partons,false); const map& showerIndexMap = theColourBasis->indexMap().at(partons); // ensure the maps are of the same size assert( mbIndexMap.size() == showerIndexMap.size() ); // Loop over both sets of partons to determine how the order between them // differs, then translate the index to the colour basis and put the // key-value pair into the map map cb2cbMap; // Get the momenta for comparison const vector& meMomenta = MBXCombPtr->matchboxME()->lastMEMomenta(); // Make sure there's the same number of momenta in both vectors assert( momenta.size() == meMomenta.size() ); bool done; // Boost the momenta to the same frame const vector momentaCM = boostToRestFrame(momenta); for ( size_t i = 0; i < momentaCM.size(); i++ ) { // The cb2cb map is intended to translate how the indices // are different in the colour basis due to the different // orderings of the partons, so it should only bother // with coloured particles. if ( partons[i]->coloured() ) { done = false; for ( size_t j = 0; j < meMomenta.size(); j++ ) { if ( !done ) { if ( compareMomentum(momentaCM[i],meMomenta[j]) ) { cb2cbMap[mbIndexMap.at(i)] = showerIndexMap.at(j); done = true; } } } } } // Make sure all momenta have been identified assert( cb2cbMap.size() == mePartonDataColoured.size() ); // Add the map to the cache theColourBasisToColourBasisMap[mePartonData] = cb2cbMap; // Get the iterator cb2cbit = theColourBasisToColourBasisMap.find(mePartonData); } // With the cb2cb index map we can create the basis vector index map const map vectorMap = theColourBasis->indexChange(mePartonDataColoured, dim,cb2cbit->second); // Prepare density operator (allocate) for set of particles prepare(mePartonData); // Check that density operator (place holder for) exist in density operator map map,matrix >::iterator dOit = theDensityOperatorMap.find(mePartonDataColoured); assert(dOit != theDensityOperatorMap.end()); // Initialize the density operator matrix& densOp = dOit->second; for(unsigned int i = 0; i < (amplitudeMap.begin()->second).size(); i++){ for(unsigned int j = 0; j < (amplitudeMap.begin()->second).size(); j++){ densOp (i,j) = 0; } } // Fill the density operator, ok to sum over helicities since the density operator // exist at the probability level CVector amplitude; for ( map,CVector>::const_iterator itHel = amplitudeMap.begin(); itHel != amplitudeMap.end(); itHel++ ) { amplitude = itHel->second; for ( unsigned int i = 0; i < amplitude.size(); i++ ) { for ( unsigned int j = 0; j < amplitude.size(); j++ ) { // vectorMap is used such that densOp is filled according to the // basis order defined by the input cPDVector partons. densOp (vectorMap.at(i),vectorMap.at(j)) += amplitude(i)*std::conj(amplitude(j)); } } } // Colour conservation check colourConservation(mePartonData); } // Update density matrix after emitting or splitting a gluon // The emissionsMap argument contains the relation of the indices in the smaller (before) // and larger basis (after). the first 3-tuple contains the indices // (emitter before, emitter after, emitted parton) // The map contains the old and new indices of all other partons (not involved) void DensityOperator::evolve(const map,Complex>& Vijk, const cPDVector& before, const cPDVector& after, const map,map >& emissionsMap, const bool splitAGluon, const bool initialGluonSplitting) { size_t dimBefore = theColourBasis->prepare( before, false ); size_t dimAfter = theColourBasis->prepare( after, false ); vector beforeColoured = theColourBasis->normalOrderMap(before); vector afterColoured = theColourBasis->normalOrderMap(after); const map,matrix >::iterator dOit = theDensityOperatorMap.find(beforeColoured); assert(dOit != theDensityOperatorMap.end()); const matrix& densOpBefore = dOit->second; prepare(after); matrix& densOpAfter = theDensityOperatorMap[afterColoured]; for(size_t i = 0; i < densOpAfter.size1(); i++){ for(size_t j = 0; j < densOpAfter.size2(); j++){ densOpAfter(i,j) = 0; } } compressed_matrix Tij; matrix TijMn (dimAfter,dimBefore); compressed_matrix Tk; matrix TijMnTkdagger (dimAfter,dimAfter); Complex V; // Compensate for sign from updated density matrix // TODO Check signs again double sign = -1.0; if ( splitAGluon ) sign = 1.0; // Loop over emitter legs ij and recoil legs k and add the contribution, // for Vijk is assumed to contain the factor 4*pi*\alpha_s/pi.pj typedef map,map > dictMap; int ij,k; // Loop over emitters for(dictMap::const_iterator ijit = emissionsMap.begin(); ijit != emissionsMap.end(); ijit++) { // get first element in 3-tuple, i.e., emitter before ij = std::get<0>(ijit->first); assert(before[ij]->coloured()); // Get rectangular matrices T_{ij} or S_{ij} taking us from the // smaller to the larger basis depending on // the involved particles before and after, the emitter index before ij, // the emitter after and the emitted parton index int i=std::get<1>(ijit->first); // emitter after int j=std::get<2>(ijit->first);// emitted parton Tij = theColourBasis->charge(before,after, ij,i,j,ijit->second).first; TijMn = prodSparseDense(Tij,densOpBefore); // Loop over spectators for(dictMap::const_iterator kit = emissionsMap.begin(); kit != emissionsMap.end(); kit++) { k = std::get<0>(kit->first); assert(before[k]->coloured()); // Standard case of gluon radiation if ( ijit != kit || splitAGluon || initialGluonSplitting ) { int k_after=std::get<1>(kit->first); // For color structure k now has role of emitter int k_emission= std::get<2>(kit->first); // Emitted parton index Tk = theColourBasis->charge(before,after, k, k_after, k_emission, kit->second).first; TijMnTkdagger = prodDenseSparse(TijMn,Tk); // sign == -1.0 if it isn't a gluon splitting into a qqbar pair V = sign*(1.0/colourNorm(before[ij]))* Vijk.at(make_pair(ij,k)); densOpAfter += V*TijMnTkdagger; } } } // Check that the density operator does not vanish assert( theColourBasis->me2(after,densOpAfter) != 0.0 ); colourConservation(after); } // The 3-tuple contains (emitter index, spectator index, emission pid) double DensityOperator::colourMatrixElementCorrection(const std::tuple& ikemission, const cPDVector& particles) { const int i = std::get<0>(ikemission); const int k = std::get<1>(ikemission); //const long emissionID = std::get<2>(ikemission); // Get the density operator // normal order particles as in ColourBasis vector particlesColoured = theColourBasis->normalOrderMap(particles); // ... and find corresponding density operator const map,matrix >::iterator particlesit = theDensityOperatorMap.find(particlesColoured); assert(particlesit != theDensityOperatorMap.end()); const matrix& densOp = particlesit->second; double Ti2 = colourNorm(particles[i]); // Emitter-spectator pair const pair ik = make_pair(i,k); // Create key for color matrix element correction map pair,pair > particlesAndLegs = make_pair(particlesColoured,ik); // Result double res = 0; // Check if it has already been calculated // TODO: move this check earlier (we only need particlesColoured to check if it has been // calculated). // Check for color matrix element associated with key, and calculate if not done const map,pair >,double >::const_iterator corrit = theCorrelatorMap.find(particlesAndLegs); if ( corrit == theCorrelatorMap.end() ) { double corrME2 = theColourBasis->colourCorrelatedME2(ik,particles,densOp); double me2 = theColourBasis->me2(particles,densOp); res = -(1/Ti2)*corrME2/me2; if ( particles[i]->id() == ParticleID::g ) res *= 2.; theCorrelatorMap.insert(make_pair(particlesAndLegs,res)); } else { res = corrit->second; } return res; } double DensityOperator::colourNorm(const cPDPtr particle) { if ( particle->id() == ParticleID::g ) { return Nc; //is 3.0 for Nc = 3, TR = 1/2 } else if ( particle->iColour() == PDT::Colour3 || particle->iColour() == PDT::Colour3bar ) { return TR*(Nc*Nc-1.)/Nc; // is 4.0/3.0 for Nc = 3, TR = 1/2 } else { throw Exception() << "Colour matrix element corrections only work " << "on quark and gluon legs. " << Exception::runerror; } } void DensityOperator::colourConservation(const cPDVector& particles) { // To contain (emitter, spectator, emission pid) std::tuple ikemission; // Normal order particles as defined in ColourBasis const vector particlesColoured = theColourBasis->normalOrderMap(particles); vector sum(particlesColoured.size(),0.0); size_t iterm = 0; // To compensate for the CMEC having a 1/(1+\delta(i is gluon)) factor // in the splitting kernel double gluonFactor = 1.0; // Loop over "emitters" to check color conservation for for ( size_t i = 0; i < particles.size(); i++ ) { if ( particles[i]->coloured() ) { for ( size_t k = 0; k < particles.size(); k++ ) { if ( particles[k]->coloured() && i != k ) { ikemission = std::make_tuple(i,k,ParticleID::g); if ( particles[i]->id() != ParticleID::g ) { gluonFactor = 1.0; } else { gluonFactor = 1./2.; } sum[iterm] += gluonFactor*colourMatrixElementCorrection(ikemission,particles); } } iterm++; } } for ( size_t i = 0; i < sum.size(); i++ ) assert( std::abs(sum[i]-1.0) < pow(10.0,-10.0)); } matrix DensityOperator::prodSparseDense(const compressed_matrix& Tij, const matrix& Mn){ // Dimension before emission size_t dimBefore = Tij.size2(); // Dimension after emission size_t dimAfter = Tij.size1(); //Check matrix dimensions assert( dimBefore == Mn.size1() ); // Allocate memory for the matrix matrix TijMn (dimAfter,dimBefore,0); // Use iterators for the compressed matrix to iterate only over the non-zero // elements. size_t ii; size_t jj; for ( compressed_matrix::const_iterator1 it1 = Tij.begin1(); it1 != Tij.end1(); it1++ ) { for ( compressed_matrix::const_iterator2 it2 = it1.begin(); it2 != it1.end(); it2++ ) { ii = it2.index1(); jj = it2.index2(); for ( size_t kk = 0; kk < dimBefore; kk++ ) { // *it2 is Tij(ii,jj) TijMn(ii,kk) += (*it2)*Mn(jj,kk); } } } return TijMn; } matrix DensityOperator::prodDenseSparse(const matrix& TijMn, const compressed_matrix& Tk){ // The compressed matrix comes from the charge method, do not transpose yet // Dimension after emission size_t dimAfter = Tk.size1();//Since this method returns TijMn*Tk^\dagger //Check matrix dimensions assert( TijMn.size2() == Tk.size2() ); // Allocate memory for the matrix matrix TijMnTkdagger (dimAfter,dimAfter,0); size_t jj; size_t kk; for ( compressed_matrix::const_iterator1 it1 = Tk.begin1(); it1 != Tk.end1(); it1++ ) { for ( compressed_matrix::const_iterator2 it2 = it1.begin(); it2 != it1.end(); it2++ ) { jj = it2.index2();//transposing Tk jj is index2(), not index1() kk = it2.index1();//transposing Tk for ( size_t ii = 0; ii < dimAfter; ii++ ) { // *it2 is Tk(kk,jj) = trans(Tk)(jj,kk) TijMnTkdagger(ii,kk) += TijMn(ii,jj)*(*it2); } } } return TijMnTkdagger; } vector DensityOperator::boostToRestFrame(const vector& momenta) { // We need 2 initial particles assert(momenta.size() >= 2); // The boosted vectors vector vboosted = momenta; // The boost should be to the rest frame of the initial particles Boost b = (momenta[0] + momenta[1]).findBoostToCM(); // Boost all of the vectors for ( size_t i = 0; i < momenta.size(); i++ ) { vboosted[i].boost(b); } return vboosted; } bool DensityOperator::compareMomentum(const Lorentz5Momentum& p, const Lorentz5Momentum& q) { bool equal = true; // Compares two momentum vectors p and q, if they are close enough (defined by the // double eps below) it returns true. This is sufficient to distinguish particles // in the matrix element as we are not interested in the matrix element for extremely // collinear radiation (better described by the parton shower). // Create the difference of the two vectors const Lorentz5Momentum l = p - q; // A relevant size that is guaranteed to be larger than 0 const Energy2 e2 = p.t()*p.t(); // Size of the difference that would be considered equal const double eps = pow(10.,-15.); if ( l.x()*l.x()/e2 > eps ) equal = false; if ( l.y()*l.y()/e2 > eps ) equal = false; if ( l.z()*l.z()/e2 > eps ) equal = false; if ( l.t()*l.t()/e2 > eps ) equal = false; return equal; } IBPtr DensityOperator::clone() const { return new_ptr(*this); } IBPtr DensityOperator::fullclone() const { return new_ptr(*this); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void DensityOperator::persistentOutput(PersistentOStream &) const { // *** ATTENTION *** os << ; // Add all member variable which should be written persistently here. } void DensityOperator::persistentInput(PersistentIStream &, int) { // *** ATTENTION *** is >> ; // Add all member variable which should be read persistently here. } // *** 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 describeHerwigDensityOperator("Herwig::DensityOperator", "DensityOperator.so"); void DensityOperator::Init() { static ClassDocumentation documentation ("There is no documentation for the DensityOperator class"); } diff --git a/MatrixElement/Matchbox/Utility/DensityOperator.h b/MatrixElement/Matchbox/Utility/DensityOperator.h --- a/MatrixElement/Matchbox/Utility/DensityOperator.h +++ b/MatrixElement/Matchbox/Utility/DensityOperator.h @@ -1,248 +1,240 @@ // -*- C++ -*- // // DensityOperator.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 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef Herwig_DensityOperator_H #define Herwig_DensityOperator_H // // This is the declaration of the DensityOperator class. // #include "ThePEG/Handlers/HandlerBase.h" #include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h" #include #include #include namespace Herwig { using namespace ThePEG; typedef boost::numeric::ublas::vector CVector; /** * Here is the documentation of the DensityOperator class. * * @see \ref DensityOperatorInterfaces "The interfaces" * defined for DensityOperator. */ class DensityOperator: public HandlerBase { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ DensityOperator(); - /** - * The destructor. - */ - virtual ~DensityOperator(); - //@} - public: /** * Clears theDensityOperatorMap. */ void clear(); /** * Prepare for the given sub process. */ void prepare(const cPDVector&); /** * Fill the density operator for the given hard subprocess, summing over all * helicity configurations. */ void fill(const Ptr::ptr, const cPDVector&, const vector& momenta); /** * Evolve the density operator, by * M_{n+1} = -\sum_{i,k}{-4*pi*alpha_s/Ti2*V_{ij,k} T_{i,n}M_nT_{k,n}^\dag}, * see arXiv:1206.0180 eq. (5), note that the pi*pj factor is assumed to be * included in V_{ij,k}. */ void evolve(const map,Complex>& Vijk, const cPDVector& before, const cPDVector& after, const map,map >& emissionsMap, const bool splitAGluon, const bool initialGluonSplitting); /** * Calculate the colour matrix element correction. * -(1+delta(i,gluon))/Ti^2 Tr(Sn+1 Ti Mn Tk^dagger)/Tr(Sn Mn) * where the bracket in front compensates for the gluon symmetry factor, * Ti^2 is C_f or C_a, Sn+1 is the matrix of scalar products, and * Ti is the radiation matrix. * The first arg contains (emitter index, spectator index, emission pid) * */ double colourMatrixElementCorrection(const std::tuple& ikemission, const cPDVector& particles); /** * Checking colour conservation for the colour matrix element corrections. */ void colourConservation(const cPDVector& particles); /** * Get the colour basis. */ Ptr::tptr colourBasis() { return theColourBasis; } /** * Get the colour basis. */ const Ptr::tptr colourBasis() const { return theColourBasis; } /** * Set the colour basis. */ void colourBasis(Ptr::ptr ptr) { theColourBasis = ptr; } /** * Get the correlator map. */ const map,pair >,double>& correlatorMap() const { return theCorrelatorMap; } /** @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: /** * Number of colours used in colourNorm. */ double Nc; /** * QCD vertex normalization. */ double TR; /** * Normalization of colour charges \mathbf{T}_{ij}^2. */ double colourNorm(const cPDPtr particle); /** * Fast evaluation of Tij*Mn, where a Tij is the matrix from ColourBasis::charge, * which is a sparse matrix, and Mn is the density operator, a dense matrix. * */ matrix prodSparseDense(const compressed_matrix&, const matrix&); /** * Fast evaluation of TijMn*Tkdagger, where a TijMn is the result from the method * prodSparseDense, a dense matrix, and Tkdagger is the transponse conjugate of * the matrix from ColourBasis::charge, a sparse matrix. * */ matrix prodDenseSparse(const matrix&, const compressed_matrix&); /** * Boosts a vector of momenta to the rest frame of the initial pair * of particles (the first 2 elements of the argument vector). Returns * the boosted vectors */ vector boostToRestFrame(const vector& momenta); /** * Boosts a vector of momenta to the rest frame of the initial pair */ bool compareMomentum(const Lorentz5Momentum& p, const Lorentz5Momentum& q); /** * Mapping of colour structures to density operator matrices. * */ map,matrix > theDensityOperatorMap; /** * Mapping of colour structures and legs to colour correlators. */ map,pair >,double> theCorrelatorMap; /** * A map from the hard subprocess particles to a map of amplitude colour * basis order to the normal ordered colour basis. */ map > theColourBasisToColourBasisMap; /** * Colour basis used. */ Ptr::ptr theColourBasis; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ DensityOperator & operator=(const DensityOperator &) = delete; }; } #endif /* Herwig_DensityOperator_H */ diff --git a/MatrixElement/Matchbox/Utility/MatchboxFactoryMatcher.cc b/MatrixElement/Matchbox/Utility/MatchboxFactoryMatcher.cc --- a/MatrixElement/Matchbox/Utility/MatchboxFactoryMatcher.cc +++ b/MatrixElement/Matchbox/Utility/MatchboxFactoryMatcher.cc @@ -1,102 +1,100 @@ // -*- C++ -*- // // MatchboxFactoryMatcher.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 MatchboxFactoryMatcher class. // #include "MatchboxFactoryMatcher.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; MatchboxFactoryMatcher::MatchboxFactoryMatcher() : theGroup("") {} -MatchboxFactoryMatcher::~MatchboxFactoryMatcher() {} - IBPtr MatchboxFactoryMatcher::clone() const { return new_ptr(*this); } IBPtr MatchboxFactoryMatcher::fullclone() const { return new_ptr(*this); } PMPtr MatchboxFactoryMatcher::pmclone() const { return new_ptr(*this); } bool MatchboxFactoryMatcher::check(const ParticleData & data) const { return theIds.find(data.id()) != theIds.end(); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void MatchboxFactoryMatcher::persistentOutput(PersistentOStream & os) const { os << theGroup << theIds; } void MatchboxFactoryMatcher::persistentInput(PersistentIStream & is, int) { is >> theGroup >> theIds; } void MatchboxFactoryMatcher::doinit() { if ( !MatchboxFactory::isMatchboxRun() ) return; if ( !MatchboxFactory::currentFactory() ) throw Exception() << "MatchboxFactoryMatcher::doinit(): No factory object is available in the matcher '" << name() << "'" << Exception::runerror; map::const_iterator grp = MatchboxFactory::currentFactory()->particleGroups().find(theGroup); if ( grp == MatchboxFactory::currentFactory()->particleGroups().end() ) throw Exception() << "MatchboxFactoryMatcher::doinit(): Particle group '" << theGroup << "' not defined in factory object '" << MatchboxFactory::currentFactory()->name() << "'" << Exception::runerror; theIds.clear(); for ( PDVector::const_iterator p = grp->second.begin(); p != grp->second.end(); ++p ) theIds.insert((**p).id()); MatcherBase::doinit(); } // *** 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 describeHerwigMatchboxFactoryMatcher("Herwig::MatchboxFactoryMatcher", "Herwig.so"); void MatchboxFactoryMatcher::Init() { static ClassDocumentation documentation ("MatchboxFactoryMatcher matches particles according to MatchboxFactory particle groups"); static Parameter interfaceGroup ("Group", "Set the group name to match.", &MatchboxFactoryMatcher::theGroup, "", false, false); } diff --git a/MatrixElement/Matchbox/Utility/MatchboxFactoryMatcher.h b/MatrixElement/Matchbox/Utility/MatchboxFactoryMatcher.h --- a/MatrixElement/Matchbox/Utility/MatchboxFactoryMatcher.h +++ b/MatrixElement/Matchbox/Utility/MatchboxFactoryMatcher.h @@ -1,147 +1,139 @@ // -*- C++ -*- // // MatchboxFactoryMatcher.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_MatchboxFactoryMatcher_H #define Herwig_MatchboxFactoryMatcher_H // // This is the declaration of the MatchboxFactoryMatcher class. // #include "ThePEG/PDT/MatcherBase.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief MatchboxFactoryMatcher matches particles according to * MatchboxFatory particle groups * * @see \ref MatchboxFactoryMatcherInterfaces "The interfaces" * defined for MatchboxFactoryMatcher. */ class MatchboxFactoryMatcher: public ThePEG::MatcherBase { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ MatchboxFactoryMatcher(); - /** - * The destructor. - */ - virtual ~MatchboxFactoryMatcher(); - //@} - public: /** * Check if a particle type meets the criteria. */ virtual bool check(const ParticleData &) const; /** * Specialized clone method for MatcherBase used by the * Repository. A sub class must make sure that also the MatcherBase * object corresponding to the complex conjugate of this is cloned. */ virtual PMPtr pmclone() const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // 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(); //@} private: /** * The particle group to be matched */ string theGroup; /** * The set of particle ids to be matched */ set theIds; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MatchboxFactoryMatcher & operator=(const MatchboxFactoryMatcher &) = delete; }; } #endif /* Herwig_MatchboxFactoryMatcher_H */ diff --git a/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.cc b/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.cc --- a/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.cc +++ b/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.cc @@ -1,83 +1,81 @@ // -*- C++ -*- // // MatchboxScaleChoice.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 MatchboxScaleChoice class. // #include "MatchboxScaleChoice.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; MatchboxScaleChoice::MatchboxScaleChoice() : theFixedScale(ZERO), theFixedQEDScale(ZERO) {} -MatchboxScaleChoice::~MatchboxScaleChoice() {} - IBPtr MatchboxScaleChoice::clone() const { return new_ptr(*this); } IBPtr MatchboxScaleChoice::fullclone() const { return new_ptr(*this); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void MatchboxScaleChoice::persistentOutput(PersistentOStream & os) const { os << theLastXComb << ounit(theFixedScale,GeV) << ounit(theFixedQEDScale,GeV); } void MatchboxScaleChoice::persistentInput(PersistentIStream & is, int) { is >> theLastXComb >> iunit(theFixedScale,GeV) >> iunit(theFixedQEDScale,GeV); } // *** Attention *** The following static variable is needed for the type // description system in ThePEG. Please check that the template arguments // are correct (the class and its base class), and that the constructor // arguments are correct (the class name and the name of the dynamically // loadable library where the class implementation can be found). DescribeClass describeHerwigMatchboxScaleChoice("Herwig::MatchboxScaleChoice", "Herwig.so"); void MatchboxScaleChoice::Init() { static ClassDocumentation documentation ("MatchboxScaleChoice is the base class for scale choices " "within Matchbox."); static Parameter interfaceFixedScale ("FixedScale", "Set a fixed scale.", &MatchboxScaleChoice::theFixedScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Parameter interfaceFixedQEDScale ("FixedQEDScale", "Set a fixed QED scale.", &MatchboxScaleChoice::theFixedQEDScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); } diff --git a/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h b/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h --- a/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h +++ b/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h @@ -1,169 +1,161 @@ // -*- C++ -*- // // MatchboxScaleChoice.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_MatchboxScaleChoice_H #define Herwig_MatchboxScaleChoice_H // // This is the declaration of the MatchboxScaleChoice class. // #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Handlers/LastXCombInfo.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief MatchboxScaleChoice is the base class for scale choices * within Matchbox. * */ class MatchboxScaleChoice: public HandlerBase, public LastXCombInfo { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ MatchboxScaleChoice(); - /** - * The destructor. - */ - virtual ~MatchboxScaleChoice(); - //@} - public: /** * Clone this scale choice. */ Ptr::ptr cloneMe() const { return dynamic_ptr_cast::ptr>(clone()); } /** * Set the XComb object. */ virtual void setXComb(tStdXCombPtr xc) { theLastXComb = xc; } /** * Return the renormalization scale. This default version returns * shat. */ virtual Energy2 renormalizationScale() const { return theFixedScale == ZERO ? lastSHat() : sqr(theFixedScale); } /** * Return the factorization scale. This default version returns * shat. */ virtual Energy2 factorizationScale() const { return theFixedScale == ZERO ? lastSHat() : sqr(theFixedScale); } /** * Return the QED renormalization scale. This default version returns * the Z mass squared. */ virtual Energy2 renormalizationScaleQED() const { if ( theFixedQEDScale != ZERO ) return sqr(theFixedQEDScale); Energy mZ = getParticleData(ParticleID::Z0)->hardProcessMass(); return mZ*mZ; } /** * Return the shower hard scale. This default implementation returns the * factorization scale. */ virtual Energy2 showerScale() const { return factorizationScale(); } 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: /** * A fixed scale choice. If zero, shat will be used. */ Energy theFixedScale; /** * A fixed QED scale choice. If zero, shat will be used. */ Energy theFixedQEDScale; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MatchboxScaleChoice & operator=(const MatchboxScaleChoice &) = delete; }; } #endif /* Herwig_MatchboxScaleChoice_H */ diff --git a/MatrixElement/Matchbox/Utility/MatchboxXComb.cc b/MatrixElement/Matchbox/Utility/MatchboxXComb.cc --- a/MatrixElement/Matchbox/Utility/MatchboxXComb.cc +++ b/MatrixElement/Matchbox/Utility/MatchboxXComb.cc @@ -1,78 +1,76 @@ // -*- C++ -*- // // MatchboxXComb.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 MatchboxXComb class. // #include "MatchboxXComb.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; MatchboxXComb::MatchboxXComb() : StandardXComb() { flushCaches(); } -MatchboxXComb::~MatchboxXComb() {} - MatchboxXComb::MatchboxXComb(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, tMEPtr newME, const DiagramVector & newDiagrams, bool mir, tStdXCombPtr newHead) : StandardXComb(newMaxEnergy, inc, newEventHandler,newSubProcessHandler, newExtractor, newCKKW, newPartonBins, newCuts, newME, newDiagrams, mir, newHead), MatchboxXCombData(newME) { flushCaches(); } MatchboxXComb::MatchboxXComb(tStdXCombPtr newHead, const PBPair & newPartonBins, tMEPtr newME, const DiagramVector & newDiagrams) : StandardXComb(newHead, newPartonBins, newME, newDiagrams), MatchboxXCombData(newME) { flushCaches(); } void MatchboxXComb::clean() { StandardXComb::clean(); flushCaches(); } void MatchboxXComb::persistentOutput(PersistentOStream & os) const { MatchboxXCombData::persistentOutput(os); } void MatchboxXComb::persistentInput(PersistentIStream & is, int version) { MatchboxXCombData::persistentInput(is,version); } // *** 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 describeHerwigMatchboxXComb("Herwig::MatchboxXComb", "Herwig.so"); void MatchboxXComb::Init() { static ClassDocumentation documentation ("MatchboxXComb extents StandardXComb."); } diff --git a/MatrixElement/Matchbox/Utility/MatchboxXComb.h b/MatrixElement/Matchbox/Utility/MatchboxXComb.h --- a/MatrixElement/Matchbox/Utility/MatchboxXComb.h +++ b/MatrixElement/Matchbox/Utility/MatchboxXComb.h @@ -1,108 +1,102 @@ // -*- C++ -*- // // MatchboxXComb.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_MatchboxXComb_H #define Herwig_MatchboxXComb_H // // This is the declaration of the MatchboxXComb class. // #include "ThePEG/Handlers/StandardXComb.h" #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXCombData.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief Matchbox extensions to StandardXComb */ class MatchboxXComb: public StandardXComb, public MatchboxXCombData { public: /** @name Standard constructors and destructors. */ //@{ /** * Standard constructor. */ MatchboxXComb(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, tMEPtr newME, const DiagramVector & newDiagrams, bool mir, tStdXCombPtr newHead = tStdXCombPtr()); /** * Constructor given a head xcomb. */ MatchboxXComb(tStdXCombPtr newHead, const PBPair & newPartonBins, tMEPtr newME, const DiagramVector & newDiagrams); /** * Default constructor. */ MatchboxXComb(); - - /** - * Destructor. - */ - virtual ~MatchboxXComb(); - //@} public: /** * Reset all saved data about last generated phasespace point; */ virtual void clean(); 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(); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MatchboxXComb & operator=(const MatchboxXComb &) = delete; }; } #endif /* Herwig_MatchboxXComb_H */ diff --git a/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.cc b/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.cc --- a/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.cc +++ b/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.cc @@ -1,70 +1,68 @@ // -*- C++ -*- // // MatchboxXCombGroup.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 MatchboxXCombGroup class. // #include "MatchboxXCombGroup.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; MatchboxXCombGroup::MatchboxXCombGroup() : StdXCombGroup() { flushCaches(); } -MatchboxXCombGroup::~MatchboxXCombGroup() {} - MatchboxXCombGroup::MatchboxXCombGroup(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, tMEGroupPtr newME, const DiagramVector & newDiagrams, bool mir, tStdXCombPtr newHead) : StdXCombGroup(newMaxEnergy, inc, newEventHandler, newSubProcessHandler, newExtractor, newCKKW, newPartonBins, newCuts, newME, newDiagrams, mir, newHead), MatchboxXCombData(newME) { flushCaches(); } void MatchboxXCombGroup::clean() { StdXCombGroup::clean(); flushCaches(); } void MatchboxXCombGroup::persistentOutput(PersistentOStream & os) const { MatchboxXCombData::persistentOutput(os); } void MatchboxXCombGroup::persistentInput(PersistentIStream & is, int version) { MatchboxXCombData::persistentInput(is,version); } // *** 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 describeHerwigMatchboxXCombGroup("Herwig::MatchboxXCombGroup", "Herwig.so"); void MatchboxXCombGroup::Init() { static ClassDocumentation documentation ("MatchboxXCombGroup extents StdXCombGroup."); } diff --git a/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.h b/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.h --- a/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.h +++ b/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.h @@ -1,102 +1,96 @@ // -*- C++ -*- // // MatchboxXCombGroup.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_MatchboxXCombGroup_H #define Herwig_MatchboxXCombGroup_H // // This is the declaration of the MatchboxXCombGroup class. // #include "ThePEG/Handlers/StdXCombGroup.h" #include "ThePEG/MatrixElement/MEGroup.h" #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXCombData.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief Matchbox extensions to StandardXComb */ class MatchboxXCombGroup: public StdXCombGroup, public MatchboxXCombData { public: /** @name Standard constructors and destructors. */ //@{ /** * Standard constructor. */ MatchboxXCombGroup(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, tMEGroupPtr newME, const DiagramVector & newDiagrams, bool mir, tStdXCombPtr newHead = tStdXCombPtr()); /** * Default constructor. */ MatchboxXCombGroup(); - - /** - * Destructor. - */ - virtual ~MatchboxXCombGroup(); - //@} public: /** * Reset all saved data about last generated phasespace point; */ virtual void clean(); 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(); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MatchboxXCombGroup & operator=(const MatchboxXCombGroup &) = delete; }; } #endif /* Herwig_MatchboxXCombGroup_H */ diff --git a/MatrixElement/Matchbox/Utility/SimpleColourBasis.cc b/MatrixElement/Matchbox/Utility/SimpleColourBasis.cc --- a/MatrixElement/Matchbox/Utility/SimpleColourBasis.cc +++ b/MatrixElement/Matchbox/Utility/SimpleColourBasis.cc @@ -1,414 +1,409 @@ // -*- C++ -*- // // SimpleColourBasis.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 SimpleColourBasis class. // #include "SimpleColourBasis.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/StandardModel/StandardModelBase.h" - #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; -SimpleColourBasis::SimpleColourBasis() {} - -SimpleColourBasis::~SimpleColourBasis() {} - IBPtr SimpleColourBasis::clone() const { return new_ptr(*this); } IBPtr SimpleColourBasis::fullclone() const { return new_ptr(*this); } size_t SimpleColourBasis::prepareBasis(const vector& basis) { if ( id33bar.empty() ) makeIds(); if ( basis == id88 || basis == id33bar || basis == id33bar8 ) return 1; if ( basis == id888 || basis == id33bar88 || basis == id33bar33bar ) return 2; if ( basis == id8888 ) return 6; throw Exception() << "SimpleColourBasis::prepareBasis(): Cannot handle colour configuration" << Exception::runerror; return 0; } double SimpleColourBasis::scalarProduct(size_t a, size_t b, const vector& abBasis) const { if ( id33bar.empty() ) makeIds(); double Nc = SM().Nc(); double Nc2 = sqr(Nc); double Nc3 = Nc*Nc2; double Nc4 = sqr(Nc2); double Nc6 = Nc2*Nc4; if ( a > b ) swap(a,b); if ( !largeN() ) { if ( abBasis == id88 ) { return ( Nc2 - 1. )/4.; } if ( abBasis == id33bar ) { return Nc; } if ( abBasis == id888 ) { if ( a == b ) return ( Nc4 - 3.*Nc2 + 2. )/(8.*Nc); return -( Nc2 - 1. )/(4.*Nc); } if ( abBasis == id33bar8 ) { return ( Nc2 - 1. )/2.; } if ( abBasis == id8888 ) { if ( a == b ) return ( Nc6 - 4.*Nc4 + 6.*Nc2 - 3. )/(16.*Nc2); if ( ( a == 0 && b == 1 ) || ( a == 2 && b == 3 ) || ( a == 4 && b == 5 ) ) return ( Nc4 + 2.*Nc2 - 3. )/(16.*Nc2); return -( Nc2 - 4. + 3./Nc2 )/16.; } if ( abBasis == id33bar88 ) { if ( a == b ) return ( Nc4 - 2.*Nc2 + 1 )/(4.*Nc); return -( Nc2 - 1. )/(4.*Nc); } if ( abBasis == id33bar33bar ) { if ( a == b ) return Nc2; return Nc; } } else { if ( a != b ) return 0.; if ( abBasis == id88 ) { return Nc2/4.; } if ( abBasis == id33bar ) { return Nc; } if ( abBasis == id888 ) { return Nc3/8.; } if ( abBasis == id33bar8 ) { return Nc2/2.; } if ( abBasis == id8888 ) { return Nc4/16.; } if ( abBasis == id33bar88 ) { return Nc3/4.; } if ( abBasis == id33bar33bar ) { return Nc2; } } throw Exception() << "SimpleColourBasis::scalarProduct(): Cannot handle colour configuration" << Exception::runerror; } double SimpleColourBasis::tMatrixElement(size_t i, size_t a, size_t b, const vector&, const vector& bBasis, #ifndef NDEBUG size_t k, size_t l, #else size_t , size_t , #endif const map& dict) const { // Check indices k and l assert( k == i ); assert( l == bBasis.size() ); // Check that dict is the standardMap assert( dict.size()+1 == bBasis.size() ); map::const_iterator tmp; for ( size_t ii = 0; ii < bBasis.size(); ii++ ) if ( ii != i ) { tmp = dict.find(ii); assert( tmp != dict.end() ); assert( tmp->second == ii ); } if ( id33bar.empty() ) makeIds(); if ( bBasis == id88 ) { if ( i == 0 ) return a == 0 ? -1. : 1.; else return a == 0 ? 1. : -1.; } if ( bBasis == id33bar ) { return i == 0 ? 1. : -1.; } if ( bBasis == id888 ) { if ( i == 0 ) { if ( a == 3 && b == 0 ) return 1.; if ( a == 0 && b == 0 ) return -1.; if ( a == 1 && b == 1 ) return 1.; if ( a == 2 && b == 1 ) return -1.; } if ( i == 1 ) { if ( a == 4 && b == 0 ) return 1.; if ( a == 3 && b == 0 ) return -1.; if ( a == 2 && b == 1 ) return 1.; if ( a == 5 && b == 1 ) return -1.; } if ( i == 2 ) { if ( a == 0 && b == 0 ) return 1.; if ( a == 4 && b == 0 ) return -1.; if ( a == 5 && b == 1 ) return 1.; if ( a == 1 && b == 1 ) return -1.; } return 0.; } if ( bBasis == id33bar8 ) { if ( i == 0 ) return a == 1 ? 1. : 0.; if ( i == 1 ) return a == 0 ? -1. : 0.; if ( i == 2 ) return a == 0 ? 1. : -1.; } throw Exception() << "SimpleColourBasis::tMatrixElement(): Cannot handle colour configuration" << Exception::runerror; return 0.; } bool SimpleColourBasis::colourConnected(const cPDVector& sub, const vector& basis, const pair& i, const pair& j, size_t a) const { if ( id33bar.empty() ) makeIds(); // translate process to basis ids map >::const_iterator trans = indexMap().find(sub); assert(trans != indexMap().end()); int idColoured = i.second ? j.first : i.first; idColoured = trans->second.find(idColoured)->second; int idAntiColoured = i.second ? i.first : j.first; idAntiColoured = trans->second.find(idAntiColoured)->second; if ( basis == id88 ) { return ( idColoured == 0 && idAntiColoured == 1 ) || ( idColoured == 1 && idAntiColoured == 0 ); } if ( basis == id33bar ) { return idColoured == 0 && idAntiColoured == 1; } if ( basis == id888 ) { if ( a == 0 ) return ( idColoured == 0 && idAntiColoured == 1 ) || ( idColoured == 1 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 0 ); if ( a == 1 ) return ( idColoured == 0 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 1 ) || ( idColoured == 1 && idAntiColoured == 0 ); } if ( basis == id33bar8 ) { return ( idColoured == 0 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 1 ); } if ( basis == id8888 ) { if ( a == 0 ) return ( idColoured == 0 && idAntiColoured == 1 ) || ( idColoured == 1 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 3 ) || ( idColoured == 3 && idAntiColoured == 0 ); if ( a == 1 ) return ( idColoured == 0 && idAntiColoured == 3 ) || ( idColoured == 3 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 1 ) || ( idColoured == 1 && idAntiColoured == 0 ); if ( a == 2 ) return ( idColoured == 0 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 1 ) || ( idColoured == 1 && idAntiColoured == 3 ) || ( idColoured == 3 && idAntiColoured == 0 ); if ( a == 3 ) return ( idColoured == 0 && idAntiColoured == 3 ) || ( idColoured == 3 && idAntiColoured == 1 ) || ( idColoured == 1 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 0 ); if ( a == 4 ) return ( idColoured == 0 && idAntiColoured == 1 ) || ( idColoured == 1 && idAntiColoured == 3 ) || ( idColoured == 3 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 0 ); if ( a == 5 ) return ( idColoured == 0 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 3 ) || ( idColoured == 3 && idAntiColoured == 1 ) || ( idColoured == 1 && idAntiColoured == 0 ); } if ( basis == id33bar88 ) { if ( a == 0 ) return ( idColoured == 0 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 3 ) || ( idColoured == 3 && idAntiColoured == 1 ); if ( a == 1 ) return ( idColoured == 0 && idAntiColoured == 3 ) || ( idColoured == 3 && idAntiColoured == 2 ) || ( idColoured == 2 && idAntiColoured == 1 ); } if ( basis == id33bar33bar ) { if ( a == 0 ) return ( idColoured == 0 && idAntiColoured == 1 ) || ( idColoured == 2 && idAntiColoured == 3 ); if ( a == 1 ) return ( idColoured == 0 && idAntiColoured == 3 ) || ( idColoured == 2 && idAntiColoured == 1 ); } return false; } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void SimpleColourBasis::makeIds() const { id88.push_back(PDT::Colour8); id88.push_back(PDT::Colour8); id33bar.push_back(PDT::Colour3); id33bar.push_back(PDT::Colour3bar); id888.push_back(PDT::Colour8); id888.push_back(PDT::Colour8); id888.push_back(PDT::Colour8); id33bar8.push_back(PDT::Colour3); id33bar8.push_back(PDT::Colour3bar); id33bar8.push_back(PDT::Colour8); id8888.push_back(PDT::Colour8); id8888.push_back(PDT::Colour8); id8888.push_back(PDT::Colour8); id8888.push_back(PDT::Colour8); id33bar88.push_back(PDT::Colour3); id33bar88.push_back(PDT::Colour3bar); id33bar88.push_back(PDT::Colour8); id33bar88.push_back(PDT::Colour8); id33bar33bar.push_back(PDT::Colour3); id33bar33bar.push_back(PDT::Colour3bar); id33bar33bar.push_back(PDT::Colour3); id33bar33bar.push_back(PDT::Colour3bar); } void SimpleColourBasis::persistentOutput(PersistentOStream &) const {} void SimpleColourBasis::persistentInput(PersistentIStream &, int) {} // *** 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 describeHerwigSimpleColourBasis("Herwig::SimpleColourBasis", "Herwig.so"); void SimpleColourBasis::Init() { static ClassDocumentation documentation ("SimpleColourBasis implements the colour algebra needed for " "electroweak boson and electroweak boson + jet production at NLO. It mainly " "serves as an example for the general ColourBasis interface."); } diff --git a/MatrixElement/Matchbox/Utility/SimpleColourBasis.h b/MatrixElement/Matchbox/Utility/SimpleColourBasis.h --- a/MatrixElement/Matchbox/Utility/SimpleColourBasis.h +++ b/MatrixElement/Matchbox/Utility/SimpleColourBasis.h @@ -1,209 +1,194 @@ // -*- C++ -*- // // SimpleColourBasis.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_SimpleColourBasis_H #define Herwig_SimpleColourBasis_H // // This is the declaration of the SimpleColourBasis class. // #include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief SimpleColourBasis implements the colour algebra needed for * electroweak boson and electroweak boson + jet production at NLO. It * mainly serves as an example for the general ColourBasis interface. * */ class SimpleColourBasis: public ColourBasis { public: - /** @name Standard constructors and destructors. */ - //@{ - /** - * The default constructor. - */ - SimpleColourBasis(); - - /** - * The destructor. - */ - virtual ~SimpleColourBasis(); - //@} - -public: - /** * Prepare the basis for the normal ordered legs and return the * dimensionality of the basis. */ virtual size_t prepareBasis(const vector&); /** * Return the scalar product of basis tensors labelled a and b in * the basis used for the given normal ordered legs. */ virtual double scalarProduct(size_t a, size_t b, const vector& abBasis) const; /** * Return the matrix element of a colour charge * between basis tensors a and b, with * respect to aBasis and bBasis */ virtual double tMatrixElement(size_t i, size_t a, size_t b, const vector& aBasis, const vector& bBasis, size_t k, size_t l, const map& dict ) const; virtual double sMatrixElement(size_t, size_t, size_t, const vector&, const vector&, size_t, size_t, const map& ) const { assert( 0 == 1 ); return 0; } /** * Return true, if this colour basis supports gluon splittings. */ virtual bool canSplitGluons() const { return false; } /** * Return true, if a large-N colour connection exists for the * given external legs and basis tensor. */ virtual bool colourConnected(const cPDVector&, const vector&, const pair&, const pair&, size_t) const; /** * Return true, if the colour basis is capable of assigning colour * flows. */ virtual bool haveColourFlows() const { return true; } /** * Create ids for bases */ void makeIds() const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * id for 88 */ mutable vector id88; /** * id for 33bar */ mutable vector id33bar; /** * id for 888 */ mutable vector id888; /** * id for 33bar8 */ mutable vector id33bar8; /** * id for 8888 */ mutable vector id8888; /** * id for 33bar88 */ mutable vector id33bar88; /** * id for 33bar33bar */ mutable vector id33bar33bar; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ SimpleColourBasis & operator=(const SimpleColourBasis &) = delete; }; } #endif /* Herwig_SimpleColourBasis_H */ diff --git a/MatrixElement/Matchbox/Utility/SimpleColourBasis2.cc b/MatrixElement/Matchbox/Utility/SimpleColourBasis2.cc --- a/MatrixElement/Matchbox/Utility/SimpleColourBasis2.cc +++ b/MatrixElement/Matchbox/Utility/SimpleColourBasis2.cc @@ -1,2882 +1,2877 @@ // -*- C++ -*- // // SimpleColourBasis2.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 SimpleColourBasis2 class. // #include "SimpleColourBasis2.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/StandardModel/StandardModelBase.h" - #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; -SimpleColourBasis2::SimpleColourBasis2() {} - -SimpleColourBasis2::~SimpleColourBasis2() {} - IBPtr SimpleColourBasis2::clone() const { return new_ptr(*this); } IBPtr SimpleColourBasis2::fullclone() const { return new_ptr(*this); } size_t SimpleColourBasis2::prepareBasis(const vector& basis) { if ( id33bar.empty() ) makeIds(); if ( basis == id88 ) { return 1; } if ( basis == id33bar ) { return 1; } if ( basis == id888 ) { return 2; } if ( basis == id33bar8 ) { return 1; } if ( basis == id8888 ) { return 9; } if ( basis == id33bar88 ) { return 3; } if ( basis == id33bar33bar ) { return 2; } if ( basis == id88888 ) { return 44; } if ( basis == id33bar888 ) { return 11; } if ( basis == id33bar33bar8 ) { return 4; } throw Exception() << "SimpleColourBasis2::prepareBasis(): Cannot handle colour configuration" << Exception::runerror; return 0; } double SimpleColourBasis2::scalarProduct(size_t a, size_t b, const vector& basis) const { if ( id33bar.empty() ) makeIds(); double Nc = SM().Nc(); double Nc2 = sqr(Nc); double Nc3 = Nc*Nc2; double Nc4 = sqr(Nc2); double Nc5 = Nc*Nc4; double Nc6 = Nc2*Nc4; double Nc8 = Nc2*Nc6; if ( a > b ) swap(a,b); if ( !largeN() ) { if ( basis == id88 ) { if( a == 0 && b == 0 ) return (-1 + Nc2)/4.; } if ( basis == id33bar ) { if( a == 0 && b == 0 ) return Nc; } if ( basis == id888 ) { if( ( a == 0 && b == 0 ) || ( a == 1 && b == 1 ) ) return (2 - 3*Nc2 + Nc4)/(8.*Nc); if( ( a == 0 ) && ( b == 1 ) ) return -(-1 + Nc2)/(4.*Nc); } if ( basis == id33bar8 ) { if( a == 0 && b == 0 ) return (-1 + Nc2)/2.; } if ( basis == id8888 ) { if( ( a == 0 && b == 0 ) || ( a == 1 && b == 1 ) || ( a == 2 && b == 2 ) ) return sqr(-1 + Nc2)/16.; if( ( a == 0 && b == 1 ) || ( a == 0 && b == 2 ) || ( a == 1 && b == 2 ) ) return (-1 + Nc2)/16.; if( ( a == 0 && b == 3 ) || ( a == 0 && b == 5 ) || ( a == 0 && b == 7 ) || ( a == 0 && b == 8 ) || ( a == 1 && b == 4 ) || ( a == 1 && b == 5 ) || ( a == 1 && b == 6 ) || ( a == 1 && b == 7 ) || ( a == 2 && b == 3 ) || ( a == 2 && b == 4 ) || ( a == 2 && b == 6 ) || ( a == 2 && b == 8 ) ) return sqr(-1 + Nc2)/(16.*Nc); if( ( a == 0 && b == 4 ) || ( a == 0 && b == 6 ) || ( a == 1 && b == 3 ) || ( a == 1 && b == 8 ) || ( a == 2 && b == 5 ) || ( a == 2 && b == 7 ) ) return -(-1 + Nc2)/(16.*Nc); if( ( a == 3 && b == 3 ) || ( a == 4 && b == 4 ) || ( a == 5 && b == 5 ) || ( a == 6 && b == 6 ) || ( a == 7 && b == 7 ) || ( a == 8 && b == 8 ) ) return (-3 + 6*Nc2 - 4*Nc4 + Nc6)/(16.*Nc2); if( ( a == 3 && b == 4 ) || ( a == 3 && b == 5 ) || ( a == 3 && b == 6 ) || ( a == 3 && b == 7 ) || ( a == 4 && b == 5 ) || ( a == 4 && b == 7 ) || ( a == 4 && b == 8 ) || ( a == 5 && b == 6 ) || ( a == 5 && b == 8 ) || ( a == 6 && b == 7 ) || ( a == 6 && b == 8 ) || ( a == 7 && b == 8 ) ) return (4 - 3/Nc2 - Nc2)/16.; if( ( a == 3 && b == 8 ) || ( a == 4 && b == 6 ) || ( a == 5 && b == 7 ) ) return (-3 + 2*Nc2 + Nc4)/(16.*Nc2); } if ( basis == id33bar88 ) { if( a == 0 && b == 0 ) return (Nc*(-1 + Nc2))/4.; if( ( a == 0 && b == 1 ) || ( a == 0 && b == 2 ) ) return (-1 + Nc2)/4.; if( ( a == 1 && b == 1 ) || ( a == 2 && b == 2 ) ) return sqr(-1 + Nc2)/(4.*Nc); if( a == 1 && b == 2 ) return -(-1 + Nc2)/(4.*Nc); } if ( basis == id33bar33bar ) { if( ( a == 0 && b == 0 ) || ( a == 1 && b == 1 ) ) return Nc2; if( a == 0 && b == 1 ) return Nc; } if ( basis == id88888 ) { if( ( a == 0 && b == 0 ) || ( a == 1 && b == 1 ) || ( a == 2 && b == 2 ) || ( a == 3 && b == 3 ) || ( a == 4 && b == 4 ) || ( a == 5 && b == 5 ) || ( a == 6 && b == 6 ) || ( a == 7 && b == 7 ) || ( a == 8 && b == 8 ) || ( a == 9 && b == 9 ) || ( a == 10 && b == 10 ) || ( a == 11 && b == 11 ) || ( a == 12 && b == 12 ) || ( a == 13 && b == 13 ) || ( a == 14 && b == 14 ) || ( a == 15 && b == 15 ) || ( a == 16 && b == 16 ) || ( a == 17 && b == 17 ) || ( a == 18 && b == 18 ) || ( a == 19 && b == 19 ) ) return ((-2 + Nc2)*sqr(-1 + Nc2))/(32.*Nc); if( ( a == 0 && b == 1 ) || ( a == 0 && b == 2 ) || ( a == 0 && b == 7 ) || ( a == 0 && b == 10 ) || ( a == 0 && b == 12 ) || ( a == 0 && b == 13 ) || ( a == 1 && b == 2 ) || ( a == 1 && b == 4 ) || ( a == 1 && b == 11 ) || ( a == 1 && b == 14 ) || ( a == 1 && b == 15 ) || ( a == 2 && b == 5 ) || ( a == 2 && b == 8 ) || ( a == 2 && b == 16 ) || ( a == 2 && b == 17 ) || ( a == 3 && b == 4 ) || ( a == 3 && b == 5 ) || ( a == 3 && b == 6 ) || ( a == 3 && b == 9 ) || ( a == 3 && b == 14 ) || ( a == 3 && b == 16 ) || ( a == 4 && b == 5 ) || ( a == 4 && b == 11 ) || ( a == 4 && b == 12 ) || ( a == 4 && b == 18 ) || ( a == 5 && b == 8 ) || ( a == 5 && b == 13 ) || ( a == 5 && b == 19 ) || ( a == 6 && b == 7 ) || ( a == 6 && b == 8 ) || ( a == 6 && b == 9 ) || ( a == 6 && b == 12 ) || ( a == 6 && b == 17 ) || ( a == 7 && b == 8 ) || ( a == 7 && b == 10 ) || ( a == 7 && b == 14 ) || ( a == 7 && b == 19 ) || ( a == 8 && b == 15 ) || ( a == 8 && b == 18 ) || ( a == 9 && b == 10 ) || ( a == 9 && b == 11 ) || ( a == 9 && b == 13 ) || ( a == 9 && b == 15 ) || ( a == 10 && b == 11 ) || ( a == 10 && b == 16 ) || ( a == 10 && b == 18 ) || ( a == 11 && b == 17 ) || ( a == 11 && b == 19 ) || ( a == 12 && b == 13 ) || ( a == 12 && b == 17 ) || ( a == 12 && b == 18 ) || ( a == 13 && b == 15 ) || ( a == 13 && b == 19 ) || ( a == 14 && b == 15 ) || ( a == 14 && b == 16 ) || ( a == 14 && b == 19 ) || ( a == 15 && b == 18 ) || ( a == 16 && b == 17 ) || ( a == 16 && b == 18 ) || ( a == 17 && b == 19 ) ) return (2 - 3*Nc2 + Nc4)/(32.*Nc); if( ( a == 0 && b == 3 ) || ( a == 1 && b == 6 ) || ( a == 2 && b == 9 ) || ( a == 4 && b == 7 ) || ( a == 5 && b == 10 ) || ( a == 8 && b == 11 ) || ( a == 12 && b == 14 ) || ( a == 13 && b == 16 ) || ( a == 15 && b == 17 ) || ( a == 18 && b == 19 ) ) return -sqr(-1 + Nc2)/(16.*Nc); if( ( a == 0 && b == 4 ) || ( a == 0 && b == 5 ) || ( a == 0 && b == 6 ) || ( a == 0 && b == 9 ) || ( a == 0 && b == 14 ) || ( a == 0 && b == 16 ) || ( a == 1 && b == 3 ) || ( a == 1 && b == 7 ) || ( a == 1 && b == 8 ) || ( a == 1 && b == 9 ) || ( a == 1 && b == 12 ) || ( a == 1 && b == 17 ) || ( a == 2 && b == 3 ) || ( a == 2 && b == 6 ) || ( a == 2 && b == 10 ) || ( a == 2 && b == 11 ) || ( a == 2 && b == 13 ) || ( a == 2 && b == 15 ) || ( a == 3 && b == 7 ) || ( a == 3 && b == 10 ) || ( a == 3 && b == 12 ) || ( a == 3 && b == 13 ) || ( a == 4 && b == 6 ) || ( a == 4 && b == 8 ) || ( a == 4 && b == 10 ) || ( a == 4 && b == 14 ) || ( a == 4 && b == 19 ) || ( a == 5 && b == 7 ) || ( a == 5 && b == 9 ) || ( a == 5 && b == 11 ) || ( a == 5 && b == 16 ) || ( a == 5 && b == 18 ) || ( a == 6 && b == 11 ) || ( a == 6 && b == 14 ) || ( a == 6 && b == 15 ) || ( a == 7 && b == 11 ) || ( a == 7 && b == 12 ) || ( a == 7 && b == 18 ) || ( a == 8 && b == 9 ) || ( a == 8 && b == 10 ) || ( a == 8 && b == 17 ) || ( a == 8 && b == 19 ) || ( a == 9 && b == 16 ) || ( a == 9 && b == 17 ) || ( a == 10 && b == 13 ) || ( a == 10 && b == 19 ) || ( a == 11 && b == 15 ) || ( a == 11 && b == 18 ) || ( a == 12 && b == 15 ) || ( a == 12 && b == 16 ) || ( a == 12 && b == 19 ) || ( a == 13 && b == 14 ) || ( a == 13 && b == 17 ) || ( a == 13 && b == 18 ) || ( a == 14 && b == 17 ) || ( a == 14 && b == 18 ) || ( a == 15 && b == 16 ) || ( a == 15 && b == 19 ) || ( a == 16 && b == 19 ) || ( a == 17 && b == 18 ) ) return -(-1 + Nc2)/(16.*Nc); if( ( a == 0 && b == 8 ) || ( a == 0 && b == 11 ) || ( a == 0 && b == 15 ) || ( a == 0 && b == 17 ) || ( a == 0 && b == 18 ) || ( a == 0 && b == 19 ) || ( a == 1 && b == 5 ) || ( a == 1 && b == 10 ) || ( a == 1 && b == 13 ) || ( a == 1 && b == 16 ) || ( a == 1 && b == 18 ) || ( a == 1 && b == 19 ) || ( a == 2 && b == 4 ) || ( a == 2 && b == 7 ) || ( a == 2 && b == 12 ) || ( a == 2 && b == 14 ) || ( a == 2 && b == 18 ) || ( a == 2 && b == 19 ) || ( a == 3 && b == 8 ) || ( a == 3 && b == 11 ) || ( a == 3 && b == 15 ) || ( a == 3 && b == 17 ) || ( a == 3 && b == 18 ) || ( a == 3 && b == 19 ) || ( a == 4 && b == 9 ) || ( a == 4 && b == 13 ) || ( a == 4 && b == 15 ) || ( a == 4 && b == 16 ) || ( a == 4 && b == 17 ) || ( a == 5 && b == 6 ) || ( a == 5 && b == 12 ) || ( a == 5 && b == 14 ) || ( a == 5 && b == 15 ) || ( a == 5 && b == 17 ) || ( a == 6 && b == 10 ) || ( a == 6 && b == 13 ) || ( a == 6 && b == 16 ) || ( a == 6 && b == 18 ) || ( a == 6 && b == 19 ) || ( a == 7 && b == 9 ) || ( a == 7 && b == 13 ) || ( a == 7 && b == 15 ) || ( a == 7 && b == 16 ) || ( a == 7 && b == 17 ) || ( a == 8 && b == 12 ) || ( a == 8 && b == 13 ) || ( a == 8 && b == 14 ) || ( a == 8 && b == 16 ) || ( a == 9 && b == 12 ) || ( a == 9 && b == 14 ) || ( a == 9 && b == 18 ) || ( a == 9 && b == 19 ) || ( a == 10 && b == 12 ) || ( a == 10 && b == 14 ) || ( a == 10 && b == 15 ) || ( a == 10 && b == 17 ) || ( a == 11 && b == 12 ) || ( a == 11 && b == 13 ) || ( a == 11 && b == 14 ) || ( a == 11 && b == 16 ) ) return 0; if( ( a == 0 && b == 20 ) || ( a == 0 && b == 21 ) || ( a == 0 && b == 23 ) || ( a == 0 && b == 25 ) || ( a == 0 && b == 36 ) || ( a == 0 && b == 42 ) || ( a == 1 && b == 21 ) || ( a == 1 && b == 22 ) || ( a == 1 && b == 23 ) || ( a == 1 && b == 24 ) || ( a == 1 && b == 30 ) || ( a == 1 && b == 40 ) || ( a == 2 && b == 20 ) || ( a == 2 && b == 22 ) || ( a == 2 && b == 24 ) || ( a == 2 && b == 25 ) || ( a == 2 && b == 28 ) || ( a == 2 && b == 34 ) || ( a == 3 && b == 26 ) || ( a == 3 && b == 27 ) || ( a == 3 && b == 29 ) || ( a == 3 && b == 31 ) || ( a == 3 && b == 37 ) || ( a == 3 && b == 43 ) || ( a == 4 && b == 24 ) || ( a == 4 && b == 27 ) || ( a == 4 && b == 28 ) || ( a == 4 && b == 29 ) || ( a == 4 && b == 30 ) || ( a == 4 && b == 38 ) || ( a == 5 && b == 22 ) || ( a == 5 && b == 26 ) || ( a == 5 && b == 28 ) || ( a == 5 && b == 30 ) || ( a == 5 && b == 31 ) || ( a == 5 && b == 32 ) || ( a == 6 && b == 31 ) || ( a == 6 && b == 32 ) || ( a == 6 && b == 33 ) || ( a == 6 && b == 35 ) || ( a == 6 && b == 37 ) || ( a == 6 && b == 41 ) || ( a == 7 && b == 25 ) || ( a == 7 && b == 33 ) || ( a == 7 && b == 34 ) || ( a == 7 && b == 35 ) || ( a == 7 && b == 36 ) || ( a == 7 && b == 39 ) || ( a == 8 && b == 20 ) || ( a == 8 && b == 26 ) || ( a == 8 && b == 32 ) || ( a == 8 && b == 34 ) || ( a == 8 && b == 36 ) || ( a == 8 && b == 37 ) || ( a == 9 && b == 29 ) || ( a == 9 && b == 35 ) || ( a == 9 && b == 38 ) || ( a == 9 && b == 39 ) || ( a == 9 && b == 41 ) || ( a == 9 && b == 43 ) || ( a == 10 && b == 23 ) || ( a == 10 && b == 33 ) || ( a == 10 && b == 39 ) || ( a == 10 && b == 40 ) || ( a == 10 && b == 41 ) || ( a == 10 && b == 42 ) || ( a == 11 && b == 21 ) || ( a == 11 && b == 27 ) || ( a == 11 && b == 38 ) || ( a == 11 && b == 40 ) || ( a == 11 && b == 42 ) || ( a == 11 && b == 43 ) || ( a == 12 && b == 20 ) || ( a == 12 && b == 28 ) || ( a == 12 && b == 32 ) || ( a == 12 && b == 38 ) || ( a == 12 && b == 41 ) || ( a == 12 && b == 42 ) || ( a == 13 && b == 21 ) || ( a == 13 && b == 30 ) || ( a == 13 && b == 32 ) || ( a == 13 && b == 35 ) || ( a == 13 && b == 36 ) || ( a == 13 && b == 38 ) || ( a == 14 && b == 22 ) || ( a == 14 && b == 26 ) || ( a == 14 && b == 34 ) || ( a == 14 && b == 39 ) || ( a == 14 && b == 40 ) || ( a == 14 && b == 43 ) || ( a == 15 && b == 23 ) || ( a == 15 && b == 26 ) || ( a == 15 && b == 29 ) || ( a == 15 && b == 30 ) || ( a == 15 && b == 36 ) || ( a == 15 && b == 39 ) || ( a == 16 && b == 24 ) || ( a == 16 && b == 27 ) || ( a == 16 && b == 33 ) || ( a == 16 && b == 34 ) || ( a == 16 && b == 37 ) || ( a == 16 && b == 40 ) || ( a == 17 && b == 25 ) || ( a == 17 && b == 27 ) || ( a == 17 && b == 28 ) || ( a == 17 && b == 31 ) || ( a == 17 && b == 33 ) || ( a == 17 && b == 42 ) || ( a == 18 && b == 20 ) || ( a == 18 && b == 23 ) || ( a == 18 && b == 24 ) || ( a == 18 && b == 29 ) || ( a == 18 && b == 37 ) || ( a == 18 && b == 41 ) || ( a == 19 && b == 21 ) || ( a == 19 && b == 22 ) || ( a == 19 && b == 25 ) || ( a == 19 && b == 31 ) || ( a == 19 && b == 35 ) || ( a == 19 && b == 43 ) ) return ((-2 + Nc2)*sqr(-1 + Nc2))/(32.*Nc2); if( ( a == 0 && b == 22 ) || ( a == 0 && b == 24 ) || ( a == 0 && b == 32 ) || ( a == 0 && b == 33 ) || ( a == 0 && b == 38 ) || ( a == 0 && b == 39 ) || ( a == 1 && b == 20 ) || ( a == 1 && b == 25 ) || ( a == 1 && b == 26 ) || ( a == 1 && b == 27 ) || ( a == 1 && b == 38 ) || ( a == 1 && b == 39 ) || ( a == 2 && b == 21 ) || ( a == 2 && b == 23 ) || ( a == 2 && b == 26 ) || ( a == 2 && b == 27 ) || ( a == 2 && b == 32 ) || ( a == 2 && b == 33 ) || ( a == 3 && b == 28 ) || ( a == 3 && b == 30 ) || ( a == 3 && b == 34 ) || ( a == 3 && b == 35 ) || ( a == 3 && b == 40 ) || ( a == 3 && b == 41 ) || ( a == 4 && b == 20 ) || ( a == 4 && b == 21 ) || ( a == 4 && b == 26 ) || ( a == 4 && b == 31 ) || ( a == 4 && b == 40 ) || ( a == 4 && b == 41 ) || ( a == 5 && b == 20 ) || ( a == 5 && b == 21 ) || ( a == 5 && b == 27 ) || ( a == 5 && b == 29 ) || ( a == 5 && b == 34 ) || ( a == 5 && b == 35 ) || ( a == 6 && b == 28 ) || ( a == 6 && b == 29 ) || ( a == 6 && b == 34 ) || ( a == 6 && b == 36 ) || ( a == 6 && b == 42 ) || ( a == 6 && b == 43 ) || ( a == 7 && b == 22 ) || ( a == 7 && b == 23 ) || ( a == 7 && b == 32 ) || ( a == 7 && b == 37 ) || ( a == 7 && b == 42 ) || ( a == 7 && b == 43 ) || ( a == 8 && b == 22 ) || ( a == 8 && b == 23 ) || ( a == 8 && b == 28 ) || ( a == 8 && b == 29 ) || ( a == 8 && b == 33 ) || ( a == 8 && b == 35 ) || ( a == 9 && b == 30 ) || ( a == 9 && b == 31 ) || ( a == 9 && b == 36 ) || ( a == 9 && b == 37 ) || ( a == 9 && b == 40 ) || ( a == 9 && b == 42 ) || ( a == 10 && b == 24 ) || ( a == 10 && b == 25 ) || ( a == 10 && b == 36 ) || ( a == 10 && b == 37 ) || ( a == 10 && b == 38 ) || ( a == 10 && b == 43 ) || ( a == 11 && b == 24 ) || ( a == 11 && b == 25 ) || ( a == 11 && b == 30 ) || ( a == 11 && b == 31 ) || ( a == 11 && b == 39 ) || ( a == 11 && b == 41 ) || ( a == 12 && b == 21 ) || ( a == 12 && b == 24 ) || ( a == 12 && b == 29 ) || ( a == 12 && b == 31 ) || ( a == 12 && b == 33 ) || ( a == 12 && b == 36 ) || ( a == 13 && b == 20 ) || ( a == 13 && b == 22 ) || ( a == 13 && b == 29 ) || ( a == 13 && b == 31 ) || ( a == 13 && b == 39 ) || ( a == 13 && b == 42 ) || ( a == 14 && b == 23 ) || ( a == 14 && b == 25 ) || ( a == 14 && b == 27 ) || ( a == 14 && b == 30 ) || ( a == 14 && b == 35 ) || ( a == 14 && b == 37 ) || ( a == 15 && b == 20 ) || ( a == 15 && b == 22 ) || ( a == 15 && b == 35 ) || ( a == 15 && b == 37 ) || ( a == 15 && b == 38 ) || ( a == 15 && b == 40 ) || ( a == 16 && b == 23 ) || ( a == 16 && b == 25 ) || ( a == 16 && b == 26 ) || ( a == 16 && b == 28 ) || ( a == 16 && b == 41 ) || ( a == 16 && b == 43 ) || ( a == 17 && b == 21 ) || ( a == 17 && b == 24 ) || ( a == 17 && b == 32 ) || ( a == 17 && b == 34 ) || ( a == 17 && b == 41 ) || ( a == 17 && b == 43 ) || ( a == 18 && b == 26 ) || ( a == 18 && b == 28 ) || ( a == 18 && b == 33 ) || ( a == 18 && b == 36 ) || ( a == 18 && b == 38 ) || ( a == 18 && b == 40 ) || ( a == 19 && b == 27 ) || ( a == 19 && b == 30 ) || ( a == 19 && b == 32 ) || ( a == 19 && b == 34 ) || ( a == 19 && b == 39 ) || ( a == 19 && b == 42 ) ) return -(2 - 3*Nc2 + Nc4)/(32.*Nc2); if( ( a == 0 && b == 26 ) || ( a == 0 && b == 27 ) || ( a == 0 && b == 29 ) || ( a == 0 && b == 31 ) || ( a == 0 && b == 37 ) || ( a == 0 && b == 43 ) || ( a == 1 && b == 31 ) || ( a == 1 && b == 32 ) || ( a == 1 && b == 33 ) || ( a == 1 && b == 35 ) || ( a == 1 && b == 37 ) || ( a == 1 && b == 41 ) || ( a == 2 && b == 29 ) || ( a == 2 && b == 35 ) || ( a == 2 && b == 38 ) || ( a == 2 && b == 39 ) || ( a == 2 && b == 41 ) || ( a == 2 && b == 43 ) || ( a == 3 && b == 20 ) || ( a == 3 && b == 21 ) || ( a == 3 && b == 23 ) || ( a == 3 && b == 25 ) || ( a == 3 && b == 36 ) || ( a == 3 && b == 42 ) || ( a == 4 && b == 25 ) || ( a == 4 && b == 33 ) || ( a == 4 && b == 34 ) || ( a == 4 && b == 35 ) || ( a == 4 && b == 36 ) || ( a == 4 && b == 39 ) || ( a == 5 && b == 23 ) || ( a == 5 && b == 33 ) || ( a == 5 && b == 39 ) || ( a == 5 && b == 40 ) || ( a == 5 && b == 41 ) || ( a == 5 && b == 42 ) || ( a == 6 && b == 21 ) || ( a == 6 && b == 22 ) || ( a == 6 && b == 23 ) || ( a == 6 && b == 24 ) || ( a == 6 && b == 30 ) || ( a == 6 && b == 40 ) || ( a == 7 && b == 24 ) || ( a == 7 && b == 27 ) || ( a == 7 && b == 28 ) || ( a == 7 && b == 29 ) || ( a == 7 && b == 30 ) || ( a == 7 && b == 38 ) || ( a == 8 && b == 21 ) || ( a == 8 && b == 27 ) || ( a == 8 && b == 38 ) || ( a == 8 && b == 40 ) || ( a == 8 && b == 42 ) || ( a == 8 && b == 43 ) || ( a == 9 && b == 20 ) || ( a == 9 && b == 22 ) || ( a == 9 && b == 24 ) || ( a == 9 && b == 25 ) || ( a == 9 && b == 28 ) || ( a == 9 && b == 34 ) || ( a == 10 && b == 22 ) || ( a == 10 && b == 26 ) || ( a == 10 && b == 28 ) || ( a == 10 && b == 30 ) || ( a == 10 && b == 31 ) || ( a == 10 && b == 32 ) || ( a == 11 && b == 20 ) || ( a == 11 && b == 26 ) || ( a == 11 && b == 32 ) || ( a == 11 && b == 34 ) || ( a == 11 && b == 36 ) || ( a == 11 && b == 37 ) || ( a == 12 && b == 22 ) || ( a == 12 && b == 26 ) || ( a == 12 && b == 34 ) || ( a == 12 && b == 39 ) || ( a == 12 && b == 40 ) || ( a == 12 && b == 43 ) || ( a == 13 && b == 24 ) || ( a == 13 && b == 27 ) || ( a == 13 && b == 33 ) || ( a == 13 && b == 34 ) || ( a == 13 && b == 37 ) || ( a == 13 && b == 40 ) || ( a == 14 && b == 20 ) || ( a == 14 && b == 28 ) || ( a == 14 && b == 32 ) || ( a == 14 && b == 38 ) || ( a == 14 && b == 41 ) || ( a == 14 && b == 42 ) || ( a == 15 && b == 25 ) || ( a == 15 && b == 27 ) || ( a == 15 && b == 28 ) || ( a == 15 && b == 31 ) || ( a == 15 && b == 33 ) || ( a == 15 && b == 42 ) || ( a == 16 && b == 21 ) || ( a == 16 && b == 30 ) || ( a == 16 && b == 32 ) || ( a == 16 && b == 35 ) || ( a == 16 && b == 36 ) || ( a == 16 && b == 38 ) || ( a == 17 && b == 23 ) || ( a == 17 && b == 26 ) || ( a == 17 && b == 29 ) || ( a == 17 && b == 30 ) || ( a == 17 && b == 36 ) || ( a == 17 && b == 39 ) || ( a == 18 && b == 21 ) || ( a == 18 && b == 22 ) || ( a == 18 && b == 25 ) || ( a == 18 && b == 31 ) || ( a == 18 && b == 35 ) || ( a == 18 && b == 43 ) || ( a == 19 && b == 20 ) || ( a == 19 && b == 23 ) || ( a == 19 && b == 24 ) || ( a == 19 && b == 29 ) || ( a == 19 && b == 37 ) || ( a == 19 && b == 41 ) ) return -sqr(-1 + Nc2)/(16.*Nc2); if( ( a == 0 && b == 28 ) || ( a == 0 && b == 30 ) || ( a == 0 && b == 34 ) || ( a == 0 && b == 35 ) || ( a == 0 && b == 40 ) || ( a == 0 && b == 41 ) || ( a == 1 && b == 28 ) || ( a == 1 && b == 29 ) || ( a == 1 && b == 34 ) || ( a == 1 && b == 36 ) || ( a == 1 && b == 42 ) || ( a == 1 && b == 43 ) || ( a == 2 && b == 30 ) || ( a == 2 && b == 31 ) || ( a == 2 && b == 36 ) || ( a == 2 && b == 37 ) || ( a == 2 && b == 40 ) || ( a == 2 && b == 42 ) || ( a == 3 && b == 22 ) || ( a == 3 && b == 24 ) || ( a == 3 && b == 32 ) || ( a == 3 && b == 33 ) || ( a == 3 && b == 38 ) || ( a == 3 && b == 39 ) || ( a == 4 && b == 22 ) || ( a == 4 && b == 23 ) || ( a == 4 && b == 32 ) || ( a == 4 && b == 37 ) || ( a == 4 && b == 42 ) || ( a == 4 && b == 43 ) || ( a == 5 && b == 24 ) || ( a == 5 && b == 25 ) || ( a == 5 && b == 36 ) || ( a == 5 && b == 37 ) || ( a == 5 && b == 38 ) || ( a == 5 && b == 43 ) || ( a == 6 && b == 20 ) || ( a == 6 && b == 25 ) || ( a == 6 && b == 26 ) || ( a == 6 && b == 27 ) || ( a == 6 && b == 38 ) || ( a == 6 && b == 39 ) || ( a == 7 && b == 20 ) || ( a == 7 && b == 21 ) || ( a == 7 && b == 26 ) || ( a == 7 && b == 31 ) || ( a == 7 && b == 40 ) || ( a == 7 && b == 41 ) || ( a == 8 && b == 24 ) || ( a == 8 && b == 25 ) || ( a == 8 && b == 30 ) || ( a == 8 && b == 31 ) || ( a == 8 && b == 39 ) || ( a == 8 && b == 41 ) || ( a == 9 && b == 21 ) || ( a == 9 && b == 23 ) || ( a == 9 && b == 26 ) || ( a == 9 && b == 27 ) || ( a == 9 && b == 32 ) || ( a == 9 && b == 33 ) || ( a == 10 && b == 20 ) || ( a == 10 && b == 21 ) || ( a == 10 && b == 27 ) || ( a == 10 && b == 29 ) || ( a == 10 && b == 34 ) || ( a == 10 && b == 35 ) || ( a == 11 && b == 22 ) || ( a == 11 && b == 23 ) || ( a == 11 && b == 28 ) || ( a == 11 && b == 29 ) || ( a == 11 && b == 33 ) || ( a == 11 && b == 35 ) || ( a == 12 && b == 23 ) || ( a == 12 && b == 25 ) || ( a == 12 && b == 27 ) || ( a == 12 && b == 30 ) || ( a == 12 && b == 35 ) || ( a == 12 && b == 37 ) || ( a == 13 && b == 23 ) || ( a == 13 && b == 25 ) || ( a == 13 && b == 26 ) || ( a == 13 && b == 28 ) || ( a == 13 && b == 41 ) || ( a == 13 && b == 43 ) || ( a == 14 && b == 21 ) || ( a == 14 && b == 24 ) || ( a == 14 && b == 29 ) || ( a == 14 && b == 31 ) || ( a == 14 && b == 33 ) || ( a == 14 && b == 36 ) || ( a == 15 && b == 21 ) || ( a == 15 && b == 24 ) || ( a == 15 && b == 32 ) || ( a == 15 && b == 34 ) || ( a == 15 && b == 41 ) || ( a == 15 && b == 43 ) || ( a == 16 && b == 20 ) || ( a == 16 && b == 22 ) || ( a == 16 && b == 29 ) || ( a == 16 && b == 31 ) || ( a == 16 && b == 39 ) || ( a == 16 && b == 42 ) || ( a == 17 && b == 20 ) || ( a == 17 && b == 22 ) || ( a == 17 && b == 35 ) || ( a == 17 && b == 37 ) || ( a == 17 && b == 38 ) || ( a == 17 && b == 40 ) || ( a == 18 && b == 27 ) || ( a == 18 && b == 30 ) || ( a == 18 && b == 32 ) || ( a == 18 && b == 34 ) || ( a == 18 && b == 39 ) || ( a == 18 && b == 42 ) || ( a == 19 && b == 26 ) || ( a == 19 && b == 28 ) || ( a == 19 && b == 33 ) || ( a == 19 && b == 36 ) || ( a == 19 && b == 38 ) || ( a == 19 && b == 40 ) ) return (1 - 1/Nc2)/16.; if( ( a == 20 && b == 20 ) || ( a == 21 && b == 21 ) || ( a == 22 && b == 22 ) || ( a == 23 && b == 23 ) || ( a == 24 && b == 24 ) || ( a == 25 && b == 25 ) || ( a == 26 && b == 26 ) || ( a == 27 && b == 27 ) || ( a == 28 && b == 28 ) || ( a == 29 && b == 29 ) || ( a == 30 && b == 30 ) || ( a == 31 && b == 31 ) || ( a == 32 && b == 32 ) || ( a == 33 && b == 33 ) || ( a == 34 && b == 34 ) || ( a == 35 && b == 35 ) || ( a == 36 && b == 36 ) || ( a == 37 && b == 37 ) || ( a == 38 && b == 38 ) || ( a == 39 && b == 39 ) || ( a == 40 && b == 40 ) || ( a == 41 && b == 41 ) || ( a == 42 && b == 42 ) || ( a == 43 && b == 43 ) ) return (4 - 10*Nc2 + 10*Nc4 - 5*Nc6 + Nc8)/(32.*Nc3); if( ( a == 20 && b == 21 ) || ( a == 20 && b == 22 ) || ( a == 20 && b == 26 ) || ( a == 20 && b == 29 ) || ( a == 20 && b == 38 ) || ( a == 21 && b == 24 ) || ( a == 21 && b == 27 ) || ( a == 21 && b == 31 ) || ( a == 21 && b == 32 ) || ( a == 22 && b == 23 ) || ( a == 22 && b == 32 ) || ( a == 22 && b == 35 ) || ( a == 22 && b == 39 ) || ( a == 23 && b == 25 ) || ( a == 23 && b == 26 ) || ( a == 23 && b == 33 ) || ( a == 23 && b == 37 ) || ( a == 24 && b == 25 ) || ( a == 24 && b == 33 ) || ( a == 24 && b == 38 ) || ( a == 24 && b == 41 ) || ( a == 25 && b == 27 ) || ( a == 25 && b == 39 ) || ( a == 25 && b == 43 ) || ( a == 26 && b == 27 ) || ( a == 26 && b == 28 ) || ( a == 26 && b == 40 ) || ( a == 27 && b == 30 ) || ( a == 27 && b == 34 ) || ( a == 28 && b == 29 ) || ( a == 28 && b == 33 ) || ( a == 28 && b == 34 ) || ( a == 28 && b == 41 ) || ( a == 29 && b == 31 ) || ( a == 29 && b == 35 ) || ( a == 29 && b == 36 ) || ( a == 30 && b == 31 ) || ( a == 30 && b == 35 ) || ( a == 30 && b == 39 ) || ( a == 30 && b == 40 ) || ( a == 31 && b == 41 ) || ( a == 31 && b == 42 ) || ( a == 32 && b == 33 ) || ( a == 32 && b == 34 ) || ( a == 32 && b == 42 ) || ( a == 33 && b == 36 ) || ( a == 34 && b == 35 ) || ( a == 34 && b == 43 ) || ( a == 35 && b == 37 ) || ( a == 36 && b == 37 ) || ( a == 36 && b == 38 ) || ( a == 36 && b == 42 ) || ( a == 37 && b == 40 ) || ( a == 37 && b == 43 ) || ( a == 38 && b == 39 ) || ( a == 38 && b == 40 ) || ( a == 39 && b == 42 ) || ( a == 40 && b == 41 ) || ( a == 41 && b == 43 ) || ( a == 42 && b == 43 ) ) return -(-4 + 7*Nc2 - 4*Nc4 + Nc6)/(32.*Nc3); if( ( a == 20 && b == 23 ) || ( a == 20 && b == 24 ) || ( a == 20 && b == 28 ) || ( a == 20 && b == 32 ) || ( a == 20 && b == 36 ) || ( a == 21 && b == 22 ) || ( a == 21 && b == 25 ) || ( a == 21 && b == 30 ) || ( a == 21 && b == 38 ) || ( a == 21 && b == 42 ) || ( a == 22 && b == 25 ) || ( a == 22 && b == 26 ) || ( a == 22 && b == 30 ) || ( a == 22 && b == 34 ) || ( a == 23 && b == 24 ) || ( a == 23 && b == 36 ) || ( a == 23 && b == 39 ) || ( a == 23 && b == 40 ) || ( a == 24 && b == 27 ) || ( a == 24 && b == 28 ) || ( a == 24 && b == 40 ) || ( a == 25 && b == 33 ) || ( a == 25 && b == 34 ) || ( a == 25 && b == 42 ) || ( a == 26 && b == 29 ) || ( a == 26 && b == 30 ) || ( a == 26 && b == 34 ) || ( a == 26 && b == 37 ) || ( a == 27 && b == 28 ) || ( a == 27 && b == 31 ) || ( a == 27 && b == 40 ) || ( a == 27 && b == 43 ) || ( a == 28 && b == 31 ) || ( a == 28 && b == 32 ) || ( a == 29 && b == 30 ) || ( a == 29 && b == 37 ) || ( a == 29 && b == 38 ) || ( a == 29 && b == 41 ) || ( a == 30 && b == 38 ) || ( a == 31 && b == 32 ) || ( a == 31 && b == 35 ) || ( a == 31 && b == 43 ) || ( a == 32 && b == 35 ) || ( a == 32 && b == 36 ) || ( a == 33 && b == 34 ) || ( a == 33 && b == 37 ) || ( a == 33 && b == 41 ) || ( a == 33 && b == 42 ) || ( a == 34 && b == 37 ) || ( a == 35 && b == 36 ) || ( a == 35 && b == 39 ) || ( a == 35 && b == 43 ) || ( a == 36 && b == 39 ) || ( a == 37 && b == 41 ) || ( a == 38 && b == 41 ) || ( a == 38 && b == 42 ) || ( a == 39 && b == 40 ) || ( a == 39 && b == 43 ) || ( a == 40 && b == 43 ) || ( a == 41 && b == 42 ) ) return (2 - 3*Nc2 + Nc4)/(16.*Nc3); if( ( a == 20 && b == 25 ) || ( a == 20 && b == 34 ) || ( a == 20 && b == 37 ) || ( a == 20 && b == 41 ) || ( a == 20 && b == 42 ) || ( a == 21 && b == 23 ) || ( a == 21 && b == 35 ) || ( a == 21 && b == 36 ) || ( a == 21 && b == 40 ) || ( a == 21 && b == 43 ) || ( a == 22 && b == 24 ) || ( a == 22 && b == 28 ) || ( a == 22 && b == 31 ) || ( a == 22 && b == 40 ) || ( a == 22 && b == 43 ) || ( a == 23 && b == 29 ) || ( a == 23 && b == 30 ) || ( a == 23 && b == 41 ) || ( a == 23 && b == 42 ) || ( a == 24 && b == 29 ) || ( a == 24 && b == 30 ) || ( a == 24 && b == 34 ) || ( a == 24 && b == 37 ) || ( a == 25 && b == 28 ) || ( a == 25 && b == 31 ) || ( a == 25 && b == 35 ) || ( a == 25 && b == 36 ) || ( a == 26 && b == 31 ) || ( a == 26 && b == 32 ) || ( a == 26 && b == 36 ) || ( a == 26 && b == 39 ) || ( a == 26 && b == 43 ) || ( a == 27 && b == 29 ) || ( a == 27 && b == 33 ) || ( a == 27 && b == 37 ) || ( a == 27 && b == 38 ) || ( a == 27 && b == 42 ) || ( a == 28 && b == 30 ) || ( a == 28 && b == 38 ) || ( a == 28 && b == 42 ) || ( a == 29 && b == 39 ) || ( a == 29 && b == 43 ) || ( a == 30 && b == 32 ) || ( a == 30 && b == 36 ) || ( a == 31 && b == 33 ) || ( a == 31 && b == 37 ) || ( a == 32 && b == 37 ) || ( a == 32 && b == 38 ) || ( a == 32 && b == 41 ) || ( a == 33 && b == 35 ) || ( a == 33 && b == 39 ) || ( a == 33 && b == 40 ) || ( a == 34 && b == 36 ) || ( a == 34 && b == 39 ) || ( a == 34 && b == 40 ) || ( a == 35 && b == 38 ) || ( a == 35 && b == 41 ) || ( a == 38 && b == 43 ) || ( a == 39 && b == 41 ) || ( a == 40 && b == 42 ) ) return (4 - 3*Nc2 - 2*Nc4 + Nc6)/(32.*Nc3); if( ( a == 20 && b == 27 ) || ( a == 20 && b == 31 ) || ( a == 20 && b == 35 ) || ( a == 20 && b == 39 ) || ( a == 20 && b == 40 ) || ( a == 21 && b == 26 ) || ( a == 21 && b == 29 ) || ( a == 21 && b == 33 ) || ( a == 21 && b == 34 ) || ( a == 21 && b == 41 ) || ( a == 22 && b == 29 ) || ( a == 22 && b == 33 ) || ( a == 22 && b == 37 ) || ( a == 22 && b == 38 ) || ( a == 22 && b == 42 ) || ( a == 23 && b == 27 ) || ( a == 23 && b == 28 ) || ( a == 23 && b == 32 ) || ( a == 23 && b == 35 ) || ( a == 23 && b == 43 ) || ( a == 24 && b == 31 ) || ( a == 24 && b == 32 ) || ( a == 24 && b == 36 ) || ( a == 24 && b == 39 ) || ( a == 24 && b == 43 ) || ( a == 25 && b == 26 ) || ( a == 25 && b == 30 ) || ( a == 25 && b == 37 ) || ( a == 25 && b == 38 ) || ( a == 25 && b == 41 ) || ( a == 26 && b == 33 ) || ( a == 26 && b == 38 ) || ( a == 26 && b == 41 ) || ( a == 27 && b == 32 ) || ( a == 27 && b == 35 ) || ( a == 27 && b == 39 ) || ( a == 28 && b == 35 ) || ( a == 28 && b == 36 ) || ( a == 28 && b == 40 ) || ( a == 28 && b == 43 ) || ( a == 29 && b == 33 ) || ( a == 29 && b == 34 ) || ( a == 29 && b == 42 ) || ( a == 30 && b == 34 ) || ( a == 30 && b == 37 ) || ( a == 30 && b == 41 ) || ( a == 30 && b == 42 ) || ( a == 31 && b == 36 ) || ( a == 31 && b == 39 ) || ( a == 31 && b == 40 ) || ( a == 32 && b == 39 ) || ( a == 32 && b == 43 ) || ( a == 33 && b == 38 ) || ( a == 34 && b == 41 ) || ( a == 34 && b == 42 ) || ( a == 35 && b == 40 ) || ( a == 36 && b == 40 ) || ( a == 36 && b == 43 ) || ( a == 37 && b == 38 ) || ( a == 37 && b == 42 ) ) return -(-1 + Nc2)/(8.*Nc3); if( ( a == 20 && b == 30 ) || ( a == 20 && b == 33 ) || ( a == 21 && b == 28 ) || ( a == 21 && b == 39 ) || ( a == 22 && b == 27 ) || ( a == 22 && b == 36 ) || ( a == 23 && b == 34 ) || ( a == 23 && b == 38 ) || ( a == 24 && b == 26 ) || ( a == 24 && b == 42 ) || ( a == 25 && b == 32 ) || ( a == 25 && b == 40 ) || ( a == 26 && b == 35 ) || ( a == 27 && b == 41 ) || ( a == 28 && b == 37 ) || ( a == 29 && b == 32 ) || ( a == 29 && b == 40 ) || ( a == 30 && b == 43 ) || ( a == 31 && b == 34 ) || ( a == 31 && b == 38 ) || ( a == 33 && b == 43 ) || ( a == 35 && b == 42 ) || ( a == 36 && b == 41 ) || ( a == 37 && b == 39 ) ) return (4 - 5*Nc2 + Nc4)/(32.*Nc3); if( ( a == 20 && b == 43 ) || ( a == 21 && b == 37 ) || ( a == 22 && b == 41 ) || ( a == 23 && b == 31 ) || ( a == 24 && b == 35 ) || ( a == 25 && b == 29 ) || ( a == 26 && b == 42 ) || ( a == 27 && b == 36 ) || ( a == 28 && b == 39 ) || ( a == 30 && b == 33 ) || ( a == 32 && b == 40 ) || ( a == 34 && b == 38 ) ) return -(-1 + Nc4)/(8.*Nc3); } if ( basis == id33bar888 ) { if( ( a == 0 && b == 0 ) || ( a == 1 && b == 1 ) ) return (2 - 3*Nc2 + Nc4)/8.; if( a == 0 && b == 1 ) return (1 - Nc2)/4.; if( ( a == 0 && b == 2 ) || ( a == 0 && b == 3 ) || ( a == 0 && b == 4 ) || ( a == 1 && b == 2 ) || ( a == 1 && b == 3 ) || ( a == 1 && b == 4 ) ) return 0; if( ( a == 0 && b == 5 ) || ( a == 0 && b == 8 ) || ( a == 0 && b == 9 ) || ( a == 1 && b == 6 ) || ( a == 1 && b == 7 ) || ( a == 1 && b == 10 ) ) return (2 - 3*Nc2 + Nc4)/(8.*Nc); if( ( a == 0 && b == 6 ) || ( a == 0 && b == 7 ) || ( a == 0 && b == 10 ) || ( a == 1 && b == 5 ) || ( a == 1 && b == 8 ) || ( a == 1 && b == 9 ) ) return -(-1 + Nc2)/(4.*Nc); if( ( a == 2 && b == 2 ) || ( a == 3 && b == 3 ) || ( a == 4 && b == 4 ) ) return sqr(-1 + Nc2)/8.; if( ( a == 2 && b == 3 ) || ( a == 2 && b == 4 ) || ( a == 3 && b == 4 ) ) return (-1 + Nc2)/8.; if( ( a == 2 && b == 5 ) || ( a == 2 && b == 6 ) || ( a == 2 && b == 8 ) || ( a == 2 && b == 10 ) || ( a == 3 && b == 6 ) || ( a == 3 && b == 7 ) || ( a == 3 && b == 8 ) || ( a == 3 && b == 9 ) || ( a == 4 && b == 5 ) || ( a == 4 && b == 7 ) || ( a == 4 && b == 9 ) || ( a == 4 && b == 10 ) ) return sqr(-1 + Nc2)/(8.*Nc); if( ( a == 2 && b == 7 ) || ( a == 2 && b == 9 ) || ( a == 3 && b == 5 ) || ( a == 3 && b == 10 ) || ( a == 4 && b == 6 ) || ( a == 4 && b == 8 ) ) return -(-1 + Nc2)/(8.*Nc); if( ( a == 5 && b == 5 ) || ( a == 6 && b == 6 ) || ( a == 7 && b == 7 ) || ( a == 8 && b == 8 ) || ( a == 9 && b == 9 ) || ( a == 10 && b == 10 ) ) return pow(-1 + Nc2,3.)/(8.*Nc2); if( ( a == 5 && b == 6 ) || ( a == 5 && b == 7 ) || ( a == 6 && b == 9 ) || ( a == 7 && b == 8 ) || ( a == 8 && b == 10 ) || ( a == 9 && b == 10 ) ) return -sqr(-1 + Nc2)/(8.*Nc2); if( ( a == 5 && b == 8 ) || ( a == 5 && b == 9 ) || ( a == 6 && b == 7 ) || ( a == 6 && b == 10 ) || ( a == 7 && b == 10 ) || ( a == 8 && b == 9 ) ) return 0.125 - 1/(8.*Nc2); if( ( a == 5 && b == 10 ) || ( a == 6 && b == 8 ) || ( a == 7 && b == 9 ) ) return (-1 + Nc4)/(8.*Nc2); } if ( basis == id33bar33bar8 ) { if( ( a == 0 && b == 0 ) || ( a == 1 && b == 1 ) || ( a == 2 && b == 2 ) || ( a == 3 && b == 3 ) ) return (Nc*(-1 + Nc2))/2.; if( ( a == 0 && b == 1 ) || ( a == 0 && b == 3 ) || ( a == 1 && b == 2 ) || ( a == 2 && b == 3 ) ) return (-1 + Nc2)/2.; if( ( a == 0 && b == 2 ) || ( a == 1 && b == 3 ) ) return 0; } } else { if ( a != b ) return 0.; if ( basis == id88 ) { return Nc2/4.; } if ( basis == id33bar ) { return Nc; } if ( basis == id888 ) { return Nc3/8.; } if ( basis == id33bar8 ) { return Nc2/2.; } if ( basis == id8888 ) { return Nc4/16.; } if ( basis == id33bar88 ) { return Nc3/4.; } if ( basis == id33bar33bar ) { return Nc2; } if ( basis == id88888 ) { return Nc5/32.; } if ( basis == id33bar888 ) { return Nc4/8.; } if ( basis == id33bar33bar8 ) { return Nc3/2.; } } throw Exception() << "SimpleColourBasis2::scalarProduct(): Cannot handle colour configuration" << Exception::runerror; } double SimpleColourBasis2::tMatrixElement(size_t i, size_t a, size_t b, const vector&, const vector& basis, #ifndef NDEBUG size_t k, size_t l, #else size_t , size_t , #endif const map& dict) const { // Check indices k and l assert( k == i ); assert( l == basis.size() ); // Check that dict is the standardMap assert( dict.size()+1 == basis.size() ); map::const_iterator tmp; for ( size_t ii = 0; ii < basis.size(); ii++ ) if ( ii != i ) { tmp = dict.find(ii); assert( tmp != dict.end() ); assert( tmp->second == ii ); } if ( id33bar.empty() ) makeIds(); if ( basis == id88 ) { if(i == 0 && a == 0 && b == 0) return -1; if(i == 0 && a == 1 && b == 0) return 1; if(i == 1 && a == 0 && b == 0) return 1; if(i == 1 && a == 1 && b == 0) return -1; return 0.; } if ( basis == id33bar ) { if(i == 0 && a == 0 && b == 0) return 1; if(i == 1 && a == 0 && b == 0) return -1; return 0.; } if ( basis == id888 ) { if(i == 0 && a == 3 && b == 0) return -1; if(i == 0 && a == 5 && b == 1) return -1; if(i == 0 && a == 7 && b == 0) return 1; if(i == 0 && a == 8 && b == 1) return 1; if(i == 1 && a == 4 && b == 0) return 1; if(i == 1 && a == 5 && b == 1) return 1; if(i == 1 && a == 6 && b == 1) return -1; if(i == 1 && a == 7 && b == 0) return -1; if(i == 2 && a == 3 && b == 0) return 1; if(i == 2 && a == 4 && b == 0) return -1; if(i == 2 && a == 6 && b == 1) return 1; if(i == 2 && a == 8 && b == 1) return -1; return 0.; } if ( basis == id33bar8 ) { if(i == 0 && a == 2 && b == 0) return 1; if(i == 1 && a == 1 && b == 0) return -1; if(i == 2 && a == 1 && b == 0) return 1; if(i == 2 && a == 2 && b == 0) return -1; return 0.; } if ( basis == id8888 ) { if(i == 0 && a == 2 && b == 2) return -1; if(i == 0 && a == 5 && b == 1) return -1; if(i == 0 && a == 8 && b == 0) return -1; if(i == 0 && a == 9 && b == 2) return 1; if(i == 0 && a == 10 && b == 1) return 1; if(i == 0 && a == 11 && b == 0) return 1; if(i == 0 && a == 20 && b == 3) return -1; if(i == 0 && a == 22 && b == 4) return -1; if(i == 0 && a == 26 && b == 5) return -1; if(i == 0 && a == 28 && b == 6) return -1; if(i == 0 && a == 32 && b == 7) return -1; if(i == 0 && a == 34 && b == 8) return -1; if(i == 0 && a == 38 && b == 3) return 1; if(i == 0 && a == 39 && b == 4) return 1; if(i == 0 && a == 40 && b == 5) return 1; if(i == 0 && a == 41 && b == 6) return 1; if(i == 0 && a == 42 && b == 7) return 1; if(i == 0 && a == 43 && b == 8) return 1; if(i == 1 && a == 2 && b == 2) return 1; if(i == 1 && a == 9 && b == 2) return -1; if(i == 1 && a == 13 && b == 0) return -1; if(i == 1 && a == 15 && b == 1) return -1; if(i == 1 && a == 16 && b == 0) return 1; if(i == 1 && a == 17 && b == 1) return 1; if(i == 1 && a == 24 && b == 3) return 1; if(i == 1 && a == 25 && b == 4) return 1; if(i == 1 && a == 27 && b == 5) return 1; if(i == 1 && a == 28 && b == 6) return 1; if(i == 1 && a == 29 && b == 6) return -1; if(i == 1 && a == 30 && b == 5) return -1; if(i == 1 && a == 33 && b == 7) return 1; if(i == 1 && a == 34 && b == 8) return 1; if(i == 1 && a == 35 && b == 8) return -1; if(i == 1 && a == 36 && b == 7) return -1; if(i == 1 && a == 38 && b == 3) return -1; if(i == 1 && a == 39 && b == 4) return -1; if(i == 2 && a == 5 && b == 1) return 1; if(i == 2 && a == 10 && b == 1) return -1; if(i == 2 && a == 13 && b == 0) return 1; if(i == 2 && a == 16 && b == 0) return -1; if(i == 2 && a == 18 && b == 2) return -1; if(i == 2 && a == 19 && b == 2) return 1; if(i == 2 && a == 21 && b == 3) return 1; if(i == 2 && a == 22 && b == 4) return 1; if(i == 2 && a == 23 && b == 4) return -1; if(i == 2 && a == 24 && b == 3) return -1; if(i == 2 && a == 30 && b == 5) return 1; if(i == 2 && a == 31 && b == 6) return 1; if(i == 2 && a == 32 && b == 7) return 1; if(i == 2 && a == 33 && b == 7) return -1; if(i == 2 && a == 35 && b == 8) return 1; if(i == 2 && a == 37 && b == 8) return -1; if(i == 2 && a == 40 && b == 5) return -1; if(i == 2 && a == 41 && b == 6) return -1; if(i == 3 && a == 8 && b == 0) return 1; if(i == 3 && a == 11 && b == 0) return -1; if(i == 3 && a == 15 && b == 1) return 1; if(i == 3 && a == 17 && b == 1) return -1; if(i == 3 && a == 18 && b == 2) return 1; if(i == 3 && a == 19 && b == 2) return -1; if(i == 3 && a == 20 && b == 3) return 1; if(i == 3 && a == 21 && b == 3) return -1; if(i == 3 && a == 23 && b == 4) return 1; if(i == 3 && a == 25 && b == 4) return -1; if(i == 3 && a == 26 && b == 5) return 1; if(i == 3 && a == 27 && b == 5) return -1; if(i == 3 && a == 29 && b == 6) return 1; if(i == 3 && a == 31 && b == 6) return -1; if(i == 3 && a == 36 && b == 7) return 1; if(i == 3 && a == 37 && b == 8) return 1; if(i == 3 && a == 42 && b == 7) return -1; if(i == 3 && a == 43 && b == 8) return -1; return 0.; } if ( basis == id33bar88 ) { if(i == 0 && a == 4 && b == 0) return 1; if(i == 0 && a == 9 && b == 1) return 1; if(i == 0 && a == 10 && b == 2) return 1; if(i == 1 && a == 4 && b == 0) return -1; if(i == 1 && a == 5 && b == 1) return -1; if(i == 1 && a == 7 && b == 2) return -1; if(i == 2 && a == 0 && b == 0) return -1; if(i == 2 && a == 1 && b == 0) return 1; if(i == 2 && a == 6 && b == 1) return 1; if(i == 2 && a == 7 && b == 2) return 1; if(i == 2 && a == 8 && b == 2) return -1; if(i == 2 && a == 9 && b == 1) return -1; if(i == 3 && a == 0 && b == 0) return 1; if(i == 3 && a == 1 && b == 0) return -1; if(i == 3 && a == 5 && b == 1) return 1; if(i == 3 && a == 6 && b == 1) return -1; if(i == 3 && a == 8 && b == 2) return 1; if(i == 3 && a == 10 && b == 2) return -1; return 0.; } if ( basis == id33bar33bar ) { if(i == 0 && a == 0 && b == 0) return 1; if(i == 0 && a == 1 && b == 1) return 1; if(i == 1 && a == 1 && b == 1) return -1; if(i == 1 && a == 2 && b == 0) return -1; if(i == 2 && a == 2 && b == 0) return 1; if(i == 2 && a == 3 && b == 1) return 1; if(i == 3 && a == 0 && b == 0) return -1; if(i == 3 && a == 3 && b == 1) return -1; return 0.; } throw Exception() << "SimpleColourBasis2::tMatrixElement(): Cannot handle colour configuration" << Exception::runerror; return 0.; } bool SimpleColourBasis2::colourConnected(const cPDVector& sub, const vector& basis, const pair& i, const pair& j, size_t a) const { if ( id33bar.empty() ) makeIds(); // translate process to basis ids map >::const_iterator trans = indexMap().find(sub); assert(trans != indexMap().end()); int idColoured = i.second ? j.first : i.first; idColoured = trans->second.find(idColoured)->second; int idAntiColoured = i.second ? i.first : j.first; idAntiColoured = trans->second.find(idAntiColoured)->second; if ( basis == id88 ) { return a == 0 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0)); } if ( basis == id33bar ) { return a == 0 && (idColoured == 0 && idAntiColoured == 1); } if ( basis == id888 ) { return (a == 0 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 1 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))); } if ( basis == id33bar8 ) { return a == 0 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1)); } if ( basis == id8888 ) { return (a == 0 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1))) || (a == 1 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1))) || (a == 2 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2))) || (a == 3 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 4 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 5 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 6 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))) || (a == 7 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 8 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))); } if ( basis == id33bar88 ) { return (a == 0 && ((idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 0 && idAntiColoured == 1))) || (a == 1 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 2 && idAntiColoured == 3))) || (a == 2 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 3 && idAntiColoured == 2))); } if ( basis == id33bar33bar ) { return (a == 0 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 2 && idAntiColoured == 1))) || (a == 1 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 2 && idAntiColoured == 3))); } if ( basis == id88888 ) { return (a == 0 && ((idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 1 && ((idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 2 && ((idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0))) || (a == 3 && ((idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))) || (a == 4 && ((idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 5 && ((idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0))) || (a == 6 && ((idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))) || (a == 7 && ((idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 8 && ((idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0))) || (a == 9 && ((idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))) || (a == 10 && ((idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 11 && ((idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 12 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1))) || (a == 13 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1))) || (a == 14 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1))) || (a == 15 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1))) || (a == 16 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0) || (idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1))) || (a == 17 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0) || (idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1))) || (a == 18 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2))) || (a == 19 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2))) || (a == 20 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0))) || (a == 21 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 22 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0))) || (a == 23 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 24 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 25 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 26 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0))) || (a == 27 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 28 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0))) || (a == 29 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))) || (a == 30 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 31 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))) || (a == 32 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0))) || (a == 33 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 34 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 0))) || (a == 35 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))) || (a == 36 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 37 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))) || (a == 38 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 39 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 40 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 0))) || (a == 41 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))) || (a == 42 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 0))) || (a == 43 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 1 && idAntiColoured == 0))); } if ( basis == id33bar888 ) { return (a == 0 && ((idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 0 && idAntiColoured == 1))) || (a == 1 && ((idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 0 && idAntiColoured == 1))) || (a == 2 && ((idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 0 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 1))) || (a == 3 && ((idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 0 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 1))) || (a == 4 && ((idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1))) || (a == 5 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 2 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 4))) || (a == 6 && ((idColoured == 0 && idAntiColoured == 2) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3))) || (a == 7 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 3 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 4))) || (a == 8 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 3 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 2))) || (a == 9 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 3 && idAntiColoured == 1) || (idColoured == 4 && idAntiColoured == 2) || (idColoured == 2 && idAntiColoured == 3))) || (a == 10 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 2 && idAntiColoured == 1) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 3 && idAntiColoured == 2))); } if ( basis == id33bar33bar8 ) { return (a == 0 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3) || (idColoured == 2 && idAntiColoured == 1))) || (a == 1 && ((idColoured == 0 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1) || (idColoured == 2 && idAntiColoured == 3))) || (a == 2 && ((idColoured == 0 && idAntiColoured == 3) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 1))) || (a == 3 && ((idColoured == 0 && idAntiColoured == 1) || (idColoured == 2 && idAntiColoured == 4) || (idColoured == 4 && idAntiColoured == 3))); } throw Exception() << "SimpleColourBasis2::colourConnected(): Cannot handle colour configuration" << Exception::runerror; return false; } map > > SimpleColourBasis2::basisList(const vector& basis) const { if ( id33bar.empty() ) makeIds(); map > > blist; vector > structures; vector structure; if ( basis == id88 ) { structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); blist[0] = structures; return blist; } if ( basis == id33bar ) { structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); blist[0] = structures; return blist; } if ( basis == id888 ) { structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(2); structures.push_back(structure); blist[0] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[1] = structures; return blist; } if ( basis == id33bar8 ) { structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[0] = structures; return blist; } if ( basis == id8888 ) { structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(1); structure.push_back(2); structures.push_back(structure); blist[0] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structures.push_back(structure); structure.clear(); structure.push_back(1); structure.push_back(3); structures.push_back(structure); blist[1] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); structure.clear(); structure.push_back(2); structure.push_back(3); structures.push_back(structure); blist[2] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(2); structure.push_back(3); structures.push_back(structure); blist[3] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(3); structure.push_back(2); structures.push_back(structure); blist[4] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(1); structure.push_back(3); structures.push_back(structure); blist[5] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(3); structure.push_back(1); structures.push_back(structure); blist[6] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(1); structure.push_back(2); structures.push_back(structure); blist[7] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[8] = structures; return blist; } if ( basis == id33bar88 ) { structures.clear(); structure.clear(); structure.push_back(2); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); blist[0] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(3); structure.push_back(1); structures.push_back(structure); blist[1] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[2] = structures; return blist; } if ( basis == id33bar33bar ) { structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[0] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); structure.clear(); structure.push_back(2); structure.push_back(3); structures.push_back(structure); blist[1] = structures; return blist; } if ( basis == id88888 ) { structures.clear(); structure.clear(); structure.push_back(3); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(2); structures.push_back(structure); blist[0] = structures; structures.clear(); structure.clear(); structure.push_back(2); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(3); structures.push_back(structure); blist[1] = structures; structures.clear(); structure.clear(); structure.push_back(2); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(4); structures.push_back(structure); blist[2] = structures; structures.clear(); structure.clear(); structure.push_back(3); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[3] = structures; structures.clear(); structure.clear(); structure.push_back(1); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(3); structures.push_back(structure); blist[4] = structures; structures.clear(); structure.clear(); structure.push_back(1); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(4); structures.push_back(structure); blist[5] = structures; structures.clear(); structure.clear(); structure.push_back(2); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(1); structures.push_back(structure); blist[6] = structures; structures.clear(); structure.clear(); structure.push_back(1); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(2); structures.push_back(structure); blist[7] = structures; structures.clear(); structure.clear(); structure.push_back(1); structure.push_back(2); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(4); structures.push_back(structure); blist[8] = structures; structures.clear(); structure.clear(); structure.push_back(2); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(1); structures.push_back(structure); blist[9] = structures; structures.clear(); structure.clear(); structure.push_back(1); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(2); structures.push_back(structure); blist[10] = structures; structures.clear(); structure.clear(); structure.push_back(1); structure.push_back(2); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(3); structures.push_back(structure); blist[11] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(1); structure.push_back(2); structure.push_back(3); structures.push_back(structure); blist[12] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(1); structure.push_back(2); structure.push_back(4); structures.push_back(structure); blist[13] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(1); structure.push_back(3); structure.push_back(2); structures.push_back(structure); blist[14] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structures.push_back(structure); structure.clear(); structure.push_back(1); structure.push_back(3); structure.push_back(4); structures.push_back(structure); blist[15] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(1); structure.push_back(4); structure.push_back(2); structures.push_back(structure); blist[16] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structures.push_back(structure); structure.clear(); structure.push_back(1); structure.push_back(4); structure.push_back(3); structures.push_back(structure); blist[17] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); structure.clear(); structure.push_back(2); structure.push_back(3); structure.push_back(4); structures.push_back(structure); blist[18] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); structure.clear(); structure.push_back(2); structure.push_back(4); structure.push_back(3); structures.push_back(structure); blist[19] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(2); structure.push_back(3); structure.push_back(4); structures.push_back(structure); blist[20] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(2); structure.push_back(4); structure.push_back(3); structures.push_back(structure); blist[21] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(3); structure.push_back(2); structure.push_back(4); structures.push_back(structure); blist[22] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(3); structure.push_back(4); structure.push_back(2); structures.push_back(structure); blist[23] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(4); structure.push_back(2); structure.push_back(3); structures.push_back(structure); blist[24] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structure.push_back(4); structure.push_back(3); structure.push_back(2); structures.push_back(structure); blist[25] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(1); structure.push_back(3); structure.push_back(4); structures.push_back(structure); blist[26] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(1); structure.push_back(4); structure.push_back(3); structures.push_back(structure); blist[27] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(3); structure.push_back(1); structure.push_back(4); structures.push_back(structure); blist[28] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(3); structure.push_back(4); structure.push_back(1); structures.push_back(structure); blist[29] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(4); structure.push_back(1); structure.push_back(3); structures.push_back(structure); blist[30] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(4); structure.push_back(3); structure.push_back(1); structures.push_back(structure); blist[31] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(1); structure.push_back(2); structure.push_back(4); structures.push_back(structure); blist[32] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(1); structure.push_back(4); structure.push_back(2); structures.push_back(structure); blist[33] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(2); structure.push_back(1); structure.push_back(4); structures.push_back(structure); blist[34] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(2); structure.push_back(4); structure.push_back(1); structures.push_back(structure); blist[35] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(4); structure.push_back(1); structure.push_back(2); structures.push_back(structure); blist[36] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(4); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[37] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(1); structure.push_back(2); structure.push_back(3); structures.push_back(structure); blist[38] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(1); structure.push_back(3); structure.push_back(2); structures.push_back(structure); blist[39] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(2); structure.push_back(1); structure.push_back(3); structures.push_back(structure); blist[40] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(2); structure.push_back(3); structure.push_back(1); structures.push_back(structure); blist[41] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(3); structure.push_back(1); structure.push_back(2); structures.push_back(structure); blist[42] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(3); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[43] = structures; return blist; } if ( basis == id33bar888 ) { structures.clear(); structure.clear(); structure.push_back(2); structure.push_back(3); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); blist[0] = structures; structures.clear(); structure.clear(); structure.push_back(2); structure.push_back(4); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); blist[1] = structures; structures.clear(); structure.clear(); structure.push_back(3); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[2] = structures; structures.clear(); structure.clear(); structure.push_back(2); structure.push_back(4); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(1); structures.push_back(structure); blist[3] = structures; structures.clear(); structure.clear(); structure.push_back(2); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(1); structures.push_back(structure); blist[4] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(3); structure.push_back(4); structure.push_back(1); structures.push_back(structure); blist[5] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(2); structure.push_back(4); structure.push_back(3); structure.push_back(1); structures.push_back(structure); blist[6] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(2); structure.push_back(4); structure.push_back(1); structures.push_back(structure); blist[7] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structure.push_back(4); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[8] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(2); structure.push_back(3); structure.push_back(1); structures.push_back(structure); blist[9] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(3); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[10] = structures; return blist; } if ( basis == id33bar33bar8 ) { structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(2); structure.push_back(1); structures.push_back(structure); blist[0] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(4); structure.push_back(1); structures.push_back(structure); structure.clear(); structure.push_back(2); structure.push_back(3); structures.push_back(structure); blist[1] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(3); structures.push_back(structure); structure.clear(); structure.push_back(2); structure.push_back(4); structure.push_back(1); structures.push_back(structure); blist[2] = structures; structures.clear(); structure.clear(); structure.push_back(0); structure.push_back(1); structures.push_back(structure); structure.clear(); structure.push_back(2); structure.push_back(4); structure.push_back(3); structures.push_back(structure); blist[3] = structures; return blist; } throw Exception() << "SimpleColourBasis2::basisList(): Cannot handle colour configuration" << Exception::runerror; return blist; } void SimpleColourBasis2::makeIds() const { id88 = vector(2,PDT::Colour8); id33bar.push_back(PDT::Colour3); id33bar.push_back(PDT::Colour3bar); id888 = vector(3,PDT::Colour8); id33bar8 = id33bar; id33bar8.push_back(PDT::Colour8); id8888 = vector(4,PDT::Colour8); id33bar88 = id33bar8; id33bar88.push_back(PDT::Colour8); id33bar33bar = id33bar; id33bar33bar.push_back(PDT::Colour3); id33bar33bar.push_back(PDT::Colour3bar); id88888 = vector(5,PDT::Colour8); id33bar888 = id33bar88; id33bar888.push_back(PDT::Colour8); id33bar33bar8 = id33bar33bar; id33bar33bar8.push_back(PDT::Colour8); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void SimpleColourBasis2::persistentOutput(PersistentOStream &) const {} void SimpleColourBasis2::persistentInput(PersistentIStream &, int) {} // *** 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 describeHerwigSimpleColourBasis2("Herwig::SimpleColourBasis2", "Herwig.so"); void SimpleColourBasis2::Init() { static ClassDocumentation documentation ("SimpleColourBasis2 implements the colour algebra needed for " "processes with four coloured legs at NLO. It mainly " "serves as an example for the general ColourBasis interface."); } diff --git a/MatrixElement/Matchbox/Utility/SimpleColourBasis2.h b/MatrixElement/Matchbox/Utility/SimpleColourBasis2.h --- a/MatrixElement/Matchbox/Utility/SimpleColourBasis2.h +++ b/MatrixElement/Matchbox/Utility/SimpleColourBasis2.h @@ -1,233 +1,218 @@ // -*- C++ -*- // // SimpleColourBasis2.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_SimpleColourBasis2_H #define Herwig_SimpleColourBasis2_H // // This is the declaration of the SimpleColourBasis2 class. // #include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief SimpleColourBasis2 implements the colour algebra needed for * processes with four coloured legs at NLO. It mainly serves as an * example for the general ColourBasis interface. * */ class SimpleColourBasis2: public ColourBasis { public: - /** @name Standard constructors and destructors. */ - //@{ - /** - * The default constructor. - */ - SimpleColourBasis2(); - - /** - * The destructor. - */ - virtual ~SimpleColourBasis2(); - //@} - -public: - /** * Prepare the basis for the normal ordered legs and return the * dimensionality of the basis. */ virtual size_t prepareBasis(const vector&); /** * Return the scalar product of basis tensors labelled a and b in * the basis used for the given normal ordered legs. */ virtual double scalarProduct(size_t a, size_t b, const vector& abBasis) const; /** * Return the matrix element of a colour charge * between basis tensors a and b, with * respect to aBasis and bBasis */ virtual double tMatrixElement(size_t i, size_t a, size_t b, const vector& aBasis, const vector& bBasis, size_t k, size_t l, const map& dict) const; /* * Temporary to make it compile */ virtual double sMatrixElement(size_t, size_t, size_t, const vector&, const vector&, size_t, size_t, const map& ) const { assert( 0 == 1 ); return 0; } /** * Return true, if this colour basis supports gluon splittings. */ virtual bool canSplitGluons() const { return false; } /** * Return true, if a large-N colour connection exists for the * given external legs and basis tensor. */ virtual bool colourConnected(const cPDVector&, const vector&, const pair&, const pair&, size_t) const; /** * Return true, if the colour basis is capable of assigning colour * flows. */ virtual bool haveColourFlows() const { return true; } /** * Create ids for bases */ void makeIds() const; /** * Return a map of basis tensor indices to vectors identifying a * certain ordering corresponding to the given colour structure. May * not be supported by all colour basis implementations. */ virtual map > > basisList(const vector&) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * id for 88 */ mutable vector id88; /** * id for 33bar */ mutable vector id33bar; /** * id for 888 */ mutable vector id888; /** * id for 33bar8 */ mutable vector id33bar8; /** * id for 8888 */ mutable vector id8888; /** * id for 33bar88 */ mutable vector id33bar88; /** * id for 33bar33bar */ mutable vector id33bar33bar; /** * id for 88888 */ mutable vector id88888; /** * id for 33bar888 */ mutable vector id33bar888; /** * id for 33bar33bar8 */ mutable vector id33bar33bar8; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ SimpleColourBasis2 & operator=(const SimpleColourBasis2 &) = delete; }; } #endif /* Herwig_SimpleColourBasis2_H */ diff --git a/MatrixElement/Matchbox/Utility/Tree2toNGenerator.cc b/MatrixElement/Matchbox/Utility/Tree2toNGenerator.cc --- a/MatrixElement/Matchbox/Utility/Tree2toNGenerator.cc +++ b/MatrixElement/Matchbox/Utility/Tree2toNGenerator.cc @@ -1,490 +1,488 @@ // -*- C++ -*- // // Tree2toNGenerator.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 Tree2toNGenerator class. // #include "Tree2toNGenerator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Utilities/StringUtils.h" using namespace Herwig; Tree2toNGenerator::Tree2toNGenerator() : maxOrderGs(0), maxOrderGem(0), prepared(false) {} -Tree2toNGenerator::~Tree2toNGenerator() {} - IBPtr Tree2toNGenerator::clone() const { return new_ptr(*this); } IBPtr Tree2toNGenerator::fullclone() const { return new_ptr(*this); } vector::ptr> Tree2toNGenerator:: generate(const PDVector& legs, unsigned int orderInGs, unsigned int orderInGem) { vector::ptr> res; list > prog = clusterAll(legs,orderInGs,orderInGem); int count = 1; for ( auto & d : prog ) { assert(d.size() == 1); Tree2toNDiagram diag = d.front().generate(count); if ( !spaceLikeAllowed.empty() ) { map counts; for ( int k = 1; k < diag.nSpace()-1; ++k ) counts[diag.allPartons()[k]] += 1; for ( auto & m : spaceLikeAllowed ) { m.reset(); for ( auto const & c : counts ) m.add(c.first,c.second); } bool failed = false; for ( auto const & m : spaceLikeAllowed) { if ( !m.check() ) { failed = true; break; } } if ( failed ) continue; } if ( !timeLikeAllowed.empty() ) { map counts; int all = diag.allPartons().size(); for ( int k = diag.nSpace(); k < all; ++k ) { if ( diag.children(k).first < 0 ) continue; counts[diag.allPartons()[k]] += 1; } for ( auto & m : timeLikeAllowed ) { m.reset(); for ( auto const & c : counts ) m.add(c.first,c.second); } bool failed = false; for ( auto const & m : timeLikeAllowed ) { if ( !m.check() ) { failed = true; break; } } if ( failed ) continue; } bool internalVeto = false; set external; int nex = diag.partons().size(); for ( int i = 0; i < nex; ++i ) { external.insert(diag.diagramId(i)); } int n = diag.allPartons().size(); for ( int i = 0; i < n; ++i ) { if ( external.find(i) != external.end() ) continue; if ( find(excludeInternal().begin(), excludeInternal().end(), diag.allPartons()[i]) != excludeInternal().end() ) { internalVeto = true; break; } } if ( internalVeto ) continue; bool gotit = false; for ( auto const & d : res) { map checkPermutation; if ( diag.isSame(d,checkPermutation) ) { gotit = true; for ( auto const & p : checkPermutation ) if ( p.first != p.second ) gotit = false; if ( gotit ) break; } } if ( !gotit ) { res.push_back(new_ptr(diag)); ++count; } } return res; } list > Tree2toNGenerator:: cluster(const vector& children, unsigned int orderInGs, unsigned int orderInGem) const { list > res; bool externalCluster = children[1].externalId != -1; if ( children.size() == 3 ) { for ( auto const & v : theVertices ) { if ( v->getNpoint() != 3 ) continue; if ( find(theExcludeVertices.begin(), theExcludeVertices.end(), v) != theExcludeVertices.end() ) continue; bool noMatch = v->orderInGs() != int(orderInGs) || v->orderInGem() != int(orderInGem) || !v->isIncoming(children[0].parent); long idij = children[0].parent->id(); long idi = children[2].parent->id(); long idj = children[1].parent->id(); if ( externalCluster && children[1].parent->CC() ) idj = -idj; if ( children[0].parent->CC() ) idij = -idij; if ( !externalCluster ) noMatch |= !v->isOutgoing(children[1].parent) || !v->isOutgoing(children[2].parent); else noMatch |= !v->isIncoming(children[1].parent) || !v->isOutgoing(children[2].parent); noMatch |= !( v->allowed(idij,idi,idj) || v->allowed(idj,idij,idi) || v->allowed(idi,idj,idij) || v->allowed(idij,idj,idi) || v->allowed(idi,idij,idj) || v->allowed(idj,idi,idij) ); if ( noMatch ) continue; Vertex last; last.spacelike = true; last.parent = children[0].parent; last.externalId = 0; last.children.push_back(children[1]); last.children.push_back(children[2]); res.push_back(vector(1,last)); // only one possible break; } return res; } // spacelike clusterings (cluster on second one) for ( size_t i = 2; i < children.size(); ++i ) { for ( auto const & v : theVertices ) { if ( v->getNpoint() != 3 ) continue; if ( find(theExcludeVertices.begin(), theExcludeVertices.end(), v) != theExcludeVertices.end() ) continue; bool noMatch = false; noMatch |= v->orderInGs() != int(orderInGs) || v->orderInGem() != int(orderInGem); if ( !externalCluster ) noMatch |= !v->isOutgoing(children[1].parent) || !v->isOutgoing(children[i].parent); else noMatch |= !v->isIncoming(children[1].parent) || !v->isOutgoing(children[i].parent); if ( noMatch ) continue; long idi = children[i].parent->id(); long idj = children[1].parent->id(); if ( externalCluster && children[1].parent->CC() ) idj = -idj; for ( set::const_iterator pij = v->outgoing().begin(); pij != v->outgoing().end() ; ++pij ) { long idij = (**pij).id(); if ( v->allowed(idij,idi,idj) || v->allowed(idj,idij,idi) || v->allowed(idi,idj,idij) || v->allowed(idij,idj,idi) || v->allowed(idi,idij,idj) || v->allowed(idj,idi,idij) ) { PDPtr dij = (**pij).CC() ? (**pij).CC() : *pij; vector cled; for ( size_t k = 0; k < children.size(); ++k ) { if ( k != 1 && k != i ) cled.push_back(children[k]); if ( k == 1 ) { Vertex merge; merge.children.push_back(children[1]); merge.children.push_back(children[i]); merge.parent = dij; merge.spacelike = true; cled.push_back(merge); } if ( k == i ) continue; } res.push_back(cled); } } } } // timelike clusterings for ( size_t i = 2; i < children.size(); ++i ) { for ( size_t j = i+1; j < children.size(); ++j ) { for ( auto const & v : theVertices ) { if ( v->getNpoint() != 3 ) continue; if ( find(theExcludeVertices.begin(), theExcludeVertices.end(), v) != theExcludeVertices.end() ) continue; if ( v->orderInGs() != int(orderInGs) || v->orderInGem() != int(orderInGem) || !v->isOutgoing(children[i].parent) || !v->isOutgoing(children[j].parent) ) continue; long idi = children[i].parent->id(); long idj = children[j].parent->id(); for ( set::const_iterator pij = v->outgoing().begin(); pij != v->outgoing().end() ; ++pij ) { long idij = (**pij).id(); if ( v->allowed(idij,idi,idj) || v->allowed(idj,idij,idi) || v->allowed(idi,idj,idij) || v->allowed(idij,idj,idi) || v->allowed(idi,idij,idj) || v->allowed(idj,idi,idij) ) { PDPtr dij = (**pij).CC() ? (**pij).CC() : *pij; vector cled; for ( size_t k = 0; k < children.size(); ++k ) { if ( k != i && k != j ) cled.push_back(children[k]); if ( k == i ) { Vertex merge; merge.children.push_back(children[i]); merge.children.push_back(children[j]); merge.parent = dij; merge.spacelike = false; cled.push_back(merge); } if ( k == j ) continue; } res.push_back(cled); } } } } } return res; } list > Tree2toNGenerator:: clusterAll(const list >& current, unsigned int orderInGs, unsigned int orderInGem) const { list > res; for ( list >::const_iterator c = current.begin(); c != current.end(); ++c ) { if ( c->size() == 1 ) { if ( orderInGs == 0 && orderInGem == 0 ) res.push_back(*c); continue; } for ( unsigned int gs = 0; gs <= maxOrderGs; ++gs ) for ( unsigned int gem = 0; gem <= maxOrderGem; ++gem ) { if ( gs == 0 && gem == 0 ) continue; if ( gs > orderInGs || gem > orderInGem ) continue; list > next = cluster(*c,gs,gem); if ( next.empty() ) continue; list > cled = clusterAll(next,orderInGs-gs,orderInGem-gem); copy(cled.begin(),cled.end(),back_inserter(res)); } } return res; } list > Tree2toNGenerator:: clusterAll(const PDVector& external, unsigned int orderInGs, unsigned int orderInGem) { if ( !prepared ) { for ( auto & v : theVertices ) { if ( find(theExcludeVertices.begin(), theExcludeVertices.end(), v) != theExcludeVertices.end() ) continue; v->init(); maxOrderGs = max(maxOrderGs,v->orderInGs()); maxOrderGem = max(maxOrderGem,v->orderInGem()); } for ( auto & m : spaceLikeAllowed) m.rebind(this); for ( auto & m : timeLikeAllowed ) m.rebind(this); prepared = true; } vector legs; for ( unsigned int k = 0; k < external.size(); ++k ) { Vertex v; v.parent = external[k]; v.externalId = k; v.spacelike = k < 2; legs.push_back(v); } list > firstlegs; firstlegs.push_back(legs); return clusterAll(firstlegs,orderInGs,orderInGem); } string Tree2toNGenerator::doSpaceLikeRange(string range) { if ( theRestrictLines.empty() ) return "No particle data specified to restrict internal lines."; vector bounds = StringUtils::split(range); if ( bounds.empty() || bounds.size() > 2 ) return "Need to specify a minimum, or a minimum and maximum number of internal lines."; pair irange(0,-1); istringstream in1(bounds[0]); in1 >> irange.first; if ( bounds.size() == 2 ) { istringstream in2(bounds[1]); in2 >> irange.second; } else { irange.second = irange.first; } if ( irange.second >= 0 && irange.first > irange.second ) return "invalid range specified"; spaceLikeAllowed.push_back(LineMatcher(theRestrictLines,irange)); return ""; } string Tree2toNGenerator::doTimeLikeRange(string range) { if ( theRestrictLines.empty() ) return "No particle data specified to restrict internal lines."; vector bounds = StringUtils::split(range); if ( bounds.empty() || bounds.size() > 2 ) return "Need to specify a minimum, or a minimum and maximum number of internal lines."; pair irange(0,-1); istringstream in1(bounds[0]); in1 >> irange.first; if ( bounds.size() == 2 ) { istringstream in2(bounds[1]); in2 >> irange.second; } else { irange.second = irange.first; } if ( irange.second >= 0 && irange.first > irange.second ) return "invalid range specified"; timeLikeAllowed.push_back(LineMatcher(theRestrictLines,irange)); return ""; } string Tree2toNGenerator::doClearRestrictLines(string) { theRestrictLines.clear(); return ""; } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void Tree2toNGenerator::persistentOutput(PersistentOStream & os) const { os << theVertices << theExcludeInternal << maxOrderGs << maxOrderGem << prepared << theExcludeVertices << spaceLikeAllowed << timeLikeAllowed; } void Tree2toNGenerator::persistentInput(PersistentIStream & is, int) { is >> theVertices >> theExcludeInternal >> maxOrderGs >> maxOrderGem >> prepared >> theExcludeVertices >> spaceLikeAllowed >> timeLikeAllowed; } // *** 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 describeHerwigTree2toNGenerator("Herwig::Tree2toNGenerator", "Herwig.so"); void Tree2toNGenerator::Init() { static ClassDocumentation documentation ("Generate Tree2toNDiagrams for a given process."); static RefVector interfaceVertices ("Vertices", "All vertices to consider.", &Tree2toNGenerator::theVertices, -1, false, false, true, false, false); static RefVector interfaceExcludeVertices ("ExcludeVertices", "The vertices to exclude.", &Tree2toNGenerator::theExcludeVertices, -1, false, false, true, false, false); static RefVector interfaceExcludeInternal ("ExcludeInternal", "Particles to be exluded from becoming internal lines.", &Tree2toNGenerator::theExcludeInternal, -1, false, false, true, false, false); static RefVector interfaceRestrictLines ("RestrictLines", "Particles to be exluded from becoming internal lines.", &Tree2toNGenerator::theRestrictLines, -1, false, false, true, false, false); static Command interfaceSpaceLikeRange ("SpaceLikeRange", "Limit the number of spacelike occurences of the specified particle.", &Tree2toNGenerator::doSpaceLikeRange, false); static Command interfaceTimeLikeRange ("TimeLikeRange", "Limit the number of timelike occurences of the specified particle.", &Tree2toNGenerator::doTimeLikeRange, false); static Command interfaceClearRestrictLines ("ClearRestrictLines", "Clear the container of lines to be considered for restrictions.", &Tree2toNGenerator::doClearRestrictLines, false); } diff --git a/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h b/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h --- a/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h +++ b/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h @@ -1,438 +1,430 @@ // -*- C++ -*- // // Tree2toNGenerator.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_Tree2toNGenerator_H #define Herwig_Tree2toNGenerator_H // // This is the declaration of the Tree2toNGenerator class. // #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/Helicity/Vertex/VertexBase.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" namespace Herwig { using namespace ThePEG; /** * \ingroup Matchbox * \author Simon Platzer * * \brief Generate Tree2toNDiagrams for a given process. * * @see \ref Tree2toNGeneratorInterfaces "The interfaces" * defined for Tree2toNGenerator. */ class Tree2toNGenerator: public HandlerBase { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ Tree2toNGenerator(); - /** - * The destructor. - */ - virtual ~Tree2toNGenerator(); - //@} - public: /** * Generate all diagrams for the given process. */ vector::ptr> generate(const PDVector&, unsigned int orderInGs, unsigned int orderInGem); typedef vector::ptr> VertexVector; /** * Access the vertices */ VertexVector& vertices() { return theVertices; } /** * Return the vertices */ const VertexVector& vertices() const { return theVertices; } /** * Access the particles to be excluded from internal lines */ PDVector& excludeInternal() { return theExcludeInternal; } /** * Return the particles to be excluded from internal lines */ const PDVector& excludeInternal() const { return theExcludeInternal; } public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); public: /** * A node in internally used trees. */ struct Vertex { /** * The outgoing particles. If this is a spacelike node, the first * child is considered the next spacelike (or second incoming) * line. If children are empty, this is an external line. */ vector children; /** * The incoming line at this node. */ PDPtr parent; /** * True, if this is spacelike node. */ bool spacelike; /** * The external leg id. */ int externalId; /** * The parent diagram id. */ int parentId; /** * The default constructor. */ Vertex() : spacelike(false), externalId(-1), parentId(-1) {} /** * Debug printout. */ void print(ostream& os, const string& prefix = "") const { os << prefix << parent->PDGName() << "[" << (spacelike ? "s" : "t") << "] ("; if ( externalId < 0 ) os << "x)\n"; else os << externalId << ")\n"; if ( !children.empty() ) { os << prefix << "|__\n"; children[0].print(os,prefix + "| "); os << prefix << "|__\n"; children[1].print(os,prefix + "| "); } } /** * Count the number of spacelike lines */ int nspace() const { if ( children.empty() ) return 1; int ret = 1; ret += children[0].nspace(); return ret; } /** * Update diagram returning a map of external ids to diagram id * parents. */ void update(Tree2toNDiagram& diag, map >& outgoing, int& lastUsed) { if ( externalId == 0 ) { assert(lastUsed==0); ++lastUsed; diag.operator,(parent); children[0].parentId = lastUsed; children[1].parentId = lastUsed; children[0].update(diag,outgoing,lastUsed); children[1].update(diag,outgoing,lastUsed); for ( map >::iterator out = outgoing.begin(); out != outgoing.end(); ++out ) { diag.operator,(out->second.first); diag.operator,(out->second.second); } return; } if ( spacelike ) { ++lastUsed; diag.operator,(parent); if ( externalId == 1 ) return; children[0].parentId = lastUsed; children[1].parentId = lastUsed; children[0].update(diag,outgoing,lastUsed); children[1].update(diag,outgoing,lastUsed); return; } if ( children.empty() ) { outgoing[externalId] = make_pair(parentId,parent); return; } diag.operator,(parentId); diag.operator,(parent); ++lastUsed; children[0].parentId = lastUsed; children[1].parentId = lastUsed; children[0].update(diag,outgoing,lastUsed); children[1].update(diag,outgoing,lastUsed); } /** * Generate a diagram of given id. */ Tree2toNDiagram generate(int id) { int nsp = nspace(); Tree2toNDiagram res(nsp); int diagid = 0; map > out; update(res,out,diagid); res.operator,(-id); return res; } }; /** * For the given set of trees determine all allowed clusterings. */ list > cluster(const vector& children, unsigned int orderInGs, unsigned int orderInGem) const; /** * For the given set of outgoing lines cluster recursively. */ list > clusterAll(const list >& current, unsigned int orderInGs, unsigned int orderInGem) const; /** * For the given set of outgoing lines cluster recursively. */ list > clusterAll(const PDVector& external, unsigned int orderInGs, unsigned int orderInGem); /** * Helper for topology restrictions */ struct LineMatcher { /** * The group of lines to be considered */ set particles; /** * The range allowed */ pair range; /** * The current count */ int count; /** * Default constructor */ LineMatcher() : range(0,0), count(0) {} /** * Construct given particles and a range */ LineMatcher(const PDVector& p, const pair& r) : range(r), count(0) { copy(p.begin(),p.end(),inserter(particles,particles.begin())); } /** * Rebind the particle data pointers */ void rebind(Tree2toNGenerator* g) { set oldp = particles; particles.clear(); for ( set::const_iterator p = oldp.begin(); p != oldp.end(); ++p ) particles.insert(g->getParticleData((**p).id())); } /** * Reset this matcher */ void reset() { count = 0; } /** * Count the given multiplicity */ void add(tcPDPtr p, int n) { if ( particles.find(p) == particles.end() ) return; count += n; } /** * Ceck if restrictions are met */ bool check() const { return count >= range.first && count <= range.second; } }; 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 vertices to be used. */ VertexVector theVertices; /** * The particles to be excluded from internal lines */ PDVector theExcludeInternal; /** * Maximum order in gs to consider. */ unsigned int maxOrderGs; /** * Maximum order in gem to consider. */ unsigned int maxOrderGem; /** * Wether or not the generator has been prepared */ bool prepared; /** * The vertices to be excluded. */ VertexVector theExcludeVertices; /** * Minimal and maximal occurences of spacelike internal lines */ vector spaceLikeAllowed; /** * Minimal and maximal occurences of timelike internal lines */ vector timeLikeAllowed; /** * The next particle for which internal lines need to be restricted */ PDVector theRestrictLines; /** * Command to set an allowed range of spacelike internal lines */ string doSpaceLikeRange(string); /** * Command to set an allowed range of timelike internal lines */ string doTimeLikeRange(string); /** * Command to clear the restrict lines container */ string doClearRestrictLines(string); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ Tree2toNGenerator & operator=(const Tree2toNGenerator &) = delete; }; inline PersistentOStream& operator<<(PersistentOStream& os, const Tree2toNGenerator::LineMatcher& m) { os << m.particles << m.range << m.count; return os; } inline PersistentIStream& operator>>(PersistentIStream& is, Tree2toNGenerator::LineMatcher& m) { is >> m.particles >> m.range >> m.count; return is; } } #endif /* Herwig_Tree2toNGenerator_H */