diff --git a/Decay/Radiation/QEDRadiationHandler.cc b/Decay/Radiation/QEDRadiationHandler.cc --- a/Decay/Radiation/QEDRadiationHandler.cc +++ b/Decay/Radiation/QEDRadiationHandler.cc @@ -1,171 +1,171 @@ // -*- C++ -*- // // QEDRadiationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the QEDRadiationHandler class. // #include "QEDRadiationHandler.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "DecayRadiationGenerator.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Repository/EventGenerator.h" #include "Herwig/Decay/DecayIntegrator.h" #include "Herwig/Decay/DecayPhaseSpaceMode.h" #include "ThePEG/PDT/DecayMode.h" using namespace Herwig; namespace { /** * A struct to order the particles in the same way as in the DecayMode's */ struct ParticleOrdering { /** * Operator for the ordering * @param p1 The first ParticleData object * @param p2 The second ParticleData object */ - bool operator()(cPDPtr p1, cPDPtr p2) { + bool operator()(cPDPtr p1, cPDPtr 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() ); } }; /** * A set of ParticleData objects ordered as for the DecayMode's */ typedef multiset OrderedParticles; } QEDRadiationHandler::QEDRadiationHandler() { // only include electroweak gauge bosons _decayingParticles.push_back( 22); _decayingParticles.push_back( 23); _decayingParticles.push_back( 24); _decayingParticles.push_back(-24); // only include the charged leptons _decayProducts.push_back( 11); _decayProducts.push_back( 13); _decayProducts.push_back( 15); _decayProducts.push_back(-11); _decayProducts.push_back(-13); _decayProducts.push_back(-15); } void QEDRadiationHandler:: handle(EventHandler & eh, const tPVector & tagged, const Hint &) { // find the potential decaying particles to be considered set parents; for(unsigned int ix=0;ixid(); if(find(_decayProducts.begin(),_decayProducts.end(),id)!=_decayProducts.end()) { PPtr par=tagged[ix]->parents()[0]; id=par->id(); if(tagged[ix]->parents()[0]->mass()>ZERO&& find(_decayingParticles.begin(),_decayingParticles.end(),id)!= _decayingParticles.end()) parents.insert(par); } } set::const_iterator sit; StepPtr step=eh.currentStep(); // loop over parents for(sit=parents.begin();sit!=parents.end();++sit) { // extract children ParticleVector children=(**sit).children(); // store number of children unsigned int initsize = children.size(); // boost the decay to the parent rest frame Boost boost = - (**sit).momentum().boostVector(); (**sit).deepBoost(boost); // construct the tag for the decay mode to find the decayer OrderedParticles products; for(unsigned int ix=0;ixdataPtr()); string tag = (**sit).dataPtr()->name() + "->"; unsigned int iprod=0; for(OrderedParticles::const_iterator it = products.begin(); it != products.end(); ++it) { ++iprod; tag += (**it).name(); if(iprod != initsize) tag += ","; } tag += ";"; tDMPtr dm = generator()->findDecayMode(tag); tDecayIntegratorPtr decayer; if(dm) { tDecayerPtr dtemp = dm->decayer(); decayer = dynamic_ptr_cast(dtemp); bool cc; tPDVector ctemp; for(unsigned int ix=0;ix(children[ix]->dataPtr())); unsigned int imode = decayer->modeNumber(cc,(**sit).dataPtr(),ctemp); decayer->me2(imode,**sit,children,DecayIntegrator::Initialize); } // generate photons ParticleVector newchildren = _generator->generatePhotons(**sit,children,decayer); // if photons produced add as children and to step for(unsigned int ix=initsize;ixaddDecayProduct(newchildren[ix]); } (**sit).deepBoost(-boost); } } void QEDRadiationHandler::persistentOutput(PersistentOStream & os) const { os << _generator << _decayingParticles << _decayProducts; } void QEDRadiationHandler::persistentInput(PersistentIStream & is, int) { is >> _generator >> _decayingParticles >> _decayProducts; } ClassDescription QEDRadiationHandler::initQEDRadiationHandler; // Definition of the static class description member. void QEDRadiationHandler::Init() { static ClassDocumentation documentation ("The QEDRadiationHandler class is designed to be used as a PostSubProcessHandler" "so that the same approach as for radiation in decays can be used for resonances" "produced as part of the hard process"); static Reference interfaceRadiationGenerator ("RadiationGenerator", "Reference to the DecayRadiationGenerator", &QEDRadiationHandler::_generator, false, false, true, false, false); static ParVector interfaceDecayingParticles ("DecayingParticles", "List of PDF codes of the particles which should have radiation" " generated for them.", &QEDRadiationHandler::_decayingParticles, -1, long(24), 0, 0, false, false, Interface::nolimits); static ParVector interfaceDecayProducts ("DecayProducts", "List of PDG codes of the particles which should be present" " as decay products for the radiation to be generated.", &QEDRadiationHandler::_decayProducts, -1, long(11), 0, 0, false, false, Interface::nolimits); } 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,556 @@ // -*- C++ -*- // // ModelGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the 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()); } } } if( parent->massGenerator() ) { Energy minMass = minimumMass(parent); Energy offShellNess = howOffShell_*parent->width(); if(minMass>parent->mass()-offShellNess) { offShellNess = parent->mass()-minMass; } parent->widthCut(offShellNess); 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) { + 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/PrototypeVertex.h b/Models/General/PrototypeVertex.h --- a/Models/General/PrototypeVertex.h +++ b/Models/General/PrototypeVertex.h @@ -1,477 +1,477 @@ // -*- C++ -*- // // PrototypeVertex.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_PrototypeVertex_H #define HERWIG_PrototypeVertex_H #include #include "ThePEG/Helicity/Vertex/VertexBase.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/EnumIO.h" #include "NBodyDecayConstructorBase.fh" // // This is the declaration of the PrototypeVertex class. // namespace Herwig { using namespace ThePEG; using Helicity::VertexBasePtr; class PrototypeVertex; ThePEG_DECLARE_POINTERS(Herwig::PrototypeVertex,PrototypeVertexPtr); /** Pair of int,double */ typedef pair CFPair; /** * A struct to order the particles in the same way as in the DecayModes */ struct ParticleOrdering { /** * Operator for the ordering * @param p1 The first ParticleData object * @param p2 The second ParticleData object */ - bool operator() (tcPDPtr p1, tcPDPtr p2) { + 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() ); } }; /** * A struct to order the particles in the same way as in the DecayMode's */ struct VertexOrdering { /** * Operator for the ordering * @param p1 The first ParticleData object * @param p2 The second ParticleData object */ bool operator() (const pair< tPDPtr, PrototypeVertexPtr > & p1, const pair< tPDPtr, PrototypeVertexPtr > & p2) const { return abs(p1.first->id()) > abs(p2.first->id()) || ( abs(p1.first->id()) == abs(p2.first->id()) && p1.first->id() > p2.first->id() ) || ( p1.first->id() == p2.first->id() && p1.first->fullName() > p2.first->fullName() ); } }; typedef multiset,VertexOrdering > OrderedVertices; /** * A set of ParticleData objects ordered as for the DecayMode's */ typedef multiset OrderedParticles; /** * Storage of a potenital n-body decay */ class PrototypeVertex : public Base { public: /** * Default Constructor */ PrototypeVertex() : npart(0), possibleOnShell(false) {} /** * Constructor */ PrototypeVertex(tPDPtr in, OrderedVertices out, VertexBasePtr v, int n) : incoming(in), outgoing(out), vertex(v), npart(n), possibleOnShell(false) {} /** * Incoming particle */ tPDPtr incoming; /** * Outgoing particles */ OrderedVertices outgoing; /** * The vertex for the interaction */ VertexBasePtr vertex; /** * The parent of the vertex */ tPrototypeVertexPtr parent; /** * Number of particles */ unsigned int npart; /** * Outgoing particles */ mutable OrderedParticles outPart; /** * Can have on-shell intermediates */ bool possibleOnShell; /** * Increment the number of particles */ void incrementN(int in) { npart += in; if(parent) parent->incrementN(in); } /** * Mass of the incoming particle */ Energy incomingMass() { return incoming->mass(); } /** * Total mass of all the outgoing particles */ Energy outgoingMass() { Energy mass(ZERO); for(OrderedVertices::const_iterator it = outgoing.begin(); it!=outgoing.end();++it) { mass += it->second ? it->second->outgoingMass() : it->first->mass(); } return mass; } /** * Total constituent mass of all the outgoing particles */ Energy outgoingConstituentMass() { Energy mass(ZERO); for(OrderedVertices::const_iterator it = outgoing.begin(); it!=outgoing.end();++it) { mass += it->second ? it->second->outgoingConstituentMass() : it->first->constituentMass(); } return mass; } /** * Check the external particles */ bool checkExternal(bool first=true) { if(outPart.empty()) setOutgoing(); if(first&&outPart.find(incoming)!=outPart.end()) return false; bool output = true; for(OrderedVertices::const_iterator it = outgoing.begin(); it!=outgoing.end();++it) { if(it->second&& !it->second->checkExternal(false)) output = false; } return output; } /** * Set the outgoing particles */ void setOutgoing() const { assert(outPart.empty()); for(OrderedVertices::const_iterator it = outgoing.begin(); it!=outgoing.end();++it) { if(it->second) { it->second->setOutgoing(); outPart.insert(it->second->outPart.begin(), it->second->outPart.end()); } else outPart.insert(it->first); } } /** * Are there potential on-shell intermediates? */ bool canBeOnShell(unsigned int opt,Energy maxMass,bool first); /** * Check if same external particles */ bool sameDecay(const PrototypeVertex & x) const; /** * Create a \f$1\to2\f$ prototype */ static void createPrototypes(tPDPtr inpart, VertexBasePtr vertex, std::stack & prototypes, NBodyDecayConstructorBasePtr decayCon); /** * Expand the prototypes by adding more legs */ static void expandPrototypes(PrototypeVertexPtr proto, VertexBasePtr vertex, std::stack & prototypes, const set & excluded, NBodyDecayConstructorBasePtr decayCon); /** * Copy the whole structure with a new branching */ static PrototypeVertexPtr replicateTree(PrototypeVertexPtr parent, PrototypeVertexPtr oldChild, PrototypeVertexPtr & newChild); }; /** * Output to a stream */ inline ostream & operator<<(ostream & os, const PrototypeVertex & diag) { os << diag.incoming->PDGName() << " -> "; bool seq=false; for(OrderedVertices::const_iterator it = diag.outgoing.begin(); it!=diag.outgoing.end();++it) { os << it->first->PDGName() << " "; if(it->second) seq = true; } os << " decays via " << diag.vertex->fullName() << " in a " << diag.npart << "-body decay\n"; if(!seq) return os; os << "Followed by\n"; for(OrderedVertices::const_iterator it = diag.outgoing.begin(); it!=diag.outgoing.end();++it) { if(it->second) os << *it->second; } return os; } /** * Test whether two diagrams are identical. */ inline bool operator==(const PrototypeVertex & x, const PrototypeVertex & y) { if(x.incoming != y.incoming) return false; if(x.vertex != y.vertex) return false; if(x.npart != y.npart) return false; if(x.outgoing.empty()&&y.outgoing.empty()) return true; if(x.outgoing.size() != y.outgoing.size()) return false; OrderedVertices::const_iterator xt = x.outgoing.begin(); OrderedVertices::const_iterator yt = y.outgoing.begin(); while(xt!=x.outgoing.end()) { if(xt->first != yt->first) return false; // special for identical particles OrderedVertices::const_iterator lxt = x.outgoing.lower_bound(*xt); OrderedVertices::const_iterator uxt = x.outgoing.upper_bound(*xt); --uxt; // just one particle if(lxt==uxt) { if(xt->second && yt->second) { if(*(xt->second)==*(yt->second)) { ++xt; ++yt; continue; } else return false; } else if(xt->second || yt->second) return false; ++xt; ++yt; } // identical particles else { ++uxt; OrderedVertices::const_iterator lyt = y.outgoing.lower_bound(*xt); OrderedVertices::const_iterator uyt = y.outgoing.upper_bound(*xt); unsigned int nx=0; for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {++nx;} unsigned int ny=0; for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {++ny;} if(nx!=ny) return false; vector matched(ny,false); for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) { bool found = false; unsigned int iy=0; for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) { if(matched[iy]) { ++iy; continue; } if( (!ixt->second &&!iyt->second) || ( ixt->second&&iyt->second && *(ixt->second)==*(iyt->second)) ) { matched[iy] = true; found = true; break; } ++iy; } if(!found) return false; } xt=uxt; yt=uyt; } } return true; } /** * A simple vertex for the N-body diagram */ struct NBVertex { /** * Constructor taking a prototype vertex as the arguments */ NBVertex(PrototypeVertexPtr proto = PrototypeVertexPtr() ); /** * Incoming particle */ tPDPtr incoming; /** * Outgoing particles */ mutable OrderedParticles outgoing; /** * The vertices */ list > vertices; /** * The vertex */ VertexBasePtr vertex; }; /** * The NBDiagram struct contains information about a \f$1\to n\f$ decay * that has been automatically generated. */ struct NBDiagram : public NBVertex { /** * Constructor taking a prototype vertex as the arguments*/ NBDiagram(PrototypeVertexPtr proto=PrototypeVertexPtr()); /** * The type of channel */ vector channelType; /** Store colour flow at \f$N_c=3\f$ information */ mutable vector colourFlow; /** Store colour flow at \f$N_c=\infty\f$ information */ mutable vector largeNcColourFlow; }; /** * Output operator to allow the structure to be persistently written * @param os The output stream * @param x The NBVertex */ inline PersistentOStream & operator<<(PersistentOStream & os, const NBVertex & x) { os << x.incoming << x.outgoing << x.vertices << x.vertex; return os; } /** * Input operator to allow persistently written data to be read in * @param is The input stream * @param x The NBVertex */ inline PersistentIStream & operator>>(PersistentIStream & is, NBVertex & x) { is >> x.incoming >> x.outgoing >> x.vertices >> x.vertex; return is; } /** * Output operator to allow the structure to be persistently written * @param os The output stream * @param x The NBDiagram */ inline PersistentOStream & operator<<(PersistentOStream & os, const NBDiagram & x) { os << x.incoming << x.channelType << x.outgoing << x.vertices << x.vertex << x.colourFlow << x.largeNcColourFlow; return os; } /** * Input operator to allow persistently written data to be read in * @param is The input stream * @param x The NBDiagram */ inline PersistentIStream & operator>>(PersistentIStream & is, NBDiagram & x) { is >> x.incoming >> x.channelType >> x.outgoing >> x.vertices >> x.vertex >> x.colourFlow >> x.largeNcColourFlow; return is; } /** * Output a NBVertex to a stream */ inline ostream & operator<<(ostream & os, const NBVertex & vertex) { os << vertex.incoming->PDGName() << " -> "; bool seq=false; for(list >::const_iterator it=vertex.vertices.begin(); it!=vertex.vertices.end();++it) { os << it->first->PDGName() << " "; if(it->second.incoming) seq = true; } os << "via vertex " << vertex.vertex->fullName() << "\n"; if(!seq) return os; os << "Followed by\n"; for(list >::const_iterator it=vertex.vertices.begin(); it!=vertex.vertices.end();++it) { if(it->second.incoming) os << it->second; } return os; } /** * Output a NBDiagram to a stream */ inline ostream & operator<<(ostream & os, const NBDiagram & diag) { os << diag.incoming->PDGName() << " -> "; for(OrderedParticles::const_iterator it=diag.outgoing.begin(); it!=diag.outgoing.end();++it) { os << (**it).PDGName() << " "; } os << " has order "; for(unsigned int ix=0;ixPDGName() << " -> "; bool seq=false; for(list >::const_iterator it=diag.vertices.begin(); it!=diag.vertices.end();++it) { os << it->first->PDGName() << " "; if(it->second.incoming) seq = true; } os << "via vertex " << diag.vertex->fullName() << "\n"; if(!seq) return os; os << "Followed by\n"; for(list >::const_iterator it=diag.vertices.begin(); it!=diag.vertices.end();++it) { if(it->second.incoming) os << it->second; } return os; } } #endif /* HERWIG_PrototypeVertex_H */ diff --git a/Shower/ShowerHandler.cc b/Shower/ShowerHandler.cc --- a/Shower/ShowerHandler.cc +++ b/Shower/ShowerHandler.cc @@ -1,1121 +1,1121 @@ // -*- C++ -*- // // ShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ShowerHandler class. // #include "ShowerHandler.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/PDF/PartonExtractor.h" #include "ThePEG/PDF/PartonBinInstance.h" #include "Herwig/PDT/StandardMatchers.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/Utilities/StringUtils.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Repository/EventGenerator.h" #include "Herwig/Utilities/EnumParticles.h" #include "Herwig/PDF/MPIPDF.h" #include "Herwig/PDF/MinBiasPDF.h" #include "ThePEG/Handlers/EventHandler.h" #include "Herwig/Shower/Core/Base/ShowerTree.h" #include "Herwig/PDF/HwRemDecayer.h" #include #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Decay/DecayIntegrator.h" #include "Herwig/Decay/DecayPhaseSpaceMode.h" using namespace Herwig; DescribeClass describeShowerHandler ("Herwig::ShowerHandler","HwShower.so"); ShowerHandler::~ShowerHandler() {} tShowerHandlerPtr ShowerHandler::currentHandler_ = tShowerHandlerPtr(); void ShowerHandler::doinit() { CascadeHandler::doinit(); // copy particles to decay before showering from input vector to the // set used in the simulation if ( particlesDecayInShower_.empty() ) { for(unsigned int ix=0;ixunrestrictedPhasespace() && restrictPhasespace() ) { generator()->log() << "ShowerApproximation warning: The scale profile chosen requires an unrestricted phase space,\n" << "however, the phase space was set to be restricted. Will switch to unrestricted phase space.\n" << flush; restrictPhasespace_ = false; } } } IBPtr ShowerHandler::clone() const { return new_ptr(*this); } IBPtr ShowerHandler::fullclone() const { return new_ptr(*this); } ShowerHandler::ShowerHandler() : maxtry_(10),maxtryMPI_(10),maxtryDP_(10),maxtryDecay_(100), factorizationScaleFactor_(1.0), renormalizationScaleFactor_(1.0), hardScaleFactor_(1.0), restrictPhasespace_(true), maxPtIsMuF_(false), pdfFreezingScale_(2.5*GeV), doFSR_(true), doISR_(true), splitHardProcess_(true), includeSpaceTime_(false), vMin_(0.1*GeV2), reweight_(1.0) { inputparticlesDecayInShower_.push_back( 6 ); // top inputparticlesDecayInShower_.push_back( 23 ); // Z0 inputparticlesDecayInShower_.push_back( 24 ); // W+/- inputparticlesDecayInShower_.push_back( 25 ); // h0 } void ShowerHandler::doinitrun(){ CascadeHandler::doinitrun(); //can't use isMPIOn here, because the EventHandler is not set at that stage if(MPIHandler_) { MPIHandler_->initialize(); if(MPIHandler_->softInt()) remDec_->initSoftInteractions(MPIHandler_->Ptmin(), MPIHandler_->beta()); } ShowerTree::_vmin2 = vMin_; ShowerTree::_spaceTime = includeSpaceTime_; } void ShowerHandler::dofinish() { CascadeHandler::dofinish(); if(MPIHandler_) MPIHandler_->finalize(); } void ShowerHandler::persistentOutput(PersistentOStream & os) const { os << remDec_ << ounit(pdfFreezingScale_,GeV) << maxtry_ << maxtryMPI_ << maxtryDP_ << maxtryDecay_ << inputparticlesDecayInShower_ << particlesDecayInShower_ << MPIHandler_ << PDFA_ << PDFB_ << PDFARemnant_ << PDFBRemnant_ << includeSpaceTime_ << ounit(vMin_,GeV2) << factorizationScaleFactor_ << renormalizationScaleFactor_ << hardScaleFactor_ << restrictPhasespace_ << maxPtIsMuF_ << hardScaleProfile_ << showerVariations_ << doFSR_ << doISR_ << splitHardProcess_ << useConstituentMasses_; } void ShowerHandler::persistentInput(PersistentIStream & is, int) { is >> remDec_ >> iunit(pdfFreezingScale_,GeV) >> maxtry_ >> maxtryMPI_ >> maxtryDP_ >> maxtryDecay_ >> inputparticlesDecayInShower_ >> particlesDecayInShower_ >> MPIHandler_ >> PDFA_ >> PDFB_ >> PDFARemnant_ >> PDFBRemnant_ >> includeSpaceTime_ >> iunit(vMin_,GeV2) >> factorizationScaleFactor_ >> renormalizationScaleFactor_ >> hardScaleFactor_ >> restrictPhasespace_ >> maxPtIsMuF_ >> hardScaleProfile_ >> showerVariations_ >> doFSR_ >> doISR_ >> splitHardProcess_ >> useConstituentMasses_; } void ShowerHandler::Init() { static ClassDocumentation documentation ("Main driver class for the showering."); static Reference interfaceRemDecayer("RemDecayer", "A reference to the Remnant Decayer object", &Herwig::ShowerHandler::remDec_, false, false, true, false); static Parameter interfacePDFFreezingScale ("PDFFreezingScale", "The PDF freezing scale", &ShowerHandler::pdfFreezingScale_, GeV, 2.5*GeV, 2.0*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfaceMaxTry ("MaxTry", "The maximum number of attempts for the main showering loop", &ShowerHandler::maxtry_, 10, 1, 100, false, false, Interface::limited); static Parameter interfaceMaxTryMPI ("MaxTryMPI", "The maximum number of regeneration attempts for an additional scattering", &ShowerHandler::maxtryMPI_, 10, 0, 100, false, false, Interface::limited); static Parameter interfaceMaxTryDP ("MaxTryDP", "The maximum number of regeneration attempts for an additional hard scattering", &ShowerHandler::maxtryDP_, 10, 0, 100, false, false, Interface::limited); static ParVector interfaceDecayInShower ("DecayInShower", "PDG codes of the particles to be decayed in the shower", &ShowerHandler::inputparticlesDecayInShower_, -1, 0l, -10000000l, 10000000l, false, false, Interface::limited); static Reference interfaceMPIHandler ("MPIHandler", "The object that administers all additional scatterings.", &ShowerHandler::MPIHandler_, false, false, true, true); static Reference interfacePDFA ("PDFA", "The PDF for beam particle A. Overrides the particle's own PDF setting." "By default used for both the shower and forced splitting in the remnant", &ShowerHandler::PDFA_, false, false, true, true, false); static Reference interfacePDFB ("PDFB", "The PDF for beam particle B. Overrides the particle's own PDF setting." "By default used for both the shower and forced splitting in the remnant", &ShowerHandler::PDFB_, false, false, true, true, false); static Reference interfacePDFARemnant ("PDFARemnant", "The PDF for beam particle A used to generate forced splittings of the remnant." " This overrides both the particle's own PDF setting and the value set by PDFA if used.", &ShowerHandler::PDFARemnant_, false, false, true, true, false); static Reference interfacePDFBRemnant ("PDFBRemnant", "The PDF for beam particle B used to generate forced splittings of the remnant." " This overrides both the particle's own PDF setting and the value set by PDFB if used.", &ShowerHandler::PDFBRemnant_, false, false, true, true, false); static Switch interfaceIncludeSpaceTime ("IncludeSpaceTime", "Whether to include the model for the calculation of space-time distances", &ShowerHandler::includeSpaceTime_, false, false, false); static SwitchOption interfaceIncludeSpaceTimeYes (interfaceIncludeSpaceTime, "Yes", "Include the model", true); static SwitchOption interfaceIncludeSpaceTimeNo (interfaceIncludeSpaceTime, "No", "Only include the displacement from the particle-s lifetime for decaying particles", false); static Parameter interfaceMinimumVirtuality ("MinimumVirtuality", "The minimum virtuality for the space-time model", &ShowerHandler::vMin_, GeV2, 0.1*GeV2, 0.0*GeV2, 1000.0*GeV2, false, false, Interface::limited); static Parameter interfaceFactorizationScaleFactor ("FactorizationScaleFactor", "The factorization scale factor.", &ShowerHandler::factorizationScaleFactor_, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceRenormalizationScaleFactor ("RenormalizationScaleFactor", "The renormalization scale factor.", &ShowerHandler::renormalizationScaleFactor_, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceHardScaleFactor ("HardScaleFactor", "The hard scale factor.", &ShowerHandler::hardScaleFactor_, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceMaxTryDecay ("MaxTryDecay", "The maximum number of attempts to generate a decay", &ShowerHandler::maxtryDecay_, 200, 10, 0, false, false, Interface::lowerlim); static Reference interfaceHardScaleProfile ("HardScaleProfile", "The hard scale profile to use.", &ShowerHandler::hardScaleProfile_, false, false, true, true, false); static Switch interfaceMaxPtIsMuF ("MaxPtIsMuF", "", &ShowerHandler::maxPtIsMuF_, false, false, false); static SwitchOption interfaceMaxPtIsMuFYes (interfaceMaxPtIsMuF, "Yes", "", true); static SwitchOption interfaceMaxPtIsMuFNo (interfaceMaxPtIsMuF, "No", "", false); static Switch interfaceRestrictPhasespace ("RestrictPhasespace", "Switch on or off phasespace restrictions", &ShowerHandler::restrictPhasespace_, true, false, false); static SwitchOption interfaceRestrictPhasespaceYes (interfaceRestrictPhasespace, "Yes", "Perform phasespace restrictions", true); static SwitchOption interfaceRestrictPhasespaceNo (interfaceRestrictPhasespace, "No", "Do not perform phasespace restrictions", false); static Command interfaceAddVariation ("AddVariation", "Add a shower variation.", &ShowerHandler::doAddVariation, false); static Switch interfaceDoFSR ("DoFSR", "Switch on or off final state radiation.", &ShowerHandler::doFSR_, true, false, false); static SwitchOption interfaceDoFSRYes (interfaceDoFSR, "Yes", "Switch on final state radiation.", true); static SwitchOption interfaceDoFSRNo (interfaceDoFSR, "No", "Switch off final state radiation.", false); static Switch interfaceDoISR ("DoISR", "Switch on or off initial state radiation.", &ShowerHandler::doISR_, true, false, false); static SwitchOption interfaceDoISRYes (interfaceDoISR, "Yes", "Switch on initial state radiation.", true); static SwitchOption interfaceDoISRNo (interfaceDoISR, "No", "Switch off initial state radiation.", false); static Switch interfaceSplitHardProcess ("SplitHardProcess", "Whether or not to try and split the hard process into production and decay processes", &ShowerHandler::splitHardProcess_, true, false, false); static SwitchOption interfaceSplitHardProcessYes (interfaceSplitHardProcess, "Yes", "Split the hard process", true); static SwitchOption interfaceSplitHardProcessNo (interfaceSplitHardProcess, "No", "Don't split the hard process", false); static Switch interfaceUseConstituentMasses ("UseConstituentMasses", "Whether or not to use constituent masses for the reconstruction of the particle after showering.", &ShowerHandler::useConstituentMasses_, true, false, false); static SwitchOption interfaceUseConstituentMassesYes (interfaceUseConstituentMasses, "Yes", "Use constituent masses.", true); static SwitchOption interfaceUseConstituentMassesNo (interfaceUseConstituentMasses, "No", "Don't use constituent masses.", false); } Energy ShowerHandler::hardScale() const { assert(false); } void ShowerHandler::cascade() { useMe(); // Initialise the weights in the event object // so that any variations are output regardless of // whether showering occurs for the given event initializeWeights(); // get the PDF's from ThePEG (if locally overridden use the local versions) tcPDFPtr first = PDFA_ ? tcPDFPtr(PDFA_) : firstPDF().pdf(); tcPDFPtr second = PDFB_ ? tcPDFPtr(PDFB_) : secondPDF().pdf(); resetPDFs(make_pair(first,second)); // set the PDFs for the remnant if( ! rempdfs_.first) rempdfs_.first = PDFARemnant_ ? PDFPtr(PDFARemnant_) : const_ptr_cast(first); if( ! rempdfs_.second) rempdfs_.second = PDFBRemnant_ ? PDFPtr(PDFBRemnant_) : const_ptr_cast(second); // get the incoming partons tPPair incomingPartons = eventHandler()->currentCollision()->primarySubProcess()->incoming(); // and the parton bins PBIPair incomingBins = make_pair(lastExtractor()->partonBinInstance(incomingPartons.first), lastExtractor()->partonBinInstance(incomingPartons.second)); // and the incoming hadrons tPPair incomingHadrons = eventHandler()->currentCollision()->incoming(); remnantDecayer()->setHadronContent(incomingHadrons); // check if incoming hadron == incoming parton // and get the incoming hadron if exists or parton otherwise incoming_ = make_pair(incomingBins.first ? incomingBins.first ->particle() : incomingPartons.first, incomingBins.second ? incomingBins.second->particle() : incomingPartons.second); // check the collision is of the beam particles // and if not boost collision to the right frame // i.e. the hadron-hadron CMF of the collision bool btotal(false); LorentzRotation rtotal; if(incoming_.first != incomingHadrons.first || incoming_.second != incomingHadrons.second ) { btotal = true; boostCollision(false); } // set the current ShowerHandler setCurrentHandler(); // first shower the hard process try { SubProPtr sub = eventHandler()->currentCollision()->primarySubProcess(); incomingPartons = cascade(sub,lastXCombPtr()); } catch(ShowerTriesVeto &veto){ throw Exception() << "Failed to generate the shower after " << veto.tries << " attempts in ShowerHandler::cascade()" << Exception::eventerror; } if(showerHardProcessVeto()) throw Veto(); // if a non-hadron collision return (both incoming non-hadronic) if( ( !incomingBins.first|| !isResolvedHadron(incomingBins.first ->particle()))&& ( !incomingBins.second|| !isResolvedHadron(incomingBins.second->particle()))) { // boost back to lab if needed if(btotal) boostCollision(true); // perform the reweighting for the hard process shower combineWeights(); // unset the current ShowerHandler unSetCurrentHandler(); return; } // get the remnants for hadronic collision pair remnants(getRemnants(incomingBins)); // set the starting scale of the forced splitting to the PDF freezing scale remnantDecayer()->initialize(remnants, incoming_, *currentStep(), pdfFreezingScale()); // do the first forcedSplitting try { remnantDecayer()->doSplit(incomingPartons, make_pair(rempdfs_.first,rempdfs_.second), true); } catch (ExtraScatterVeto) { throw Exception() << "Remnant extraction failed in " << "ShowerHandler::cascade() from primary interaction" << Exception::eventerror; } // perform the reweighting for the hard process shower combineWeights(); // if no MPI return if( !isMPIOn() ) { remnantDecayer()->finalize(); // boost back to lab if needed if(btotal) boostCollision(true); // unset the current ShowerHandler unSetCurrentHandler(); return; } // generate the multiple scatters use modified pdf's now: setMPIPDFs(); // additional "hard" processes unsigned int tries(0); // This is the loop over additional hard scatters (most of the time // only one, but who knows...) for(unsigned int i=1; i <= getMPIHandler()->additionalHardProcs(); i++){ //counter for regeneration unsigned int multSecond = 0; // generate the additional scatters while( multSecond < getMPIHandler()->multiplicity(i) ) { // generate the hard scatter tStdXCombPtr lastXC = getMPIHandler()->generate(i); SubProPtr sub = lastXC->construct(); // add to the Step newStep()->addSubProcess(sub); // increment the counters tries++; multSecond++; if(tries == maxtryDP_) throw Exception() << "Failed to establish the requested number " << "of additional hard processes. If this error " << "occurs often, your selection of additional " << "scatter is probably unphysical" << Exception::eventerror; // Generate the shower. If not possible veto the event try { incomingPartons = cascade(sub,lastXC); } catch(ShowerTriesVeto &veto){ throw Exception() << "Failed to generate the shower of " << "a secondary hard process after " << veto.tries << " attempts in Evolver::showerHardProcess()" << Exception::eventerror; } try { // do the forcedSplitting remnantDecayer()->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false); } catch(ExtraScatterVeto){ //remove all particles associated with the subprocess newStep()->removeParticle(incomingPartons.first); newStep()->removeParticle(incomingPartons.second); //remove the subprocess from the list newStep()->removeSubProcess(sub); //regenerate the scattering multSecond--; continue; } // connect with the remnants but don't set Remnant colour, // because that causes problems due to the multiple colour lines. if ( !remnants.first ->extract(incomingPartons.first , false) || !remnants.second->extract(incomingPartons.second, false) ) throw Exception() << "Remnant extraction failed in " << "ShowerHandler::cascade() for additional scatter" << Exception::runerror; } // perform the reweighting for the additional hard scatter shower combineWeights(); } // the underlying event processes unsigned int ptveto(1), veto(0); unsigned int max(getMPIHandler()->multiplicity()); for(unsigned int i=0; i maxtryMPI_) break; //generate PSpoint tStdXCombPtr lastXC = getMPIHandler()->generate(); SubProPtr sub = lastXC->construct(); //If Algorithm=1 additional scatters of the signal type // with pt > ptmin have to be vetoed //with probability 1/(m+1), where m is the number of occurances in this event if( getMPIHandler()->Algorithm() == 1 ){ //get the pT Energy pt = sub->outgoing().front()->momentum().perp(); if(pt > getMPIHandler()->PtForVeto() && UseRandom::rnd() < 1./(ptveto+1) ){ ptveto++; i--; continue; } } // add to the SubProcess to the step newStep()->addSubProcess(sub); // Run the Shower. If not possible veto the scattering try { incomingPartons = cascade(sub,lastXC); } // discard this extra scattering, but try the next one catch(ShowerTriesVeto) { newStep()->removeSubProcess(sub); //regenerate the scattering veto++; i--; continue; } try{ //do the forcedSplitting remnantDecayer()->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false); } catch (ExtraScatterVeto) { //remove all particles associated with the subprocess newStep()->removeParticle(incomingPartons.first); newStep()->removeParticle(incomingPartons.second); //remove the subprocess from the list newStep()->removeSubProcess(sub); //regenerate the scattering veto++; i--; continue; } //connect with the remnants but don't set Remnant colour, //because that causes problems due to the multiple colour lines. if ( !remnants.first ->extract(incomingPartons.first , false) || !remnants.second->extract(incomingPartons.second, false) ) throw Exception() << "Remnant extraction failed in " << "ShowerHandler::cascade() for MPI hard scattering" << Exception::runerror; //reset veto counter veto = 0; // perform the reweighting for the MPI process shower combineWeights(); } // finalize the remnants remnantDecayer()->finalize(getMPIHandler()->colourDisrupt(), getMPIHandler()->softMultiplicity()); // boost back to lab if needed if(btotal) boostCollision(true); // unset the current ShowerHandler unSetCurrentHandler(); getMPIHandler()->clean(); resetPDFs(make_pair(first,second)); } void ShowerHandler::initializeWeights() { if ( !showerVariations().empty() ) { tEventPtr event = eventHandler()->currentEvent(); for ( map::const_iterator var = showerVariations().begin(); var != showerVariations().end(); ++var ) { // Check that this is behaving as intended //map::iterator wi = event->optionalWeights().find(var->first); //assert(wi == event->optionalWeights().end() ); event->optionalWeights()[var->first] = 1.0; currentWeights_[var->first] = 1.0; } } reweight_ = 1.0; } void ShowerHandler::resetWeights() { for ( map::iterator w = currentWeights_.begin(); w != currentWeights_.end(); ++w ) { w->second = 1.0; } reweight_ = 1.0; } void ShowerHandler::combineWeights() { tEventPtr event = eventHandler()->currentEvent(); for ( map::const_iterator w = currentWeights_.begin(); w != currentWeights_.end(); ++w ) { map::iterator ew = event->optionalWeights().find(w->first); if ( ew != event->optionalWeights().end() ) ew->second *= w->second; else { assert(false && "Weight name unknown."); //event->optionalWeights()[w->first] = w->second; } } if ( reweight_ != 1.0 ) { Ptr::tptr eh = dynamic_ptr_cast::tptr>(eventHandler()); if ( !eh ) { throw Exception() << "ShowerHandler::combineWeights() : Cross section reweighting " << "through the shower is currently only available with standard " << "event generators" << Exception::runerror; } eh->reweight(reweight_); } } string ShowerHandler::doAddVariation(string in) { if ( in.empty() ) return "expecting a name and a variation specification"; string name = StringUtils::car(in); ShowerVariation var; string res = var.fromInFile(StringUtils::cdr(in)); if ( res.empty() ) { if ( !var.firstInteraction && !var.secondaryInteractions ) { // TODO what about decay showers? return "variation does not apply to any shower"; } if ( var.renormalizationScaleFactor == 1.0 && var.factorizationScaleFactor == 1.0 ) { return "variation does not vary anything"; } /* Repository::clog() << "adding a variation with tag '" << name << "' using\nxir = " << var.renormalizationScaleFactor << " xif = " << var.factorizationScaleFactor << "\napplying to:\n" << "first interaction = " << var.firstInteraction << " " << "secondary interactions = " << var.secondaryInteractions << "\n" << flush; */ showerVariations()[name] = var; } return res; } tPPair ShowerHandler::cascade(tSubProPtr, XCPtr) { assert(false); } ShowerHandler::RemPair ShowerHandler::getRemnants(PBIPair incomingBins) { RemPair remnants; // first beam particle if(incomingBins.first&&!incomingBins.first->remnants().empty()) { remnants.first = dynamic_ptr_cast(incomingBins.first->remnants()[0] ); if(remnants.first) { ParticleVector children=remnants.first->children(); for(unsigned int ix=0;ixdataPtr()==remnants.first->dataPtr()) remnants.first = dynamic_ptr_cast(children[ix]); } //remove existing colour lines from the remnants if(remnants.first->colourLine()) remnants.first->colourLine()->removeColoured(remnants.first); if(remnants.first->antiColourLine()) remnants.first->antiColourLine()->removeAntiColoured(remnants.first); } } // seconnd beam particle if(incomingBins.second&&!incomingBins. second->remnants().empty()) { remnants.second = dynamic_ptr_cast(incomingBins.second->remnants()[0] ); if(remnants.second) { ParticleVector children=remnants.second->children(); for(unsigned int ix=0;ixdataPtr()==remnants.second->dataPtr()) remnants.second = dynamic_ptr_cast(children[ix]); } //remove existing colour lines from the remnants if(remnants.second->colourLine()) remnants.second->colourLine()->removeColoured(remnants.second); if(remnants.second->antiColourLine()) remnants.second->antiColourLine()->removeAntiColoured(remnants.second); } } assert(remnants.first || remnants.second); return remnants; } namespace { void addChildren(tPPtr in,set & particles) { particles.insert(in); for(unsigned int ix=0;ixchildren().size();++ix) addChildren(in->children()[ix],particles); } } void ShowerHandler::boostCollision(bool boost) { // calculate boost from lab to rest if(!boost) { Lorentz5Momentum ptotal=incoming_.first ->momentum()+incoming_.second->momentum(); boost_ = LorentzRotation(-ptotal.boostVector()); Axis axis((boost_*incoming_.first ->momentum()).vect().unit()); if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); boost_.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } } // first call performs the boost and second inverse // get the particles to be boosted set particles; addChildren(incoming_.first,particles); addChildren(incoming_.second,particles); // apply the boost for(set::const_iterator cit=particles.begin(); cit!=particles.end();++cit) { (*cit)->transform(boost_); } if(!boost) boost_.invert(); } void ShowerHandler::setMPIPDFs() { if ( !mpipdfs_.first ) { // first have to check for MinBiasPDF tcMinBiasPDFPtr first = dynamic_ptr_cast(firstPDF().pdf()); if(first) mpipdfs_.first = new_ptr(MPIPDF(first->originalPDF())); else mpipdfs_.first = new_ptr(MPIPDF(firstPDF().pdf())); } if ( !mpipdfs_.second ) { tcMinBiasPDFPtr second = dynamic_ptr_cast(secondPDF().pdf()); if(second) mpipdfs_.second = new_ptr(MPIPDF(second->originalPDF())); else mpipdfs_.second = new_ptr(MPIPDF(secondPDF().pdf())); } if( !remmpipdfs_.first ) { tcMinBiasPDFPtr first = dynamic_ptr_cast(rempdfs_.first); if(first) remmpipdfs_.first = new_ptr(MPIPDF(first->originalPDF())); else remmpipdfs_.first = new_ptr(MPIPDF(rempdfs_.first)); } if( !remmpipdfs_.second ) { tcMinBiasPDFPtr second = dynamic_ptr_cast(rempdfs_.second); if(second) remmpipdfs_.second = new_ptr(MPIPDF(second->originalPDF())); else remmpipdfs_.second = new_ptr(MPIPDF(rempdfs_.second)); } // reset the PDFs stored in the base class resetPDFs(mpipdfs_); } bool ShowerHandler::isResolvedHadron(tPPtr particle) { if(!HadronMatcher::Check(particle->data())) return false; for(unsigned int ix=0;ixchildren().size();++ix) { if(particle->children()[ix]->id()==ParticleID::Remnant) return true; } return false; } namespace { bool decayProduct(tSubProPtr subProcess, tPPtr particle) { // must be time-like and not incoming if(particle->momentum().m2()<=ZERO|| particle == subProcess->incoming().first|| particle == subProcess->incoming().second) return false; // if only 1 outgoing and this is it if(subProcess->outgoing().size()==1 && subProcess->outgoing()[0]==particle) return true; // must not be the s-channel intermediate otherwise if(find(subProcess->incoming().first->children().begin(), subProcess->incoming().first->children().end(),particle)!= subProcess->incoming().first->children().end()&& find(subProcess->incoming().second->children().begin(), subProcess->incoming().second->children().end(),particle)!= subProcess->incoming().second->children().end()&& subProcess->incoming().first ->children().size()==1&& subProcess->incoming().second->children().size()==1) return false; // if non-coloured this is enough if(!particle->dataPtr()->coloured()) return true; // if coloured must be unstable if(particle->dataPtr()->stable()) return false; // must not have same particle type as a child int id = particle->id(); for(unsigned int ix=0;ixchildren().size();++ix) if(particle->children()[ix]->id()==id) return false; // otherwise its a decaying particle return true; } PPtr findParent(PPtr original, bool & isHard, set outgoingset, tSubProPtr subProcess) { PPtr parent=original; isHard |=(outgoingset.find(original) != outgoingset.end()); if(!original->parents().empty()) { PPtr orig=original->parents()[0]; if(decayProduct(subProcess,orig)) parent=findParent(orig,isHard,outgoingset,subProcess); } return parent; } } void ShowerHandler::findDecayProducts(PPtr in,PerturbativeProcessPtr hard, DecayProcessMap & decay) const { ParticleVector children=in->children(); for(ParticleVector::const_iterator it=children.begin(); it!=children.end();++it) { // if decayed or should be decayed in shower make the PerturbaitveProcess bool radiates = false; if(!(**it).children().empty()) { // remove d,u,s,c,b quarks and leptons other than on-shell taus if( StandardQCDPartonMatcher::Check((**it).id()) || ( LeptonMatcher::Check((**it).id()) && !(abs((**it).id())==ParticleID::tauminus && abs((**it).mass()-(**it).dataPtr()->mass())id()==(**it).id()) { foundParticle = true; } else if((**it).children()[iy]->id()==ParticleID::g || (**it).children()[iy]->id()==ParticleID::gamma) { foundGauge = true; } } radiates = foundParticle && foundGauge; } } if(radiates) { findDecayProducts(*it,hard,decay); } else if(!(**it).children().empty()|| (decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) { createDecayProcess(*it,hard,decay); } else { hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr())); } } } void ShowerHandler::splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard, DecayProcessMap & decay) const { // temporary storage of the particles set hardParticles; // tagged particles in a set set outgoingset(tagged.begin(),tagged.end()); bool isHard=false; // loop over the tagged particles for (tParticleVector::const_iterator taggedP = tagged.begin(); taggedP != tagged.end(); ++taggedP) { // skip remnants if (eventHandler()->currentCollision()&& eventHandler()->currentCollision()->isRemnant(*taggedP)) continue; // find the parent and whether its a decaying particle bool isDecayProd=false; // check if hard isHard |=(outgoingset.find(*taggedP) != outgoingset.end()); if(splitHardProcess_) { tPPtr parent = *taggedP; // check if from s channel decaying colourless particle while(parent&&!parent->parents().empty()&&!isDecayProd) { parent = parent->parents()[0]; if(parent == subProcess_->incoming().first || parent == subProcess_->incoming().second ) break; isDecayProd = decayProduct(subProcess_,parent); } if (isDecayProd) hardParticles.insert(findParent(parent,isHard,outgoingset,subProcess_)); } if (!isDecayProd) hardParticles.insert(*taggedP); } // there must be something to shower if(hardParticles.empty()) throw Exception() << "No particles to shower in " << "ShowerHandler::splitHardProcess()" << Exception::eventerror; // must be a hard process if(!isHard) throw Exception() << "Starting on decay not yet implemented in " << "ShowerHandler::splitHardProcess()" << Exception::runerror; // create the hard process hard = new_ptr(PerturbativeProcess()); // incoming particles hard->incoming().push_back(make_pair(subProcess_->incoming().first ,PerturbativeProcessPtr())); hard->incoming().push_back(make_pair(subProcess_->incoming().second,PerturbativeProcessPtr())); // outgoing particles for(set::const_iterator it=hardParticles.begin();it!=hardParticles.end();++it) { // if decayed or should be decayed in shower make the tree PPtr orig = *it; bool radiates = false; if(!orig->children().empty()) { // remove d,u,s,c,b quarks and leptons other than on-shell taus if( StandardQCDPartonMatcher::Check(orig->id()) || ( LeptonMatcher::Check(orig->id()) && !(abs(orig->id())==ParticleID::tauminus && abs(orig->mass()-orig->dataPtr()->mass())children().size();++iy) { if(orig->children()[iy]->id()==orig->id()) { foundParticle = true; } else if(orig->children()[iy]->id()==ParticleID::g || orig->children()[iy]->id()==ParticleID::gamma) { foundGauge = true; } } radiates = foundParticle && foundGauge; } } if(radiates) { findDecayProducts(orig,hard,decay); } else if(!(**it).children().empty()|| (decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) { createDecayProcess(*it,hard,decay); } else { hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr())); } } } void ShowerHandler::createDecayProcess(PPtr in,PerturbativeProcessPtr hard, DecayProcessMap & decay) const { // there must be an incoming particle assert(in); // create the new process and connect with the parent PerturbativeProcessPtr newDecay=new_ptr(PerturbativeProcess()); newDecay->incoming().push_back(make_pair(in,hard)); Energy width=in->dataPtr()->generateWidth(in->mass()); decay.insert(make_pair(width,newDecay)); hard->outgoing().push_back(make_pair(in,newDecay)); // we need to deal with the decay products if decayed ParticleVector children = in->children(); if(!children.empty()) { for(ParticleVector::const_iterator it = children.begin(); it!= children.end(); ++it) { // if decayed or should be decayed in shower make the tree in->abandonChild(*it); bool radiates = false; if(!(**it).children().empty()) { if(StandardQCDPartonMatcher::Check((**it).id())|| (LeptonMatcher::Check((**it).id())&& !(abs((**it).id())==ParticleID::tauminus && abs((**it).mass()-(**it).dataPtr()->mass())id()==(**it).id()) { foundParticle = true; } else if((**it).children()[iy]->id()==ParticleID::g || (**it).children()[iy]->id()==ParticleID::gamma) { foundGauge = true; } } radiates = foundParticle && foundGauge; } // finally assume all non-decaying particles are in this class // pr 27/11/15 not sure about this bit // if(!radiates) { // radiates = !decaysInShower((**it).id()); // } } if(radiates) { findDecayProducts(*it,newDecay,decay); } else if(!(**it).children().empty()|| (decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) { createDecayProcess(*it,newDecay,decay); } else { newDecay->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr())); } } } } tDMPtr ShowerHandler::decay(PerturbativeProcessPtr process, DecayProcessMap & decayMap, bool radPhotons ) const { PPtr parent = process->incoming()[0].first; assert(parent); if(parent->spinInfo()) parent->spinInfo()->decay(true); unsigned int ntry = 0; ParticleVector children; tDMPtr dm = DMPtr(); while (true) { // exit if fails if (++ntry>=maxtryDecay_) throw Exception() << "Failed to perform decay in ShowerHandler::decay()" << " after " << maxtryDecay_ << " attempts for " << parent->PDGName() << Exception::eventerror; // select decay mode dm = parent->data().selectMode(*parent); if(!dm) throw Exception() << "Failed to select decay mode in ShowerHandler::decay()" << "for " << parent->PDGName() << Exception::eventerror; if(!dm->decayer()) throw Exception() << "No Decayer for selected decay mode " << " in ShowerHandler::decay()" << Exception::runerror; // start of try block try { children = dm->decayer()->decay(*dm, *parent); // if no children have another go if(children.empty()) continue; if(radPhotons){ // generate radiation in the decay tDecayIntegratorPtr hwdec=dynamic_ptr_cast(dm->decayer()); if (hwdec && hwdec->canGeneratePhotons()) children = hwdec->generatePhotons(*parent,children); } // set up parent parent->decayMode(dm); // add children for (unsigned int i = 0, N = children.size(); i < N; ++i ) { children[i]->setLabVertex(parent->labDecayVertex()); //parent->addChild(children[i]); } // if succeeded break out of loop break; } catch(Veto) { } } assert(!children.empty()); for(ParticleVector::const_iterator it = children.begin(); it!= children.end(); ++it) { if(!(**it).children().empty()|| (decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) { createDecayProcess(*it,process,decayMap); } else { process->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr())); } } return dm; } // Note: The tag must be constructed from an ordered particle container. tDMPtr ShowerHandler::findDecayMode(const string & tag) const { static map cache; map::const_iterator pos = cache.find(tag); if ( pos != cache.end() ) return pos->second; tDMPtr dm = CurrentGenerator::current().findDecayMode(tag); cache[tag] = dm; return dm; } /** * Operator for the particle ordering * @param p1 The first ParticleData object * @param p2 The second ParticleData object */ -bool ShowerHandler::ParticleOrdering::operator() (tcPDPtr p1, tcPDPtr p2) { +bool ShowerHandler::ParticleOrdering::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() ); } diff --git a/Shower/ShowerHandler.h b/Shower/ShowerHandler.h --- a/Shower/ShowerHandler.h +++ b/Shower/ShowerHandler.h @@ -1,848 +1,848 @@ // -*- C++ -*- // // ShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ShowerHandler_H #define HERWIG_ShowerHandler_H // // This is the declaration of the ShowerHandler class. // #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Handlers/CascadeHandler.h" #include "ShowerVariation.h" #include "Herwig/PDF/HwRemDecayer.fh" #include "ThePEG/EventRecord/RemnantParticle.fh" #include "UEBase.h" #include "PerturbativeProcess.h" #include "Herwig/MatrixElement/Matchbox/Matching/HardScaleProfile.h" #include "ShowerHandler.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This class is the main driver of the shower: it is responsible for * the proper handling of all other specific collaborating classes * and for the storing of the produced particles in the event record. * * @see \ref ShowerHandlerInterfaces "The interfaces" * * @see ThePEG::CascadeHandler * @see MPIHandler * @see HwRemDecayer */ class ShowerHandler: public CascadeHandler { public: /** * Typedef for a pair of ThePEG::RemnantParticle pointers. */ typedef pair RemPair; public: /** * Default constructor */ ShowerHandler(); /** * Destructor */ virtual ~ShowerHandler(); public: /** * The main method which manages the multiple interactions and starts * the shower by calling cascade(sub, lastXC). */ virtual void cascade(); /** * pointer to "this", the current ShowerHandler. */ static const tShowerHandlerPtr currentHandler() { assert(currentHandler_); return currentHandler_; } public: /** * Hook to allow vetoing of event after showering hard sub-process * as in e.g. MLM merging. */ virtual bool showerHardProcessVeto() const { return false; } /** * Return true, if this cascade handler will perform reshuffling from hard * process masses. */ virtual bool isReshuffling() const { return true; } /** * Return true, if this cascade handler will put the final state * particles to their constituent mass. If false the nominal mass is used. */ virtual bool retConstituentMasses() const { return useConstituentMasses_; } /** * Return true, if the shower handler can generate a truncated * shower for POWHEG style events generated using Matchbox */ virtual bool canHandleMatchboxTrunc() const { return false; } /** * Get the PDF freezing scale */ Energy pdfFreezingScale() const { return pdfFreezingScale_; } /** * Get the local PDFs. */ PDFPtr getPDFA() const {return PDFA_;} /** * Get the local PDFs. */ PDFPtr getPDFB() const {return PDFB_;} /** * Return true if currently the primary subprocess is showered. */ bool firstInteraction() const { if (!eventHandler()->currentCollision())return true; return ( subProcess_ == eventHandler()->currentCollision()->primarySubProcess() ); } /** * Return the remnant decayer. */ tHwRemDecPtr remnantDecayer() const { return remDec_; } /** * Split the hard process into production and decays * @param tagged The tagged particles from the StepHandler * @param hard The hard perturbative process * @param decay The decay particles */ void splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard, DecayProcessMap & decay) const; /** * Information if the Showerhandler splits the hard process. */ bool doesSplitHardProcess()const {return splitHardProcess_;} /** * Decay a particle. * radPhotons switches the generation of photon * radiation on/off. * Required for Dipole Shower but not QTilde Shower. */ tDMPtr decay(PerturbativeProcessPtr, DecayProcessMap & decay, bool radPhotons = false) const; /** * Cached lookup of decay modes. * Generator::findDecayMode() is not efficient. */ tDMPtr findDecayMode(const string & tag) const; /** * A struct to order the particles in the same way as in the DecayMode's */ struct ParticleOrdering { - bool operator() (tcPDPtr p1, tcPDPtr p2); + bool operator() (tcPDPtr p1, tcPDPtr p2) const; }; /** * A container for ordered particles required * for constructing tags for decay mode lookup. */ typedef multiset OrderedParticles; public: /** * @name Switches for initial- and final-state radiation */ //@{ /** * Switch for any radiation */ bool doRadiation() const {return doFSR_ || doISR_;} /** * Switch on or off final state radiation. */ bool doFSR() const { return doFSR_;} /** * Switch on or off initial state radiation. */ bool doISR() const { return doISR_;} //@} public: /** * @name Switches for scales */ //@{ /** * Return true if maximum pt should be deduced from the factorization scale */ bool hardScaleIsMuF() const { return maxPtIsMuF_; } /** * The factorization scale factor. */ double factorizationScaleFactor() const { return factorizationScaleFactor_; } /** * The renormalization scale factor. */ double renFac() const { return renormalizationScaleFactor_; } /** * The factorization scale factor. */ double facFac() const { return factorizationScaleFactor_; } /** * The renormalization scale factor. */ double renormalizationScaleFactor() const { return renormalizationScaleFactor_; } /** * The scale factor for the hard scale */ double hardScaleFactor() const { return hardScaleFactor_; } /** * Return true, if the phase space restrictions of the dipole shower should * be applied. */ bool restrictPhasespace() const { return restrictPhasespace_; } /** * Return profile scales */ Ptr::tptr profileScales() const { return hardScaleProfile_; } /** * Return the relevant hard scale to be used in the profile scales */ virtual Energy hardScale() const; /** * Return information about shower phase space choices */ virtual int showerPhaseSpaceOption() const { assert(false && "not implemented in general"); return -1; } //@} public: /** * Access the shower variations */ map& showerVariations() { return showerVariations_; } /** * Return the shower variations */ const map& showerVariations() const { return showerVariations_; } /** * Access the current Weights */ map& currentWeights() { return currentWeights_; } /** * Return the current Weights */ const map& currentWeights() const { return currentWeights_; } /** * Change the current reweighting factor */ void reweight(double w) { reweight_ = w; } /** * Return the current reweighting factor */ double reweight() const { return reweight_; } public: /** * struct that is used to catch exceptions which are thrown * due to energy conservation issues of additional scatters */ struct ExtraScatterVeto {}; /** * struct that is used to catch exceptions which are thrown * due to fact that the Shower has been invoked more than * a defined threshold on a certain configuration */ struct ShowerTriesVeto { /** variable to store the number of attempts */ const int tries; /** constructor */ ShowerTriesVeto(int t) : tries(t) {} }; 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 Functions to perform the cascade */ //@{ /** * The main method which manages the showering of a subprocess. */ virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb); /** * Set up for the cascade */ void prepareCascade(tSubProPtr sub) { current_ = currentStep(); subProcess_ = sub; } /** * Boost all the particles in the collision so that the collision always occurs * in the rest frame with the incoming particles along the z axis */ void boostCollision(bool boost); //@} protected: /** * Set/unset the current shower handler */ //@{ /** * Set the current handler */ void setCurrentHandler() { currentHandler_ = tShowerHandlerPtr(this); } /** * Unset the current handler */ void unSetCurrentHandler() { currentHandler_ = tShowerHandlerPtr(); } //@} protected: /** * @name Members relating to the underlying event and MPI */ //@{ /** * Return true if multiple parton interactions are switched on * and can be used for this beam setup. */ bool isMPIOn() const { return MPIHandler_ && MPIHandler_->beamOK(); } /** * Access function for the MPIHandler, it should only be called after * checking with isMPIOn. */ tUEBasePtr getMPIHandler() const { assert(MPIHandler_); return MPIHandler_; } /** * Is a beam particle where hadronic structure is resolved */ bool isResolvedHadron(tPPtr); /** * Get the remnants from the ThePEG::PartonBinInstance es and * do some checks. */ RemPair getRemnants(PBIPair incbins); /** * Reset the PDF's after the hard collision has been showered */ void setMPIPDFs(); //@} public: /** * Check if a particle decays in the shower * @param id The PDG code for the particle */ bool decaysInShower(long id) const { return ( particlesDecayInShower_.find( abs(id) ) != particlesDecayInShower_.end() ); } protected: /** * Members to handle splitting up of hard process and decays */ //@{ /** * Find decay products from the hard process and create decay processes * @param parent The parent particle * @param hard The hard process * @param decay The decay processes */ void findDecayProducts(PPtr parent, PerturbativeProcessPtr hard, DecayProcessMap & decay) const; /** * Find decay products from the hard process and create decay processes * @param parent The parent particle * @param hard The parent hard process * @param decay The decay processes */ void createDecayProcess(PPtr parent,PerturbativeProcessPtr hard, DecayProcessMap & decay) const; //@} /** * @name Functions to return information relevant to the process being showered */ //@{ /** * Return the currently used SubProcess. */ tSubProPtr currentSubProcess() const { assert(subProcess_); return subProcess_; } /** * Access to the incoming beam particles */ tPPair incomingBeams() const { return incoming_; } //@} protected: /** * Weight handling for shower variations */ //@ /** * Combine the variation weights which have been encountered */ void combineWeights(); /** * Initialise the weights in currentEvent() */ void initializeWeights(); /** * Reset the current weights */ void resetWeights(); //@} protected: /** * Return the maximum number of attempts for showering * a given subprocess. */ unsigned int maxtry() const { return maxtry_; } protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ShowerHandler & operator=(const ShowerHandler &) = delete; private: /** * pointer to "this", the current ShowerHandler. */ static tShowerHandlerPtr currentHandler_; /** * a MPIHandler to administer the creation of several (semihard) * partonic interactions. */ UEBasePtr MPIHandler_; /** * Pointer to the HwRemDecayer */ HwRemDecPtr remDec_; private: /** * Maximum tries for various stages of the showering process */ //@{ /** * Maximum number of attempts for the * main showering loop */ unsigned int maxtry_; /** * Maximum number of attempts for the regeneration of an additional * scattering, before the number of scatters is reduced. */ unsigned int maxtryMPI_; /** * Maximum number of attempts for the regeneration of an additional * hard scattering, before this event is vetoed. */ unsigned int maxtryDP_; /** * Maximum number of attempts to generate a decay */ unsigned int maxtryDecay_; //@} private: /** * Factors for the various scales */ //@{ /** * The factorization scale factor. */ double factorizationScaleFactor_; /** * The renormalization scale factor. */ double renormalizationScaleFactor_; /** * The scale factor for the hard scale */ double hardScaleFactor_; /** * True, if the phase space restrictions of the dipole shower should * be applied. */ bool restrictPhasespace_; /** * True if maximum pt should be deduced from the factorization scale */ bool maxPtIsMuF_; /** * The profile scales */ Ptr::ptr hardScaleProfile_; //@} private: /** * Storage of information about the current event */ //@{ /** * The incoming beam particles for the current collision */ tPPair incoming_; /** * Boost to get back to the lab */ LorentzRotation boost_; /** * Const pointer to the currently handeled ThePEG::SubProcess */ tSubProPtr subProcess_; /** * Const pointer to the current step */ tcStepPtr current_; //@} private: /** * PDFs to be used for the various stages and related parameters */ //@{ /** * The PDF freezing scale */ Energy pdfFreezingScale_; /** * PDFs to be used for the various stages and related parameters */ //@{ /** * The PDF for beam particle A. Overrides the particle's own PDF setting. */ PDFPtr PDFA_; /** * The PDF for beam particle B. Overrides the particle's own PDF setting. */ PDFPtr PDFB_; /** * The PDF for beam particle A for remnant splitting. Overrides the particle's own PDF setting. */ PDFPtr PDFARemnant_; /** * The PDF for beam particle B for remnant splitting. Overrides the particle's own PDF setting. */ PDFPtr PDFBRemnant_; /** * The MPI PDF's to be used for secondary scatters. */ pair mpipdfs_; /** * The MPI PDF's to be used for secondary scatters. */ pair rempdfs_; /** * The MPI PDF's to be used for secondary scatters. */ pair remmpipdfs_; //@} private: /** * @name Parameters for initial- and final-state radiation */ //@{ /** * Switch on or off final state radiation. */ bool doFSR_; /** * Switch on or off initial state radiation. */ bool doISR_; //@} private: /** * @name Parameters for particle decays */ //@{ /** * Whether or not to split into hard and decay trees */ bool splitHardProcess_; /** * PDG codes of the particles which decay during showering * this is fast storage for use during running */ set particlesDecayInShower_; /** * PDG codes of the particles which decay during showering * this is a vector that is interfaced so they can be changed */ vector inputparticlesDecayInShower_; //@} private: /** * Parameters for the space-time model */ //@{ /** * Whether or not to include spa-cetime distances in the shower */ bool includeSpaceTime_; /** * The minimum virtuality for the space-time model */ Energy2 vMin_; //@} private: /** * Parameters for the constituent mass treatment. */ //@{ // True if shower should return constituent masses. bool useConstituentMasses_=true; //@} private: /** * Parameters relevant for reweight and variations */ //@{ /** * The shower variations */ map showerVariations_; /** * Command to add a shower variation */ string doAddVariation(string); /** * A reweighting factor applied by the showering */ double reweight_; /** * The shower variation weights */ map currentWeights_; //@} }; } #endif /* HERWIG_ShowerHandler_H */