diff --git a/Models/General/DecayConstructor.cc b/Models/General/DecayConstructor.cc --- a/Models/General/DecayConstructor.cc +++ b/Models/General/DecayConstructor.cc @@ -1,166 +1,166 @@ // -*- C++ -*- // // DecayConstructor.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 DecayConstructor class. // #include "DecayConstructor.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Repository/BaseRepository.h" #include using namespace Herwig; using namespace ThePEG; IBPtr DecayConstructor::clone() const { return new_ptr(*this); } IBPtr DecayConstructor::fullclone() const { return new_ptr(*this); } void DecayConstructor::persistentOutput(PersistentOStream & os) const { os << NBodyDecayConstructors_ << QEDGenerator_; } void DecayConstructor::persistentInput(PersistentIStream & is, int) { is >> NBodyDecayConstructors_ >> QEDGenerator_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigDecayConstructor("Herwig::DecayConstructor", "Herwig.so"); void DecayConstructor::Init() { static ClassDocumentation documentation ("There is no documentation for the TwoBodyDecayConstructor class"); static RefVector interfaceNBodyDecayConstructors ("NBodyDecayConstructors", "Vector of references to NBodyDecayConstructors", &DecayConstructor::NBodyDecayConstructors_, -1, false, false, true, false, false); static ParVector interfaceDisableModes ("DisableModes", "A list of decay modes to disable", &DecayConstructor::_disableDMTags, -1, string(""), string(""), string(""), false, false, Interface::nolimits); static Reference interfaceQEDGenerator ("QEDGenerator", "Object to generate QED radiation in particle decays", &DecayConstructor::QEDGenerator_, false, false, true, true, false); } /** A helper function for for_each to sort the decay mode tags into the * standard order. */ namespace { void adjustFSOrder(string & tag) { string::size_type sep = tag.find(">"); string head = tag.substr(0, sep + 1); string products = tag.substr(sep + 1); OrderedParticles finalstate; bool loopbreak(true); while ( loopbreak ) { sep = products.find(","); string child; if( sep != string::npos ) { child = products.substr(0, sep); products = products.substr(sep + 1); } else { child = string(products.begin(), products.end() - 1); loopbreak = false; } PDPtr p = BaseRepository::GetObject (string("/Herwig/Particles/" + child)); if( p ) finalstate.insert(p); } if( finalstate.empty() ) return; tag = head; OrderedParticles::const_iterator iend = finalstate.end(); OrderedParticles::size_type count(0), npr(finalstate.size()); for( OrderedParticles::const_iterator it = finalstate.begin(); it != iend; ++it ) { tag += (**it).name(); if( ++count != npr ) tag += string(","); } tag += string(";"); } } namespace { /// Helper function for sorting by number of outgoing lines inline bool orderNBodyConstructors(tNBodyDecayConstructorBasePtr a, tNBodyDecayConstructorBasePtr b) { return a->numBodies() < b->numBodies(); } } void DecayConstructor::doinit() { Interfaced::doinit(); //Need to check that the stored decay mode tags have the //products in the standard order for_each( _disableDMTags.begin(), _disableDMTags.end(), adjustFSOrder ); sort(NBodyDecayConstructors_.begin(), NBodyDecayConstructors_.end(), orderNBodyConstructors); } void DecayConstructor::createDecayers(const PDVector & particles, double minBR) { _minBR = minBR; if ( particles.empty() || NBodyDecayConstructors_.empty() ) return; // turn the vector into a set to avoid duplicates - set particleSet(particles.begin(),particles.end()); + set particleSet(particles.begin(),particles.end()); // remove any antiparticles for(set::iterator it=particleSet.begin();it!=particleSet.end();++it) { PDPtr cc = (**it).CC(); if(!cc) continue; set::iterator ic = particleSet.find(cc); if(ic!=particleSet.end()) particleSet.erase(ic); } // set the decay list in the NBodyDecayConstructors typedef vector::iterator NBDecayIterator; NBDecayIterator it = NBodyDecayConstructors_.begin(); NBDecayIterator iend = NBodyDecayConstructors_.end(); for( ; it != iend; ++it ) { (**it).init(); (**it).decayConstructor(this); (**it).DecayList(particleSet); } } bool DecayConstructor::disableDecayMode(string tag) const { if( _disableDMTags.empty() ) return false; vector::const_iterator dit = _disableDMTags.begin(); vector::const_iterator dend = _disableDMTags.end(); bool disable(false); for( ; dit != dend; ++dit ) { if( *dit == tag ) { disable = true; break; } } return disable; } diff --git a/Models/General/FourBodyDecayConstructor.cc b/Models/General/FourBodyDecayConstructor.cc --- a/Models/General/FourBodyDecayConstructor.cc +++ b/Models/General/FourBodyDecayConstructor.cc @@ -1,259 +1,259 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the FourBodyDecayConstructor class. // #include "FourBodyDecayConstructor.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/RefVector.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/PDT/DecayMode.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "Herwig/Decay/General/GeneralFourBodyDecayer.h" #include "DecayConstructor.h" using namespace Herwig; FourBodyDecayConstructor::~FourBodyDecayConstructor() {} IBPtr FourBodyDecayConstructor::clone() const { return new_ptr(*this); } IBPtr FourBodyDecayConstructor::fullclone() const { return new_ptr(*this); } void FourBodyDecayConstructor::persistentOutput(PersistentOStream & os) const { os << interOpt_ << widthOpt_ << particles_; } void FourBodyDecayConstructor::persistentInput(PersistentIStream & is, int) { is >> interOpt_ >> widthOpt_ >> particles_; } DescribeClass describeFourBodyDecayConstructor("Herwig::FourBodyDecayConstructor","Herwig.so"); void FourBodyDecayConstructor::Init() { static ClassDocumentation documentation ("The FourBodyDecayConstructor class implements a small number" " of 4-body decays in general models"); static Switch interfaceWidthOption ("WidthOption", "Option for the treatment of the widths of the intermediates", &FourBodyDecayConstructor::widthOpt_, 1, false, false); static SwitchOption interfaceWidthOptionFixed (interfaceWidthOption, "Fixed", "Use fixed widths", 1); static SwitchOption interfaceWidthOptionRunning (interfaceWidthOption, "Running", "Use running widths", 2); static SwitchOption interfaceWidthOptionZero (interfaceWidthOption, "Zero", "Set the widths to zero", 3); static Switch interfaceIntermediateOption ("IntermediateOption", "Option for the inclusion of intermediates in the event", &FourBodyDecayConstructor::interOpt_, 0, false, false); static SwitchOption interfaceIntermediateOptionAlways (interfaceIntermediateOption, "Always", "Always include the intermediates", 1); static SwitchOption interfaceIntermediateOptionNever (interfaceIntermediateOption, "Never", "Never include the intermediates", 2); static SwitchOption interfaceIntermediateOptionOnlyIfOnShell (interfaceIntermediateOption, "OnlyIfOnShell", "Only if there are on-shell diagrams", 0); static RefVector interfaceParticles ("Particles", "Particles to override the choice in the DecayConstructor for 4-body decays," " if empty the defaults from the DecayConstructor are used.", &FourBodyDecayConstructor::particles_, -1, false, false, true, true, false); static Switch interfaceParticleType ("ParticleType", "Which types of particles to calculate four body decay modes for", &FourBodyDecayConstructor::particleType_, false, false, false); static SwitchOption interfaceParticleTypeStable (interfaceParticleType, "Stable", "Only calculate four-body decays in no 2/3 body modes", false); static SwitchOption interfaceParticleTypeAll (interfaceParticleType, "All", "Calculate 4-body modes for all particles", true); } -void FourBodyDecayConstructor::DecayList(const set & particles) { +void FourBodyDecayConstructor::DecayList(const set & particles) { if( particles.empty() ) return; - set new_particles; - for(set::const_iterator it=particles.begin();it!=particles.end();++it) { + set new_particles; + for(set::const_iterator it=particles.begin();it!=particles.end();++it) { if(!particles_.empty() && find(particles_.begin(),particles_.end(),*it)==particles_.end()) continue; if(!(**it).stable()&&!particleType_) continue; new_particles.insert(*it); } if(!new_particles.empty()) NBodyDecayConstructorBase::DecayList(new_particles); } void FourBodyDecayConstructor:: createDecayMode(vector & diagrams, bool possibleOnShell, double symfac) { // some basic checks for the modes we are interested in // only looking at scalars if(diagrams[0].incoming->iSpin()!=PDT::Spin0) return; // which decay to 4 fermions unsigned int nferm=0; for(OrderedParticles::const_iterator it=diagrams[0].outgoing.begin(); it!=diagrams[0].outgoing.end();++it) { if((**it).iSpin()==PDT::Spin1Half) ++nferm; } if(nferm!=4) return; // check for on-shell intermediates bool inter = interOpt_ == 1 || (interOpt_ == 0 && possibleOnShell); // incoming particle tPDPtr inpart = diagrams[0].incoming; // outgoing particles OrderedParticles outgoing=diagrams[0].outgoing; // incoming particle is now unstable inpart->stable(false); // construct the tag for the decay mode string tag = inpart->name() + "->"; for(OrderedParticles::const_iterator it = outgoing.begin(); it != outgoing.end(); ++it) { if(it!=outgoing.begin()) tag += ","; tag += (**it).name(); } tag += ";"; tDMPtr dm = generator()->findDecayMode(tag); // create mode if needed if( createDecayModes() && (!dm || inpart->id() == ParticleID::h0) ) { // create the decayer GeneralFourBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac); if(!decayer) { if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for " << tag << " so mode not created\n"; return; } // create the decay mode tDMPtr ndm = generator()->preinitCreateDecayMode(tag); if(ndm) { string test = generator()->preinitInterface(ndm, "Decayer", "set", decayer->fullName()); generator()->preinitInterface(ndm, "Active", "set", "Yes"); Energy width = decayer->partialWidth(inpart,outgoing); setBranchingRatio(ndm, width); } else throw NBodyDecayConstructorError() << "FourBodyDecayConstructor::createDecayMode - Needed to create " << "new decaymode but one could not be created for the tag " << tag << Exception::warning; } // otherwise else if (dm && (dm->decayer()->fullName()).find("Mambo") != string::npos) { // create the decayer GeneralFourBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac); if(!decayer) { if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for " << tag << " so mode not created\n"; return; } generator()->preinitInterface(dm, "Decayer", "set", decayer->fullName()); } //update CC mode if it exists if( inpart->CC() ) inpart->CC()->synchronize(); } GeneralFourBodyDecayerPtr FourBodyDecayConstructor::createDecayer(vector & diagrams, bool inter, double symfac) const { if(diagrams.empty()) return GeneralFourBodyDecayerPtr(); // extract the external particles for the process PDPtr incoming = diagrams[0].incoming; // outgoing particles vector outgoing(diagrams[0].outgoing.begin(), diagrams[0].outgoing.end()); // get the name for the object string objectname ("/Herwig/Decays/"); string classname = DecayerClassName(incoming, diagrams[0].outgoing, objectname); if(classname=="") return GeneralFourBodyDecayerPtr(); // create the object GeneralFourBodyDecayerPtr decayer = dynamic_ptr_cast (generator()->preinitCreate(classname, objectname)); // set up the decayer and return if doesn't work if(!decayer->setDecayInfo(incoming,outgoing,diagrams,symfac)) return GeneralFourBodyDecayerPtr(); // set decayer options from base class setDecayerInterfaces(objectname); // set the width option ostringstream value; value << widthOpt_; generator()->preinitInterface(objectname, "WidthOption", "set", value.str()); // set the intermediates option ostringstream value2; value2 << inter; generator()->preinitInterface(objectname, "GenerateIntermediates", "set", value2.str()); // initialize the decayer decayer->init(); // return the decayer return decayer; } string FourBodyDecayConstructor::DecayerClassName(tcPDPtr incoming, const OrderedParticles & outgoing, string & objname) const { string classname("Herwig::"); // spins of the outgoing particles unsigned int ns(0),nf(0),nv(0); objname += incoming->PDGName() + "2"; for(OrderedParticles::const_iterator it=outgoing.begin(); it!=outgoing.end();++it) { if ((**it).iSpin()==PDT::Spin0 ) ++ns; else if((**it).iSpin()==PDT::Spin1Half) ++nf; else if((**it).iSpin()==PDT::Spin1 ) ++nv; objname += (**it).PDGName(); } objname += "Decayer"; if(incoming->iSpin()==PDT::Spin0) { if(nf==4) classname += "StoFFFFDecayer"; else classname = ""; } else { classname=""; } return classname; } diff --git a/Models/General/FourBodyDecayConstructor.h b/Models/General/FourBodyDecayConstructor.h --- a/Models/General/FourBodyDecayConstructor.h +++ b/Models/General/FourBodyDecayConstructor.h @@ -1,150 +1,150 @@ // -*- C++ -*- #ifndef THEPEG_FourBodyDecayConstructor_H #define THEPEG_FourBodyDecayConstructor_H // // This is the declaration of the FourBodyDecayConstructor class. // #include "NBodyDecayConstructorBase.h" #include "ThePEG/Helicity/Vertex/VertexBase.h" #include "Herwig/Decay/General/GeneralFourBodyDecayer.fh" #include "PrototypeVertex.h" namespace Herwig { using namespace ThePEG; using Helicity::VertexBasePtr; /** * Here is the documentation of the FourBodyDecayConstructor class. * * @see \ref FourBodyDecayConstructorInterfaces "The interfaces" * defined for FourBodyDecayConstructor. */ class FourBodyDecayConstructor: public NBodyDecayConstructorBase { public: /** * The default constructor. */ FourBodyDecayConstructor() : interOpt_(0), widthOpt_(1), particleType_(false) {} /** * Destructor */ ~FourBodyDecayConstructor(); /** * Function used to determine allowed decaymodes, to be implemented * in derived class. * @param particles vector of ParticleData pointers containing * particles in model */ - virtual void DecayList(const set & particles); + virtual void DecayList(const set & particles); /** * Number of outgoing lines. Required for correct ordering. */ virtual unsigned int numBodies() const {return 4;} /** * Create a decay mode */ void createDecayMode(vector &,bool,double); /** * Create the decayer * @param diagrams The diagrams for the decay * @param inter Option for intermediates */ GeneralFourBodyDecayerPtr createDecayer(vector & diagrams, bool inter, double symfac) const; /** * Contruct the classname and object name for the Decayer * @param incoming The incoming particle * @param outgoing The decay products * @param objname a string containing the default path of the Decayer object */ string DecayerClassName(tcPDPtr incoming, const OrderedParticles & outgoing, string & objname) 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; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ FourBodyDecayConstructor & operator=(const FourBodyDecayConstructor &) = delete; private: /** * Option for the inclusion of intermediates */ unsigned int interOpt_; /** * How to treat the widths of the intermediate particles */ unsigned int widthOpt_; /** * Particles to override the default list */ vector particles_; /** * Types of particles */ bool particleType_; }; } #endif /* THEPEG_FourBodyDecayConstructor_H */ diff --git a/Models/General/ModelGenerator.cc b/Models/General/ModelGenerator.cc --- a/Models/General/ModelGenerator.cc +++ b/Models/General/ModelGenerator.cc @@ -1,556 +1,557 @@ // -*- C++ -*- // // ModelGenerator.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 ModelGenerator class. // #include "ModelGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Repository/CurrentGenerator.h" #include "BSMWidthGenerator.h" #include "Herwig/PDT/GenericMassGenerator.h" #include "Herwig/Decay/DecayIntegrator.h" #include "ThePEG/Repository/BaseRepository.h" using namespace Herwig; IBPtr ModelGenerator::clone() const { return new_ptr(*this); } IBPtr ModelGenerator::fullclone() const { return new_ptr(*this); } void ModelGenerator::persistentOutput(PersistentOStream & os) const { os << hardProcessConstructors_ << _theDecayConstructor << particles_ << offshell_ << Offsel_ << BRnorm_ << twoBodyOnly_ << howOffShell_ << Npoints_ << Iorder_ << BWshape_ << brMin_ << decayOutput_; } void ModelGenerator::persistentInput(PersistentIStream & is, int) { is >> hardProcessConstructors_ >> _theDecayConstructor >> particles_ >> offshell_ >> Offsel_ >> BRnorm_ >> twoBodyOnly_ >> howOffShell_ >> Npoints_ >> Iorder_ >> BWshape_ >> brMin_ >> decayOutput_; } // Static variable needed for the type description system in ThePEG. DescribeClass describeThePEGModelGenerator("Herwig::ModelGenerator", "Herwig.so"); void ModelGenerator::Init() { static ClassDocumentation documentation ("This class controls the the use of BSM physics.", "BSM physics was produced using the algorithm of " "\\cite{Gigg:2007cr,Gigg:2008yc}", "\\bibitem{Gigg:2007cr} M.~Gigg and P.~Richardson, \n" "Eur.\\ Phys.\\ J.\\ C {\\bf 51} (2007) 989.\n" "%%CITATION = EPHJA,C51,989;%%\n" " %\\cite{Gigg:2008yc}\n" "\\bibitem{Gigg:2008yc}\n" " M.~A.~Gigg and P.~Richardson,\n" " %``Simulation of Finite Width Effects in Physics Beyond the Standard Model,''\n" " arXiv:0805.3037 [hep-ph].\n" " %%CITATION = ARXIV:0805.3037;%%\n" ); static RefVector interfaceHardProcessConstructors ("HardProcessConstructors", "The objects to construct hard processes", &ModelGenerator::hardProcessConstructors_, -1, false, false, true, false, false); static Reference interfaceDecayConstructor ("DecayConstructor", "Pointer to DecayConstructor helper class", &ModelGenerator::_theDecayConstructor, false, false, true, false); static RefVector interfaceModelParticles ("DecayParticles", "ParticleData pointers to the particles requiring spin correlation " "decayers. If decay modes do not exist they will also be created.", &ModelGenerator::particles_, -1, false, false, true, false); static RefVector interfaceOffshell ("Offshell", "The particles to treat as off-shell", &ModelGenerator::offshell_, -1, false, false, true, false); static Switch interfaceWhichOffshell ("WhichOffshell", "A switch to determine which particles to create mass and width " "generators for.", &ModelGenerator::Offsel_, 0, false, false); static SwitchOption interfaceWhichOffshellSelected (interfaceWhichOffshell, "Selected", "Only create mass and width generators for the particles specified", 0); static SwitchOption interfaceWhichOffshellAll (interfaceWhichOffshell, "All", "Treat all particles specified in the DecayParticles " "list as off-shell", 1); static Switch interfaceBRNormalize ("BRNormalize", "Whether to normalize the partial widths to BR*total width for an " "on-shell particle", &ModelGenerator::BRnorm_, true, false, false); static SwitchOption interfaceBRNormalizeNormalize (interfaceBRNormalize, "Yes", "Normalize the partial widths", true); static SwitchOption interfaceBRNormalizeNoNormalize (interfaceBRNormalize, "No", "Do not normalize the partial widths", false); static Parameter interfacePoints ("InterpolationPoints", "Number of points to use for interpolation tables when needed", &ModelGenerator::Npoints_, 10, 5, 1000, false, false, true); static Parameter interfaceInterpolationOrder ("InterpolationOrder", "The interpolation order for the tables", &ModelGenerator::Iorder_, 1, 1, 5, false, false, Interface::limited); static Switch interfaceBreitWignerShape ("BreitWignerShape", "Controls the shape of the mass distribution generated", &ModelGenerator::BWshape_, 0, false, false); static SwitchOption interfaceBreitWignerShapeDefault (interfaceBreitWignerShape, "Default", "Running width with q in numerator and denominator width factor", 0); static SwitchOption interfaceBreitWignerShapeFixedWidth (interfaceBreitWignerShape, "FixedWidth", "Use a fixed width", 1); static SwitchOption interfaceBreitWignerShapeNoq (interfaceBreitWignerShape, "Noq", "Use M rather than q in the numerator and denominator width factor", 2); static SwitchOption interfaceBreitWignerShapeNoNumerator (interfaceBreitWignerShape, "NoNumerator", "Neglect the numerator factors", 3); static Switch interfaceTwoBodyOnly ("TwoBodyOnly", "Whether to use only two-body or all modes in the running width calculation", &ModelGenerator::twoBodyOnly_, false, false, false); static SwitchOption interfaceTwoBodyOnlyYes (interfaceTwoBodyOnly, "Yes", "Only use two-body modes", true); static SwitchOption interfaceTwoBodyOnlyNo (interfaceTwoBodyOnly, "No", "Use all modes", false); static Parameter interfaceMinimumBR ("MinimumBR", "The minimum branching fraction to include", &ModelGenerator::brMin_, 1e-6, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceDecayOutput ("DecayOutput", "Option to control the output of the decay mode information", &ModelGenerator::decayOutput_, 1, false, false); static SwitchOption interfaceDecayOutputNone (interfaceDecayOutput, "None", "No output", 0); static SwitchOption interfaceDecayOutputPlain (interfaceDecayOutput, "Plain", "Default plain text output", 1); static SwitchOption interfaceDecayOutputSLHA (interfaceDecayOutput, "SLHA", "Output in the Susy Les Houches Accord format", 2); static Parameter interfaceMinimumWidthFraction ("MinimumWidthFraction", "Minimum fraction of the particle's mass the width can be" " for the off-shell treatment.", &ModelGenerator::minWidth_, 1e-6, 1e-15, 1., false, false, Interface::limited); static Parameter interfaceHowMuchOffShell ("HowMuchOffShell", "The multiple of the particle's width by which it is allowed to be off-shell", &ModelGenerator::howOffShell_, 5., 0.0, 100., false, false, Interface::limited); } namespace { + /// Helper function for sorting by mass inline bool massIsLess(tcPDPtr a, tcPDPtr b) { return a->mass() < b->mass(); } // Helper function to find minimum possible mass of a particle inline Energy minimumMass(tcPDPtr parent) { Energy output(Constants::MaxEnergy); for(set::const_iterator dit = parent->decayModes().begin(); dit != parent->decayModes().end(); ++dit) { Energy outMass(ZERO); for(unsigned int ix=0;ix<(**dit).orderedProducts().size();++ix) { outMass += (**dit).orderedProducts()[ix]->massMin(); } output = min(output,outMass); } return output; } } void ModelGenerator::doinit() { useMe(); Interfaced::doinit(); // make sure the model is initialized Ptr::pointer model = dynamic_ptr_cast::pointer>(generator()->standardModel()); model->init(); // and the vertices for(size_t iv = 0; iv < model->numberOfVertices(); ++iv) model->vertex(iv)->init(); // uniq and sort DecayParticles list by mass set tmp(particles_.begin(),particles_.end()); particles_.assign(tmp.begin(),tmp.end()); sort(particles_.begin(),particles_.end(),massIsLess); //create decayers and decaymodes (if necessary) if( _theDecayConstructor ) { _theDecayConstructor->init(); _theDecayConstructor->createDecayers(particles_,brMin_); } // write out decays with spin correlations ostream & os = CurrentGenerator::current().misc(); ofstream ofs; if ( decayOutput_ > 1 ) { string filename = CurrentGenerator::current().filename() + "-BR.spc"; ofs.open(filename.c_str()); } if(decayOutput_!=0) { if(decayOutput_==1) { os << "# The decay modes listed below will have spin\n" << "# correlations included when they are generated.\n#\n#"; } else { ofs << "# Herwig decay tables in SUSY Les Houches accord format\n"; ofs << "Block DCINFO # Program information\n"; ofs << "1 Herwig # Decay Calculator\n"; ofs << "2 " << generator()->strategy()->versionstring() << " # Version number\n"; } } //create mass and width generators for the requested particles set offShell; if( Offsel_ == 0 ) offShell = set(offshell_.begin() ,offshell_.end() ); else offShell = set(particles_.begin(),particles_.end()); for(PDVector::iterator pit = particles_.begin(); pit != particles_.end(); ++pit) { tPDPtr parent = *pit; // Check decays for ones where quarks cannot be put on constituent // mass-shell checkDecays(parent); parent->reset(); // Now switch off the modes if needed for(DecaySet::const_iterator it=parent->decayModes().begin(); it!=parent->decayModes().end();++it) { if( _theDecayConstructor->disableDecayMode((**it).tag()) ) { DMPtr mode = *it; generator()->preinitInterface(mode, "Active", "set", "No"); if(mode->CC()) { DMPtr CCmode = mode->CC(); generator()->preinitInterface(CCmode, "Active", "set", "No"); } } } parent->update(); if( parent->CC() ) parent->CC()->synchronize(); if( parent->decaySelector().empty() ) { parent->stable(true); parent->width(ZERO); parent->widthCut(ZERO); parent->massGenerator(tGenericMassGeneratorPtr()); parent->widthGenerator(tGenericWidthGeneratorPtr()); } else { if(parent->mass()*minWidth_>parent->width()) { parent->massGenerator(tGenericMassGeneratorPtr()); parent->widthGenerator(tGenericWidthGeneratorPtr()); } else { if( offShell.find(*pit) != offShell.end() ) { createWidthGenerator(*pit); } else { parent->massGenerator(tGenericMassGeneratorPtr()); parent->widthGenerator(tGenericWidthGeneratorPtr()); } } } // set the offshellness Energy minMass = minimumMass(parent); Energy offShellNess = howOffShell_*parent->width(); if(minMass>parent->mass()-offShellNess) { offShellNess = parent->mass()-minMass; } parent->widthCut(offShellNess); if( parent->massGenerator() ) { parent->massGenerator()->reset(); if(decayOutput_==1) os << "# " <PDGName() << " will be considered off-shell.\n#\n"; } if( parent->widthGenerator() ) parent->widthGenerator()->reset(); } // loop again to initialise mass and width generators // switch off modes and write output for(PDVector::iterator pit = particles_.begin(); pit != particles_.end(); ++pit) { tPDPtr parent = *pit; if(parent->widthGenerator()) parent->widthGenerator()->init(); if(parent->massGenerator()) parent->massGenerator()->init(); // output the modes if needed if( !parent->decaySelector().empty() ) { if ( decayOutput_ == 2 ) writeDecayModes(ofs, parent); else writeDecayModes(os, parent); } } //Now construct hard processes given that we know which //objects have running widths for(unsigned int ix=0;ixinit(); hardProcessConstructors_[ix]->constructDiagrams(); } } void ModelGenerator::checkDecays(PDPtr parent) { if( parent->stable() ) { if(parent->coloured()) cerr << "Warning: No decays for coloured particle " << parent->PDGName() << "\n\n" << "have been calcluated in BSM model.\n" << "This may cause problems in the hadronization phase.\n" << "You may have forgotten to switch on the decay mode calculation using\n" << " set TwoBodyDC:CreateDecayModes Yes\n" << " set ThreeBodyDC:CreateDecayModes Yes\n" << " set WeakDecayConstructor:CreateDecayModes Yes\n" << "or the decays of this particle are missing from your\n" << "input spectrum and decay file in the SLHA format.\n\n"; return; } DecaySet::iterator dit = parent->decayModes().begin(); DecaySet::iterator dend = parent->decayModes().end(); Energy oldwidth(parent->width()), newwidth(ZERO); bool rescalebrat(false); double brsum(0.); for(; dit != dend; ++dit ) { if( !(**dit).on() ) continue; Energy release((**dit).parent()->mass()); tPDVector::const_iterator pit = (**dit).orderedProducts().begin(); tPDVector::const_iterator pend =(**dit).orderedProducts().end(); for( ; pit != pend; ++pit ) { release -= (**pit).constituentMass(); } if( (**dit).brat() < brMin_ || release < ZERO ) { if( release < ZERO ) cerr << "Warning: The shower cannot be generated using this decay " << (**dit).tag() << " because it is too close to threshold.\nIt " << "will be switched off and the branching fractions of the " << "remaining modes rescaled.\n"; rescalebrat = true; generator()->preinitInterface(*dit, "Active", "set", "No"); generator()->preinitInterface(*dit, "BranchingRatio", "set", "0.0"); DecayIntegratorPtr decayer = dynamic_ptr_cast((**dit).decayer()); if(decayer) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); } } else { brsum += (**dit).brat(); newwidth += (**dit).brat()*oldwidth; } } // if no modes left set stable if(newwidth==ZERO) { parent->stable(true); parent->width(ZERO); parent->widthCut(ZERO); parent->massGenerator(tGenericMassGeneratorPtr()); parent->widthGenerator(tGenericWidthGeneratorPtr()); } // otherwise rescale if needed else if( ( rescalebrat || abs(brsum - 1.) > 1e-12 ) && !parent->decayModes().empty()) { dit = parent->decayModes().begin(); dend = parent->decayModes().end(); double factor = oldwidth/newwidth; brsum = 0.; for( ; dit != dend; ++dit ) { if( !(**dit).on() ) continue; double newbrat = ((**dit).brat())*factor; brsum += newbrat; ostringstream brf; brf << setprecision(13) << newbrat; generator()->preinitInterface(*dit, "BranchingRatio", "set", brf.str()); } parent->width(newwidth); if( newwidth > ZERO ) parent->cTau(hbarc/newwidth); } } namespace { struct DecayModeOrdering { bool operator()(tcDMPtr m1, tcDMPtr m2) const { if(m1->brat()!=m2->brat()) { return m1->brat()>m2->brat(); } else { if(m1->products().size()==m2->products().size()) { ParticleMSet::const_iterator it1=m1->products().begin(); ParticleMSet::const_iterator it2=m2->products().begin(); do { if((**it1).id()!=(**it2).id()) { return (**it1).id()>(**it2).id(); } ++it1; ++it2; } while(it1!=m1->products().end()&& it2!=m2->products().end()); assert(false); } else return m1->products().size()products().size(); } return false; } }; } void ModelGenerator::writeDecayModes(ostream & os, tcPDPtr parent) const { if(decayOutput_==0) return; set modes(parent->decayModes().begin(), parent->decayModes().end()); if(decayOutput_==1) { os << " Parent: " << parent->PDGName() << " Mass (GeV): " << parent->mass()/GeV << " Total Width (GeV): " << parent->width()/GeV << endl; os << std::left << std::setw(40) << '#' << std::left << std::setw(20) << "Partial Width/GeV" << std::left << std::setw(20) << "BR" << "Yes/No\n"; for(set::iterator dit=modes.begin(); dit!=modes.end();++dit) os << std::left << std::setw(40) << (**dit).tag() << std::left << std::setw(20) << (**dit).brat()*parent->width()/GeV << std::left << std::setw(20) << (**dit).brat() << ((**dit).on() ? "Yes" : "No" ) << '\n'; os << "#\n#"; } else if(decayOutput_==2) { os << "# \t PDG \t Width\n"; os << "DECAY\t" << parent->id() << "\t" << parent->width()/GeV << "\t # " << parent->PDGName() << "\n"; for(set::iterator dit=modes.begin(); dit!=modes.end();++dit) { os << "\t" << std::left << std::setw(10) << (**dit).brat() << "\t" << (**dit).orderedProducts().size() << "\t"; for(unsigned int ix=0;ix<(**dit).orderedProducts().size();++ix) os << std::right << std::setw(10) << (**dit).orderedProducts()[ix]->id() ; for(unsigned int ix=(**dit).orderedProducts().size();ix<4;++ix) os << "\t"; os << "# " << (**dit).tag() << "\n"; } } } void ModelGenerator::createWidthGenerator(tPDPtr p) { string wn = p->fullName() + string("-WGen"); string mn = p->fullName() + string("-MGen"); GenericMassGeneratorPtr mgen = dynamic_ptr_cast (generator()->preinitCreate("Herwig::GenericMassGenerator", mn)); BSMWidthGeneratorPtr wgen = dynamic_ptr_cast (generator()->preinitCreate("Herwig::BSMWidthGenerator", wn)); //set the particle interface mgen->particle(p); wgen->particle(p); //set the generator interfaces in the ParticleData object generator()->preinitInterface(p, "Mass_generator","set", mn); generator()->preinitInterface(p, "Width_generator","set", wn); //allow the branching fraction of this particle type to vary p->variableRatio(true); if( p->CC() ) p->CC()->variableRatio(true); //initialize the generators generator()->preinitInterface(mgen, "Initialize", "set", "Yes"); generator()->preinitInterface(wgen, "Initialize", "set", "Yes"); string norm = BRnorm_ ? "Yes" : "No"; generator()->preinitInterface(wgen, "BRNormalize", "set", norm); string twob = twoBodyOnly_ ? "Yes" : "No"; generator()->preinitInterface(wgen, "TwoBodyOnly", "set", twob); ostringstream os; os << Npoints_; generator()->preinitInterface(wgen, "Points", "set", os.str()); os.str(""); os << Iorder_; generator()->preinitInterface(wgen, "InterpolationOrder", "set", os.str()); os.str(""); os << BWshape_; generator()->preinitInterface(mgen, "BreitWignerShape", "set", os.str()); } diff --git a/Models/General/NBodyDecayConstructorBase.cc b/Models/General/NBodyDecayConstructorBase.cc --- a/Models/General/NBodyDecayConstructorBase.cc +++ b/Models/General/NBodyDecayConstructorBase.cc @@ -1,682 +1,684 @@ // -*- C++ -*- // // NBodyDecayConstructorBase.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 NBodyDecayConstructorBase class. // #include "NBodyDecayConstructorBase.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Repository/Repository.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "DecayConstructor.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; using namespace ThePEG; void NBodyDecayConstructorBase::persistentOutput(PersistentOStream & os ) const { os << init_ << iteration_ << points_ << info_ << decayConstructor_ << createModes_ << minReleaseFraction_ << maxBoson_ << maxList_ << removeOnShell_ << excludeEffective_ << includeTopOnShell_ << excludedVerticesVector_ << excludedVerticesSet_ << excludedParticlesVector_ << excludedParticlesSet_ << removeFlavourChangingVertices_ << removeSmallVertices_ << minVertexNorm_; } void NBodyDecayConstructorBase::persistentInput(PersistentIStream & is , int) { is >> init_ >> iteration_ >> points_ >> info_ >> decayConstructor_ >> createModes_ >> minReleaseFraction_ >> maxBoson_ >> maxList_ >> removeOnShell_ >> excludeEffective_ >> includeTopOnShell_ >> excludedVerticesVector_ >> excludedVerticesSet_ >> excludedParticlesVector_ >> excludedParticlesSet_ >> removeFlavourChangingVertices_ >> removeSmallVertices_ >> minVertexNorm_; } // Static variable needed for the type description system in ThePEG. DescribeAbstractClass describeThePEGNBodyDecayConstructorBase("Herwig::NBodyDecayConstructorBase", "Herwig.so"); // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractClass describeHerwigNBodyDecayConstructorBase("Herwig::NBodyDecayConstructorBase", "Herwig.so"); void NBodyDecayConstructorBase::Init() { static ClassDocumentation documentation ("The NBodyDecayConstructorBase class is the base class for the automatic" "construction of the decay modes"); static Switch interfaceInitializeDecayers ("InitializeDecayers", "Initialize new decayers", &NBodyDecayConstructorBase::init_, true, false, false); static SwitchOption interfaceInitializeDecayersInitializeDecayersYes (interfaceInitializeDecayers, "Yes", "Initialize new decayers to find max weights", true); static SwitchOption interfaceInitializeDecayersNo (interfaceInitializeDecayers, "No", "Use supplied weights for integration", false); static Parameter interfaceInitIteration ("InitIteration", "Number of iterations to optimise integration weights", &NBodyDecayConstructorBase::iteration_, 1, 0, 10, false, false, true); static Parameter interfaceInitPoints ("InitPoints", "Number of points to generate when optimising integration", &NBodyDecayConstructorBase::points_, 1000, 100, 100000000, false, false, true); static Switch interfaceOutputInfo ("OutputInfo", "Whether to output information about the decayers", &NBodyDecayConstructorBase::info_, false, false, false); static SwitchOption interfaceOutputInfoNo (interfaceOutputInfo, "No", "Do not output information regarding the created decayers", false); static SwitchOption interfaceOutputInfoYes (interfaceOutputInfo, "Yes", "Output information regarding the decayers", true); static Switch interfaceCreateDecayModes ("CreateDecayModes", "Whether to create the ThePEG::DecayMode objects as well as the decayers", &NBodyDecayConstructorBase::createModes_, true, false, false); static SwitchOption interfaceCreateDecayModesYes (interfaceCreateDecayModes, "Yes", "Create the ThePEG::DecayMode objects", true); static SwitchOption interfaceCreateDecayModesNo (interfaceCreateDecayModes, "No", "Only create the Decayer objects", false); static Switch interfaceRemoveOnShell ("RemoveOnShell", "Remove on-shell diagrams as should be treated as a sequence of 1->2 decays", &NBodyDecayConstructorBase::removeOnShell_, 1, false, false); static SwitchOption interfaceRemoveOnShellYes (interfaceRemoveOnShell, "Yes", "Remove the diagrams if neither the production of decay or the intermediate" " can happen", 1); static SwitchOption interfaceRemoveOnShellNo (interfaceRemoveOnShell, "No", "Never remove the intermediate", 0); static SwitchOption interfaceRemoveOnShellProduction (interfaceRemoveOnShell, "Production", "Remove the diagram if the on-shell production of the intermediate is allowed", 2); static RefVector interfaceExcludedVertices ("ExcludedVertices", "Vertices which are not included in the three-body decayers", &NBodyDecayConstructorBase::excludedVerticesVector_, -1, false, false, true, true, false); static RefVector interfaceExcludedIntermediates ("ExcludedIntermediates", "Excluded intermediate particles", &NBodyDecayConstructorBase::excludedParticlesVector_, -1, false, false, true, true, false); static Switch interfaceExcludeEffectiveVertices ("ExcludeEffectiveVertices", "Exclude effectice vertices", &NBodyDecayConstructorBase::excludeEffective_, true, false, false); static SwitchOption interfaceExcludeEffectiveVerticesYes (interfaceExcludeEffectiveVertices, "Yes", "Exclude the effective vertices", true); static SwitchOption interfaceExcludeEffectiveVerticesNo (interfaceExcludeEffectiveVertices, "No", "Don't exclude the effective vertices", false); static Parameter interfaceMinReleaseFraction ("MinReleaseFraction", "The minimum energy release for a three-body decay, as a " "fraction of the parent mass.", &NBodyDecayConstructorBase::minReleaseFraction_, 1e-3, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceMaximumGaugeBosons ("MaximumGaugeBosons", "Maximum number of electroweak gauge bosons" " to be produced as decay products", &NBodyDecayConstructorBase::maxBoson_, 1, false, false); static SwitchOption interfaceMaximumGaugeBosonsNone (interfaceMaximumGaugeBosons, "None", "Produce no W/Zs", 0); static SwitchOption interfaceMaximumGaugeBosonsSingle (interfaceMaximumGaugeBosons, "Single", "Produce at most one W/Zs", 1); static SwitchOption interfaceMaximumGaugeBosonsDouble (interfaceMaximumGaugeBosons, "Double", "Produce at most two W/Zs", 2); static SwitchOption interfaceMaximumGaugeBosonsTriple (interfaceMaximumGaugeBosons, "Triple", "Produce at most three W/Zs", 3); static Switch interfaceMaximumNewParticles ("MaximumNewParticles", "Maximum number of particles from the list of " "decaying particles to be allowed as decay products", &NBodyDecayConstructorBase::maxList_, 0, false, false); static SwitchOption interfaceMaximumNewParticlesNone (interfaceMaximumNewParticles, "None", "No particles from the list", 0); static SwitchOption interfaceMaximumNewParticlesOne (interfaceMaximumNewParticles, "One", "A single particle from the list", 1); static SwitchOption interfaceMaximumNewParticlesTwo (interfaceMaximumNewParticles, "Two", "Two particles from the list", 2); static SwitchOption interfaceMaximumNewParticlesThree (interfaceMaximumNewParticles, "Three", "Three particles from the list", 3); static SwitchOption interfaceMaximumNewParticlesFour (interfaceMaximumNewParticles, "Four", "Four particles from the list", 4); static Switch interfaceIncludeOnShellTop ("IncludeOnShellTop", "Include the on-shell diagrams involving t -> bW", &NBodyDecayConstructorBase::includeTopOnShell_, false, false, false); static SwitchOption interfaceIncludeOnShellTopYes (interfaceIncludeOnShellTop, "Yes", "Inlude them", true); static SwitchOption interfaceIncludeOnShellTopNo (interfaceIncludeOnShellTop, "No", "Don't include them", true); static Switch interfaceRemoveSmallVertices ("RemoveSmallVertices", "Remove vertices with norm() below minVertexNorm", &NBodyDecayConstructorBase::removeSmallVertices_, false, false, false); static SwitchOption interfaceRemoveSmallVerticesYes (interfaceRemoveSmallVertices, "Yes", "Remove them", true); static SwitchOption interfaceRemoveSmallVerticesNo (interfaceRemoveSmallVertices, "No", "Don't remove them", false); static Parameter interfaceMinVertexNorm ("MinVertexNorm", "Minimum allowed value of the notm() of the vertex if removing small vertices", &NBodyDecayConstructorBase::minVertexNorm_, 1e-8, 1e-300, 1., false, false, Interface::limited); static Switch interfaceRemoveFlavourChangingVertices ("RemoveFlavourChangingVertices", "Remove flavour changing interactions with the photon and gluon", &NBodyDecayConstructorBase::removeFlavourChangingVertices_, false, false, false); static SwitchOption interfaceRemoveFlavourChangingVerticesYes (interfaceRemoveFlavourChangingVertices, "Yes", "Remove them", true); static SwitchOption interfaceRemoveFlavourChangingVerticesNo (interfaceRemoveFlavourChangingVertices, "No", "Don't remove them", false); } void NBodyDecayConstructorBase::setBranchingRatio(tDMPtr dm, Energy pwidth) { // if zero width just set BR to zero if(pwidth==ZERO) { generator()->preinitInterface(dm, "BranchingRatio","set", "0."); generator()->preinitInterface(dm, "OnOff","set", "Off"); return; } // Need width and branching ratios for all currently created decay modes PDPtr parent = const_ptr_cast(dm->parent()); DecaySet modes = parent->decayModes(); unsigned int nmodes=0; for( auto dm : modes ) { if(dm->on()) ++nmodes; } if( nmodes==0 ) return; double dmbrat(0.); if( nmodes == 1 ) { parent->width(pwidth); if( pwidth > ZERO ) parent->cTau(hbarc/pwidth); dmbrat = 1.; } else { Energy currentwidth(parent->width()); Energy newWidth(currentwidth + pwidth); parent->width(newWidth); if( newWidth > ZERO ) parent->cTau(hbarc/newWidth); //need to reweight current branching fractions if there are any double factor = newWidth > ZERO ? double(currentwidth/newWidth) : 0.; for(DecaySet::const_iterator dit = modes.begin(); dit != modes.end(); ++dit) { if( **dit == *dm || !(**dit).on() ) continue; double newbrat = (**dit).brat()*factor; ostringstream brf; brf << setprecision(13)<< newbrat; generator()->preinitInterface(*dit, "BranchingRatio", "set", brf.str()); } //set brat for current mode dmbrat = newWidth > ZERO ? double(pwidth/newWidth) : 0.; } ostringstream br; br << setprecision(13) << dmbrat; generator()->preinitInterface(dm, "BranchingRatio", "set", br.str()); } void NBodyDecayConstructorBase::setDecayerInterfaces(string fullname) const { if( initialize() ) { ostringstream value; value << initialize(); generator()->preinitInterface(fullname, "Initialize", "set", value.str()); value.str(""); value << iteration(); generator()->preinitInterface(fullname, "Iteration", "set", value.str()); value.str(""); value << points(); generator()->preinitInterface(fullname, "Points", "set", value.str()); } // QED stuff if needed if(decayConstructor()->QEDGenerator()) generator()->preinitInterface(fullname, "PhotonGenerator", "set", decayConstructor()->QEDGenerator()->fullName()); string outputmodes; if( info() ) outputmodes = string("Output"); else outputmodes = string("NoOutput"); generator()->preinitInterface(fullname, "OutputModes", "set", outputmodes); } void NBodyDecayConstructorBase::doinit() { Interfaced::doinit(); excludedVerticesSet_ = set(excludedVerticesVector_.begin(), excludedVerticesVector_.end()); excludedParticlesSet_ = set(excludedParticlesVector_.begin(), excludedParticlesVector_.end()); if(removeOnShell_==0&&numBodies()>2) generator()->log() << "Warning: Including diagrams with on-shell " << "intermediates in " << numBodies() << "-body BSM decays, this" << " can lead to double counting and is not" << " recommended unless you really know what you are doing\n" << "This can be switched off using\n set " << fullName() << ":RemoveOnShell Yes\n"; } namespace { void constructIdenticalSwaps(unsigned int depth, vector > identical, vector channelType, list > & swaps) { if(depth==0) { unsigned int size = identical.size(); for(unsigned ix=0;ix newType=channelType; swap(newType[identical[depth][ix]],newType[identical[depth][iy]]); // at bottom of chain if(depth+1==identical.size()) { swaps.push_back(newType); } else { constructIdenticalSwaps(depth+1,identical,newType,swaps); } } } } void identicalFromSameDecay(unsigned int & loc, const NBVertex & vertex, vector > & sameDecay) { list >::const_iterator it = vertex.vertices.begin(); while(it!=vertex.vertices.end()) { if(it->second.incoming) { identicalFromSameDecay(loc,it->second,sameDecay); ++it; continue; } ++loc; long id = it->first->id(); ++it; if(it == vertex.vertices.end()) break; if(it->second.incoming) continue; if(it->first->id()!=id) continue; sameDecay.push_back(vector()); sameDecay.back().push_back(loc-1); while(it != vertex.vertices.end() && !it->second.incoming && it->first->id()==id) { ++loc; ++it; sameDecay.back().push_back(loc-1); } }; } } -void NBodyDecayConstructorBase::DecayList(const set & particles) { +void NBodyDecayConstructorBase::DecayList(const set & particles) { if( particles.empty() ) return; // cast the StandardModel to the Hw++ one to get the vertices tHwSMPtr model = dynamic_ptr_cast(generator()->standardModel()); unsigned int nv(model->numberOfVertices()); // loop over the particles and create the modes - for(set::const_iterator ip=particles.begin(); + for(set::const_iterator ip=particles.begin(); ip!=particles.end();++ip) { // get the decaying particle tPDPtr parent = *ip; if ( Debug::level > 0 ) Repository::cout() << "Constructing N-body decays for " << parent->PDGName() << '\n'; // first create prototype 1->2 decays std::stack prototypes; for(unsigned int iv = 0; iv < nv; ++iv) { VertexBasePtr vertex = model->vertex(iv); if(excluded(vertex)) continue; PrototypeVertex::createPrototypes(parent, vertex, prototypes,this); } // now expand them into potential decay modes list > modes; while(!prototypes.empty()) { // get the first prototype from the stack PrototypeVertexPtr proto = prototypes.top(); prototypes.pop(); // multiplcity too low if(proto->npartvertex(iv); if(excluded(vertex) || proto->npart+vertex->getNpoint()>numBodies()+2) continue; PrototypeVertex::expandPrototypes(proto,vertex,prototypes, excludedParticlesSet_,this); } } // multiplcity too high disgard else if(proto->npart>numBodies()) { continue; } // right multiplicity else { // check it's kinematical allowed for physical masses if( proto->incomingMass() < proto->outgoingMass() ) continue; // and for constituent masses Energy outgoingMass = proto->outgoingConstituentMass(); if( proto->incomingMass() < proto->outgoingConstituentMass() ) continue; // remove processes which aren't kinematically allowed within // the release fraction if( proto->incomingMass() - outgoingMass <= minReleaseFraction_ * proto->incomingMass() ) continue; // check the external particles if(!proto->checkExternal()) continue; // check if first piece on-shell bool onShell = proto->canBeOnShell(removeOnShell(),proto->incomingMass(),true); // special treatment for three-body top decays if(onShell) { if(includeTopOnShell_ && abs(proto->incoming->id())==ParticleID::t && proto->npart==3) { unsigned int nprimary=0; bool foundW=false, foundQ=false; for(OrderedVertices::const_iterator it = proto->outgoing.begin(); it!=proto->outgoing.end();++it) { if(abs(it->first->id())==ParticleID::Wplus) foundW = true; if(abs(it->first->id())==ParticleID::b || abs(it->first->id())==ParticleID::s || abs(it->first->id())==ParticleID::d) foundQ = true; ++nprimary; } if(!(foundW&&foundQ&&nprimary==2)) continue; } else continue; } // check if should be added to an existing decaymode bool added = false; for(list >::iterator it=modes.begin(); it!=modes.end();++it) { // is the same ? if(!(*it)[0]->sameDecay(*proto)) continue; // it is the same added = true; // check if it is a duplicate bool already = false; for(unsigned int iz = 0; iz < (*it).size(); ++iz) { if( *(*it)[iz] == *proto) { already = true; break; } } if(!already) (*it).push_back(proto); break; } if(!added) modes.push_back(vector(1,proto)); } } // now look at the decay modes for(list >::iterator mit = modes.begin(); mit!=modes.end();++mit) { // count the number of gauge bosons and particles from the list unsigned int nlist(0),nbos(0); for(OrderedParticles::const_iterator it=(*mit)[0]->outPart.begin(); it!=(*mit)[0]->outPart.end();++it) { if(abs((**it).id()) == ParticleID::Wplus || abs((**it).id()) == ParticleID::Z0) ++nbos; if(particles.find(*it)!=particles.end()) ++nlist; if((**it).CC() && particles.find((**it).CC())!=particles.end()) ++nlist; } // if too many ignore the mode if(nbos > maxBoson_ || nlist > maxList_) continue; // translate the prototypes into diagrams vector newDiagrams; double symfac(1.); bool possibleOnShell=false; for(unsigned int ix=0;ix<(*mit).size();++ix) { symfac = 1.; possibleOnShell |= (*mit)[ix]->possibleOnShell; NBDiagram templateDiagram = NBDiagram((*mit)[ix]); // extract the ordering vector order(templateDiagram.channelType.size()); for(unsigned int iz=0;iz > identical; OrderedParticles::const_iterator it=templateDiagram.outgoing.begin(); unsigned int iloc=0; while(it!=templateDiagram.outgoing.end()) { OrderedParticles::const_iterator lt = templateDiagram.outgoing.lower_bound(*it); OrderedParticles::const_iterator ut = templateDiagram.outgoing.upper_bound(*it); unsigned int nx=0; for(OrderedParticles::const_iterator jt=lt;jt!=ut;++jt) {++nx;} if(nx==1) { ++it; ++iloc; } else { symfac *= factorial(nx); identical.push_back(vector()); for(OrderedParticles::const_iterator jt=lt;jt!=ut;++jt) { identical.back().push_back(order[iloc]); ++iloc; } it = ut; } } // that's it if there outgoing are unqiue if(identical.empty()) { newDiagrams.push_back(templateDiagram); continue; } // find any identical particles which shouldn't be swapped unsigned int loc=0; vector > sameDecay; identicalFromSameDecay(loc,templateDiagram,sameDecay); // compute the swaps list > swaps; constructIdenticalSwaps(0,identical,templateDiagram.channelType,swaps); // special if identical from same decay if(!sameDecay.empty()) { for(vector >::const_iterator st=sameDecay.begin(); st!=sameDecay.end();++st) { for(list >::iterator it=swaps.begin(); it!=swaps.end();++it) { for(unsigned int iy=0;iy<(*st).size();++iy) { for(unsigned int iz=iy+1;iz<(*st).size();++iz) { if((*it)[(*st)[iy]]>(*it)[(*st)[iz]]) swap((*it)[(*st)[iy]],(*it)[(*st)[iz]]); } } } } } // remove any dupliciates for(list >::iterator it=swaps.begin(); it!=swaps.end();++it) { for(list >::iterator jt=it; jt!=swaps.end();++jt) { if(it==jt) continue; bool different=false; for(unsigned int iz=0;iz<(*it).size();++iz) { if((*it)[iz]!=(*jt)[iz]) { different = true; break; } } if(!different) { jt = swaps.erase(jt); --jt; } } } // special for identical decay products if(templateDiagram.vertices.begin()->second.incoming) { if( templateDiagram.vertices.begin() ->first == (++templateDiagram.vertices.begin())->first) { if(*( (*mit)[ix]->outgoing.begin() ->second) == *((++(*mit)[ix]->outgoing.begin())->second)) { for(list >::iterator it=swaps.begin(); it!=swaps.end();++it) { for(list >::iterator jt=it; jt!=swaps.end();++jt) { if(it==jt) continue; if((*it)[0]==(*jt)[2]&&(*it)[1]==(*jt)[3]) { jt = swaps.erase(jt); --jt; } } } } } } for(list >::iterator it=swaps.begin(); it!=swaps.end();++it) { newDiagrams.push_back(templateDiagram); newDiagrams.back().channelType = *it; } } // create the decay if( Debug::level > 1 ) { generator()->log() << "Mode: "; generator()->log() << (*mit)[0]->incoming->PDGName() << " -> "; for(OrderedParticles::const_iterator it=(*mit)[0]->outPart.begin(); it!=(*mit)[0]->outPart.end();++it) generator()->log() << (**it).PDGName() << " "; generator()->log() << "\n"; generator()->log() << "There are " << (*mit).size() << " diagrams\n"; for(unsigned int iy=0;iylog() << "Diagram: " << iy << "\n"; generator()->log() << newDiagrams[iy] << "\n"; } } createDecayMode(newDiagrams,possibleOnShell,symfac); } } } void NBodyDecayConstructorBase::createDecayMode(vector &, bool, double) { + assert(false); throw Exception() << "In NBodyDecayConstructorBase::createDecayMode() which" - << " should have be overridden in an inheriting class" + << " should have be overridden in an inheriting class " + << fullName() << Exception::abortnow; } diff --git a/Models/General/NBodyDecayConstructorBase.h b/Models/General/NBodyDecayConstructorBase.h --- a/Models/General/NBodyDecayConstructorBase.h +++ b/Models/General/NBodyDecayConstructorBase.h @@ -1,340 +1,348 @@ // -*- C++ -*- // // NBodyDecayConstructorBase.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_NBodyDecayConstructorBase_H #define HERWIG_NBodyDecayConstructorBase_H // // This is the declaration of the NBodyDecayConstructorBase class. // #include "ThePEG/Interface/Interfaced.h" #include "ThePEG/Utilities/Exception.h" #include "ThePEG/PDT/ParticleData.h" #include "NBodyDecayConstructorBase.fh" #include "PrototypeVertex.h" #include "DecayConstructor.fh" namespace Herwig { using namespace ThePEG; /** * This is the base class for NBodyDecayConstructors. An N-body * decay constructor should inherit from this and implement the * DecayList virtual funtcion to create the decays and decayers. * * @see \ref NBodyDecayConstructorBaseInterfaces "The interfaces" * defined for NBodyDecayConstructor. */ class NBodyDecayConstructorBase: public Interfaced { public: + struct MassOrdering { + bool operator()(PDPtr p1, PDPtr p2) const { + return p1->mass() < p2->mass() || (p1->mass()==p2->mass() && p1->id()>p2->id()); + } + }; + +public: + /** * The default constructor. */ NBodyDecayConstructorBase() : init_(true),iteration_(1), points_(1000), info_(false), createModes_(true), removeOnShell_(1), excludeEffective_(true), minReleaseFraction_(1e-3), maxBoson_(1), maxList_(1), includeTopOnShell_(false ), removeFlavourChangingVertices_(false), removeSmallVertices_(false), minVertexNorm_(1e-8) {} /** * Function used to determine allowed decaymodes, to be implemented * in derived class. * @param particles vector of ParticleData pointers containing * particles in model */ - virtual void DecayList(const set & particles); + virtual void DecayList(const set & particles); /** * Number of outgoing lines. Required for correct ordering. */ virtual unsigned int numBodies() const = 0; /** * Set the pointer to the DecayConstrcutor */ void decayConstructor(tDecayConstructorPtr d) { decayConstructor_ = d; } /** * Remove flavour changing vertices ? */ bool removeFlavourChangingVertices() const { return removeFlavourChangingVertices_; } /** * Remove small vertices ? */ bool removeSmallVertices() const { return removeSmallVertices_; } /** * Minimum norm for vertex removal */ double minVertexNorm() const { return minVertexNorm_; } protected: /** * Method to set up the decay mode, should be overidden in inheriting class */ virtual void createDecayMode(vector & mode, bool possibleOnShell, double symfac); /** * Set the branching ratio of this mode. This requires * calculating a new width for the decaying particle and reweighting * the current branching fractions. * @param dm The decaymode for which to set the branching ratio * @param pwidth The calculated width of the mode */ void setBranchingRatio(tDMPtr dm, Energy pwidth); /** * Set the interfaces of the decayers depending on the flags stored. * @param name Fullname of the decayer in the EventGenerator * including the path */ void setDecayerInterfaces(string name) const; /** * Whether to initialize decayers or not */ bool initialize() const { return init_; } /** * Number of iterations if initializing (default 1) */ int iteration() const { return iteration_; } /** * Number of points to do in initialization */ int points() const { return points_; } /** * Whether to output information on the decayers */ bool info() const { return info_; } /** * Whether to create the DecayModes as well as the Decayer objects */ bool createDecayModes() const { return createModes_; } /** * Maximum number of electroweak gauge bosons */ unsigned int maximumGaugeBosons() const { return maxBoson_;} /** * Maximum number of particles from the list whose decays we are calculating */ unsigned int maximumList() const { return maxList_;} /** * Minimum energy release fraction */ double minimumReleaseFraction() const {return minReleaseFraction_;} /** * Get the pointer to the DecayConstructor object */ tDecayConstructorPtr decayConstructor() const { return decayConstructor_; } /** * Option for on-shell particles */ unsigned int removeOnShell() const { return removeOnShell_; } /** * Check if a vertex is excluded */ bool excluded(VertexBasePtr vertex) const { // skip an effective vertex if( excludeEffective_ && vertex->orderInAllCouplings() != int(vertex->getNpoint())-2) return true; // check if explicitly forbidden return excludedVerticesSet_.find(vertex)!=excludedVerticesSet_.end(); } 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 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 assignment operator is private and must never be called. * In fact, it should not even be implemented. */ NBodyDecayConstructorBase & operator=(const NBodyDecayConstructorBase &) = delete; private: /** * Whether to initialize decayers or not */ bool init_; /** * Number of iterations if initializing (default 1) */ int iteration_; /** * Number of points to do in initialization */ int points_; /** * Whether to output information on the decayers */ bool info_; /** * Whether to create the DecayModes as well as the Decayer objects */ bool createModes_; /** * Whether or not to remove on-shell diagrams */ unsigned int removeOnShell_; /** * Excluded Vertices */ vector excludedVerticesVector_; /** * Excluded Vertices */ set excludedVerticesSet_; /** * Excluded Particles */ vector excludedParticlesVector_; /** * Excluded Particles */ set excludedParticlesSet_; /** * Whether or not to exclude effective vertices */ bool excludeEffective_; /** * A pointer to the DecayConstructor object */ tDecayConstructorPtr decayConstructor_; /** * The minimum energy release for a three-body decay as a * fraction of the parent mass */ double minReleaseFraction_; /** * Maximum number of EW gauge bosons */ unsigned int maxBoson_; /** * Maximum number of particles from the decaying particle list */ unsigned int maxList_; /** * Include on-shell for \f$t\to b W\f$ */ bool includeTopOnShell_; /** * Remove flavour changing vertices ? */ bool removeFlavourChangingVertices_; /** * Remove small vertices ? */ bool removeSmallVertices_; /** * Minimum norm for vertex removal */ double minVertexNorm_; }; /** An Exception class that can be used by all inheriting classes to * indicate a setup problem. */ class NBodyDecayConstructorError : public Exception { public: NBodyDecayConstructorError() : Exception() {} NBodyDecayConstructorError(const string & str, Severity sev) : Exception(str,sev) {} }; } #endif /* HERWIG_NBodyDecayConstructorBase_H */ diff --git a/Models/General/ThreeBodyDecayConstructor.cc b/Models/General/ThreeBodyDecayConstructor.cc --- a/Models/General/ThreeBodyDecayConstructor.cc +++ b/Models/General/ThreeBodyDecayConstructor.cc @@ -1,377 +1,377 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the ThreeBodyDecayConstructor class. // #include "ThreeBodyDecayConstructor.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "Herwig/Decay/General/GeneralThreeBodyDecayer.h" #include "Herwig/PDT/ThreeBodyAllOnCalculator.h" #include "ThePEG/PDT/StandardMatchers.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Utilities/Throw.h" #include "DecayConstructor.h" #include "WeakCurrentDecayConstructor.h" using namespace Herwig; IBPtr ThreeBodyDecayConstructor::clone() const { return new_ptr(*this); } IBPtr ThreeBodyDecayConstructor::fullclone() const { return new_ptr(*this); } void ThreeBodyDecayConstructor::persistentOutput(PersistentOStream & os) const { os << interOpt_ << widthOpt_ << intOpt_ << relErr_ << includeIntermediatePhotons_; } void ThreeBodyDecayConstructor::persistentInput(PersistentIStream & is, int) { is >> interOpt_ >> widthOpt_ >> intOpt_ >> relErr_ >> includeIntermediatePhotons_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigThreeBodyDecayConstructor("Herwig::ThreeBodyDecayConstructor", "Herwig.so"); void ThreeBodyDecayConstructor::Init() { static ClassDocumentation documentation ("The ThreeBodyDecayConstructor class constructs the three body decay modes"); static Switch interfaceIncludeIntermediatePhotons ("IncludeIntermediatePhotons", "Whether or not ot allow intermediate photons", &ThreeBodyDecayConstructor::includeIntermediatePhotons_, false, false, false); static SwitchOption interfaceIncludeIntermediatePhotonsYes (interfaceIncludeIntermediatePhotons, "Yes", "Include them", true); static SwitchOption interfaceIncludeIntermediatePhotonsNo (interfaceIncludeIntermediatePhotons, "No", "Don't include them", false); static Switch interfaceWidthOption ("WidthOption", "Option for the treatment of the widths of the intermediates", &ThreeBodyDecayConstructor::widthOpt_, 1, false, false); static SwitchOption interfaceWidthOptionFixed (interfaceWidthOption, "Fixed", "Use fixed widths", 1); static SwitchOption interfaceWidthOptionRunning (interfaceWidthOption, "Running", "Use running widths", 2); static SwitchOption interfaceWidthOptionZero (interfaceWidthOption, "Zero", "Set the widths to zero", 3); static Switch interfaceIntermediateOption ("IntermediateOption", "Option for the inclusion of intermediates in the event", &ThreeBodyDecayConstructor::interOpt_, 0, false, false); static SwitchOption interfaceIntermediateOptionAlways (interfaceIntermediateOption, "Always", "Always include the intermediates", 1); static SwitchOption interfaceIntermediateOptionNever (interfaceIntermediateOption, "Never", "Never include the intermediates", 2); static SwitchOption interfaceIntermediateOptionOnlyIfOnShell (interfaceIntermediateOption, "OnlyIfOnShell", "Only if there are on-shell diagrams", 0); static Switch interfaceIntegrationOption ("IntegrationOption", "Option for the integration of the partial width", &ThreeBodyDecayConstructor::intOpt_, 1, false, false); static SwitchOption interfaceIntegrationOptionAllPoles (interfaceIntegrationOption, "AllPoles", "Include all potential poles", 0); static SwitchOption interfaceIntegrationOptionShallowestPole (interfaceIntegrationOption, "ShallowestPole", "Only include the shallowest pole", 1); static Parameter interfaceRelativeError ("RelativeError", "The relative error for the GQ integration", &ThreeBodyDecayConstructor::relErr_, 1e-2, 1e-10, 1., false, false, Interface::limited); } -void ThreeBodyDecayConstructor::DecayList(const set & particles) { +void ThreeBodyDecayConstructor::DecayList(const set & particles) { if( particles.empty() ) return; // special for weak decays for(unsigned int ix=0;ixdecayConstructors().size();++ix) { Ptr::pointer weak = dynamic_ptr_cast::pointer > (decayConstructor()->decayConstructors()[ix]); if(!weak) continue; weakMassCut_ = max(weakMassCut_,weak->massCut()); } NBodyDecayConstructorBase::DecayList(particles); } GeneralThreeBodyDecayerPtr ThreeBodyDecayConstructor:: createDecayer(vector & diagrams, bool inter, double symfac) const { if(diagrams.empty()) return GeneralThreeBodyDecayerPtr(); // extract the external particles for the process PDPtr incoming = getParticleData(diagrams[0].incoming); // outgoing particles OrderedParticles outgoing; outgoing.insert(getParticleData(diagrams[0].outgoing )); outgoing.insert(getParticleData(diagrams[0].outgoingPair.first )); outgoing.insert(getParticleData(diagrams[0].outgoingPair.second)); // get the object name string objectname ("/Herwig/Decays/"); string classname = DecayerClassName(incoming, outgoing, objectname); if(classname=="") return GeneralThreeBodyDecayerPtr(); // create the object GeneralThreeBodyDecayerPtr decayer = dynamic_ptr_cast (generator()->preinitCreate(classname, objectname)); // set up the decayer and return if doesn't work vector outVector(outgoing.begin(),outgoing.end()); if(!decayer->setDecayInfo(incoming,outVector,diagrams,symfac)) return GeneralThreeBodyDecayerPtr(); // set decayer options from base class setDecayerInterfaces(objectname); // options for partial width integration ostringstream value; value << intOpt_; generator()->preinitInterface(objectname, "PartialWidthIntegration", "set", value.str()); value.str(""); value << relErr_; generator()->preinitInterface(objectname, "RelativeError", "set", value.str()); // set the width option value.str(""); value << widthOpt_; generator()->preinitInterface(objectname, "WidthOption", "set", value.str()); // set the intermediates option value.str(""); value << inter; generator()->preinitInterface(objectname, "GenerateIntermediates", "set", value.str()); // initialize the decayer decayer->init(); // return the decayer return decayer; } string ThreeBodyDecayConstructor:: DecayerClassName(tcPDPtr incoming, const OrderedParticles & outgoing, string & objname) const { string classname("Herwig::"); // spins of the outgoing particles unsigned int ns(0),nf(0),nv(0); objname += incoming->PDGName() + "2"; for(OrderedParticles::const_iterator it=outgoing.begin(); it!=outgoing.end();++it) { if ((**it).iSpin()==PDT::Spin0 ) ++ns; else if((**it).iSpin()==PDT::Spin1Half) ++nf; else if((**it).iSpin()==PDT::Spin1 ) ++nv; objname += (**it).PDGName(); } objname += "Decayer"; if(incoming->iSpin()==PDT::Spin0) { if(ns==1&&nf==2) classname += "StoSFFDecayer"; else if(nf==2&&nv==1) classname += "StoFFVDecayer"; else classname = ""; } else if(incoming->iSpin()==PDT::Spin1Half) { if(nf==3) classname += "FtoFFFDecayer"; else if(nf==1&&nv==2) classname += "FtoFVVDecayer"; else classname = ""; } else if(incoming->iSpin()==PDT::Spin1) { if(nf==2&&nv==1) classname += "VtoFFVDecayer"; else classname = ""; } else { classname=""; } return classname; } void ThreeBodyDecayConstructor:: createDecayMode(vector & mode, bool possibleOnShell, double symfac) { // convert the diagrams from the general to the three body structure vector diagrams; for(unsigned int iy=0;iyZERO && diagrams.back().intermediate && abs(diagrams.back().intermediate->id())==ParticleID::Wplus) { Energy deltaM = getParticleData(diagrams.back().incoming)->mass() - getParticleData(diagrams.back().outgoing)->mass(); if(deltaMid())==ParticleID::gamma) diagrams.pop_back(); } if(diagrams.empty()) return; // check if possible on-shell internal particles bool inter = interOpt_ == 1 || (interOpt_ == 0 && possibleOnShell); // incoming particle tPDPtr inpart = getParticleData(diagrams[0].incoming); // outgoing particles OrderedParticles outgoing; outgoing.insert(getParticleData(diagrams[0].outgoing)); outgoing.insert(getParticleData(diagrams[0].outgoingPair.first )); outgoing.insert(getParticleData(diagrams[0].outgoingPair.second)); // sort out ordering and labeling of diagrams vector outVector(outgoing.begin(),outgoing.end()); for(unsigned int ix=0;ixid() != diagrams[ix].outgoingPair.first) || ( diagrams[ix].channelType == TBDiagram::channel13 && outVector[0]->id() != diagrams[ix].outgoingPair.first) || ( diagrams[ix].channelType == TBDiagram::channel12 && outVector[0]->id() != diagrams[ix].outgoingPair.first) ) swap(diagrams[ix].outgoingPair.first, diagrams[ix].outgoingPair.second); } // construct the tag for the decay mode string tag = inpart->name() + "->"; unsigned int iprod=0; for(OrderedParticles::const_iterator it = outgoing.begin(); it != outgoing.end(); ++it) { ++iprod; tag += (**it).name(); if(iprod != 3) tag += ","; } tag += ";"; tDMPtr dm = generator()->findDecayMode(tag); // create mode if needed if( createDecayModes() && (!dm || inpart->id() == ParticleID::h0) ) { // create the decayer GeneralThreeBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac); if(!decayer) { if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for " << tag << " so mode not created\n"; return; } OrderedParticles::const_iterator pit=outgoing.begin(); tPDPtr pa = *pit; ++pit; tPDPtr pb = *pit; ++pit; tPDPtr pc = *pit; Energy width = decayer->partialWidth(make_pair(inpart,inpart->mass()), make_pair(pa,pa->mass()) , make_pair(pb,pb->mass()) , make_pair(pc,pc->mass())); if ( Debug::level > 1 ) generator()->log() << "Partial width is: " << width / GeV << " GeV\n"; if(width==ZERO) { if ( Debug::level > 1 ) generator()->log() << "Partial width for " << tag << " zero so mode not created \n"; generator()->preinitRemove(decayer); return; } tDMPtr ndm = generator()->preinitCreateDecayMode(tag); if(!ndm) throw NBodyDecayConstructorError() << "ThreeBodyDecayConstructor::createDecayMode - Needed to create " << "new decaymode but one could not be created for the tag " << tag << Exception::warning; generator()->preinitInterface(ndm, "Decayer", "set", decayer->fullName()); generator()->preinitInterface(ndm, "Active", "set", "Yes"); if(!ndm->decayer()) { generator()->log() << "Can't set the decayer for " << tag << " so mode not created \n"; return; } setBranchingRatio(ndm, width); if(ndm->brat()minimumBR()) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); } // incoming particle is now unstable inpart->stable(false); if(Debug::level > 1 ) { generator()->log() << "Calculated partial width for mode " << tag << " is " << width/GeV << "\n"; } } else if( dm ) { if(dm->brat()minimumBR()) { return; } if((dm->decayer()->fullName()).find("Mambo") != string::npos) { // create the decayer GeneralThreeBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac); if(decayer) { generator()->preinitInterface(dm, "Decayer", "set", decayer->fullName()); // check not zero OrderedParticles::const_iterator pit=outgoing.begin(); tPDPtr pa = *pit; ++pit; tPDPtr pb = *pit; ++pit; tPDPtr pc = *pit; Energy width = decayer->partialWidth(make_pair(inpart,inpart->mass()), make_pair(pa,pa->mass()) , make_pair(pb,pb->mass()) , make_pair(pc,pc->mass())); if(width/(dm->brat()*inpart->width())<1e-10) { string message = "Herwig calculation of the partial width for the decay mode " + inpart->PDGName() + " -> " + pa->PDGName() + " " + pb->PDGName() + " " + pc->PDGName() + " is zero.\n This will cause problems with the calculation of" + " spin correlations.\n It is probably due to inconsistent parameters" + " and decay modes being passed to Herwig via the SLHA file.\n" + " Zeroing the branching ratio for this mode."; setBranchingRatio(dm,ZERO); generator()->logWarning(NBodyDecayConstructorError(message,Exception::warning)); } } // incoming particle is now unstable inpart->stable(false); } } //update CC mode if it exists if( inpart->CC() ) inpart->CC()->synchronize(); } diff --git a/Models/General/ThreeBodyDecayConstructor.h b/Models/General/ThreeBodyDecayConstructor.h --- a/Models/General/ThreeBodyDecayConstructor.h +++ b/Models/General/ThreeBodyDecayConstructor.h @@ -1,167 +1,167 @@ // -*- C++ -*- #ifndef HERWIG_ThreeBodyDecayConstructor_H #define HERWIG_ThreeBodyDecayConstructor_H // // This is the declaration of the ThreeBodyDecayConstructor class. // #include "NBodyDecayConstructorBase.h" #include "ThePEG/Helicity/Vertex/VertexBase.h" #include "TBDiagram.h" #include "PrototypeVertex.h" #include "Herwig/Decay/General/GeneralThreeBodyDecayer.fh" namespace Herwig { using namespace ThePEG; using Helicity::VertexBasePtr; /** * The ThreeBodyDecayConstructor class inherits from the dummy base class * NBodyDecayConstructorBase and implements the necessary functions in * order to create the 3 body decaymodes for a given set of vertices * stored in a Model class. * * @see \ref ThreeBodyDecayConstructorInterfaces "The interfaces" * defined for ThreeBodyDecayConstructor. * @see NBodyDecayConstructor */ class ThreeBodyDecayConstructor: public NBodyDecayConstructorBase { public: /** * The default constructor. */ ThreeBodyDecayConstructor() : includeIntermediatePhotons_(false), interOpt_(0), widthOpt_(1), weakMassCut_(-GeV), intOpt_(1), relErr_(1e-2) {} /** * Function used to determine allowed decaymodes, to be implemented * in derived class. *@param part vector of ParticleData pointers containing particles in model */ - virtual void DecayList(const set & part); + virtual void DecayList(const set & part); /** * Number of outgoing lines. Required for correct ordering. */ virtual unsigned int numBodies() const { return 3; } protected: /** * Create the decayer * @param diagrams The diagrams for the decay * @param inter Option for intermediates */ GeneralThreeBodyDecayerPtr createDecayer(vector & diagrams, bool inter,double symfac) const; /** * Contruct the classname and object name for the Decayer * @param incoming The incoming particle * @param outgoing The decay products * @param objname a string containing the default path of the Decayer object */ string DecayerClassName(tcPDPtr incoming, const OrderedParticles & outgoing, string & objname) const; /** * Create the DecayMode from the diagrams * @param diagrams The diagrams * @param inter Option for intermediates */ virtual void createDecayMode(vector & mode, bool possibleOnShell, double symfac); 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; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ThreeBodyDecayConstructor & operator=(const ThreeBodyDecayConstructor &) = delete; private: /** * Option for intermediate photons */ bool includeIntermediatePhotons_; /** * Option for the inclusion of intermediates */ unsigned int interOpt_; /** * How to treat the widths of the intermediate particles */ unsigned int widthOpt_; /** * Cut off or decays via the weak current */ Energy weakMassCut_; /** * Option for the integration to get the partial width */ unsigned int intOpt_; /** * Relative error for partial width integration */ double relErr_; }; } #endif /* HERWIG_ThreeBodyDecayConstructor_H */ diff --git a/Models/General/TwoBodyDecayConstructor.cc b/Models/General/TwoBodyDecayConstructor.cc --- a/Models/General/TwoBodyDecayConstructor.cc +++ b/Models/General/TwoBodyDecayConstructor.cc @@ -1,446 +1,446 @@ // -*- C++ -*- // // TwoBodyDecayConstructor.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 TwoBodyDecayConstructor class. // #include "TwoBodyDecayConstructor.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "Herwig/Decay/General/GeneralTwoBodyDecayer.h" #include "Herwig/Models/StandardModel/StandardModel.h" #include "ThePEG/PDT/EnumParticles.h" #include "DecayConstructor.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/Utilities/EnumIO.h" #include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractFFSVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractVVSVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractVSSVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractVVTVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractFFTVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractSSTVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractSSSVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractVVVVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractRFSVertex.fh" #include "ThePEG/Helicity/Vertex/AbstractRFVVertex.fh" #include "VectorCurrentDecayConstructor.h" using namespace Herwig; using ThePEG::Helicity::VertexBasePtr; IBPtr TwoBodyDecayConstructor::clone() const { return new_ptr(*this); } IBPtr TwoBodyDecayConstructor::fullclone() const { return new_ptr(*this); } void TwoBodyDecayConstructor::persistentOutput(PersistentOStream & os) const { os << alphaQCD_ << alphaQED_ << oenum(inter_); } void TwoBodyDecayConstructor::persistentInput(PersistentIStream & is, int) { is >> alphaQCD_ >> alphaQED_>> ienum(inter_); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigTwoBodyDecayConstructor("Herwig::TwoBodyDecayConstructor", "Herwig.so"); void TwoBodyDecayConstructor::Init() { static ClassDocumentation documentation ("The TwoBodyDecayConstructor implements to creation of 2 body decaymodes " "and decayers that do not already exist for the given set of vertices."); static Reference interfaceShowerAlphaQCD ("AlphaQCD", "The coupling for QCD corrections", &TwoBodyDecayConstructor::alphaQCD_, false, false, true, false, false); static Reference interfaceShowerAlphaQED ("AlphaQED", "The coupling for QED corrections", &TwoBodyDecayConstructor::alphaQED_, false, false, true, false, false); static Switch interfaceInteractions ("Interactions", "which interactions to include for the hard corrections", &TwoBodyDecayConstructor::inter_, ShowerInteraction::QCD, false, false); static SwitchOption interfaceInteractionsQCD (interfaceInteractions, "QCD", "QCD Only", ShowerInteraction::QCD); static SwitchOption interfaceInteractionsQED (interfaceInteractions, "QED", "QED only", ShowerInteraction::QED); static SwitchOption interfaceInteractionsQCDandQED (interfaceInteractions, "QCDandQED", "Both QCD and QED", ShowerInteraction::Both); } -void TwoBodyDecayConstructor::DecayList(const set & particles) { +void TwoBodyDecayConstructor::DecayList(const set & particles) { // special for weak decays for(unsigned int ix=0;ixdecayConstructors().size();++ix) { Ptr::pointer weak = dynamic_ptr_cast::pointer > (decayConstructor()->decayConstructors()[ix]); if(!weak) continue; weakMassCut_ = max(weakMassCut_,weak->massCut()); } if( particles.empty() ) return; tHwSMPtr model = dynamic_ptr_cast(generator()->standardModel()); unsigned int nv(model->numberOfVertices()); - for(set::const_iterator ip=particles.begin(); + for(set::const_iterator ip=particles.begin(); ip!=particles.end();++ip) { tPDPtr parent = *ip; if ( Debug::level > 0 ) Repository::cout() << "Constructing 2-body decays for " << parent->PDGName() << '\n'; multiset decays; for(unsigned int iv = 0; iv < nv; ++iv) { if(excluded(model->vertex(iv)) || model->vertex(iv)->getNpoint()>3) continue; for(unsigned int il = 0; il < 3; ++il) createModes(parent, model->vertex(iv), il,decays); } if( !decays.empty() ) createDecayMode(decays); } } void TwoBodyDecayConstructor:: createModes(tPDPtr inpart, VertexBasePtr vertex, unsigned int list, multiset & modes) { if( !vertex->isIncoming(inpart) || vertex->getNpoint() != 3 ) return; Energy m1(inpart->mass()); tPDPtr ccpart = inpart->CC() ? inpart->CC() : inpart; long id = ccpart->id(); tPDVector decaylist = vertex->search(list, ccpart); tPDVector::size_type nd = decaylist.size(); for( tPDVector::size_type i = 0; i < nd; i += 3 ) { tPDPtr pa(decaylist[i]), pb(decaylist[i + 1]), pc(decaylist[i + 2]); if( pb->id() == id ) swap(pa, pb); if( pc->id() == id ) swap(pa, pc); //allowed on-shell decay? if( m1 <= pb->mass() + pc->mass() ) continue; // double counting with current decays? if(inpart->iSpin()==PDT::Spin1 && inpart->iCharge()==0 && pb->id()==-pc->id() && abs(pb->id())<=3 && inpart->mass() <= weakMassCut_ ) { continue; } //vertices are defined with all particles incoming modes.insert( TwoBodyDecay(inpart,pb, pc, vertex) ); } } GeneralTwoBodyDecayerPtr TwoBodyDecayConstructor::createDecayer(TwoBodyDecay decay, vector vertices) { string name; using namespace Helicity::VertexType; PDT::Spin in = decay.parent_->iSpin(); // PDT::Spin out1 = decay.children_.first ->iSpin(); PDT::Spin out2 = decay.children_.second->iSpin(); switch(decay.vertex_->getName()) { case FFV : if(in == PDT::Spin1Half) { name = "FFVDecayer"; if(out2==PDT::Spin1Half) swap(decay.children_.first,decay.children_.second); } else { name = "VFFDecayer"; } break; case FFS : if(in == PDT::Spin1Half) { name = "FFSDecayer"; if(out2==PDT::Spin1Half) swap(decay.children_.first,decay.children_.second); } else { name = "SFFDecayer"; } break; case VVS : if(in == PDT::Spin1) { name = "VVSDecayer"; if(out2==PDT::Spin1) swap(decay.children_.first,decay.children_.second); } else { name = "SVVDecayer"; } break; case VSS : if(in == PDT::Spin1) { name = "VSSDecayer"; } else { name = "SSVDecayer"; if(out2==PDT::Spin0) swap(decay.children_.first,decay.children_.second); } break; case VVT : name = in==PDT::Spin2 ? "TVVDecayer" : "Unknown"; break; case FFT : name = in==PDT::Spin2 ? "TFFDecayer" : "Unknown"; break; case SST : name = in==PDT::Spin2 ? "TSSDecayer" : "Unknown"; break; case SSS : name = "SSSDecayer"; break; case VVV : name = "VVVDecayer"; break; case RFS : if(in==PDT::Spin1Half) { name = "FRSDecayer"; if(out2==PDT::Spin3Half) swap(decay.children_.first,decay.children_.second); } else if(in==PDT::Spin0) { name = "SRFDecayer"; if(out2==PDT::Spin3Half) swap(decay.children_.first,decay.children_.second); } else { name = "Unknown"; } break; case RFV : if(in==PDT::Spin1Half) { name = "FRVDecayer"; if(out2==PDT::Spin3Half) swap(decay.children_.first,decay.children_.second); } else name = "Unknown"; break; default : Throw() << "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. " << "Decay is " << decay.parent_->PDGName() << " -> " << decay.children_.first ->PDGName() << " " << decay.children_.second->PDGName() << Exception::runerror; } if(name=="Unknown") Throw() << "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. " << "Decay is " << decay.parent_->PDGName() << " -> " << decay.children_.first ->PDGName() << " " << decay.children_.second->PDGName() << Exception::runerror; ostringstream fullname; fullname << "/Herwig/Decays/" << name << "_" << decay.parent_->PDGName() << "_" << decay.children_.first ->PDGName() << "_" << decay.children_.second->PDGName(); string classname = "Herwig::" + name; GeneralTwoBodyDecayerPtr decayer; decayer = dynamic_ptr_cast (generator()->preinitCreate(classname,fullname.str())); if(!decayer) Throw() << "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. " << "Decay is " << decay.parent_->PDGName() << " -> " << decay.children_.first ->PDGName() << " " << decay.children_.second->PDGName() << Exception::runerror; // set the strong coupling for radiation generator()->preinitInterface(decayer, "AlphaS" , "set", alphaQCD_->fullName()); // set the EM coupling for radiation generator()->preinitInterface(decayer, "AlphaEM", "set", alphaQED_->fullName()); // set the type of interactions for the correction if(inter_==ShowerInteraction::QCD) generator()->preinitInterface(decayer, "Interactions", "set", "QCD"); else if(inter_==ShowerInteraction::QED) generator()->preinitInterface(decayer, "Interactions", "set", "QED"); else generator()->preinitInterface(decayer, "Interactions", "set", "QCDandQED"); // get the vertices for radiation from the external legs map inRad,fourRad; vector > outRad(2); vector itemp={ShowerInteraction::QCD,ShowerInteraction::QED}; for(auto & inter : itemp) { inRad[inter] = radiationVertex(decay.parent_,inter); outRad[0][inter] = radiationVertex(decay.children_.first ,inter); outRad[1][inter] = radiationVertex(decay.children_.second,inter); // get any contributing 4 point vertices fourRad[inter] = radiationVertex(decay.parent_,inter, decay.children_); } // set info on decay decayer->setDecayInfo(decay.parent_,decay.children_,vertices, inRad,outRad,fourRad); // initialised the decayer setDecayerInterfaces(fullname.str()); decayer->init(); return decayer; } void TwoBodyDecayConstructor:: createDecayMode(multiset & decays) { tPDPtr inpart = decays.begin()->parent_; for( multiset::iterator dit = decays.begin(); dit != decays.end(); ) { TwoBodyDecay mode = *dit; // get all the moees with the same in and outgoing particles pair::iterator, multiset::iterator> range = decays.equal_range(mode); // construct the decay mode tPDPtr pb((mode).children_.first), pc((mode).children_.second); string tag = inpart->name() + "->" + pb->name() + "," + pc->name() + ";"; // Does it exist already ? tDMPtr dm = generator()->findDecayMode(tag); // find the vertices vector vertices; for ( multiset::iterator dit2 = range.first; dit2 != range.second; ++dit2) { vertices.push_back(dit2->vertex_); } dit=range.second; // now create DecayMode objects that do not already exist if( createDecayModes() && (!dm || inpart->id() == ParticleID::h0) ) { tDMPtr ndm = generator()->preinitCreateDecayMode(tag); if(ndm) { inpart->stable(false); GeneralTwoBodyDecayerPtr decayer=createDecayer(mode,vertices); if(!decayer) continue; generator()->preinitInterface(ndm, "Decayer", "set", decayer->fullName()); generator()->preinitInterface(ndm, "Active", "set", "Yes"); Energy width = decayer->partialWidth(make_pair(inpart,inpart->mass()), make_pair(pb,pb->mass()) , make_pair(pc,pc->mass())); setBranchingRatio(ndm, width); if(width==ZERO || ndm->brat()minimumBR()) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); } } else Throw() << "TwoBodyDecayConstructor::createDecayMode - Needed to create " << "new decaymode but one could not be created for the tag " << tag << Exception::warning; } else if( dm ) { if(dm->brat()minimumBR()) { continue; } if((dm->decayer()->fullName()).find("Mambo") != string::npos) { inpart->stable(false); GeneralTwoBodyDecayerPtr decayer=createDecayer(mode,vertices); if(!decayer) continue; generator()->preinitInterface(dm, "Decayer", "set", decayer->fullName()); Energy width = decayer->partialWidth(make_pair(inpart,inpart->mass()), make_pair(pb,pb->mass()) , make_pair(pc,pc->mass())); if(width/(dm->brat()*inpart->width())<1e-10) { string message = "Herwig calculation of the partial width for the decay mode " + inpart->PDGName() + " -> " + pb->PDGName() + " " + pc->PDGName() + " is zero.\n This will cause problems with the calculation of" + " spin correlations.\n It is probably due to inconsistent parameters" + " and decay modes being passed to Herwig via the SLHA file.\n" + " Zeroing the branching ratio for this mode."; setBranchingRatio(dm,ZERO); generator()->logWarning(NBodyDecayConstructorError(message,Exception::warning)); } } } } // update CC mode if it exists if( inpart->CC() ) inpart->CC()->synchronize(); } VertexBasePtr TwoBodyDecayConstructor::radiationVertex(tPDPtr particle, ShowerInteraction inter, tPDPair children) { tHwSMPtr model = dynamic_ptr_cast(generator()->standardModel()); map::iterator rit = radiationVertices_[inter].find(particle); tPDPtr cc = particle->CC() ? particle->CC() : particle; if(children==tPDPair() && rit!=radiationVertices_[inter].end()) return rit->second; unsigned int nv(model->numberOfVertices()); long bosonID = inter==ShowerInteraction::QCD ? ParticleID::g : ParticleID::gamma; tPDPtr gluon = getParticleData(bosonID); // look for radiation vertices for incoming and outgoing particles for(unsigned int iv=0;ivvertex(iv); // look for 3 point vertices if (children==tPDPair()){ if( !vertex->isIncoming(particle) || vertex->getNpoint() != 3 || !vertex->isOutgoing(particle) || !vertex->isOutgoing(gluon)) continue; for(unsigned int list=0;list<3;++list) { tPDVector decaylist = vertex->search(list, particle); for( tPDVector::size_type i = 0; i < decaylist.size(); i += 3 ) { tPDPtr pa(decaylist[i]), pb(decaylist[i + 1]), pc(decaylist[i + 2]); if( pb->id() == bosonID ) swap(pa, pb); if( pc->id() == bosonID ) swap(pa, pc); if( pb->id() != particle->id()) swap(pb, pc); if( pa->id() != bosonID) continue; if( pb != particle) continue; if( pc != cc) continue; radiationVertices_[inter][particle] = vertex; return vertex; } } } // look for 4 point vertex including a gluon else { if( !vertex->isIncoming(particle) || vertex->getNpoint()!=4 || !vertex->isOutgoing(children.first) || !vertex->isOutgoing(children.second) || !vertex->isOutgoing(gluon)) continue; for(unsigned int list=0;list<4;++list) { tPDVector decaylist = vertex->search(list, particle); for( tPDVector::size_type i = 0; i < decaylist.size(); i += 4 ) { tPDPtr pa(decaylist[i]), pb(decaylist[i+1]), pc(decaylist[i+2]), pd(decaylist[i+3]); // order so that a = g, b = parent if( pb->id() == bosonID ) swap(pa, pb); if( pc->id() == bosonID ) swap(pa, pc); if( pd->id() == bosonID ) swap(pa, pd); if( pc->id() == particle->id()) swap(pb, pc); if( pd->id() == particle->id()) swap(pb, pd); if( pa->id() != bosonID) continue; if( pb->id() != particle->id()) continue; if( !((abs(pd->id()) == abs(children. first->id()) && abs(pc->id()) == abs(children.second->id())) || (abs(pc->id()) == abs(children. first->id()) && abs(pd->id()) == abs(children.second->id())))) continue; return vertex; } } } } return VertexBasePtr(); } diff --git a/Models/General/TwoBodyDecayConstructor.h b/Models/General/TwoBodyDecayConstructor.h --- a/Models/General/TwoBodyDecayConstructor.h +++ b/Models/General/TwoBodyDecayConstructor.h @@ -1,182 +1,182 @@ // -*- C++ -*- // // TwoBodyDecayConstructor.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_TwoBodyDecayConstructor_H #define HERWIG_TwoBodyDecayConstructor_H // // This is the declaration of the TwoBodyDecayConstructor class. // #include "NBodyDecayConstructorBase.h" #include "ThePEG/Helicity/Vertex/VertexBase.h" #include "Herwig/Decay/General/GeneralTwoBodyDecayer.fh" #include "Herwig/Shower/ShowerAlpha.h" #include "TwoBodyDecay.h" namespace Herwig { using namespace ThePEG; using Helicity::VertexBasePtr; using Helicity::tVertexBasePtr; /** * The TwoBodyDecayConstructor class inherits from the dummy base class * NBodyDecayConstructorBase and implements the necessary functions in * order to create the 2 body decay modes for a given set of vertices * stored in a Model class. * * @see \ref TwoBodyDecayConstructorInterfaces "The interfaces" * defined for TwoBodyDecayConstructor. * @see NBodyDecayConstructor **/ class TwoBodyDecayConstructor: public NBodyDecayConstructorBase { public: /** * The default constructor. */ TwoBodyDecayConstructor() : inter_(ShowerInteraction::Both), weakMassCut_(-GeV) { radiationVertices_[ShowerInteraction::QCD] = map(); radiationVertices_[ShowerInteraction::QED] = map(); } /** * Function used to determine allowed decaymodes *@param part vector of ParticleData pointers containing particles in model */ - virtual void DecayList(const set & part); + virtual void DecayList(const set & part); /** * Number of outgoing lines. Required for correct ordering. */ virtual unsigned int numBodies() const { return 2; } 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; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ TwoBodyDecayConstructor & operator=(const TwoBodyDecayConstructor &) = delete; private: /** @name Functions to create decayers and decaymodes. */ //@{ /** * Function to create decays * @param inpart Incoming particle * @param vert The vertex to create decays for * @param ilist Which list to search * @param iv Row number in _theExistingDecayers member * @return A vector a decay modes */ void createModes(tPDPtr inpart, VertexBasePtr vert, unsigned int ilist, multiset & modes); /** * Function to create decayer for specific vertex * @param decay decay mode for this decay * member variable */ GeneralTwoBodyDecayerPtr createDecayer(TwoBodyDecay decay, vector ); /** * Create decay mode(s) from given part and decay modes * @param decays The vector of decay modes * @param decayer The decayer responsible for this decay */ void createDecayMode(multiset & decays); //@} /** * Get the vertex for QED/QCD radiation */ VertexBasePtr radiationVertex(tPDPtr particle,ShowerInteraction inter, tPDPair children = tPDPair ()); private: /** * Map of particles and the vertices which generate their QCD * radiation */ map > radiationVertices_; /** * Default choice for the strong coupling object for hard QCD radiation */ ShowerAlphaPtr alphaQCD_; /** * Default choice for the strong coupling object for hard QED radiation */ ShowerAlphaPtr alphaQED_; /** * Which type of corrections to the decays to include */ ShowerInteraction inter_; /** * Cut off or decays via the weak current */ Energy weakMassCut_; }; } #endif /* HERWIG_TwoBodyDecayConstructor_H */ diff --git a/Models/General/VectorCurrentDecayConstructor.cc b/Models/General/VectorCurrentDecayConstructor.cc --- a/Models/General/VectorCurrentDecayConstructor.cc +++ b/Models/General/VectorCurrentDecayConstructor.cc @@ -1,172 +1,172 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the VectorCurrentDecayConstructor class. // #include "VectorCurrentDecayConstructor.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/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Decay/General/VectorCurrentDecayer.h" namespace { struct ParticleOrdering { /** * Operator for the ordering * @param p1 The first ParticleData object * @param p2 The second ParticleData object */ bool operator() (tcPDPtr p1, tcPDPtr p2) const { return abs(p1->id()) > abs(p2->id()) || ( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) || ( p1->id() == p2->id() && p1->fullName() > p2->fullName() ); } }; } using namespace Herwig; IBPtr VectorCurrentDecayConstructor::clone() const { return new_ptr(*this); } IBPtr VectorCurrentDecayConstructor::fullclone() const { return new_ptr(*this); } void VectorCurrentDecayConstructor::persistentOutput(PersistentOStream & os) const { os << ounit(massCut_,GeV) << current_; } void VectorCurrentDecayConstructor::persistentInput(PersistentIStream & is, int) { is >> iunit(massCut_,GeV) >> current_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigVectorCurrentDecayConstructor("Herwig::VectorCurrentDecayConstructor", "Herwig.so"); void VectorCurrentDecayConstructor::Init() { static ClassDocumentation documentation ("The VectorCurrentDecayConstructor class constructs the decays of low mass vector bosons" " to hadrons via the weak current"); static RefVector interfaceCurrent ("Current", "The current for the decay mode", &VectorCurrentDecayConstructor::current_, -1, false, false, true, false, false); static Parameter interfaceMassCut ("MassCut", "The maximum mass difference for the decay", &VectorCurrentDecayConstructor::massCut_, GeV, 2.0*GeV, 1.0*GeV, 10.0*GeV, false, false, Interface::limited); } void VectorCurrentDecayConstructor::doinit() { NBodyDecayConstructorBase::doinit(); model_ = dynamic_ptr_cast::pointer>(generator()->standardModel()); } -void VectorCurrentDecayConstructor::DecayList(const set & particles) { +void VectorCurrentDecayConstructor::DecayList(const set & particles) { if( particles.empty() ) return; unsigned int nv(model_->numberOfVertices()); for(PDPtr part : particles) { if(part->iSpin()!=PDT::Spin1) continue; if(part->iCharge()!=0) continue; bool foundD(false),foundU(false),foundS(false); if(part->mass()>massCut_) continue; for(unsigned int iv = 0; iv < nv; ++iv) { VertexBasePtr vertex = model_->vertex(iv); if( !vertex->isIncoming(part) || vertex->getNpoint() != 3 ) continue; for(unsigned int iloc = 0;iloc < 3; ++iloc) { vector ext = vertex->search(iloc, part->id()); if(ext.empty()) continue; for(unsigned int ioff=0;ioffinit(); for(unsigned int imode=0;imodenumberOfModes();++imode) { // get the external particles for this mode int iq(0),ia(0); tPDVector out = current->particles(0,imode,iq,ia); current->decayModeInfo(imode,iq,ia); if(iq==2&&ia==-2) continue; // order the particles bool skip=false; for(unsigned int ix=0;ix outgoing(out.begin(),out.end()); Energy minMass(ZERO); string tag = part->PDGName() + "->"; bool first=false; int charge(0); for(tcPDPtr part : outgoing) { if(!first) first=true; else tag+=","; tag+=part->PDGName(); minMass+=part->mass(); charge+=part->iCharge(); } tag+=";"; if(minMass>part->mass()) continue; if(charge!=0) continue; // create the decayer ostringstream fullname; fullname << "/Herwig/Decays/DMMediator_" << part->PDGName(); for(tcPDPtr part : out) fullname << "_" << part->PDGName(); string classname = "Herwig::VectorCurrentDecayer"; VectorCurrentDecayerPtr decayer = dynamic_ptr_cast (generator()->preinitCreate(classname,fullname.str())); decayer->setDecayInfo(part,out,current); // // set decayer options from base class // setDecayerInterfaces(fullname.str()); // initialize the decayer decayer->init(); // calculate the width Energy pWidth = decayer->partialWidth(part,out); if(pWidth<=ZERO) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); continue; } tDMPtr ndm = generator()->preinitCreateDecayMode(tag); generator()->preinitInterface(ndm, "Decayer", "set", decayer->fullName()); part->stable(false); generator()->preinitInterface(ndm, "Active", "set", "Yes"); setBranchingRatio(ndm, pWidth); } } } } diff --git a/Models/General/VectorCurrentDecayConstructor.h b/Models/General/VectorCurrentDecayConstructor.h --- a/Models/General/VectorCurrentDecayConstructor.h +++ b/Models/General/VectorCurrentDecayConstructor.h @@ -1,131 +1,131 @@ // -*- C++ -*- #ifndef Herwig_VectorCurrentDecayConstructor_H #define Herwig_VectorCurrentDecayConstructor_H // // This is the declaration of the VectorCurrentDecayConstructor class. // #include "NBodyDecayConstructorBase.h" #include "Herwig/Decay/WeakCurrents/WeakCurrent.h" namespace Herwig { using namespace ThePEG; /** * The VectorCurrentDecayConstructor class constructs the decay of low mass vector bosons via the weak currents. * * @see \ref VectorCurrentDecayConstructorInterfaces "The interfaces" * defined for VectorCurrentDecayConstructor. */ class VectorCurrentDecayConstructor: public NBodyDecayConstructorBase { public: /** * The default constructor. */ VectorCurrentDecayConstructor() : massCut_(2.*GeV) {} /** * Function used to determine allowed decaymodes, to be implemented * in derived class. *@param part vector of ParticleData pointers containing particles in model */ - virtual void DecayList(const set & part); + virtual void DecayList(const set & part); /** * Number of outgoing lines. Required for correct ordering (do this one next-to-last) */ virtual unsigned int numBodies() const { return 999; } /** * Cut off */ Energy massCut() const { return massCut_;} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ VectorCurrentDecayConstructor & operator=(const VectorCurrentDecayConstructor &) = delete; private: /** * Model Pointer */ Ptr::pointer model_; /** * Cut-off on the mass difference */ Energy massCut_; /** * The current for the mode */ vector current_; }; } #endif /* Herwig_VectorCurrentDecayConstructor_H */ diff --git a/Models/General/WeakCurrentDecayConstructor.cc b/Models/General/WeakCurrentDecayConstructor.cc --- a/Models/General/WeakCurrentDecayConstructor.cc +++ b/Models/General/WeakCurrentDecayConstructor.cc @@ -1,295 +1,295 @@ // -*- C++ -*- // // WeakCurrentDecayConstructor.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 WeakCurrentDecayConstructor class. // #include "WeakCurrentDecayConstructor.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh" #include "DecayConstructor.h" using namespace Herwig; using ThePEG::Helicity::VertexBasePtr; IBPtr WeakCurrentDecayConstructor::clone() const { return new_ptr(*this); } IBPtr WeakCurrentDecayConstructor::fullclone() const { return new_ptr(*this); } void WeakCurrentDecayConstructor::doinit() { NBodyDecayConstructorBase::doinit(); model_ = dynamic_ptr_cast::pointer>(generator()->standardModel()); unsigned int isize=decayTags_.size(); if(isize!=_norm .size()||isize!=_current.size()) throw InitException() << "Invalid sizes for the decay mode vectors in " << " WeakCurrentDecayConstructor " << decayTags_.size() << " " << _norm.size() << " " << _current.size() << Exception::runerror; // get the particles from the tags for(unsigned int ix=0;ixinit(); particles_.push_back(vector()); string tag=decayTags_[ix]; do { string::size_type next = min(tag.find(','), tag.find(';')); particles_.back().push_back(generator()->findParticle(tag.substr(0,next))); if(!particles_.back().back()) throw Exception() << "Failed to find particle " << tag.substr(0,next) << " in DecayMode " << decayTags_[ix] << " in WeakCurrentDecayConstructor::doinit()" << Exception::runerror; if(tag[next]==';') break; tag = tag.substr(next+1); } while(true); } } void WeakCurrentDecayConstructor::persistentOutput(PersistentOStream & os) const { os << ounit(massCut_,GeV) << decayTags_ << particles_ << _norm << _current; } void WeakCurrentDecayConstructor::persistentInput(PersistentIStream & is, int) { is >> iunit(massCut_,GeV) >> decayTags_ >> particles_ >> _norm >> _current; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigWeakCurrentDecayConstructor("Herwig::WeakCurrentDecayConstructor", "Herwig.so"); void WeakCurrentDecayConstructor::Init() { static ClassDocumentation documentation ("The WeakCurrentDecayConstructor class implemets the decay of BSM particles " "to low mass hadronic states using the Weak current"); static ParVector interfaceDecayModes ("DecayModes", "The decays of the weak current", &WeakCurrentDecayConstructor::decayTags_, -1, "", "", "", false, false, Interface::nolimits); static ParVector interfaceNormalisation ("Normalisation", "The normalisation of the different modes", &WeakCurrentDecayConstructor::_norm, -1, 1.0, 0.0, 10.0, false, false, Interface::limited); static RefVector interfaceCurrent ("Current", "The current for the decay mode", &WeakCurrentDecayConstructor::_current, -1, false, false, true, false, false); static Parameter interfaceMassCut ("MassCut", "The maximum mass difference for the decay", &WeakCurrentDecayConstructor::massCut_, GeV, 5.0*GeV, 1.0*GeV, 10.0*GeV, false, false, Interface::limited); } -void WeakCurrentDecayConstructor::DecayList(const set & part) { +void WeakCurrentDecayConstructor::DecayList(const set & part) { if( part.empty() ) return; unsigned int nv(model_->numberOfVertices()); - for(set::const_iterator ip=part.begin();ip!=part.end();++ip) { + for(set::const_iterator ip=part.begin();ip!=part.end();++ip) { for(unsigned int iv = 0; iv < nv; ++iv) { for(unsigned int ilist = 0; ilist < 3; ++ilist) { vector decays = createModes(*ip, model_->vertex(iv),ilist); if(!decays.empty()) createDecayMode(decays); } } } } vector WeakCurrentDecayConstructor::createModes(const PDPtr inpart, const VertexBasePtr vert, unsigned int ilist) { int id = inpart->id(); if( !vert->isIncoming(inpart) || vert->getNpoint() != 3 ) return vector(); Energy m1(inpart->mass()); vector decaylist; decaylist = vert->search(ilist,inpart); tPDVector::size_type nd = decaylist.size(); vector decays; for( tPDVector::size_type i = 0; i < nd; i += 3 ) { tPDPtr pa(decaylist[i]), pb(decaylist.at(i + 1)), pc(decaylist.at(i + 2)); if( pb->id() == id ) swap(pa, pb); if( pc->id() == id ) swap(pa, pc); //One of the products must be a W Energy mp(ZERO); if( abs(pb->id()) == ParticleID::Wplus ) mp = pc->mass(); else if( abs(pc->id()) == ParticleID::Wplus ) mp = pb->mass(); else continue; //allowed on-shell decay and passes mass cut if( m1 >= pb->mass() + pc->mass() ) continue; if( m1 < mp ) continue; if( m1 - mp >= massCut_ ) continue; //vertices are defined with all particles incoming if( pb->CC() ) pb = pb->CC(); if( pc->CC() ) pc = pc->CC(); decays.push_back( TwoBodyDecay(inpart,pb, pc, vert) ); if(abs(decays.back().children_.second->id())!=ParticleID::Wplus) swap(decays.back().children_.first,decays.back().children_.second); assert(abs(decays.back().children_.second->id())==ParticleID::Wplus); } return decays; } GeneralCurrentDecayerPtr WeakCurrentDecayConstructor::createDecayer(PDPtr in, PDPtr out1, vector outCurrent, VertexBasePtr vertex, WeakCurrentPtr current) { string name; using namespace ThePEG::Helicity::VertexType; switch(vertex->getName()) { case FFV : name = "FFVCurrentDecayer"; break; default : ostringstream message; message << "Invalid vertex for decays of " << in->PDGName() << " -> " << out1->PDGName() << " via weak current " << vertex->fullName() << "\n"; generator()->logWarning(NBodyDecayConstructorError(message.str(), Exception::warning)); return GeneralCurrentDecayerPtr(); } ostringstream fullname; fullname << "/Herwig/Decays/" << name << "_" << in->PDGName() << "_" << out1->PDGName(); for(unsigned int ix=0;ixPDGName(); string classname = "Herwig::" + name; GeneralCurrentDecayerPtr decayer = dynamic_ptr_cast (generator()->preinitCreate(classname,fullname.str())); decayer->setDecayInfo(in,out1,outCurrent,vertex,current,massCut_); // set decayer options from base class setDecayerInterfaces(fullname.str()); // initialize the decayer decayer->init(); // return the decayer return decayer; } void WeakCurrentDecayConstructor:: createDecayMode(vector & decays) { assert(!decays.empty()); for(unsigned int ix = 0; ix < decays.size(); ++ix) { PDVector particles(3); particles[0] = decays[ix].parent_; particles[1] = decays[ix].children_.first ; bool Wplus=decays[ix].children_.second->id()==ParticleID::Wplus; for(unsigned int iy=0;iy<_current.size();++iy) { particles.resize(2); vector wprod=particles_[iy]; int icharge=0; Energy msum = particles[0]->mass()-particles[1]->mass(); for(unsigned int iz=0;iziCharge(); msum -=wprod[iz]->mass(); } if(msum<=ZERO) continue; bool cc = (Wplus&&icharge==-3)||(!Wplus&&icharge==3); OrderedParticles outgoing; outgoing.insert(particles[1]); for(unsigned int iz=0;izCC()) wprod[iz]=wprod[iz]->CC(); outgoing.insert(wprod[iz]); } // check outgoing particles initialised for(unsigned int iz=0;izinit(); // create the tag for the decay mode string tag = particles[0]->PDGName() + "->"; OrderedParticles::const_iterator it = outgoing.begin(); do { tag += (**it).name(); ++it; if(it!=outgoing.end()) tag +=","; else tag +=";"; } while(it!=outgoing.end()); // find the decay mode tDMPtr dm= generator()->findDecayMode(tag); if( !dm && createDecayModes() ) { // create the decayer GeneralCurrentDecayerPtr decayer = createDecayer(particles[0],particles[1], wprod,decays[ix].vertex_, _current[iy]); if(!decayer) continue; // calculate the width Energy pWidth = _norm[iy]*decayer->partialWidth(particles[0],particles[1],wprod); if(pWidth<=ZERO) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); continue; } tDMPtr ndm = generator()->preinitCreateDecayMode(tag); if(!ndm) throw NBodyDecayConstructorError() << "WeakCurrentDecayConstructor::createDecayMode - Needed to create " << "new decaymode but one could not be created for the tag " << tag << Exception::warning; generator()->preinitInterface(ndm, "Decayer", "set", decayer->fullName()); generator()->preinitInterface(ndm, "Active", "set", "Yes"); setBranchingRatio(ndm, pWidth); particles[0]->stable(false); if(ndm->brat()minimumBR()) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); } } else if (dm) { // create the decayer GeneralCurrentDecayerPtr decayer = createDecayer(particles[0],particles[1], wprod,decays[ix].vertex_, _current[iy]); if(!decayer) continue; generator()->preinitInterface(dm, "Decayer", "set", decayer->fullName()); particles[0]->stable(false); if(createDecayModes()) { // calculate the width Energy pWidth = _norm[iy]*decayer->partialWidth(particles[0],particles[1],wprod); if(pWidth<=ZERO) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); continue; } generator()->preinitInterface(dm, "Active", "set", "Yes"); particles[0]->width(particles[0]->width()*(1.-dm->brat())); setBranchingRatio(dm, pWidth); } if(dm->brat()minimumBR()) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); } } } } } diff --git a/Models/General/WeakCurrentDecayConstructor.h b/Models/General/WeakCurrentDecayConstructor.h --- a/Models/General/WeakCurrentDecayConstructor.h +++ b/Models/General/WeakCurrentDecayConstructor.h @@ -1,194 +1,194 @@ // -*- C++ -*- // // WeakCurrentDecayConstructor.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_WeakCurrentDecayConstructor_H #define HERWIG_WeakCurrentDecayConstructor_H // // This is the declaration of the WeakCurrentDecayConstructor class. // #include "NBodyDecayConstructorBase.h" #include "ThePEG/Helicity/Vertex/VertexBase.h" #include "Herwig/Decay/General/GeneralCurrentDecayer.fh" #include "Herwig/Models/StandardModel/StandardModel.h" #include "Herwig/Decay/WeakCurrents/WeakCurrent.h" #include "Herwig/Decay/General/GeneralCurrentDecayer.h" #include "TwoBodyDecay.h" namespace Herwig { using namespace ThePEG; /** * Here is the documentation of the WeakCurrentDecayConstructor class. * * @see \ref WeakCurrentDecayConstructorInterfaces "The interfaces" * defined for WeakCurrentDecayConstructor. */ class WeakCurrentDecayConstructor: public NBodyDecayConstructorBase { public: /** * The default constructor. */ WeakCurrentDecayConstructor() : massCut_(5.*GeV) {} /** * Function used to determine allowed decaymodes, to be implemented * in derived class. *@param part vector of ParticleData pointers containing particles in model */ - virtual void DecayList(const set & part); + virtual void DecayList(const set & part); /** * Number of outgoing lines. Required for correct ordering (do this one last) */ virtual unsigned int numBodies() const { return 1000; } /** * Cut off */ Energy massCut() const { return massCut_;} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** @name Functions to create decayers and decaymodes. */ //@{ /** * Function to create decays * @param inpart Incoming particle * @param vert The vertex to create decays for * @param ilist Which list to search * @param iv Row number in _theExistingDecayers member * @return vector of ParticleData ptrs */ vector createModes(const PDPtr inpart,const VertexBasePtr vert, unsigned int ilist); /** * Function to create decayer for specific vertex * @param vert Pointer to vertex * @param icol Integer referring to the colmun in _theExistingDecayers * @param ivert Integer referring to the row in _theExistingDecayers * member variable */ GeneralCurrentDecayerPtr createDecayer(PDPtr in, PDPtr out1, vector outCurrent, VertexBasePtr vertex, WeakCurrentPtr current); /** * Create decay mode(s) from given part and decay modes * @param inpart pointer to incoming particle * @param decays list of allowed interactions * @param decayer The decayer responsible for this decay */ void createDecayMode(vector & decays); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ WeakCurrentDecayConstructor & operator=(const WeakCurrentDecayConstructor &) = delete; private: /** * Model Pointer */ Ptr::pointer model_; /** * Cut-off on the mass difference */ Energy massCut_; /** * Tags for the modes */ vector decayTags_; /** * Particles for the mode */ vector > particles_; /** * Normalisation */ vector _norm; /** * The current for the mode */ vector _current; }; } #endif /* HERWIG_WeakCurrentDecayConstructor_H */