diff --git a/Decay/General/FtoFFFDecayer.cc b/Decay/General/FtoFFFDecayer.cc --- a/Decay/General/FtoFFFDecayer.cc +++ b/Decay/General/FtoFFFDecayer.cc @@ -1,311 +1,311 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the FtoFFFDecayer class. // #include "FtoFFFDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Decay/DecayPhaseSpaceMode.h" #include "Herwig/PDT/ThreeBodyAllOnCalculator.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include using namespace Herwig; IBPtr FtoFFFDecayer::clone() const { return new_ptr(*this); } IBPtr FtoFFFDecayer::fullclone() const { return new_ptr(*this); } void FtoFFFDecayer::persistentOutput(PersistentOStream & os) const { os << sca_ << vec_ << ten_; } void FtoFFFDecayer::persistentInput(PersistentIStream & is, int) { is >> sca_ >> vec_ >> ten_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigFtoFFFDecayer("Herwig::FtoFFFDecayer", "Herwig.so"); void FtoFFFDecayer::Init() { static ClassDocumentation documentation ("The FtoFFFDecayer class implements the general decay of a fermion to " "three fermions."); } -void FtoFFFDecayer::doinit() { - GeneralThreeBodyDecayer::doinit(); +void FtoFFFDecayer::setupDiagrams(bool kinCheck) { + GeneralThreeBodyDecayer::setupDiagrams(kinCheck); if(outgoing().empty()) return; unsigned int ndiags = getProcessInfo().size(); sca_.resize(ndiags); vec_.resize(ndiags); ten_.resize(ndiags); for(unsigned int ix = 0;ix < ndiags; ++ix) { TBDiagram current = getProcessInfo()[ix]; tcPDPtr offshell = current.intermediate; if( offshell->CC() ) offshell = offshell->CC(); if(offshell->iSpin() == PDT::Spin0) { AbstractFFSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFSVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a scalar diagram in FtoFFFDecayer::doinit()" + << "Invalid vertices for a scalar diagram in FtoFFFDecayer::setupDiagrams()" << Exception::runerror; sca_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin1) { AbstractFFVVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFVVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a vector diagram in FtoFFFDecayer::doinit()" + << "Invalid vertices for a vector diagram in FtoFFFDecayer::setupDiagrams()" << Exception::runerror; vec_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin2) { AbstractFFTVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFTVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a tensor diagram in FtoFFFDecayer::doinit()" + << "Invalid vertices for a tensor diagram in FtoFFFDecayer::setupDiagrams()" << Exception::runerror; ten_[ix] = make_pair(vert1, vert2); } } } double FtoFFFDecayer::me2(const int ichan, const Particle & inpart, const ParticleVector & decay, MEOption meopt) const { // particle or CC of particle bool cc = (*getProcessInfo().begin()).incoming != inpart.id(); // special handling or first/last call const vector > cfactors(getColourFactors()); const vector > nfactors(getLargeNcColourFactors()); const size_t ncf(numberOfFlows()); Energy2 scale(sqr(inpart.mass())); if(meopt==Initialize) { SpinorWaveFunction:: calculateWaveFunctions(inwave_.first,rho_,const_ptr_cast(&inpart), Helicity::incoming); inwave_.second.resize(2); if(inwave_.first[0].wave().Type() == SpinorType::u) { for(unsigned int ix = 0; ix < 2; ++ix) { inwave_.second[ix] = inwave_.first[ix].bar(); inwave_.second[ix].conjugate(); } } else { for(unsigned int ix = 0; ix < 2; ++ix) { inwave_.second[ix] = inwave_.first[ix].bar(); inwave_.first[ix].conjugate(); } } // fix rho if no correlations fixRho(rho_); } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()<0) SpinorWaveFunction::constructSpinInfo(inwave_.first, const_ptr_cast(&inpart), Helicity::incoming,true); else SpinorBarWaveFunction::constructSpinInfo(inwave_.second, const_ptr_cast(&inpart), Helicity::incoming,true); // outgoing particles for(unsigned int ix = 0; ix < 3; ++ix) { SpinorWaveFunction:: constructSpinInfo(outwave_[ix].first,decay[ix],Helicity::outgoing,true); } } // outgoing particles for(unsigned int ix = 0; ix < 3; ++ix) { SpinorWaveFunction:: calculateWaveFunctions(outwave_[ix].first,decay[ix],Helicity::outgoing); outwave_[ix].second.resize(2); if(outwave_[ix].first[0].wave().Type() == SpinorType::u) { for(unsigned int iy = 0; iy < 2; ++iy) { outwave_[ix].second[iy] = outwave_[ix].first[iy].bar(); outwave_[ix].first[iy].conjugate(); } } else { for(unsigned int iy = 0; iy < 2; ++iy) { outwave_[ix].second[iy] = outwave_[ix].first[iy].bar(); outwave_[ix].second[iy].conjugate(); } } } bool ferm = inpart.id()>0; vector flows(ncf, Complex(0.)),largeflows(ncf, Complex(0.)); static const unsigned int out2[3]={1,0,0},out3[3]={2,2,1}; vector mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half,PDT::Spin1Half))); vector mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half,PDT::Spin1Half))); unsigned int ihel[4]; for(ihel[0] = 0; ihel[0] < 2; ++ihel[0]) { for(ihel[1] = 0; ihel[1] < 2; ++ihel[1]) { for(ihel[2] = 0; ihel[2] < 2; ++ihel[2]) { for(ihel[3] = 0; ihel[3] < 2; ++ihel[3]) { flows = vector(ncf, Complex(0.)); largeflows = vector(ncf, Complex(0.)); unsigned int idiag=0; for(vector::const_iterator dit=getProcessInfo().begin(); dit!=getProcessInfo().end();++dit) { if(ichan>=0&&diagramMap()[ichan]!=idiag) { ++idiag; continue; } // the sign from normal ordering double sign = ferm ? 1. : -1; // outgoing wavefunction and NO sign if (dit->channelType==TBDiagram::channel23) sign *= -1.; else if(dit->channelType==TBDiagram::channel13) sign *= 1.; else if(dit->channelType==TBDiagram::channel12) sign *= -1.; else throw Exception() << "Unknown diagram type in FtoFFFDecayer::me2()" << Exception::runerror; // wavefunctions SpinorWaveFunction w0,w3; SpinorBarWaveFunction w1,w2; // incoming wavefunction if(ferm) { w0 = inwave_.first [ihel[0]]; w1 = outwave_[dit->channelType].second[ihel[dit->channelType+1]]; } else { w0 = outwave_[dit->channelType].first [ihel[dit->channelType+1]]; w1 = inwave_.second[ihel[0]]; } if(decay[out2[dit->channelType]]->id()<0&& decay[out3[dit->channelType]]->id()>0) { w2 = outwave_[out3[dit->channelType]].second[ihel[out3[dit->channelType]+1]]; w3 = outwave_[out2[dit->channelType]].first [ihel[out2[dit->channelType]+1]]; sign *= -1.; } else { w2 = outwave_[out2[dit->channelType]].second[ihel[out2[dit->channelType]+1]]; w3 = outwave_[out3[dit->channelType]].first [ihel[out3[dit->channelType]+1]]; } tcPDPtr offshell = dit->intermediate; if(cc&&offshell->CC()) offshell=offshell->CC(); Complex diag(0.); // intermediate scalar if (offshell->iSpin() == PDT::Spin0) { ScalarWaveFunction inters = sca_[idiag].first-> evaluate(scale, widthOption(), offshell, w0, w1); diag = sca_[idiag].second->evaluate(scale,w3,w2,inters); } // intermediate vector else if(offshell->iSpin() == PDT::Spin1) { VectorWaveFunction interv = vec_[idiag].first-> evaluate(scale, widthOption(), offshell, w0, w1); diag = vec_[idiag].second->evaluate(scale,w3,w2,interv); } // intermediate tensor else if(offshell->iSpin() == PDT::Spin2) { TensorWaveFunction intert = ten_[idiag].first-> evaluate(scale, widthOption(), offshell, w0, w1); diag = ten_[idiag].second->evaluate(scale,w3,w2,intert); } // unknown else throw Exception() << "Unknown intermediate in FtoFFFDecayer::me2()" << Exception::runerror; // apply NO sign diag *= sign; // matrix element for the different colour flows if(ichan<0) { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } else { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1!=colourFlow()) continue; flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1!=colourFlow()) continue; largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } ++idiag; } // now add the flows to the me2 with appropriate colour factors for(unsigned int ix = 0; ix < ncf; ++ix) { (*mes[ix])(ihel[0],ihel[1],ihel[2],ihel[3]) = flows[ix]; (*mel[ix])(ihel[0],ihel[1],ihel[2],ihel[3]) = largeflows[ix]; } } } } } double me2(0.); if(ichan<0) { vector pflows(ncf,0.); for(unsigned int ix = 0; ix < ncf; ++ix) { for(unsigned int iy = 0; iy < ncf; ++ iy) { double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real(); me2 += con; if(ix==iy) { con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real(); pflows[ix] += con; } } } double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.)); ptotal *=UseRandom::rnd(); for(unsigned int ix=0;ixcontract(*mel[iflow],rho_)).real(); } // return the matrix element squared return me2; } WidthCalculatorBasePtr FtoFFFDecayer:: threeBodyMEIntegrator(const DecayMode & ) const { vector intype; vector inmass,inwidth; vector inpow,inweights; constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights); return new_ptr(ThreeBodyAllOnCalculator (inweights,intype,inmass,inwidth,inpow,*this,0, outgoing()[0]->mass(),outgoing()[1]->mass(),outgoing()[2]->mass(), relativeError())); } diff --git a/Decay/General/FtoFFFDecayer.h b/Decay/General/FtoFFFDecayer.h --- a/Decay/General/FtoFFFDecayer.h +++ b/Decay/General/FtoFFFDecayer.h @@ -1,142 +1,137 @@ // -*- C++ -*- #ifndef HERWIG_FtoFFFDecayer_H #define HERWIG_FtoFFFDecayer_H // // This is the declaration of the FtoFFFDecayer class. // #include "GeneralThreeBodyDecayer.h" #include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h" namespace Herwig { using namespace ThePEG; /** * The FtoFFFDecayer class provides the general matrix elements for the * decay of a (anti)fermion to three (anti)fermions. * * @see \ref FtoFFFDecayerInterfaces "The interfaces" * defined for FtoFFFDecayer. */ class FtoFFFDecayer: public GeneralThreeBodyDecayer { public: /** * Return the matrix element squared for a given mode and phase-space channel * @param ichan The channel we are calculating the matrix element for. * @param part The decaying Particle. * @param decay The particles produced in the decay. * @param meopt Option for the matrix element * @return The matrix element squared for the phase-space configuration. */ virtual double me2(const int ichan, const Particle & part, const ParticleVector & decay, MEOption meopt) const; /** * Method to return an object to calculate the 3 (or higher body) partial width * @param dm The DecayMode * @return A pointer to a WidthCalculatorBase object capable of calculating the width */ virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} 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. + * Set up the diagrams etc */ - virtual void doinit(); - //@} + virtual void setupDiagrams(bool checkKinematics); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ FtoFFFDecayer & operator=(const FtoFFFDecayer &); private: /** * Store the vector of FFSVertex pairs */ vector > sca_; /** * Store the vector of FFVVertex pairs */ vector > vec_; /** * Store the vector of FFTVertex pairs */ vector > ten_; /** * Spin density matrix */ mutable RhoDMatrix rho_; /** * Spinors for incoming particle */ mutable pair,vector > inwave_; /** * Spinors for outgoing particles */ mutable pair,vector > outwave_[3]; }; } #endif /* HERWIG_FtoFFFDecayer_H */ diff --git a/Decay/General/FtoFVVDecayer.cc b/Decay/General/FtoFVVDecayer.cc --- a/Decay/General/FtoFVVDecayer.cc +++ b/Decay/General/FtoFVVDecayer.cc @@ -1,403 +1,403 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the FtoFVVDecayer class. // #include "FtoFVVDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Decay/DecayPhaseSpaceMode.h" #include "Herwig/PDT/ThreeBodyAllOnCalculator.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include using namespace Herwig; IBPtr FtoFVVDecayer::clone() const { return new_ptr(*this); } IBPtr FtoFVVDecayer::fullclone() const { return new_ptr(*this); } void FtoFVVDecayer::persistentOutput(PersistentOStream & os) const { os << sca_ << fer_ << vec_ << ten_; } void FtoFVVDecayer::persistentInput(PersistentIStream & is, int) { is >> sca_ >> fer_ >> vec_ >> ten_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigFtoFVVDecayer("Herwig::FtoFVVDecayer", "Herwig.so"); void FtoFVVDecayer::Init() { static ClassDocumentation documentation ("The FtoFVVDecayer class implements the general decay of a fermion to " "a fermion and a pair of vectors."); } -void FtoFVVDecayer::doinit() { - GeneralThreeBodyDecayer::doinit(); +void FtoFVVDecayer::setupDiagrams(bool kinCheck) { + GeneralThreeBodyDecayer::setupDiagrams(kinCheck); if(outgoing().empty()) return; unsigned int ndiags = getProcessInfo().size(); sca_.resize(ndiags); fer_.resize(ndiags); vec_.resize(ndiags); ten_.resize(ndiags); for(unsigned int ix = 0;ix < ndiags; ++ix) { TBDiagram current = getProcessInfo()[ix]; tcPDPtr offshell = current.intermediate; if(offshell->iSpin() == PDT::Spin0) { AbstractFFSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractVVSVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a scalar diagram in FtoFVVDecayer::doinit()" + << "Invalid vertices for a scalar diagram in FtoFVVDecayer::setupDiagrams()" << Exception::runerror; sca_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin1Half) { AbstractFFVVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFVVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a scalar diagram in FtoFVVDecayer::doinit()" + << "Invalid vertices for a scalar diagram in FtoFVVDecayer::setupDiagrams()" << Exception::runerror; fer_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin1) { AbstractFFVVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractVVVVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a vector diagram in FtoFVVDecayer::doinit()" + << "Invalid vertices for a vector diagram in FtoFVVDecayer::setupDiagrams()" << Exception::runerror; vec_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin2) { AbstractFFTVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractVVTVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a tensor diagram in FtoFVVDecayer::doinit()" + << "Invalid vertices for a tensor diagram in FtoFVVDecayer::setupDiagrams()" << Exception::runerror; ten_[ix] = make_pair(vert1, vert2); } } } double FtoFVVDecayer::me2(const int ichan, const Particle & inpart, const ParticleVector & decay, MEOption meopt) const { // particle or CC of particle bool cc = (*getProcessInfo().begin()).incoming != inpart.id(); // special handling or first/last call //Set up wave-functions bool ferm( inpart.id() > 0 ); if(meopt==Initialize) { if( ferm ) { SpinorWaveFunction:: calculateWaveFunctions(fwave_,rho_,const_ptr_cast(&inpart), Helicity::incoming); if( fwave_[0].wave().Type() != SpinorType::u ) fwave_[0].conjugate(); if( fwave_[1].wave().Type() != SpinorType::u ) fwave_[1].conjugate(); } else { SpinorBarWaveFunction:: calculateWaveFunctions(fbwave_, rho_, const_ptr_cast(&inpart), Helicity::incoming); if( fbwave_[0].wave().Type() != SpinorType::v ) fbwave_[0].conjugate(); if( fbwave_[1].wave().Type() != SpinorType::v ) fbwave_[1].conjugate(); } // fix rho if no correlations fixRho(rho_); } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(ferm) SpinorWaveFunction::constructSpinInfo(fwave_, const_ptr_cast(&inpart), Helicity::incoming,true); else SpinorBarWaveFunction::constructSpinInfo(fbwave_, const_ptr_cast(&inpart), Helicity::incoming,true); int ivec(-1); // outgoing particles for(int ix = 0; ix < 3; ++ix) { tPPtr p = decay[ix]; if( p->dataPtr()->iSpin() == PDT::Spin1Half ) { if( ferm ) { SpinorBarWaveFunction:: constructSpinInfo(fbwave_, p, Helicity::outgoing,true); } else { SpinorWaveFunction:: constructSpinInfo(fwave_ , p, Helicity::outgoing,true); } } else if( p->dataPtr()->iSpin() == PDT::Spin1 ) { if( ivec < 0 ) { ivec = ix; VectorWaveFunction:: constructSpinInfo(vwave_.first , p, Helicity::outgoing, true, false); } else { VectorWaveFunction:: constructSpinInfo(vwave_.second, p, Helicity::outgoing, true, false); } } } return 0.; } // outgoing, keep track of fermion and first occurrence of vector positions int isp(-1), ivec(-1); // outgoing particles pair mass = make_pair(false,false); for(int ix = 0; ix < 3; ++ix) { tPPtr p = decay[ix]; if( p->dataPtr()->iSpin() == PDT::Spin1Half ) { isp = ix; if( ferm ) { SpinorBarWaveFunction:: calculateWaveFunctions(fbwave_, p, Helicity::outgoing); if( fbwave_[0].wave().Type() != SpinorType::u ) fbwave_[0].conjugate(); if( fbwave_[1].wave().Type() != SpinorType::u ) fbwave_[1].conjugate(); } else { SpinorWaveFunction:: calculateWaveFunctions(fwave_, p, Helicity::outgoing); if( fwave_[0].wave().Type() != SpinorType::v ) fwave_[0].conjugate(); if( fwave_[1].wave().Type() != SpinorType::v ) fwave_[1].conjugate(); } } else if( p->dataPtr()->iSpin() == PDT::Spin1 ) { bool massless = p->id() == ParticleID::gamma || p->id() == ParticleID::g; if( ivec < 0 ) { ivec = ix; VectorWaveFunction:: calculateWaveFunctions(vwave_.first , p, Helicity::outgoing, massless); mass.first = massless; } else { VectorWaveFunction:: calculateWaveFunctions(vwave_.second, p, Helicity::outgoing, massless); mass.second = massless; } } } assert(isp >= 0 && ivec >= 0); Energy2 scale(sqr(inpart.mass())); const vector > cfactors(getColourFactors()); const vector > nfactors(getLargeNcColourFactors()); const size_t ncf(numberOfFlows()); vector flows(ncf, Complex(0.)), largeflows(ncf, Complex(0.)); vector mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half, (isp == 0) ? PDT::Spin1Half : PDT::Spin1, (isp == 1) ? PDT::Spin1Half : PDT::Spin1, (isp == 2) ? PDT::Spin1Half : PDT::Spin1))); vector mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half, (isp == 0) ? PDT::Spin1Half : PDT::Spin1, (isp == 1) ? PDT::Spin1Half : PDT::Spin1, (isp == 2) ? PDT::Spin1Half : PDT::Spin1))); //Helicity calculation for( unsigned int if1 = 0; if1 < 2; ++if1 ) { for( unsigned int if2 = 0; if2 < 2; ++if2 ) { for( unsigned int iv1 = 0; iv1 < 3; ++iv1 ) { if ( mass.first && iv1 == 1 ) continue; for( unsigned int iv2 = 0; iv2 < 3; ++iv2 ) { if ( mass.second && iv2 == 1 ) continue; flows = vector(ncf, Complex(0.)); largeflows = vector(ncf, Complex(0.)); unsigned int idiag(0); for(vector::const_iterator dit = getProcessInfo().begin(); dit != getProcessInfo().end(); ++dit) { // If we are selecting a particular channel if( ichan >= 0 && diagramMap()[ichan] != idiag ) { ++idiag; continue; } tcPDPtr offshell = (*dit).intermediate; if(cc&&offshell->CC()) offshell=offshell->CC(); Complex diag; if( offshell->iSpin() == PDT::Spin1Half ) { // Make sure we connect the correct particles VectorWaveFunction vw1, vw2; if( (*dit).channelType == TBDiagram::channel23 ) { vw1 = vwave_.first[iv1]; vw2 = vwave_.second[iv2]; } else if( (*dit).channelType == TBDiagram::channel12 ) { vw1 = vwave_.second[iv2]; vw2 = vwave_.first[iv1]; } else { if( ivec < isp ) { vw1 = vwave_.second[iv2]; vw2 = vwave_.first[iv1]; } else { vw1 = vwave_.first[iv1]; vw2 = vwave_.second[iv2]; } } if( ferm ) { SpinorWaveFunction inters = fer_[idiag].first->evaluate(scale, widthOption(), offshell, fwave_[if1], vw1); diag = fer_[idiag].second->evaluate(scale, inters, fbwave_[if2], vw2); } else { SpinorBarWaveFunction inters = fer_[idiag].first->evaluate(scale, widthOption(), offshell, fbwave_[if2], vw1); diag = fer_[idiag].second->evaluate(scale, fwave_[if1], inters, vw2); } } else if( offshell->iSpin() == PDT::Spin0 ) { ScalarWaveFunction inters = sca_[idiag].first->evaluate(scale, widthOption(), offshell, fwave_[if1], fbwave_[if2]); diag = sca_[idiag].second->evaluate(scale, vwave_.first[iv1], vwave_.second[iv2], inters); } else if( offshell->iSpin() == PDT::Spin1 ) { VectorWaveFunction interv = vec_[idiag].first->evaluate(scale, widthOption(), offshell, fwave_[if1], fbwave_[if2]); diag = vec_[idiag].second->evaluate(scale, vwave_.first[iv1], vwave_.second[iv2], interv); } else if( offshell->iSpin() == PDT::Spin2 ) { TensorWaveFunction intert = ten_[idiag].first->evaluate(scale, widthOption(), offshell, fwave_[if1], fbwave_[if2]); diag = ten_[idiag].second->evaluate(scale, vwave_.first[iv1], vwave_.second[iv2], intert); } else throw Exception() << "Unknown intermediate in FtoFVVDecayer::me2()" << Exception::runerror; //NO sign if( !ferm ) diag *= -1; if(ichan<0) { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } else { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1!=colourFlow()) continue; flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1!=colourFlow()) continue; largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } ++idiag; }// end diagram loop // now add the flows to the me2 unsigned int h1(if1), h2(if2); if( !ferm ) swap(h1,h2); for(unsigned int ix = 0; ix < ncf; ++ix) { if(isp == 0) { (*mes[ix])(h1, h2, iv1, iv2) = flows[ix]; (*mel[ix])(h1, h2, iv1, iv2) = largeflows[ix]; } else if(isp == 1) { (*mes[ix])(h1, iv1, h2, iv2) = flows[ix]; (*mel[ix])(h1, iv1, h2, iv2) = largeflows[ix]; } else if(isp == 2) { (*mes[ix])(h1, iv1, iv2, h2) = flows[ix]; (*mel[ix])(h1, iv1, h2, iv2) = largeflows[ix]; } } }//v2 }//v1 }//f2 }//f1 //Finally, work out me2. This depends on whether we are selecting channels //or not double me2(0.); if(ichan<0) { vector pflows(ncf,0.); for(unsigned int ix = 0; ix < ncf; ++ix) { for(unsigned int iy = 0; iy < ncf; ++ iy) { double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real(); me2 += con; if(ix==iy) { con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real(); pflows[ix] += con; } } } double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.)); ptotal *= UseRandom::rnd(); for(unsigned int ix=0;ixcontract(*mel[iflow],rho_)).real(); } // return the matrix element squared return me2; } WidthCalculatorBasePtr FtoFVVDecayer:: threeBodyMEIntegrator(const DecayMode & ) const { vector intype; vector inmass,inwidth; vector inpow,inweights; constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights); return new_ptr(ThreeBodyAllOnCalculator (inweights,intype,inmass,inwidth,inpow,*this,0, outgoing()[0]->mass(),outgoing()[1]->mass(), outgoing()[2]->mass(),relativeError())); } diff --git a/Decay/General/FtoFVVDecayer.h b/Decay/General/FtoFVVDecayer.h --- a/Decay/General/FtoFVVDecayer.h +++ b/Decay/General/FtoFVVDecayer.h @@ -1,156 +1,151 @@ // -*- C++ -*- #ifndef HERWIG_FtoFVVDecayer_H #define HERWIG_FtoFVVDecayer_H // // This is the declaration of the FtoFVVDecayer class. // #include "GeneralThreeBodyDecayer.h" #include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h" #include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h" #include "ThePEG/Helicity/Vertex/AbstractVVTVertex.h" namespace Herwig { using namespace ThePEG; /** * The FtoFVVDecayer class provides the general matrix elements for the * decay of a fermion to a fermion and two vector bosons. * * @see \ref FtoFVVDecayerInterfaces "The interfaces" * defined for FtoFVVDecayer. */ class FtoFVVDecayer: public GeneralThreeBodyDecayer { public: /** * Return the matrix element squared for a given mode and phase-space channel * @param ichan The channel we are calculating the matrix element for. * @param part The decaying Particle. * @param decay The particles produced in the decay. * @param meopt Option for the calculation of the matrix element * @return The matrix element squared for the phase-space configuration. */ virtual double me2(const int ichan, const Particle & part, const ParticleVector & decay, MEOption meopt) const; /** * Method to return an object to calculate the 3 (or higher body) partial width * @param dm The DecayMode * @return A pointer to a WidthCalculatorBase object capable of calculating the width */ virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} 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. + * Set up the diagrams etc */ - virtual void doinit(); - //@} + virtual void setupDiagrams(bool checkKinematics); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ FtoFVVDecayer & operator=(const FtoFVVDecayer &); private: /** * Store the vector of scalar intermediates */ vector > sca_; /** * Store the vector for fermion intermediates */ vector > fer_; /** * Store the vector for gauge boson intermediates */ vector > vec_; /** * Store the vector of tensor intermediates */ vector > ten_; /** * Spin density matrix */ mutable RhoDMatrix rho_; /** * Spinor wavefunctions */ mutable vector fwave_; /** * Barred spinor wavefunctions */ mutable vector fbwave_; /** * Vector wavefunctions */ mutable pair, vector > vwave_; }; } #endif /* HERWIG_FtoFVVDecayer_H */ diff --git a/Decay/General/GeneralThreeBodyDecayer.cc b/Decay/General/GeneralThreeBodyDecayer.cc --- a/Decay/General/GeneralThreeBodyDecayer.cc +++ b/Decay/General/GeneralThreeBodyDecayer.cc @@ -1,1028 +1,1052 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the GeneralThreeBodyDecayer class. // #include "GeneralThreeBodyDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Decay/DecayPhaseSpaceMode.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/PDT/ThreeBodyAllOnCalculator.h" using namespace Herwig; void GeneralThreeBodyDecayer::persistentOutput(PersistentOStream & os) const { os << incoming_ << outgoing_ << diagrams_ << diagmap_ << colour_ << colourLargeNC_ << nflow_ << widthOpt_ << refTag_ << refTagCC_ << intOpt_ << relerr_; } void GeneralThreeBodyDecayer::persistentInput(PersistentIStream & is, int) { is >> incoming_ >> outgoing_ >> diagrams_ >> diagmap_ >> colour_ >> colourLargeNC_ >> nflow_ >> widthOpt_ >> refTag_ >> refTagCC_ >> intOpt_ >> relerr_; } // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractClass describeHerwigGeneralThreeBodyDecayer("Herwig::GeneralThreeBodyDecayer", "Herwig.so"); void GeneralThreeBodyDecayer::Init() { static ClassDocumentation documentation ("The GeneralThreeBodyDecayer class is the base class for the implementation of" " all three body decays based on spin structures in Herwig."); static Switch interfaceWidthOption ("WidthOption", "Option for the treatment of the widths of the intermediates", &GeneralThreeBodyDecayer::widthOpt_, 1, false, false); static SwitchOption interfaceWidthOptionFixed (interfaceWidthOption, "Fixed", "Use fixed widths", 1); static SwitchOption interfaceWidthOptionRunning (interfaceWidthOption, "Running", "Use running widths", 2); static SwitchOption interfaceWidthOptionZero (interfaceWidthOption, "Zero", "Set the widths to zero", 3); static Switch interfacePartialWidthIntegration ("PartialWidthIntegration", "Switch to control the partial width integration", &GeneralThreeBodyDecayer::intOpt_, 0, false, false); static SwitchOption interfacePartialWidthIntegrationAllPoles (interfacePartialWidthIntegration, "AllPoles", "Include all potential poles", 0); static SwitchOption interfacePartialWidthIntegrationShallowestPole (interfacePartialWidthIntegration, "ShallowestPole", "Only include the shallowest pole", 1); static Parameter interfaceRelativeError ("RelativeError", "The relative error for the GQ integration of the partial width", &GeneralThreeBodyDecayer::relerr_, 1e-2, 1e-10, 1., false, false, Interface::limited); } ParticleVector GeneralThreeBodyDecayer::decay(const Particle & parent, const tPDVector & children) const { // return empty vector if products heavier than parent Energy mout(ZERO); for(tPDVector::const_iterator it=children.begin(); it!=children.end();++it) mout+=(**it).massMin(); if(mout>parent.mass()) return ParticleVector(); // generate the decay bool cc; int imode=modeNumber(cc,parent.dataPtr(),children); // generate the kinematics ParticleVector decay=generate(generateIntermediates(),cc,imode,parent); // make the colour connections colourConnections(parent, decay); // return the answer return decay; } int GeneralThreeBodyDecayer:: modeNumber(bool & cc, tcPDPtr in, const tPDVector & outin) const { assert( !refTag_.empty() && !refTagCC_.empty() ); // check number of outgoing particles if( outin.size() != 3 || abs(in->id()) != abs(incoming_->id()) ) return -1; OrderedParticles testmode(outin.begin(), outin.end()); OrderedParticles::const_iterator dit = testmode.begin(); string testtag(in->name() + "->"); for( unsigned int i = 1; dit != testmode.end(); ++dit, ++i) { testtag += (**dit).name(); if( i != 3 ) testtag += string(","); } if( testtag == refTag_ ) { cc = false; return 0; } else if ( testtag == refTagCC_ ) { cc = true; return 0; } else return -1; } bool GeneralThreeBodyDecayer::setDecayInfo(PDPtr incoming, vector outgoing, const vector & process, double symfac) { // set the member variables from the info supplied incoming_ = incoming; outgoing_ = outgoing; diagrams_ = process; assert( outgoing_.size() == 3 ); // Construct reference tags for testing in modeNumber function OrderedParticles refmode(outgoing_.begin(), outgoing_.end()); OrderedParticles::const_iterator dit = refmode.begin(); refTag_ = incoming_->name() + "->"; for( unsigned int i = 1; dit != refmode.end(); ++dit, ++i) { refTag_ += (**dit).name(); if( i != 3 ) refTag_ += string(","); } //CC-mode refmode.clear(); refTagCC_ = incoming_->CC() ? incoming_->CC()->name() : incoming_->name(); refTagCC_ += "->"; for( unsigned int i = 0; i < 3; ++i ) { if( outgoing_[i]->CC() ) refmode.insert( outgoing_[i]->CC() ); else refmode.insert( outgoing_[i] ); } dit = refmode.begin(); for( unsigned int i = 1; dit != refmode.end(); ++dit , ++i) { refTagCC_ += (**dit).name(); if( i != 3 ) refTagCC_ += string(","); } // check if intermeidates or only four point diagrams bool intermediates(false); for(auto diagram : diagrams_) { if(diagram.intermediate) { intermediates=true; break; } } if(!intermediates) { incoming_= PDPtr(); outgoing_.clear(); generator()->log() << "Only four body diagrams for decay " << refTag_ << " in GeneralThreeBodyDecayer::" << "setDecayInfo(), omitting decay\n"; return false; } // set the colour factors and return the answer if(setColourFactors(symfac)) return true; incoming_= PDPtr(); outgoing_.clear(); return false; } void GeneralThreeBodyDecayer::doinit() { DecayIntegrator::doinit(); if(outgoing_.empty()) return; - // create the phase space integrator - tPDVector extpart(1,incoming_); - extpart.insert(extpart.end(),outgoing_.begin(),outgoing_.end()); - // create the integration channels for the decay - DecayPhaseSpaceModePtr mode(new_ptr(DecayPhaseSpaceMode(extpart,this,true))); - DecayPhaseSpaceChannelPtr newchannel; - // create the phase-space channels for the integration - unsigned int nmode(0); - for(unsigned int ix=0;ixmass() == ZERO || - diagrams_[ix].intermediate->width() == ZERO ) { - jac = 1; - power = -2.0; - } - if(diagrams_[ix].channelType==TBDiagram::channel23) { - newchannel->addIntermediate(extpart[0],0,0.0,-1,1); - newchannel->addIntermediate(diagrams_[ix].intermediate,jac,power, 2,3); - } - else if(diagrams_[ix].channelType==TBDiagram::channel13) { - newchannel->addIntermediate(extpart[0],0,0.0,-1,2); - newchannel->addIntermediate(diagrams_[ix].intermediate,jac,power, 1,3); - } - else if(diagrams_[ix].channelType==TBDiagram::channel12) { - newchannel->addIntermediate(extpart[0],0,0.0,-1,3); - newchannel->addIntermediate(diagrams_[ix].intermediate,jac,power, 1,2); - } - diagmap_.push_back(ix); - mode->addChannel(newchannel); - ++nmode; - } - if(nmode==0) { - string mode = extpart[0]->PDGName() + "->"; - for(unsigned int ix=1;ixPDGName() + " "; - throw Exception() << "No decay channels in GeneralThreeBodyDecayer::" - << "doinit() for " << mode << "\n" << Exception::runerror; - } - // add the mode - vector wgt(nmode,1./double(nmode)); - addMode(mode,1.,wgt); + setupDiagrams(false); +} + +void GeneralThreeBodyDecayer::doinitrun() { + setupDiagrams(true); + DecayIntegrator::doinitrun(); } double GeneralThreeBodyDecayer:: threeBodyMatrixElement(const int imode, const Energy2 q2, const Energy2 s3, const Energy2 s2, const Energy2 s1, const Energy m1, const Energy m2, const Energy m3) const { // calculate the momenta of the outgoing particles Energy m0=sqrt(q2); // energies Energy eout[3] = {0.5*(q2+sqr(m1)-s1)/m0, 0.5*(q2+sqr(m2)-s2)/m0, 0.5*(q2+sqr(m3)-s3)/m0}; // magnitudes of the momenta Energy pout[3] = {sqrt(sqr(eout[0])-sqr(m1)), sqrt(sqr(eout[1])-sqr(m2)), sqrt(sqr(eout[2])-sqr(m3))}; double cos2 = 0.5*(sqr(pout[0])+sqr(pout[1])-sqr(pout[2]))/pout[0]/pout[1]; double cos3 = 0.5*(sqr(pout[0])-sqr(pout[1])+sqr(pout[2]))/pout[0]/pout[2]; double sin2 = sqrt(1.-sqr(cos2)), sin3 = sqrt(1.-sqr(cos3)); Lorentz5Momentum out[3]= {Lorentz5Momentum( ZERO , ZERO , pout[0] , eout[0] , m1), Lorentz5Momentum( pout[1]*sin2 , ZERO , -pout[1]*cos2 , eout[1] , m2), Lorentz5Momentum( -pout[2]*sin3 , ZERO , -pout[2]*cos3 , eout[2] , m3)}; // create the incoming PPtr inpart=mode(imode)->externalParticles(0)-> produceParticle(Lorentz5Momentum(sqrt(q2))); // and outgoing particles ParticleVector decay; for(unsigned int ix=1;ix<4;++ix) decay.push_back(mode(imode)->externalParticles(ix)->produceParticle(out[ix-1])); // return the matrix element return me2(-1,*inpart,decay,Initialize); } double GeneralThreeBodyDecayer::brat(const DecayMode &, const Particle & p, double oldbrat) const { ParticleVector children = p.children(); if( children.size() != 3 || !p.data().widthGenerator() ) return oldbrat; // partial width for this mode Energy scale = p.mass(); Energy pwidth = partialWidth( make_pair(p.dataPtr(), scale), make_pair(children[0]->dataPtr(), children[0]->mass()), make_pair(children[1]->dataPtr(), children[1]->mass()), make_pair(children[2]->dataPtr(), children[2]->mass()) ); Energy width = p.data().widthGenerator()->width(p.data(), scale); return pwidth/width; } Energy GeneralThreeBodyDecayer::partialWidth(PMPair inpart, PMPair outa, PMPair outb, PMPair outc) const { if(inpart.secondname() + "->"; tag += outgoing_[0]->name() + "," + outgoing_[1]->name() + "," + outgoing_[2]->name() + ";"; DMPtr dm = generator()->findDecayMode(tag); widthCalc_ = threeBodyMEIntegrator(*dm); } return widthCalc_->partialWidth(sqr(inpart.second)); } void GeneralThreeBodyDecayer:: colourConnections(const Particle & parent, const ParticleVector & out) const { // first extract the outgoing particles and intermediate PPtr inter; ParticleVector outgoing; if(!generateIntermediates()) { outgoing=out; } else { // find the diagram unsigned int idiag = diagramMap()[mode(imode())->selectedChannel()]; PPtr child; for(unsigned int ix=0;ixchildren().empty()) child = out[ix]; else inter = out[ix]; } outgoing.resize(3); switch(diagrams_[idiag].channelType) { case TBDiagram::channel23: outgoing[0] = child; outgoing[1] = inter->children()[0]; outgoing[2] = inter->children()[1]; break; case TBDiagram::channel13: outgoing[0] = inter->children()[0]; outgoing[1] = child; outgoing[2] = inter->children()[1]; break; case TBDiagram::channel12: outgoing[0] = inter->children()[0]; outgoing[1] = inter->children()[1]; outgoing[2] = child; break; default: throw Exception() << "unknown diagram type in GeneralThreeBodyDecayer::" << "colourConnections()" << Exception::runerror; } } // extract colour of the incoming and outgoing particles PDT::Colour inColour(parent.data().iColour()); vector outColour; vector singlet,octet,triplet,antitriplet; for(unsigned int ix=0;ixdata().iColour()); switch(outColour.back()) { case PDT::Colour0 : singlet.push_back(ix); break; case PDT::Colour3 : triplet.push_back(ix); break; case PDT::Colour3bar: antitriplet.push_back(ix); break; case PDT::Colour8 : octet.push_back(ix); break; default: throw Exception() << "Unknown colour for particle in GeneralThreeBodyDecayer::" << "colourConnections()" << Exception::runerror; } } // colour neutral decaying particle if ( inColour == PDT::Colour0) { // options are all neutral or triplet/antitriplet+ neutral if(singlet.size()==3) return; else if(singlet.size()==1&&triplet.size()==1&&antitriplet.size()==1) { outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]); // add intermediate if needed if(inter&&inter->coloured()) { if(inter->dataPtr()->iColour()==PDT::Colour3) outgoing[triplet[0]]->colourLine()->addColoured(inter); else if(inter->dataPtr()->iColour()==PDT::Colour3bar) outgoing[triplet[0]]->colourLine()->addAntiColoured(inter); } } else if(octet.size()==1&&triplet.size()==1&&antitriplet.size()==1) { outgoing[ triplet[0]]->antiColourNeighbour(outgoing[octet[0]]); outgoing[antitriplet[0]]-> colourNeighbour(outgoing[octet[0]]); if(inter&&inter->coloured()) { if(inter->dataPtr()->iColour()==PDT::Colour3) outgoing[antitriplet[0]]->antiColourLine()->addColoured(inter); else if(inter->dataPtr()->iColour()==PDT::Colour3bar) outgoing[ triplet[0]]-> colourLine()->addAntiColoured(inter); else if(inter->dataPtr()->iColour()==PDT::Colour8) { outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter); outgoing[ triplet[0]]-> colourLine()->addColoured(inter); } } } else if(triplet.size()==3) { tColinePtr col[3] = {ColourLine::create(outgoing[0]), ColourLine::create(outgoing[1]), ColourLine::create(outgoing[2])}; col[0]->setSourceNeighbours(col[1],col[2]); } else if(antitriplet.size()==3) { tColinePtr col[3] = {ColourLine::create(outgoing[0],true), ColourLine::create(outgoing[1],true), ColourLine::create(outgoing[2],true)}; col[0]->setSinkNeighbours(col[1],col[2]); } else { string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " " + out[1]->PDGName() + " " + out[2]->PDGName(); throw Exception() << "Unknown colour structure in GeneralThreeBodyDecayer::" << "colourConnections() for singlet decaying particle " << mode << Exception::runerror; } } // colour triplet decaying particle else if( inColour == PDT::Colour3) { if(singlet.size()==2&&triplet.size()==1) { outgoing[triplet[0]]->incomingColour(const_ptr_cast(&parent)); if(inter&&inter->coloured()) outgoing[triplet[0]]->colourLine()->addColoured(inter); } else if(antitriplet.size()==1&&triplet.size()==2) { if(colourFlow()==0) { outgoing[triplet[0]]->incomingColour(const_ptr_cast(&parent)); outgoing[antitriplet[0]]->colourNeighbour(outgoing[triplet[1]]); if(inter&&inter->coloured()) { switch (inter->dataPtr()->iColour()) { case PDT::Colour8: inter->incomingColour(const_ptr_cast(&parent)); outgoing[triplet[1]]->colourLine()->addAntiColoured(inter); break; default: string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " " + out[1]->PDGName() + " " + out[2]->PDGName(); throw Exception() << "Unknown colour for intermediate in " << "GeneralThreeBodyDecayer::" << "colourConnections() for " << "decaying colour triplet " << mode << Exception::runerror; } } } else { outgoing[triplet[1]]->incomingColour(const_ptr_cast(&parent)); outgoing[antitriplet[0]]->colourNeighbour(outgoing[triplet[0]]); if(inter&&inter->coloured()) { switch (inter->dataPtr()->iColour()) { case PDT::Colour8: inter->incomingColour(const_ptr_cast(&parent)); outgoing[triplet[0]]->colourLine()->addAntiColoured(inter); break; default: string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " " + out[1]->PDGName() + " " + out[2]->PDGName(); throw Exception() << "Unknown colour for intermediate in " << "GeneralThreeBodyDecayer::" << "colourConnections() for " << "decaying colour triplet " << mode << Exception::runerror; } } } } else if (singlet.size()==1&&triplet.size()==1&&octet.size()==1) { if(inter) { if(inter->children()[0]->dataPtr()->iColour()==PDT::Colour8 || inter->children()[1]->dataPtr()->iColour()==PDT::Colour8) { inter->incomingColour(const_ptr_cast(&parent)); outgoing[octet[0]]->incomingColour(inter); outgoing[octet[0]]->colourNeighbour(outgoing[triplet[0]]); } else { outgoing[octet[0]]->incomingColour(inter); outgoing[octet[0]]->colourNeighbour(inter); outgoing[triplet[0]]->incomingColour(inter); } } else { outgoing[octet[0]]->incomingColour(const_ptr_cast(&parent)); outgoing[octet[0]]->colourNeighbour(outgoing[triplet[0]]); } } else if (singlet.size()==1&&antitriplet.size()==2) { tColinePtr col[2] = {ColourLine::create(outgoing[antitriplet[0]],true), ColourLine::create(outgoing[antitriplet[1]],true)}; parent.colourLine()->setSinkNeighbours(col[0],col[1]); } else { string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " " + out[1]->PDGName() + " " + out[2]->PDGName(); throw Exception() << "Unknown colour structure in GeneralThreeBodyDecayer::" << "colourConnections() for triplet decaying particle " << mode << Exception::runerror; } } else if( inColour == PDT::Colour3bar) { if(singlet.size()==2&&antitriplet.size()==1) { outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast(&parent)); } else if(antitriplet.size()==2&&triplet.size()==1) { if(colourFlow()==0) { outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast(&parent)); outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[1]]); if(inter&&inter->coloured()) { switch (inter->dataPtr()->iColour()) { case PDT::Colour8: inter->incomingAntiColour(const_ptr_cast(&parent)); outgoing[antitriplet[1]]->antiColourLine()->addAntiColoured(inter); break; default: string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " " + out[1]->PDGName() + " " + out[2]->PDGName(); throw Exception() << "Unknown colour for intermediate in" << " GeneralThreeBodyDecayer::" << "colourConnections() for " << "decaying colour antitriplet " << mode << Exception::runerror; } } } else { outgoing[antitriplet[1]]->incomingAntiColour(const_ptr_cast(&parent)); outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]); if(inter&&inter->coloured()) { switch (inter->dataPtr()->iColour()) { case PDT::Colour8: inter->incomingAntiColour(const_ptr_cast(&parent)); outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter); break; default: string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " " + out[1]->PDGName() + " " + out[2]->PDGName(); throw Exception() << "Unknown colour for intermediate in " << "GeneralThreeBodyDecayer::" << "colourConnections() for " << "decaying colour antitriplet " << mode << Exception::runerror; } } } } else if (singlet.size()==1&&antitriplet.size()==1&&octet.size()==1) { if(inter) { if(inter->children()[0]->dataPtr()->iColour()==PDT::Colour8 || inter->children()[1]->dataPtr()->iColour()==PDT::Colour8) { inter->incomingColour(const_ptr_cast(&parent)); outgoing[octet[0]]->incomingAntiColour(inter); outgoing[octet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]); } else { outgoing[octet[0]]->incomingAntiColour(inter); outgoing[octet[0]]->antiColourNeighbour(inter); outgoing[antitriplet[0]]->incomingAntiColour(inter); } } else { outgoing[octet[0]]->incomingAntiColour(const_ptr_cast(&parent)); outgoing[octet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]); } } else if (singlet.size()==1&&triplet.size()==2) { tColinePtr col[2] = {ColourLine::create(outgoing[triplet[0]]), ColourLine::create(outgoing[triplet[1]])}; parent.antiColourLine()->setSourceNeighbours(col[0],col[1]); } else { string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " " + out[1]->PDGName() + " " + out[2]->PDGName(); throw Exception() << "Unknown colour structure in GeneralThreeBodyDecayer::" << "colourConnections() for anti-triplet decaying particle" << mode << Exception::runerror; } } else if( inColour == PDT::Colour8) { if(triplet.size()==1&&antitriplet.size()==1&&singlet.size()==1) { outgoing[ triplet[0]]->incomingColour (const_ptr_cast(&parent)); outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast(&parent)); if(inter&&inter->coloured()) { switch (inter->dataPtr()->iColour()) { case PDT::Colour3: outgoing[triplet[0]]->colourLine()->addColoured(inter); break; case PDT::Colour3bar: outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter); break; default: string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " " + out[1]->PDGName() + " " + out[2]->PDGName(); throw Exception() << "Unknown colour for intermediate" << " in GeneralThreeBodyDecayer::" << "colourConnections() for " << "decaying colour octet " << mode << Exception::runerror; } } } else if(triplet.size()==3) { tColinePtr col[2]; if(colourFlow()==0) { outgoing[0]->incomingColour (const_ptr_cast(&parent)); col[0] = ColourLine::create(outgoing[1]); col[1] = ColourLine::create(outgoing[2]); } else if(colourFlow()==1) { outgoing[1]->incomingColour (const_ptr_cast(&parent)); col[0] = ColourLine::create(outgoing[0]); col[1] = ColourLine::create(outgoing[2]); } else if(colourFlow()==2) { outgoing[2]->incomingColour (const_ptr_cast(&parent)); col[0] = ColourLine::create(outgoing[0]); col[1] = ColourLine::create(outgoing[1]); } else assert(false); parent.antiColourLine()->setSourceNeighbours(col[0],col[1]); } else if(antitriplet.size()==3) { tColinePtr col[2]; if(colourFlow()==0) { outgoing[0]->incomingAntiColour(const_ptr_cast(&parent)); col[0] = ColourLine::create(outgoing[1],true); col[1] = ColourLine::create(outgoing[2],true); } else if(colourFlow()==1) { outgoing[1]->incomingAntiColour(const_ptr_cast(&parent)); col[0] = ColourLine::create(outgoing[0],true); col[1] = ColourLine::create(outgoing[2],true); } else if(colourFlow()==2) { outgoing[2]->incomingAntiColour(const_ptr_cast(&parent)); col[0] = ColourLine::create(outgoing[0],true); col[1] = ColourLine::create(outgoing[1],true); } else assert(false); parent.colourLine()->setSinkNeighbours(col[0],col[1]); } else { string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " " + out[1]->PDGName() + " " + out[2]->PDGName(); throw Exception() << "Unknown colour structure in GeneralThreeBodyDecayer::" << "colourConnections() for octet decaying particle" << mode << Exception::runerror; } } } void GeneralThreeBodyDecayer:: constructIntegratorChannels(vector & intype, vector & inmass, vector & inwidth, vector & inpow, vector & inweights) const { // check if any intermediate photons bool hasPhoton=false; for(unsigned int iy=0;iyid()==ParticleID::gamma) hasPhoton = true; } // loop over channels Energy min = incoming()->mass(); int nchannel(0); pair imin[4]={make_pair(-1,-1.*GeV),make_pair(-1,-1.*GeV), make_pair(-1,-1.*GeV),make_pair(-1,-1.*GeV)}; Energy absmin = -1e20*GeV; int minType = -1; for(unsigned int iy=0;iymass()); Energy dm2(getProcessInfo()[ix].intermediate->mass()); int itype(0); if (getProcessInfo()[ix].channelType==TBDiagram::channel23) { dm1 -= outgoing()[0]->mass(); dm2 -= outgoing()[1]->mass()+outgoing()[2]->mass(); itype = 3; } else if(getProcessInfo()[ix].channelType==TBDiagram::channel13) { dm1 -= outgoing()[1]->mass(); dm2 -= outgoing()[0]->mass()+outgoing()[2]->mass(); itype = 2; } else if(getProcessInfo()[ix].channelType==TBDiagram::channel12) { dm1 -= outgoing()[2]->mass(); dm2 -= outgoing()[0]->mass()+outgoing()[1]->mass(); itype = 1; } if((dm1id()!=ParticleID::gamma) { intype.push_back(itype); inpow.push_back(0.); inmass.push_back(getProcessInfo()[ix].intermediate->mass()); inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[ix].intermediate->width()); ++nchannel; } else if(getProcessInfo()[ix].intermediate->id()==ParticleID::gamma) { intype.push_back(itype); inpow.push_back(-2.); inmass.push_back(-1.*GeV); inwidth.push_back(-1.*GeV); ++nchannel; } } // physical poles, use them and return if(nchannel>0) { inweights = vector(nchannel,1./double(nchannel)); return; } // use shallowest pole else if(intOpt_==1&&minType>0&&getProcessInfo()[imin[minType].first].intermediate->id()!=ParticleID::gamma) { intype.push_back(minType); inpow.push_back(0.); inmass.push_back(getProcessInfo()[imin[minType].first].intermediate->mass()); inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[imin[minType].first].intermediate->width()); inweights = vector(1,1.); return; } for(unsigned int ix=1;ix<4;++ix) { if(imin[ix].first>=0) { intype.push_back(ix); if(getProcessInfo()[imin[ix].first].intermediate->id()!=ParticleID::gamma) { inpow.push_back(0.); inmass.push_back(getProcessInfo()[imin[ix].first].intermediate->mass()); inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[imin[ix].first].intermediate->width()); } else { inpow.push_back(-2.); inmass.push_back(-1.*GeV); inwidth.push_back(-1.*GeV); } ++nchannel; } } inweights = vector(nchannel,1./double(nchannel)); } bool GeneralThreeBodyDecayer::setColourFactors(double symfac) { string name = incoming_->PDGName() + "->"; vector sng,trip,atrip,oct; unsigned int iloc(0); for(vector::const_iterator it = outgoing_.begin(); it != outgoing_.end();++it) { name += (**it).PDGName() + " "; if ((**it).iColour() == PDT::Colour0 ) sng.push_back(iloc) ; else if((**it).iColour() == PDT::Colour3 ) trip.push_back(iloc) ; else if((**it).iColour() == PDT::Colour3bar ) atrip.push_back(iloc); else if((**it).iColour() == PDT::Colour8 ) oct.push_back(iloc) ; ++iloc; } // colour neutral decaying particle if ( incoming_->iColour() == PDT::Colour0) { // options are all neutral or triplet/antitriplet+ neutral if(sng.size()==3) { nflow_ = 1; colour_ = vector(1,DVector(1,1.)); colourLargeNC_ = vector(1,DVector(1,1.)); } else if(sng.size()==1&&trip.size()==1&&atrip.size()==1) { nflow_ = 1; colour_ = vector(1,DVector(1,3.)); colourLargeNC_ = vector(1,DVector(1,3.)); } else if(trip.size()==1&&atrip.size()==1&&oct.size()==1) { nflow_ = 1; colour_ = vector(1,DVector(1,4.)); colourLargeNC_ = vector(1,DVector(1,4.)); } else if( trip.size() == 3 || atrip.size() == 3 ) { nflow_ = 1; colour_ = vector(1,DVector(1,6.)); colourLargeNC_ = vector(1,DVector(1,6.)); for(unsigned int ix=0;ix(1,make_pair(1,sign)); diagrams_[ix].largeNcColourFlow = vector(1,make_pair(1,sign)); } } else { generator()->log() << "Unknown colour flow structure for " << "colour neutral decay " << name << " in GeneralThreeBodyDecayer::" << "setColourFactors(), omitting decay\n"; return false; } } // colour triplet decaying particle else if( incoming_->iColour() == PDT::Colour3) { if(sng.size()==2&&trip.size()==1) { nflow_ = 1; colour_ = vector(1,DVector(1,1.)); colourLargeNC_ = vector(1,DVector(1,1.)); } else if(trip.size()==2&&atrip.size()==1) { nflow_ = 2; colour_.clear(); colour_.resize(2,DVector(2,0.)); colour_[0][0] = 3.; colour_[0][1] = 1.; colour_[1][0] = 1.; colour_[1][1] = 3.; colourLargeNC_.clear(); colourLargeNC_.resize(2,DVector(2,0.)); colourLargeNC_[0][0] = 3.; colourLargeNC_[1][1] = 3.; // sort out the contribution of the different diagrams to the colour // flows for(unsigned int ix=0;ixiColour()==PDT::Colour0) { if(diagrams_[ix].channelType==trip[0]) { diagrams_[ix]. colourFlow = vector(1,make_pair(1,1.)); diagrams_[ix].largeNcColourFlow = vector(1,make_pair(1,1.)); } else { diagrams_[ix].colourFlow = vector(1,make_pair(2,1.)); diagrams_[ix].largeNcColourFlow = vector(1,make_pair(2,1.)); } } // colour octet intermediate else if(diagrams_[ix].intermediate->iColour()==PDT::Colour8) { if(diagrams_[ix].channelType==trip[0]) { vector flow(1,make_pair(2, 0.5 )); diagrams_[ix].largeNcColourFlow = flow; flow.push_back( make_pair(1,-1./6.)); diagrams_[ix].colourFlow=flow; } else { vector flow(1,make_pair(1, 0.5 )); diagrams_[ix].largeNcColourFlow = flow; flow.push_back( make_pair(2,-1./6.)); diagrams_[ix].colourFlow=flow; } } else { generator()->log() << "Unknown colour for the intermediate in " << "triplet -> triplet triplet antitriplet in " << "GeneralThreeBodyDecayer::setColourFactors()" << " for " << name << " omitting decay\n"; return false; } } } else if(trip.size()==1&&oct.size()==1&&sng.size()==1) { nflow_ = 1; colour_ = vector(1,DVector(1,4./3.)); colourLargeNC_ = vector(1,DVector(1,4./3.)); } else if(sng.size()==1&&atrip.size()==2) { nflow_ = 1; colour_ = vector(1,DVector(1,2.)); colourLargeNC_ = vector(1,DVector(1,2.)); } else { generator()->log() << "Unknown colour structure for " << "triplet decay in " << "GeneralThreeBodyDecayer::setColourFactors()" << " for " << name << " omitting decay\n"; return false; } } // colour antitriplet decaying particle else if( incoming_->iColour() == PDT::Colour3bar) { if(sng.size()==2&&atrip.size()==1) { nflow_ = 1; colour_ = vector(1,DVector(1,1.)); colourLargeNC_ = vector(1,DVector(1,1.)); } else if(atrip.size()==2&&trip.size()==1) { nflow_ = 2; colour_.clear(); colour_.resize(2,DVector(2,0.)); colour_[0][0] = 3.; colour_[0][1] = 1.; colour_[1][0] = 1.; colour_[1][1] = 3.; colourLargeNC_.clear(); colourLargeNC_.resize(2,DVector(2,0.)); colourLargeNC_[0][0] = 3.; colourLargeNC_[1][1] = 3.; // sort out the contribution of the different diagrams to the colour // flows for(unsigned int ix=0;ixiColour()==PDT::Colour0) { if(diagrams_[ix].channelType==atrip[0]) { diagrams_[ix]. colourFlow = vector(1,make_pair(1,1.)); diagrams_[ix].largeNcColourFlow = vector(1,make_pair(1,1.)); } else { diagrams_[ix].colourFlow = vector(1,make_pair(2,1.)); diagrams_[ix].largeNcColourFlow = vector(1,make_pair(2,1.)); } } // colour octet intermediate else if(diagrams_[ix].intermediate->iColour()==PDT::Colour8) { if(diagrams_[ix].channelType==atrip[0]) { vector flow(1,make_pair(2, 0.5 )); diagrams_[ix].largeNcColourFlow = flow; flow.push_back( make_pair(1,-1./6.)); diagrams_[ix].colourFlow=flow; } else { vector flow(1,make_pair(1, 0.5 )); diagrams_[ix].largeNcColourFlow = flow; flow.push_back( make_pair(2,-1./6.)); diagrams_[ix].colourFlow=flow; } } else { generator()->log() << "Unknown colour for the intermediate in " << "antitriplet -> antitriplet antitriplet triplet in " << "GeneralThreeBodyDecayer::setColourFactors()" << " for " << name << " omitting decay\n"; return false; } } } else if(atrip.size()==1&&oct.size()==1&&sng.size()==1) { nflow_ = 1; colour_ = vector(1,DVector(1,4./3.)); colourLargeNC_ = vector(1,DVector(1,4./3.)); } else if(sng.size()==1&&trip.size()==2) { nflow_ = 1; colour_ = vector(1,DVector(1,2.)); colourLargeNC_ = vector(1,DVector(1,2.)); } else { generator()->log() << "Unknown colour antitriplet decay in " << "GeneralThreeBodyDecayer::setColourFactors()" << " for " << name << " omitting decay\n"; return false; } } // colour octet particle else if( incoming_->iColour() == PDT::Colour8) { // triplet antitriplet if(trip.size() == 1 && atrip.size() == 1 && sng.size() == 1) { nflow_ = 1; colour_ = vector(1,DVector(1,0.5)); colourLargeNC_ = vector(1,DVector(1,0.5)); } // three (anti)triplets else if(trip.size()==3||atrip.size()==3) { nflow_ = 3; colour_ = vector(3,DVector(3,0.)); colourLargeNC_ = vector(3,DVector(3,0.)); colour_[0][0] = 1.; colour_[1][1] = 1.; colour_[2][2] = 1.; colour_[0][1] = -0.5; colour_[1][0] = -0.5; colour_[0][2] = -0.5; colour_[2][0] = -0.5; colour_[1][2] = -0.5; colour_[2][1] = -0.5; colourLargeNC_ = vector(3,DVector(3,0.)); colourLargeNC_[0][0] = 1.; colourLargeNC_[1][1] = 1.; colourLargeNC_[2][2] = 1.; // sett the factors for the diagrams for(unsigned int ix=0;ixCC()) inter = inter->CC(); unsigned int io[2]={1,2}; double sign = diagrams_[ix].channelType == TBDiagram::channel13 ? -1. : 1.; for(unsigned int iy=0;iy<3;++iy) { if (iy==1) io[0]=0; else if(iy==2) io[1]=1; tPDVector decaylist = diagrams_[ix].vertices.second->search(iy, inter); if(decaylist.empty()) continue; bool found=false; for(unsigned int iz=0;izid()==diagrams_[ix].outgoingPair.first && decaylist[iz+io[1]]->id()==diagrams_[ix].outgoingPair.second) { sign *= 1.; found = true; } else if(decaylist[iz+io[0]]->id()==diagrams_[ix].outgoingPair.second && decaylist[iz+io[1]]->id()==diagrams_[ix].outgoingPair.first ) { sign *= -1.; found = true; } } if(found) { if(iy==1) sign *=-1.; break; } } diagrams_[ix]. colourFlow = vector(1,make_pair(diagrams_[ix].channelType+1,sign)); diagrams_[ix].largeNcColourFlow = vector(1,make_pair(diagrams_[ix].channelType+1,sign)); } } // unknown else { generator()->log() << "Unknown colour octet decay in " << "GeneralThreeBodyDecayer::setColourFactors()" << " for " << name << " omitting decay\n"; return false; } } else if (incoming_->iColour() == PDT::Colour6 ) { generator()->log() << "Unknown colour sextet decay in " << "GeneralThreeBodyDecayer::setColourFactors()" << " for " << name << " omitting decay\n"; return false; } else if (incoming_->iColour() == PDT::Colour6bar ) { generator()->log() << "Unknown colour anti-sextet decay in " << "GeneralThreeBodyDecayer::setColourFactors()" << " for " << name << " omitting decay\n"; return false; } assert(nflow_ != 999); for(unsigned int ix=0;ix 1 ) { generator()->log() << "Mode: " << name << " has colour factors\n"; for(unsigned int ix=0;ixlog() << colour_[ix][iy] << " "; } generator()->log() << "\n"; } for(unsigned int ix=0;ixlog() << "colour flow for diagram : " << ix; for(unsigned int iy=0;iylog() << "(" << diagrams_[ix].colourFlow[iy].first << "," << diagrams_[ix].colourFlow[iy].second << "); "; generator()->log() << "\n"; } } return true; } + +void GeneralThreeBodyDecayer::setupDiagrams(bool kinCheck) { + clearModes(); + // create the phase space integrator + tPDVector extpart(1,incoming_); + extpart.insert(extpart.end(),outgoing_.begin(),outgoing_.end()); + // create the integration channels for the decay + DecayPhaseSpaceModePtr mode(new_ptr(DecayPhaseSpaceMode(extpart,this,true))); + DecayPhaseSpaceChannelPtr newchannel; + // create the phase-space channels for the integration + unsigned int nmode(0); + unsigned int idiag(0); + for(vector::iterator it = diagrams_.begin();it!=diagrams_.end();++it) { + if(it->channelType==TBDiagram::fourPoint|| + it->channelType==TBDiagram::UNDEFINED) { + idiag+=1; + continue; + } + // create the new channel + newchannel=new_ptr(DecayPhaseSpaceChannel(mode)); + int jac = 0; + double power = 0.0; + if ( it->intermediate->mass() == ZERO || + it->intermediate->width() == ZERO ) { + jac = 1; + power = -2.0; + } + if(it->channelType==TBDiagram::channel23) { + newchannel->addIntermediate(extpart[0],0,0.0,-1,1); + newchannel->addIntermediate(it->intermediate,jac,power, 2,3); + } + else if(it->channelType==TBDiagram::channel13) { + newchannel->addIntermediate(extpart[0],0,0.0,-1,2); + newchannel->addIntermediate(it->intermediate,jac,power, 1,3); + } + else if(it->channelType==TBDiagram::channel12) { + newchannel->addIntermediate(extpart[0],0,0.0,-1,3); + newchannel->addIntermediate(it->intermediate,jac,power, 1,2); + } + newchannel->init(); + if(kinCheck&&!newchannel->checkKinematics()) { + generator()->log() << "Erasing diagram " + << *it + << "from three body decay as zero width propagator can be on-shell,\n" + << "hopefully this diagram is zero in your model, but you should check this\n"; + it = diagrams_.erase(it); + continue; + } + diagmap_.push_back(idiag); + mode->addChannel(newchannel); + ++nmode; + ++idiag; + } + if(nmode==0) { + string mode = extpart[0]->PDGName() + "->"; + for(unsigned int ix=1;ixPDGName() + " "; + throw Exception() << "No decay channels in GeneralThreeBodyDecayer::" + << "doinit() for " << mode << "\n" << Exception::runerror; + } + // // add the mode + vector wgt(nmode,1./double(nmode)); + addMode(mode,1.,wgt); +} diff --git a/Decay/General/GeneralThreeBodyDecayer.h b/Decay/General/GeneralThreeBodyDecayer.h --- a/Decay/General/GeneralThreeBodyDecayer.h +++ b/Decay/General/GeneralThreeBodyDecayer.h @@ -1,321 +1,332 @@ // -*- C++ -*- #ifndef HERWIG_GeneralThreeBodyDecayer_H #define HERWIG_GeneralThreeBodyDecayer_H // // This is the declaration of the GeneralThreeBodyDecayer class. // #include "Herwig/Decay/DecayIntegrator.h" #include "Herwig/Models/General/TBDiagram.h" #include "GeneralThreeBodyDecayer.fh" namespace Herwig { using namespace ThePEG; /** * Here is the documentation of the GeneralThreeBodyDecayer class. * * @see \ref GeneralThreeBodyDecayerInterfaces "The interfaces" * defined for GeneralThreeBodyDecayer. */ class GeneralThreeBodyDecayer: public DecayIntegrator { public: /** A ParticleData ptr and (possible) mass pair.*/ typedef pair PMPair; public: /** * The default constructor. */ GeneralThreeBodyDecayer() : nflow_(999), widthOpt_(1), refTag_(), refTagCC_(), iflow_(999), intOpt_(0), relerr_(1e-2) {} /** @name Virtual functions required by the Decayer class. */ //@{ /** * For a given decay mode and a given particle instance, perform the * decay and return the decay products. As this is the base class this * is not implemented. * @return The vector of particles produced in the decay. */ virtual ParticleVector decay(const Particle & parent, const tPDVector & children) const; /** * Which of the possible decays is required * @param cc Is this mode the charge conjugate * @param parent The decaying particle * @param children The decay products */ virtual int modeNumber(bool & cc, tcPDPtr parent,const tPDVector & children) const; /** * The matrix element to be integrated for the three-body decays as a function * of the invariant masses of pairs of the outgoing particles. * @param imode The mode for which the matrix element is needed. * @param q2 The scale, \e i.e. the mass squared of the decaying particle. * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$. * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$. * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$. * @param m1 The mass of the first outgoing particle. * @param m2 The mass of the second outgoing particle. * @param m3 The mass of the third outgoing particle. * @return The matrix element */ virtual double threeBodyMatrixElement(const int imode, const Energy2 q2, const Energy2 s3, const Energy2 s2, const Energy2 s1, const Energy m1, const Energy m2, const Energy m3) const; /** * Function to return partial Width * @param inpart The decaying particle. * @param outa First decay product. * @param outb Second decay product. * @param outc Third decay product. */ virtual Energy partialWidth(PMPair inpart, PMPair outa, PMPair outb, PMPair outc) const; /** * An overidden member to calculate a branching ratio for a certain * particle instance. * @param dm The DecayMode of the particle * @param p The particle object * @param oldbrat The branching fraction given in the DecayMode object */ virtual double brat(const DecayMode & dm, const Particle & p, double oldbrat) const; //@} /** * Set the diagrams */ bool setDecayInfo(PDPtr incoming,vector outgoing, const vector & process, double symfac); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); + + /** + * Initialize this object. Called in the run phase just before + * a run begins. + */ + virtual void doinitrun(); //@} + /** + * Set up the diagrams etc + */ + virtual void setupDiagrams(bool checkKinematics); + protected: /** * Access the TBDiagrams that store the required information * to create the diagrams */ const vector & getProcessInfo() const { return diagrams_; } /** * Incoming particle */ PDPtr incoming() const { return incoming_; } /** * Outgoing particles */ const vector & outgoing() const { return outgoing_; } /** * Number of colour flows */ unsigned int numberOfFlows() const { return nflow_; } /** * Set up the colour factors */ bool setColourFactors(double symfac); /** * Return the matrix of colour factors */ const vector & getColourFactors() const { return colour_; } /** * Return the matrix of colour factors */ const vector & getLargeNcColourFactors() const { return colourLargeNC_; } /** * Get the mapping between the phase-space channel and the diagram */ const vector & diagramMap() const { return diagmap_; } /** * Option for the handling of the widths of the intermediate particles */ unsigned int widthOption() const { return widthOpt_; } /** * Set colour connections * @param parent Parent particle * @param out Particle vector containing particles to * connect colour lines */ void colourConnections(const Particle & parent, const ParticleVector & out) const; /** * Method to construct the channels for the integrator to give the partial width * @param intype Types of the channels * @param inmass Mass for the channels * @param inwidth Width for the channels * @param inpow Power for the channels * @param inweights Weights for the channels */ void constructIntegratorChannels(vector & intype, vector & inmass, vector & inwidth, vector & inpow, vector & inweights) const; /** * Set the colour flow * @param flow The value for the colour flow */ void colourFlow(unsigned int flow) const { iflow_ = flow; } /** * Set the colour flow */ unsigned int const & colourFlow() const { return iflow_; } /** * Relative error for GQ integration */ double relativeError() const {return relerr_;} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ GeneralThreeBodyDecayer & operator=(const GeneralThreeBodyDecayer &); private: /** * Store the incoming particle */ PDPtr incoming_; /** * Outgoing particles */ vector outgoing_; /** * Store the diagrams for the decay */ vector diagrams_; /** * Map between the diagrams and the phase-space channels */ vector diagmap_; /** * Store colour factors for ME calc. */ vector colour_; /** * Store cololur factors for ME calc at large N_c */ vector colourLargeNC_; /** * The number of colourflows. */ unsigned int nflow_; /** * Reference to object to calculate the partial width */ mutable WidthCalculatorBasePtr widthCalc_; /** * Option for the treatment of the widths */ unsigned int widthOpt_; /** * Store a decay tag for this mode that can be tested when * trying to determine whether it can be generated by * this Decayer */ string refTag_; /** * Store a decay tag for the cc-mode that can be tested when * trying to determine whether it can be generated by * this Decayer */ string refTagCC_; /** * The colour flow */ mutable unsigned int iflow_; /** * Option for the construction of the gaussian integrator */ unsigned int intOpt_; /** * Relative error for GQ integration of partial width */ double relerr_; }; } #endif /* HERWIG_GeneralThreeBodyDecayer_H */ diff --git a/Decay/General/StoFFVDecayer.cc b/Decay/General/StoFFVDecayer.cc --- a/Decay/General/StoFFVDecayer.cc +++ b/Decay/General/StoFFVDecayer.cc @@ -1,389 +1,389 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the StoFFVDecayer class. // #include "StoFFVDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/PDT/ThreeBodyAllOnCalculator.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include using namespace Herwig; using namespace ThePEG; using namespace ThePEG::Helicity; IBPtr StoFFVDecayer::clone() const { return new_ptr(*this); } IBPtr StoFFVDecayer::fullclone() const { return new_ptr(*this); } void StoFFVDecayer::persistentOutput(PersistentOStream & os) const { os << sca_ << fer_ << vec_ << RSfer_ << four_; } void StoFFVDecayer::persistentInput(PersistentIStream & is, int) { is >> sca_ >> fer_ >> vec_ >> RSfer_ >> four_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigStoFFVDecayer("Herwig::StoFFVDecayer", "Herwig.so"); void StoFFVDecayer::Init() { static ClassDocumentation documentation ("The StoFFVDecayer class implements the general decay of a scalar to " "a two fermions and a vector."); } WidthCalculatorBasePtr StoFFVDecayer:: threeBodyMEIntegrator(const DecayMode & ) const { vector intype; vector inmass,inwidth; vector inpow,inweights; constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights); return new_ptr(ThreeBodyAllOnCalculator (inweights,intype,inmass,inwidth,inpow,*this,0, outgoing()[0]->mass(),outgoing()[1]->mass(), outgoing()[2]->mass(),relativeError())); } -void StoFFVDecayer::doinit() { - GeneralThreeBodyDecayer::doinit(); +void StoFFVDecayer::setupDiagrams(bool kinCheck) { + GeneralThreeBodyDecayer::setupDiagrams(kinCheck); if(outgoing().empty()) return; unsigned int ndiags = getProcessInfo().size(); sca_.resize(ndiags); fer_.resize(ndiags); RSfer_.resize(ndiags); vec_.resize(ndiags); four_.resize(ndiags); for(unsigned int ix = 0;ix < ndiags; ++ix) { TBDiagram current = getProcessInfo()[ix]; tcPDPtr offshell = current.intermediate; // four point vertex if(!offshell) { four_[ix] = dynamic_ptr_cast(current.vertices.first); continue; } if( offshell->CC() ) offshell = offshell->CC(); if(offshell->iSpin() == PDT::Spin0) { AbstractVSSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFSVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a scalar diagram in StoFFVDecayer::doinit()" + << "Invalid vertices for a scalar diagram in StoFFVDecayer::setupDiagrams()" << Exception::runerror; sca_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin1Half) { AbstractFFSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFVVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a fermion diagram in StoFFVDecayer::doinit()" + << "Invalid vertices for a fermion diagram in StoFFVDecayer::setupDiagrams()" << Exception::runerror; fer_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin1) { AbstractVVSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFVVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a vector diagram in StoFFVDecayer::doinit()" + << "Invalid vertices for a vector diagram in StoFFVDecayer::setupDiagrams()" << Exception::runerror; vec_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin3Half) { AbstractRFSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractRFVVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a RS fermion diagram in StoFFVDecayer::doinit()" + << "Invalid vertices for a RS fermion diagram in StoFFVDecayer::setupDiagrams()" << Exception::runerror; RSfer_[ix] = make_pair(vert1, vert2); } } } double StoFFVDecayer::me2(const int ichan, const Particle & inpart, const ParticleVector & decay, MEOption meopt) const { // particle or CC of particle bool cc = (*getProcessInfo().begin()).incoming != inpart.id(); // special handling or first/last call if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(rho_,const_ptr_cast(&inpart), Helicity::incoming); swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(), Helicity::incoming); // fix rho if no correlations fixRho(rho_); } if(meopt==Terminate) { ScalarWaveFunction:: constructSpinInfo(const_ptr_cast(&inpart), Helicity::incoming,true); for(unsigned int ix=0;ixdataPtr()->iSpin()==PDT::Spin1) { VectorWaveFunction::constructSpinInfo(outVector_,decay[ix], Helicity::outgoing,true,false); } else { SpinorWaveFunction:: constructSpinInfo(outspin_[ix].first,decay[ix],Helicity::outgoing,true); } } } unsigned int ivec(0); bool massless(false); for(unsigned int ix = 0; ix < decay.size();++ix) { if(decay[ix]->dataPtr()->iSpin() == PDT::Spin1) { ivec = ix; massless = decay[ivec]->mass()==ZERO; VectorWaveFunction:: calculateWaveFunctions(outVector_, decay[ix], Helicity::outgoing,massless); } else { SpinorWaveFunction:: calculateWaveFunctions(outspin_[ix].first,decay[ix],Helicity::outgoing); outspin_[ix].second.resize(2); // Need a ubar and a v spinor if(outspin_[ix].first[0].wave().Type() == SpinorType::u) { for(unsigned int iy = 0; iy < 2; ++iy) { outspin_[ix].second[iy] = outspin_[ix].first[iy].bar(); outspin_[ix].first[iy].conjugate(); } } else { for(unsigned int iy = 0; iy < 2; ++iy) { outspin_[ix].second[iy] = outspin_[ix].first[iy].bar(); outspin_[ix].second[iy].conjugate(); } } } } const vector > cfactors(getColourFactors()); const vector > nfactors(getLargeNcColourFactors()); Energy2 scale(sqr(inpart.mass())); const size_t ncf(numberOfFlows()); vector flows(ncf, Complex(0.)), largeflows(ncf, Complex(0.)); // setup the DecayMatrixElement vector mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin0, ivec == 0 ? PDT::Spin1 : PDT::Spin1Half, ivec == 1 ? PDT::Spin1 : PDT::Spin1Half, ivec == 2 ? PDT::Spin1 : PDT::Spin1Half))); vector mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin0, ivec == 0 ? PDT::Spin1 : PDT::Spin1Half, ivec == 1 ? PDT::Spin1 : PDT::Spin1Half, ivec == 2 ? PDT::Spin1 : PDT::Spin1Half))); //the channel possiblities static const unsigned int out2[3] = {1,0,0}, out3[3] = {2,2,1}; for(unsigned int s1 = 0; s1 < 2; ++s1) { for(unsigned int s2 = 0; s2 < 2; ++s2) { for(unsigned int v1 = 0; v1 < 3; ++v1) { if(massless&&v1==1) ++v1; flows = vector(ncf, Complex(0.)); largeflows = vector(ncf, Complex(0.)); unsigned int idiag(0); Complex diag; for(vector::const_iterator dit=getProcessInfo().begin(); dit!=getProcessInfo().end();++dit) { // channels if selecting if( ichan >= 0 && diagramMap()[ichan] != idiag ) { ++idiag; continue; } tcPDPtr offshell = dit->intermediate; if(offshell) { if(cc&&offshell->CC()) offshell=offshell->CC(); unsigned int o2(out2[dit->channelType]), o3(out3[dit->channelType]); double sign = (o3 < o2) ? 1. : -1.; // intermediate scalar if(offshell->iSpin() == PDT::Spin0) { ScalarWaveFunction inters = sca_[idiag].first-> evaluate(scale, widthOption(), offshell, outVector_[v1], swave_); unsigned int h1(s1),h2(s2); if(o2 > o3) swap(h1, h2); if(decay[o2]->id() < 0 && decay[o3]->id() > 0) { diag = -sign*sca_[idiag].second-> evaluate(scale,outspin_[o2].first[h1], outspin_[o3].second[h2],inters); } else { diag = sign*sca_[idiag].second-> evaluate(scale, outspin_[o3].first [h2], outspin_[o2].second[h1],inters); } } // intermediate fermion else if(offshell->iSpin() == PDT::Spin1Half) { int iferm = (decay[o2]->dataPtr()->iSpin() == PDT::Spin1Half) ? o2 : o3; unsigned int h1(s1),h2(s2); if(dit->channelType > iferm) swap(h1, h2); sign = iferm < dit->channelType ? 1. : -1.; if((decay[dit->channelType]->id() < 0 && decay[iferm]->id() > 0 ) || (decay[dit->channelType]->id()*offshell->id()>0)) { SpinorWaveFunction inters = fer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].first[h1], swave_); diag = -sign*fer_[idiag].second-> evaluate(scale,inters,outspin_[iferm].second[h2], outVector_[v1]); } else { SpinorBarWaveFunction inters = fer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].second[h1],swave_); diag = sign*fer_[idiag].second-> evaluate(scale,outspin_[iferm].first [h2],inters, outVector_[v1]); } } // intermediate vector else if(offshell->iSpin() == PDT::Spin1) { VectorWaveFunction interv = vec_[idiag].first-> evaluate(scale, widthOption(), offshell, outVector_[v1], swave_); unsigned int h1(s1),h2(s2); if(o2 > o3) swap(h1,h2); if(decay[o2]->id() < 0 && decay[o3]->id() > 0) { diag =-sign*vec_[idiag].second-> evaluate(scale, outspin_[o2].first[h1], outspin_[o3].second[h2], interv); } else { diag = sign*vec_[idiag].second-> evaluate(scale, outspin_[o3].first[h2], outspin_[o2].second[h1], interv); } } // intermediate RS fermion else if(offshell->iSpin() == PDT::Spin3Half) { int iferm = (decay[o2]->dataPtr()->iSpin() == PDT::Spin1Half) ? o2 : o3; unsigned int h1(s1),h2(s2); if(dit->channelType > iferm) swap(h1, h2); sign = iferm < dit->channelType ? 1. : -1.; if((decay[dit->channelType]->id() < 0 && decay[iferm]->id() > 0 ) || (decay[dit->channelType]->id()*offshell->id()>0)) { RSSpinorWaveFunction inters = RSfer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].first[h1], swave_); diag = -sign*RSfer_[idiag].second-> evaluate(scale,inters,outspin_[iferm].second[h2], outVector_[v1]); } else { RSSpinorBarWaveFunction inters = RSfer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].second[h1],swave_); diag = sign*RSfer_[idiag].second-> evaluate(scale,outspin_[iferm].first [h2],inters, outVector_[v1]); } } // unknown else throw Exception() << "Unknown intermediate in StoFFVDecayer::me2()" << Exception::runerror; } else { unsigned int o2 = ivec > 0 ? 0 : 1; unsigned int o3 = ivec < 2 ? 2 : 1; if(decay[o2]->id() < 0 && decay[o3]->id() > 0) { diag =-four_[idiag]-> evaluate(scale, outspin_[o2].first[s1], outspin_[o3].second[s2], outVector_[v1], swave_); } else { diag = four_[idiag]-> evaluate(scale, outspin_[o3].first[s2], outspin_[o2].second[s1], outVector_[v1], swave_); } } // matrix element for the different colour flows if(ichan < 0) { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } else { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1 != colourFlow()) continue; flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1!=colourFlow()) continue; largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } ++idiag; } //end of diagrams // now add the flows to the me2 with appropriate colour factors for(unsigned int ix = 0; ix < ncf; ++ix) { if ( ivec == 0 ) { (*mes[ix])(0, v1, s1, s2) = flows[ix]; (*mel[ix])(0, v1, s1, s2) = largeflows[ix]; } else if( ivec == 1 ) { (*mes[ix])(0, s1, v1, s2) = flows[ix]; (*mel[ix])(0, s1, v1, s2) = largeflows[ix]; } else if( ivec == 2 ) { (*mes[ix])(0, s1, s2, v1) = flows[ix]; (*mel[ix])(0, s1, s2, v1) = largeflows[ix]; } } } } } double me2(0.); if(ichan < 0) { vector pflows(ncf,0.); for(unsigned int ix = 0; ix < ncf; ++ix) { for(unsigned int iy = 0; iy < ncf; ++ iy) { double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real(); me2 += con; if(ix == iy) { con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real(); pflows[ix] += con; } } } double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.)); ptotal *= UseRandom::rnd(); for(unsigned int ix = 0;ix < pflows.size(); ++ix) { if(ptotal <= pflows[ix]) { colourFlow(ix); ME(mes[ix]); break; } ptotal -= pflows[ix]; } } else { unsigned int iflow = colourFlow(); me2 = nfactors[iflow][iflow]*(mel[iflow]->contract(*mel[iflow],rho_)).real(); } // return the matrix element squared return me2; } diff --git a/Decay/General/StoFFVDecayer.h b/Decay/General/StoFFVDecayer.h --- a/Decay/General/StoFFVDecayer.h +++ b/Decay/General/StoFFVDecayer.h @@ -1,162 +1,157 @@ // -*- C++ -*- #ifndef THEPEG_StoFFVDecayer_H #define THEPEG_StoFFVDecayer_H // // This is the declaration of the StoFFVDecayer class. // #include "GeneralThreeBodyDecayer.h" #include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h" #include "ThePEG/Helicity/Vertex/AbstractRFVVertex.h" #include "ThePEG/Helicity/Vertex/AbstractRFSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFVSVertex.h" namespace Herwig { using namespace ThePEG; /** * Here is the documentation of the StoFFVDecayer class. * * @see \ref StoFFVDecayerInterfaces "The interfaces" * defined for StoFFVDecayer. */ class StoFFVDecayer: public GeneralThreeBodyDecayer { public: /** * Return the matrix element squared for a given mode and phase-space channel * @param ichan The channel we are calculating the matrix element for. * @param part The decaying Particle. * @param decay The particles produced in the decay. * @param meopt Option for the calculation of the matrix element * @return The matrix element squared for the phase-space configuration. */ virtual double me2(const int ichan, const Particle & part, const ParticleVector & decay, MEOption meopt) const; /** * Method to return an object to calculate the 3 (or higher body) partial width * @param dm The DecayMode * @return A pointer to a WidthCalculatorBase object capable of * calculating the width */ virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} 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. + * Set up the diagrams etc */ - virtual void doinit(); - //@} + virtual void setupDiagrams(bool checkKinematics); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ StoFFVDecayer & operator=(const StoFFVDecayer &); private: /** * Store the vertices for fermion intrermediate */ vector > fer_; /** * Store the vertices for fermion intrermediate */ vector > RSfer_; /** * Store the vertices for scalar intrermediate */ vector > sca_; /** * Store the vertices for vector intrermediate */ vector > vec_; /** * Store the vertices for 4-point diagrams */ vector four_; /** * Spin density matrix */ mutable RhoDMatrix rho_; /** * Scalar wavefunction */ mutable ScalarWaveFunction swave_; /** * Vector wavefunction */ mutable vector outVector_; /** * Spinor wavefunctions */ mutable pair,vector > outspin_[3]; }; } #endif /* THEPEG_StoFFVDecayer_H */ diff --git a/Decay/General/StoSFFDecayer.cc b/Decay/General/StoSFFDecayer.cc --- a/Decay/General/StoSFFDecayer.cc +++ b/Decay/General/StoSFFDecayer.cc @@ -1,424 +1,424 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the StoSFFDecayer class. // #include "StoSFFDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/PDT/ThreeBodyAllOnCalculator.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include using namespace Herwig; using namespace ThePEG; using namespace ThePEG::Helicity; IBPtr StoSFFDecayer::clone() const { return new_ptr(*this); } IBPtr StoSFFDecayer::fullclone() const { return new_ptr(*this); } void StoSFFDecayer::persistentOutput(PersistentOStream & os) const { os << sca_ << fer_ << vec_ << ten_ << RSfer_ << four_; } void StoSFFDecayer::persistentInput(PersistentIStream & is, int) { is >> sca_ >> fer_ >> vec_ >> ten_ >> RSfer_ >> four_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigStoSFFDecayer("Herwig::StoSFFDecayer", "Herwig.so"); void StoSFFDecayer::Init() { static ClassDocumentation documentation ("The StoSFFDecayer class implements the general decay of a scalar to " "a scalar and two fermions."); } WidthCalculatorBasePtr StoSFFDecayer:: threeBodyMEIntegrator(const DecayMode & ) const { vector intype; vector inmass,inwidth; vector inpow,inweights; constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights); return new_ptr(ThreeBodyAllOnCalculator (inweights,intype,inmass,inwidth,inpow,*this,0, outgoing()[0]->mass(),outgoing()[1]->mass(),outgoing()[2]->mass(), relativeError())); } -void StoSFFDecayer::doinit() { - GeneralThreeBodyDecayer::doinit(); +void StoSFFDecayer::setupDiagrams(bool kinCheck) { + GeneralThreeBodyDecayer::setupDiagrams(kinCheck); if(outgoing().empty()) return; unsigned int ndiags = getProcessInfo().size(); sca_.resize(ndiags); fer_.resize(ndiags); RSfer_.resize(ndiags); vec_.resize(ndiags); ten_.resize(ndiags); four_.resize(ndiags); for(unsigned int ix = 0;ix < ndiags; ++ix) { TBDiagram current = getProcessInfo()[ix]; tcPDPtr offshell = current.intermediate; // four point vertex if(!offshell) { four_[ix] = dynamic_ptr_cast(current.vertices.first); continue; } if( offshell->CC() ) offshell = offshell->CC(); if(offshell->iSpin() == PDT::Spin0) { AbstractSSSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFSVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a scalar diagram in StoSFFDecayer::doinit()" + << "Invalid vertices for a scalar diagram in StoSFFDecayer::setupDiagrams()" << Exception::runerror; sca_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin1Half) { AbstractFFSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFSVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a fermion diagram in StoSFFDecayer::doinit()" + << "Invalid vertices for a fermion diagram in StoSFFDecayer::setupDiagrams()" << Exception::runerror; fer_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin1) { AbstractVSSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFVVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a vector diagram in StoSFFDecayer::doinit()" + << "Invalid vertices for a vector diagram in StoSFFDecayer::setupDiagrams()" << Exception::runerror; vec_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin2) { AbstractSSTVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFTVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a tensor diagram in StoSFFDecayer::doinit()" + << "Invalid vertices for a tensor diagram in StoSFFDecayer::setupDiagrams()" << Exception::runerror; ten_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin3Half) { AbstractRFSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractRFSVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a RS fermion diagram in StoSFFDecayer::doinit()" + << "Invalid vertices for a RS fermion diagram in StoSFFDecayer::setupDiagrams()" << Exception::runerror; RSfer_[ix] = make_pair(vert1, vert2); } } } double StoSFFDecayer::me2(const int ichan, const Particle & inpart, const ParticleVector & decay, MEOption meopt) const { // particle or CC of particle bool cc = (*getProcessInfo().begin()).incoming != inpart.id(); // special handling or first/last call if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(rho_,const_ptr_cast(&inpart), Helicity::incoming); swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(), Helicity::incoming); // fix rho if no correlations fixRho(rho_); } if(meopt==Terminate) { ScalarWaveFunction:: constructSpinInfo(const_ptr_cast(&inpart), Helicity::incoming,true); for(unsigned int ix=0;ixdataPtr()->iSpin()==PDT::Spin0) { ScalarWaveFunction::constructSpinInfo(decay[ix],Helicity::outgoing,true); } else { SpinorWaveFunction:: constructSpinInfo(outspin_[ix].first,decay[ix],Helicity::outgoing,true); } } return 0.; } // get the wavefunctions for all the particles ScalarWaveFunction outScalar; unsigned int isca(0); for(unsigned int ix=0;ixdataPtr()->iSpin()==PDT::Spin0) { isca = ix; outScalar = ScalarWaveFunction(decay[ix]->momentum(), decay[ix]->dataPtr(),Helicity::outgoing); } else { SpinorWaveFunction:: calculateWaveFunctions(outspin_[ix].first,decay[ix],Helicity::outgoing); outspin_[ix].second.resize(2); if(outspin_[ix].first[0].wave().Type() == SpinorType::u) { for(unsigned int iy = 0; iy < 2; ++iy) { outspin_[ix].second[iy] = outspin_[ix].first[iy].bar(); outspin_[ix].first[iy].conjugate(); } } else { for(unsigned int iy = 0; iy < 2; ++iy) { outspin_[ix].second[iy] = outspin_[ix].first[iy].bar(); outspin_[ix].second[iy].conjugate(); } } } } const vector > cfactors(getColourFactors()); const vector > nfactors(getLargeNcColourFactors()); Energy2 scale(sqr(inpart.mass())); const size_t ncf(numberOfFlows()); vector flows(ncf, Complex(0.)), largeflows(ncf, Complex(0.)); vector mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin0, isca==0 ? PDT::Spin0 : PDT::Spin1Half, isca==1 ? PDT::Spin0 : PDT::Spin1Half, isca==2 ? PDT::Spin0 : PDT::Spin1Half))); vector mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin0, isca == 0 ? PDT::Spin0 : PDT::Spin1Half, isca == 1 ? PDT::Spin0 : PDT::Spin1Half, isca == 2 ? PDT::Spin0 : PDT::Spin1Half))); static const unsigned int out2[3]={1,0,0},out3[3]={2,2,1}; for(unsigned int s1 = 0;s1 < 2; ++s1) { for(unsigned int s2 = 0;s2 < 2; ++s2) { flows = vector(ncf, Complex(0.)); largeflows = vector(ncf, Complex(0.)); unsigned int idiag(0); for(vector::const_iterator dit = getProcessInfo().begin(); dit != getProcessInfo().end(); ++dit) { // channels if selecting if( ichan >= 0 && diagramMap()[ichan] != idiag ) { ++idiag; continue; } tcPDPtr offshell = dit->intermediate; Complex diag; if(offshell) { if(cc&&offshell->CC()) offshell=offshell->CC(); double sign = out3[dit->channelType] < out2[dit->channelType] ? 1. : -1.; // intermediate scalar if (offshell->iSpin() == PDT::Spin0) { ScalarWaveFunction inters = sca_[idiag].first-> evaluate(scale, widthOption(), offshell, swave_, outScalar); unsigned int h1(s1),h2(s2); if(out2[dit->channelType]>out3[dit->channelType]) swap(h1,h2); if(decay[out2[dit->channelType]]->id()<0&& decay[out3[dit->channelType]]->id()>0) { diag =-sign*sca_[idiag].second-> evaluate(scale, outspin_[out2[dit->channelType]].first [h1], outspin_[out3[dit->channelType]].second[h2],inters); } else { diag = sign*sca_[idiag].second-> evaluate(scale, outspin_[out3[dit->channelType]].first [h2], outspin_[out2[dit->channelType]].second[h1],inters); } } // intermediate fermion else if(offshell->iSpin() == PDT::Spin1Half) { int iferm = decay[out2[dit->channelType]]->dataPtr()->iSpin()==PDT::Spin1Half ? out2[dit->channelType] : out3[dit->channelType]; unsigned int h1(s1),h2(s2); if(dit->channelType>iferm) swap(h1,h2); sign = ifermchannelType ? 1. : -1.; if((decay[dit->channelType]->id() < 0 &&decay[iferm]->id() > 0 ) || (decay[dit->channelType]->id()*offshell->id()>0)) { SpinorWaveFunction inters = fer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].first [h1],swave_); diag = -sign*fer_[idiag].second-> evaluate(scale,inters,outspin_[iferm].second[h2],outScalar); } else { SpinorBarWaveFunction inters = fer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].second[h1],swave_); diag = sign*fer_[idiag].second-> evaluate(scale,outspin_[iferm].first [h2],inters,outScalar); } } // intermediate vector else if(offshell->iSpin() == PDT::Spin1) { VectorWaveFunction interv = vec_[idiag].first-> evaluate(scale, widthOption(), offshell, swave_, outScalar); unsigned int h1(s1),h2(s2); if(out2[dit->channelType]>out3[dit->channelType]) swap(h1,h2); if(decay[out2[dit->channelType]]->id()<0&& decay[out3[dit->channelType]]->id()>0) { diag =-sign*vec_[idiag].second-> evaluate(scale, outspin_[out2[dit->channelType]].first [h1], outspin_[out3[dit->channelType]].second[h2],interv); } else { diag = sign*vec_[idiag].second-> evaluate(scale, outspin_[out3[dit->channelType]].first [h2], outspin_[out2[dit->channelType]].second[h1],interv); } } // intermediate tensor else if(offshell->iSpin() == PDT::Spin2) { TensorWaveFunction intert = ten_[idiag].first-> evaluate(scale, widthOption(), offshell, swave_, outScalar); unsigned int h1(s1),h2(s2); if(out2[dit->channelType]>out3[dit->channelType]) swap(h1,h2); if(decay[out2[dit->channelType]]->id()<0&& decay[out3[dit->channelType]]->id()>0) { diag =-sign*ten_[idiag].second-> evaluate(scale, outspin_[out2[dit->channelType]].first [h1], outspin_[out3[dit->channelType]].second[h2],intert); } else { diag = sign*ten_[idiag].second-> evaluate(scale, outspin_[out3[dit->channelType]].first [h2], outspin_[out2[dit->channelType]].second[h1],intert); } } // intermediate RS fermion else if(offshell->iSpin() == PDT::Spin3Half) { int iferm = decay[out2[dit->channelType]]->dataPtr()->iSpin()==PDT::Spin1Half ? out2[dit->channelType] : out3[dit->channelType]; unsigned int h1(s1),h2(s2); if(dit->channelType>iferm) swap(h1,h2); sign = ifermchannelType ? 1. : -1.; if((decay[dit->channelType]->id() < 0 &&decay[iferm]->id() > 0 ) || (decay[dit->channelType]->id()*offshell->id()>0)) { RSSpinorWaveFunction inters = RSfer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].first [h1],swave_); diag = -sign*RSfer_[idiag].second-> evaluate(scale,inters,outspin_[iferm].second[h2],outScalar); } else { RSSpinorBarWaveFunction inters = RSfer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].second[h1],swave_); diag = sign*RSfer_[idiag].second-> evaluate(scale,outspin_[iferm].first [h2],inters,outScalar); } } // unknown else throw Exception() << "Unknown intermediate in StoSFFDecayer::me2()" << Exception::runerror; } // four point diagram else { unsigned int o2 = isca > 0 ? 0 : 1; unsigned int o3 = isca < 2 ? 2 : 1; if(decay[o2]->id() < 0 && decay[o3]->id() > 0) { diag =-four_[idiag]-> evaluate(scale, outspin_[o2].first[s1], outspin_[o3].second[s2], outScalar, swave_); } else { diag = four_[idiag]-> evaluate(scale, outspin_[o3].first[s2], outspin_[o2].second[s1], outScalar, swave_); } } // matrix element for the different colour flows if(ichan < 0) { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } else { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1 != colourFlow()) continue; flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1!=colourFlow()) continue; largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } ++idiag; } for(unsigned int ix = 0; ix < ncf; ++ix) { if(isca == 0) { (*mes[ix])(0, 0, s1, s2) = flows[ix]; (*mel[ix])(0, 0, s1, s2) = largeflows[ix]; } else if(isca == 1 ) { (*mes[ix])(0, s1, 0, s2) = flows[ix]; (*mel[ix])(0, s1, 0, s2) = largeflows[ix]; } else if(isca == 2) { (*mes[ix])(0, s1,s2, 0) = flows[ix]; (*mel[ix])(0, s1,s2, 0) = largeflows[ix] ; } } } } double me2(0.); if(ichan < 0) { vector pflows(ncf,0.); for(unsigned int ix = 0; ix < ncf; ++ix) { for(unsigned int iy = 0; iy < ncf; ++ iy) { double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real(); me2 += con; if(ix == iy) { con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real(); pflows[ix] += con; } } } double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.)); ptotal *= UseRandom::rnd(); for(unsigned int ix = 0;ix < pflows.size(); ++ix) { if(ptotal <= pflows[ix]) { colourFlow(ix); ME(mes[ix]); break; } ptotal -= pflows[ix]; } } else { unsigned int iflow = colourFlow(); me2 = nfactors[iflow][iflow]*(mel[iflow]->contract(*mel[iflow],rho_)).real(); } // return the matrix element squared return me2; } diff --git a/Decay/General/StoSFFDecayer.h b/Decay/General/StoSFFDecayer.h --- a/Decay/General/StoSFFDecayer.h +++ b/Decay/General/StoSFFDecayer.h @@ -1,163 +1,158 @@ // -*- C++ -*- #ifndef THEPEG_StoSFFDecayer_H #define THEPEG_StoSFFDecayer_H // // This is the declaration of the StoSFFDecayer class. // #include "GeneralThreeBodyDecayer.h" #include "ThePEG/Helicity/Vertex/AbstractSSSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractRFSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h" #include "ThePEG/Helicity/Vertex/AbstractSSTVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFSSVertex.h" namespace Herwig { using namespace ThePEG; /** * The StoSFFDecayer class provides the general matrix element for * scalar decays into another scalar and a fermion-antifermion pair. * * @see \ref StoSFFDecayerInterfaces "The interfaces" * defined for StoSFFDecayer. */ class StoSFFDecayer: public GeneralThreeBodyDecayer { public: /** * Return the matrix element squared for a given mode and phase-space channel * @param ichan The channel we are calculating the matrix element for. * @param part The decaying Particle. * @param decay The particles produced in the decay. * @param meopt Option for the calculation of the matrix element * @return The matrix element squared for the phase-space configuration. */ virtual double me2(const int ichan, const Particle & part, const ParticleVector & decay, MEOption meopt) const; /** * Method to return an object to calculate the 3 (or higher body) partial width * @param dm The DecayMode * @return A pointer to a WidthCalculatorBase object capable of calculating the width */ virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} 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. + * Set up the diagrams etc */ - virtual void doinit(); - //@} + virtual void setupDiagrams(bool checkKinematics); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ StoSFFDecayer & operator=(const StoSFFDecayer &); private: /** * Store the vertices for scalar intermediate */ vector > sca_; /** * Store the vertices for spin-\f$\frac12\f$ fermion intermediate */ vector > fer_; /** * Store the vertices for spin-\f$\frac32\f$ fermion intermediate */ vector > RSfer_; /** * Store the vertices for vector intermediate */ vector > vec_; /** * Store the vertices for tensor intermediate */ vector > ten_; /** * Store the vertices for four point diagrams */ vector four_; /** * Spin density matrix */ mutable RhoDMatrix rho_; /** * Scalar wavefunction */ mutable ScalarWaveFunction swave_; /** * Spinor wavefunctions */ mutable pair,vector > outspin_[3]; }; } #endif /* THEPEG_StoSFFDecayer_H */ diff --git a/Decay/General/VtoFFVDecayer.cc b/Decay/General/VtoFFVDecayer.cc --- a/Decay/General/VtoFFVDecayer.cc +++ b/Decay/General/VtoFFVDecayer.cc @@ -1,371 +1,371 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the VtoFFVDecayer class. // #include "VtoFFVDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/PDT/ThreeBodyAllOnCalculator.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include using namespace Herwig; using namespace ThePEG; using namespace ThePEG::Helicity; IBPtr VtoFFVDecayer::clone() const { return new_ptr(*this); } IBPtr VtoFFVDecayer::fullclone() const { return new_ptr(*this); } void VtoFFVDecayer::persistentOutput(PersistentOStream & os) const { os << sca_ << fer_ << vec_ << ten_; } void VtoFFVDecayer::persistentInput(PersistentIStream & is, int) { is >> sca_ >> fer_ >> vec_ >> ten_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigVtoFFVDecayer("Herwig::VtoFFVDecayer", "Herwig.so"); void VtoFFVDecayer::Init() { static ClassDocumentation documentation ("The VtoFFVDecayer class implements the general three-body " "decay of a vector to a two fermions and a vector."); } WidthCalculatorBasePtr VtoFFVDecayer:: threeBodyMEIntegrator(const DecayMode & ) const { vector intype; vector inmass,inwidth; vector inpow,inweights; constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights); return new_ptr(ThreeBodyAllOnCalculator (inweights,intype,inmass,inwidth,inpow,*this,0, outgoing()[0]->mass(),outgoing()[1]->mass(), outgoing()[2]->mass(),relativeError())); } -void VtoFFVDecayer::doinit() { - GeneralThreeBodyDecayer::doinit(); +void VtoFFVDecayer::setupDiagrams(bool kinCheck) { + GeneralThreeBodyDecayer::setupDiagrams(kinCheck); if(outgoing().empty()) return; unsigned int ndiags = getProcessInfo().size(); sca_.resize(ndiags); fer_.resize(ndiags); vec_.resize(ndiags); ten_.resize(ndiags); for(unsigned int ix = 0;ix < ndiags; ++ix) { TBDiagram current = getProcessInfo()[ix]; tcPDPtr offshell = current.intermediate; if( offshell->CC() ) offshell = offshell->CC(); if(offshell->iSpin() == PDT::Spin0) { AbstractVVSVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFSVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a scalar diagram in VtoFFVDecayer::doinit()" + << "Invalid vertices for a scalar diagram in VtoFFVDecayer::setupDiagrams()" << Exception::runerror; sca_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin1Half) { AbstractFFVVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFVVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a fermion diagram in VtoFFVDecayer::doinit()" + << "Invalid vertices for a fermion diagram in VtoFFVDecayer::setupDiagrams()" << Exception::runerror; fer_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin1) { AbstractVVVVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFVVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a vector diagram in VtoFFVDecayer::doinit()" + << "Invalid vertices for a vector diagram in VtoFFVDecayer::setupDiagrams()" << Exception::runerror; vec_[ix] = make_pair(vert1, vert2); } else if(offshell->iSpin() == PDT::Spin2) { AbstractVVTVertexPtr vert1 = dynamic_ptr_cast (current.vertices.first); AbstractFFTVertexPtr vert2 = dynamic_ptr_cast (current.vertices.second); if(!vert1||!vert2) throw Exception() - << "Invalid vertices for a tensor diagram in VtoFFVDecayer::doinit()" + << "Invalid vertices for a tensor diagram in VtoFFVDecayer::setupDiagrams()" << Exception::runerror; ten_[ix] = make_pair(vert1, vert2); } } } double VtoFFVDecayer::me2(const int ichan, const Particle & inpart, const ParticleVector & decay, MEOption meopt) const { // particle or CC of particle bool cc = (*getProcessInfo().begin()).incoming != inpart.id(); // special handling or first/last call if(meopt==Initialize) { VectorWaveFunction:: calculateWaveFunctions(inVector_,rho_,const_ptr_cast(&inpart), Helicity::incoming,false); // fix rho if no correlations fixRho(rho_); } if(meopt==Terminate) { VectorWaveFunction:: constructSpinInfo(inVector_,const_ptr_cast(&inpart), Helicity::incoming,true,false); for(unsigned int ix=0;ixdataPtr()->iSpin()==PDT::Spin1) { VectorWaveFunction::constructSpinInfo(outVector_,decay[ix], Helicity::outgoing,true,false); } else { SpinorWaveFunction:: constructSpinInfo(outspin_[ix].first,decay[ix],Helicity::outgoing,true); } } } unsigned int ivec(0); bool massless = false; for(unsigned int ix = 0; ix < decay.size();++ix) { if(decay[ix]->dataPtr()->iSpin() == PDT::Spin1) { massless = decay[ix]->id()==ParticleID::g || decay[ix]->id()==ParticleID::gamma; ivec = ix; VectorWaveFunction:: calculateWaveFunctions(outVector_, decay[ix], Helicity::outgoing, massless ); } else { SpinorWaveFunction:: calculateWaveFunctions(outspin_[ix].first,decay[ix],Helicity::outgoing); outspin_[ix].second.resize(2); // Need a ubar and a v spinor if(outspin_[ix].first[0].wave().Type() == SpinorType::u) { for(unsigned int iy = 0; iy < 2; ++iy) { outspin_[ix].second[iy] = outspin_[ix].first[iy].bar(); outspin_[ix].first[iy].conjugate(); } } else { for(unsigned int iy = 0; iy < 2; ++iy) { outspin_[ix].second[iy] = outspin_[ix].first[iy].bar(); outspin_[ix].second[iy].conjugate(); } } } } const vector > cfactors(getColourFactors()); const vector > nfactors(getLargeNcColourFactors()); Energy2 scale(sqr(inpart.mass())); const size_t ncf(numberOfFlows()); vector flows(ncf, Complex(0.)), largeflows(ncf, Complex(0.)); // setup the DecayMatrixElement vector mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1, ivec == 0 ? PDT::Spin1 : PDT::Spin1Half, ivec == 1 ? PDT::Spin1 : PDT::Spin1Half, ivec == 2 ? PDT::Spin1 : PDT::Spin1Half))); vector mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1, ivec == 0 ? PDT::Spin1 : PDT::Spin1Half, ivec == 1 ? PDT::Spin1 : PDT::Spin1Half, ivec == 2 ? PDT::Spin1 : PDT::Spin1Half))); //the channel possiblities static const unsigned int out2[3] = {1,0,0}, out3[3] = {2,2,1}; for(unsigned int vi = 0; vi < 3; ++vi) { for(unsigned int s1 = 0; s1 < 2; ++s1) { for(unsigned int s2 = 0; s2 < 2; ++s2) { for(unsigned int v1 = 0; v1 < 3; ++v1) { if ( massless && v1 == 1 ) continue; flows = vector(ncf, Complex(0.)); largeflows = vector(ncf, Complex(0.)); unsigned int idiag(0); for(vector::const_iterator dit=getProcessInfo().begin(); dit!=getProcessInfo().end();++dit) { // channels if selecting if( ichan >=0 && diagramMap()[ichan] != idiag ) { ++idiag; continue; } tcPDPtr offshell = dit->intermediate; if(cc&&offshell->CC()) offshell=offshell->CC(); Complex diag; unsigned int o2(out2[dit->channelType]), o3(out3[dit->channelType]); double sign = (o3 < o2) ? 1. : -1.; // intermediate scalar if(offshell->iSpin() == PDT::Spin0) { ScalarWaveFunction inters = sca_[idiag].first-> evaluate(scale, widthOption(), offshell, outVector_[v1], inVector_[vi]); unsigned int h1(s1),h2(s2); if(o2 > o3) swap(h1, h2); if(decay[o2]->id() < 0 && decay[o3]->id() > 0) { diag = -sign*sca_[idiag].second-> evaluate(scale,outspin_[o2].first[h1], outspin_[o3].second[h2],inters); } else { diag = sign*sca_[idiag].second-> evaluate(scale, outspin_[o3].first [h2], outspin_[o2].second[h1],inters); } } // intermediate fermion else if(offshell->iSpin() == PDT::Spin1Half) { int iferm = (decay[o2]->dataPtr()->iSpin() == PDT::Spin1Half) ? o2 : o3; unsigned int h1(s1),h2(s2); if(dit->channelType > iferm) swap(h1, h2); sign = iferm < dit->channelType ? 1. : -1.; if(decay[dit->channelType]->id() < 0 && decay[iferm]->id() > 0 ) { SpinorWaveFunction inters = fer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].first[h1], inVector_[vi]); diag = -sign*fer_[idiag].second-> evaluate(scale,inters,outspin_[iferm].second[h2], outVector_[v1]); } else { SpinorBarWaveFunction inters = fer_[idiag].first-> evaluate(scale,widthOption(),offshell, outspin_[dit->channelType].second[h1],inVector_[vi]); diag = sign*fer_[idiag].second-> evaluate(scale,outspin_[iferm].first [h2],inters, outVector_[v1]); } } // intermediate vector else if(offshell->iSpin() == PDT::Spin1) { VectorWaveFunction interv = vec_[idiag].first-> evaluate(scale, widthOption(), offshell, outVector_[v1], inVector_[vi]); unsigned int h1(s1),h2(s2); if(o2 > o3) swap(h1,h2); if(decay[o2]->id() < 0 && decay[o3]->id() > 0) { diag =-sign*vec_[idiag].second-> evaluate(scale, outspin_[o2].first[h1], outspin_[o3].second[h2], interv); } else { diag = sign*vec_[idiag].second-> evaluate(scale, outspin_[o3].first[h2], outspin_[o2].second[h1], interv); } } else if(offshell->iSpin() == PDT::Spin2) { TensorWaveFunction intert = ten_[idiag].first-> evaluate(scale, widthOption(), offshell, inVector_[vi], outVector_[v1]); unsigned int h1(s1),h2(s2); if(out2[dit->channelType]>out3[dit->channelType]) swap(h1,h2); if(decay[out2[dit->channelType]]->id()<0&& decay[out3[dit->channelType]]->id()>0) { diag =-sign*ten_[idiag].second-> evaluate(scale, outspin_[out2[dit->channelType]].first [h1], outspin_[out3[dit->channelType]].second[h2],intert); } else { diag = sign*ten_[idiag].second-> evaluate(scale, outspin_[out3[dit->channelType]].first [h2], outspin_[out2[dit->channelType]].second[h1],intert); } } // unknown else throw Exception() << "Unknown intermediate in VtoFFVDecayer::me2()" << Exception::runerror; // matrix element for the different colour flows if(ichan < 0) { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } else { for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1 != colourFlow()) continue; flows[dit->colourFlow[iy].first - 1] += dit->colourFlow[iy].second * diag; } for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) { if(dit->colourFlow[iy].first - 1!=colourFlow()) continue; largeflows[dit->largeNcColourFlow[iy].first - 1] += dit->largeNcColourFlow[iy].second * diag; } } ++idiag; } //end of diagrams // now add the flows to the me2 with appropriate colour factors for(unsigned int ix = 0; ix < ncf; ++ix) { if (ivec == 0) { (*mes[ix])(vi, v1, s1, s2) = flows[ix]; (*mel[ix])(vi, v1, s1, s2) = largeflows[ix]; } else if(ivec == 1) { (*mes[ix])(vi, s1, v1, s2) = flows[ix]; (*mel[ix])(vi, s1, v1, s2) = largeflows[ix]; } else if(ivec == 2) { (*mes[ix])(vi, s1, s2, v1) = flows[ix]; (*mel[ix])(vi, s1, s2, v1) = largeflows[ix]; } } } } } } double me2(0.); if(ichan < 0) { vector pflows(ncf,0.); for(unsigned int ix = 0; ix < ncf; ++ix) { for(unsigned int iy = 0; iy < ncf; ++ iy) { double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real(); me2 += con; if(ix == iy) { con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real(); pflows[ix] += con; } } } double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.)); ptotal *= UseRandom::rnd(); for(unsigned int ix = 0;ix < pflows.size(); ++ix) { if(ptotal <= pflows[ix]) { colourFlow(ix); ME(mes[ix]); break; } ptotal -= pflows[ix]; } } else { unsigned int iflow = colourFlow(); me2 = nfactors[iflow][iflow]*(mel[iflow]->contract(*mel[iflow],rho_)).real(); } // return the matrix element squared return me2; } diff --git a/Decay/General/VtoFFVDecayer.h b/Decay/General/VtoFFVDecayer.h --- a/Decay/General/VtoFFVDecayer.h +++ b/Decay/General/VtoFFVDecayer.h @@ -1,161 +1,156 @@ // -*- C++ -*- #ifndef THEPEG_VtoFFVDecayer_H #define THEPEG_VtoFFVDecayer_H // // This is the declaration of the VtoFFVDecayer class. // #include "GeneralThreeBodyDecayer.h" #include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h" #include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h" #include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h" #include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h" #include "ThePEG/Helicity/Vertex/AbstractVVTVertex.h" namespace Herwig { using namespace ThePEG; /** * Here is the documentation of the VtoFFVDecayer class. * * @see \ref VtoFFVDecayerInterfaces "The interfaces" * defined for VtoFFVDecayer. */ class VtoFFVDecayer: public GeneralThreeBodyDecayer { public: /** * Return the matrix element squared for a given mode and phase-space channel * @param ichan The channel we are calculating the matrix element for. * @param part The decaying Particle. * @param decay The particles produced in the decay. * @param meopt Option for the matrix element * @return The matrix element squared for the phase-space configuration. */ virtual double me2(const int ichan, const Particle & part, const ParticleVector & decay, MEOption meopt) const; /** * Method to return an object to calculate the 3 (or higher body) partial width * @param dm The DecayMode * @return A pointer to a WidthCalculatorBase object capable of * calculating the width */ virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} 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. + * Set up the diagrams etc */ - virtual void doinit(); - //@} + virtual void setupDiagrams(bool checkKinematics); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ VtoFFVDecayer & operator=(const VtoFFVDecayer &); private: /** * Store the vertices for scalar intrermediate */ vector > sca_; /** * Store the vertices for fermion intrermediate */ vector > fer_; /** * Store the vertices for vector intrermediate */ vector > vec_; /** * Store the vertices for vector intrermediate */ vector > ten_; /** * Spinr density matrix */ mutable RhoDMatrix rho_; /** * Polarization vectors for the decaying particle */ mutable vector inVector_; /** * Scalar wavefunction for the decay products */ mutable ScalarWaveFunction swave_; /** * Polarization vectors for the decay products */ mutable vector outVector_; /** * Spinors for the decay products */ mutable pair,vector > outspin_[3]; }; } #endif /* THEPEG_VtoFFVDecayer_H */