diff --git a/Models/General/WeakCurrentDecayConstructor.cc b/Models/General/WeakCurrentDecayConstructor.cc --- a/Models/General/WeakCurrentDecayConstructor.cc +++ b/Models/General/WeakCurrentDecayConstructor.cc @@ -1,296 +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(); - _theModel = dynamic_ptr_cast::pointer> - (generator()->standardModel()); + 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; + os << ounit(massCut_,GeV) << decayTags_ << particles_ << _norm << _current; } void WeakCurrentDecayConstructor::persistentInput(PersistentIStream & is, int) { - is >> iunit(_masscut,GeV) >> decayTags_ >> particles_ >> _norm >> _current; + 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, + &WeakCurrentDecayConstructor::massCut_, GeV, 5.0*GeV, 1.0*GeV, 10.0*GeV, false, false, Interface::limited); } void WeakCurrentDecayConstructor::DecayList(const set & part) { if( part.empty() ) return; - unsigned int nv(_theModel->numberOfVertices()); + unsigned int nv(model_->numberOfVertices()); 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, _theModel->vertex(iv),ilist); + 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; + 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); + 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) {} + 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); /** * 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;} + 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 _theModel; + Ptr::pointer model_; /** * Cut-off on the mass difference */ - Energy _masscut; + 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 */