diff --git a/Decay/General/VectorCurrentDecayer.cc b/Decay/General/VectorCurrentDecayer.cc --- a/Decay/General/VectorCurrentDecayer.cc +++ b/Decay/General/VectorCurrentDecayer.cc @@ -1,273 +1,273 @@ // -*- C++ -*- // // VectorCurrentDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the VectorCurrentDecayer class. // #include "VectorCurrentDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include "Herwig/Models/General/BSMModel.h" using namespace Herwig; IBPtr VectorCurrentDecayer::clone() const { return new_ptr(*this); } IBPtr VectorCurrentDecayer::fullclone() const { return new_ptr(*this); } void VectorCurrentDecayer::persistentOutput(PersistentOStream & os) const { os << inpart_ << currentOut_ << current_ << mode_ << wgtloc_ << wgtmax_ << weights_ << cSMmed_; } void VectorCurrentDecayer::persistentInput(PersistentIStream & is, int) { is >> inpart_ >> currentOut_ >> current_ >> mode_ >> wgtloc_ >> wgtmax_ >> weights_ >> cSMmed_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigVectorCurrentDecayer("Herwig::VectorCurrentDecayer", "Herwig.so"); void VectorCurrentDecayer::Init() { static ClassDocumentation documentation ("The VectorCurrentDecayer class is designed for the decays of low mass vector bosons"); } void VectorCurrentDecayer::setDecayInfo(PDPtr in, const vector & outCurrent, WeakCurrentPtr current) { inpart_ = in; currentOut_ = outCurrent; current_ = current; // cast the model Ptr::ptr model = dynamic_ptr_cast::ptr>(generator()->standardModel()); bool foundU(false),foundD(false),foundS(false); // find the vertices we need and extract the couplings for(unsigned int ix = 0; ix < model->numberOfVertices(); ++ix ) { VertexBasePtr vertex = model->vertex(ix); if(vertex->getNpoint()!=3) continue; for(unsigned int iloc = 0;iloc < 3; ++iloc) { vector ext = vertex->search(iloc, in->id()); if(ext.empty()) continue; for(unsigned int ioff=0;ioffsetCoupling(sqr(in->mass()),getParticleData(1),getParticleData(-1),in); cSMmed_[0] = vertex->norm(); } else if(abs(ext[ioff])==2 && abs(ext[ioff+1])==2 && ext[ioff]==-ext[ioff+1]) { foundU = true; vertex->setCoupling(sqr(in->mass()),getParticleData(2),getParticleData(-2),in); cSMmed_[1] = vertex->norm(); } else if(abs(ext[ioff])==3 && abs(ext[ioff+1])==3 && ext[ioff]==-ext[ioff+1]) { foundS = true; vertex->setCoupling(sqr(in->mass()),getParticleData(3),getParticleData(-3),in); cSMmed_[2] = vertex->norm(); } } } } if(!foundD) { throw InitException() << "Cannot find down quark coupling in VectorCurrentDecayer::doinit()"; } if(!foundU) { throw InitException() << "Cannot find up quark coupling in VectorCurrentDecayer::doinit()"; } if(!foundS) { throw InitException() << "Cannot find strange quark coupling in VectorCurrentDecayer::doinit()"; } } int VectorCurrentDecayer::modeNumber(bool & cc, tcPDPtr parent, const tPDVector & children) const { vector id; id.push_back(parent->id()); for(unsigned int ix=0;ixid()); return modeNumber(cc,id); } void VectorCurrentDecayer::doinitrun() { current_->initrun(); DecayIntegrator::doinitrun(); } void VectorCurrentDecayer::doinit() { DecayIntegrator::doinit(); // make sure the current got initialised current_->init(); // find the mode for(unsigned int ix=0;ixnumberOfModes();++ix) { // get the external particles for this mode int iq(0),ia(0); tPDVector ptemp = current_->particles(inpart_->iCharge(),ix,iq,ia); // check this is the right mode if(ptemp.size()!=currentOut_.size()) continue; vector matched(ptemp.size(),false); bool match = true; for(unsigned int iy=0;iycreateMode(inpart_->iCharge(),tcPDPtr(),FlavourInfo(), ix,mode,0,-1,channel,inpart_->mass()); if(done) { // the maximum weight and the channel weights // the weights for the channel if(weights_.empty()) { weights_.resize(mode->channels().size(),1./(mode->channels().size())); } mode_ = ix; // special for the two body modes if(out.size()==2) { weights_.clear(); mode=new_ptr(PhaseSpaceMode(inpart_,out,1.)); } mode->maxWeight(wgtmax_); mode->setWeights(weights_); addMode(mode); } break; } } int VectorCurrentDecayer::modeNumber(bool & cc, vector id) const { // incoming particle long idtemp; tPDPtr p0=getParticleData(id[0]); idtemp = p0->CC() ? -id[0] : id[0]; if( id[0] ==inpart_->id()) cc=false; else if(idtemp==inpart_->id()) cc=true ; else return -1; vector idout; for(vector::iterator it=++id.begin();it!=id.end();++it) { idout.push_back(*it); } unsigned int icurr=current_->decayMode(idout); if(mode_==icurr) return 0; else return -1; } void VectorCurrentDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast(&part), Helicity::incoming,true,false); weakCurrent()->constructSpinInfo(ParticleVector(decay.begin(),decay.end())); } double VectorCurrentDecayer::me2(const int ichan, const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { using namespace ThePEG::Helicity; // polarization vectors for the incoming particle if(meopt==Initialize) { VectorWaveFunction::calculateWaveFunctions(vectors_,rho_, const_ptr_cast(&part), incoming,false); // fix rho if no correlations fixRho(rho_); } // work out the mapping for the hadron vector int nOut = momenta.size(); vector constants(nOut+1); vector iSpin(nOut); vector hadrons(nOut); int itemp(1); int ix(nOut); do { --ix; iSpin[ix] = outgoing[ix]->iSpin(); itemp *= iSpin[ix]; constants[ix] = itemp; hadrons[ix] = outgoing[ix]->id(); } while(ix>0); constants[nOut] = 1; Energy2 scale(sqr(part.mass())); // calculate the hadron current Energy q = part.mass(); // currents for the different flavour components vector hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero), mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate)); vector hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero), mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate)); vector hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::ssbar), mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate)); // compute the matrix element GeneralDecayMEPtr newME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,iSpin))); vector ihel(momenta.size()+1); unsigned int hI0_size = hadronI0.size(); unsigned int hI1_size = hadronI1.size(); unsigned int hss_size = hadronssbar.size(); unsigned int maxsize = max(max(hadronI0.size(),hadronI1.size()),hss_size); for(unsigned int hhel=0;hhel0;--ix) { ihel[ix]=(hhel%constants[ix-1])/constants[ix]; } for(ihel[0]=0;ihel[0]<3;++ihel[0]) { Complex amp = 0.; // work on coefficients for the I1 and I0 bits if(hI0_size != 0 ) - amp += (cSMmed_[0]+cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI0[hhel])); + amp += Complex((cSMmed_[0]+cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI0[hhel]))); if(hI1_size !=0) - amp += (cSMmed_[0]-cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI1[hhel])); + amp += Complex((cSMmed_[0]-cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI1[hhel]))); if(hss_size !=0) - amp += cSMmed_[2] /q*(vectors_[ihel[0]].wave().dot(hadronssbar[hhel])); + amp += Complex(cSMmed_[2] /q*(vectors_[ihel[0]].wave().dot(hadronssbar[hhel]))); (*newME)(ihel) = amp; } } // store the matrix element ME(newME); // return the answer double output = (ME()->contract(rho_)).real(); return output; } Energy VectorCurrentDecayer::partialWidth(tPDPtr part, vector out) { vector id; id.push_back(part->id()); for(unsigned int ix=0;ixid()); bool cc; int mode=modeNumber(cc,id); imode(mode); return initializePhaseSpaceMode(mode,true,true); } diff --git a/Decay/ScalarMeson/PScalar4FermionsDecayer.cc b/Decay/ScalarMeson/PScalar4FermionsDecayer.cc --- a/Decay/ScalarMeson/PScalar4FermionsDecayer.cc +++ b/Decay/ScalarMeson/PScalar4FermionsDecayer.cc @@ -1,398 +1,398 @@ // -*- C++ -*- // // PScalar4FermionsDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the PScalar4FermionsDecayer class. // #include "PScalar4FermionsDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Helicity/ScalarSpinInfo.h" #include "ThePEG/Helicity/FermionSpinInfo.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "ThePEG/Helicity/epsilon.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" using namespace Herwig; using namespace ThePEG::Helicity; void PScalar4FermionsDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()) { for(unsigned int ix=0;ix<_incoming.size();++ix) _maxweight[ix] = mode(ix)->maxWeight(); } } PScalar4FermionsDecayer::PScalar4FermionsDecayer() : _coupling(1,0.025159062/GeV), _incoming(1,111), _outgoing1(1,11), _outgoing2(1,11), _maxweight(1,0.000234211), _includeVMD(1,2),_VMDid(1,113), _VMDmass(1,0.7758*GeV), _VMDwidth(1,0.1503*GeV), _initsize(1) { // intermediates generateIntermediates(false); } void PScalar4FermionsDecayer::doinit() { DecayIntegrator::doinit(); // check the parameters are consistent unsigned int isize=_coupling.size(); if(isize!=_incoming.size() || isize!=_outgoing1.size() || isize!=_outgoing2.size()|| isize!=_maxweight.size() || isize!=_includeVMD.size()|| isize!=_VMDid.size() || isize!=_VMDmass.size() || isize!=_VMDwidth.size()) throw InitException() << "Inconsistent parameters in PScalar4FermionsDecayer" << Exception::abortnow; // create the integration channels for each mode tPDPtr gamma=getParticleData(ParticleID::gamma); for(unsigned int ix=0;ix<_incoming.size();++ix) { tPDPtr in = getParticleData(_incoming[ix]); tPDVector out={getParticleData( _outgoing1[ix]), getParticleData(-_outgoing1[ix]), getParticleData( _outgoing2[ix]), getParticleData(-_outgoing2[ix])}; PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,_maxweight[ix])); PhaseSpaceChannel newChannel((PhaseSpaceChannel(mode),0,gamma,0,gamma,1,1,1,2,2,3,2,4)); newChannel.setJacobian(1,PhaseSpaceChannel::PhaseSpaceResonance::Power,-1.1); newChannel.setJacobian(2,PhaseSpaceChannel::PhaseSpaceResonance::Power,-1.1); mode->addChannel(newChannel); PhaseSpaceChannel newChannel2((PhaseSpaceChannel(mode),0,gamma,0,gamma,1,3,1,2,2,1,2,4)); newChannel2.setJacobian(1,PhaseSpaceChannel::PhaseSpaceResonance::Power,-1.1); newChannel2.setJacobian(2,PhaseSpaceChannel::PhaseSpaceResonance::Power,-1.1); mode->addChannel(newChannel2); addMode(mode); } // set up the values for the VMD factor if needed (copy the default mass and width // into the array) for(unsigned ix=0;ixmass(); _VMDwidth[ix]=getParticleData(_VMDid[ix])->width(); } } } int PScalar4FermionsDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { // must be four outgoing particles if(children.size()!=4) return -1; // get the id's of the outgoing particles int id[4]={0,0,0,0}; bool done[4]={false,false,false,false}; unsigned int ix(0),iy(0); // ids of the particles int id0(parent->id()),idtemp(-1),idl1(-1),idl2(-1),idt[2]; tPDVector::const_iterator pit = children.begin(); for ( ;pit!=children.end();++pit) { id[ix]=(**pit).id(); done[ix]=false; ++ix; } // find the two lepton pairs // find the first fermion ix=0; do { if( id[ix]>0 && !done[ix] ) { done[ix]=true; idtemp=id[ix]; } ++ix; } while(ix<4&&idtemp<0); if(idtemp<0) return -1; // find its antiparticle ix=0; do { if( id[ix]==-idtemp && !done[ix] ) { done[ix]=true; idl1=idtemp; } ++ix; } while( ix<4 && idl1<0 ); if(idl1<0) return -1; // find the second particle antiparticle pair for(ix=0;ix<4;++ix) { if(!done[ix]) { idt[iy]=id[ix]; ++iy; } } if(idt[0]==-idt[1]) idl2=abs(idt[0]); if(idl2<0) return -1; // loop over the modes and see if this is one of them ix=0; int imode(-1); do { if(_incoming[ix]==id0) { if((idl1==_outgoing1[ix]&&idl2==_outgoing2[ix])|| (idl2==_outgoing1[ix]&&idl1==_outgoing2[ix])) imode=ix; } ++ix; } while(imode<0&&ix<_incoming.size()); cc=false; return imode; } void PScalar4FermionsDecayer::persistentOutput(PersistentOStream & os) const { os << ounit(_coupling,1/MeV) << _incoming << _outgoing1 << _outgoing2 << _maxweight << _includeVMD << _VMDid << ounit(_VMDmass,MeV) << ounit(_VMDwidth,MeV); } void PScalar4FermionsDecayer::persistentInput(PersistentIStream & is, int) { is >> iunit(_coupling,1/MeV) >> _incoming >> _outgoing1 >> _outgoing2 >> _maxweight >> _includeVMD >> _VMDid >> iunit(_VMDmass,MeV) >> iunit(_VMDwidth,MeV); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigPScalar4FermionsDecayer("Herwig::PScalar4FermionsDecayer", "HwSMDecay.so"); void PScalar4FermionsDecayer::Init() { static ClassDocumentation documentation ("The PScalar4FermionsDecayer class is designed for the decay" " of a pseudosclar meson to four fermions. It is intended for the decay of" "the pion to two electron-positron pairs."); static ParVector interfaceIncoming ("Incoming", "The PDG code for the incoming particle", &PScalar4FermionsDecayer::_incoming, 0, 0, 0, -10000, 10000, false, false, true); static ParVector interfaceOutcoming1 ("Outgoing1", "The PDG code for the first outgoing fermion", &PScalar4FermionsDecayer::_outgoing1, 0, 0, 0, -10000, 10000, false, false, true); static ParVector interfaceOutcoming2 ("Outgoing2", "The PDG code for the second outgoing fermion", &PScalar4FermionsDecayer::_outgoing2, 0, 0, 0, -10000, 10000, false, false, true); static ParVector interfaceCoupling ("Coupling", "The coupling for the decay mode", &PScalar4FermionsDecayer::_coupling, 1/MeV, 0, ZERO, -10000/MeV, 10000/MeV, false, false, true); static ParVector interfaceMaxWeight ("MaxWeight", "The maximum weight for the decay mode", &PScalar4FermionsDecayer::_maxweight, 0, 0, 0, -10000, 10000, false, false, true); static ParVector interfaceIncludeVMD ("IncludeVMD", "There are three options for 0 the VMD factor is not included, for 1 the factor " "is included using the default mass and width of the particle specified by" " VMDID, and for 2 the factor is included using the mass and width specified" " by VMDwidth and VMDmass.", &PScalar4FermionsDecayer::_includeVMD, 0, 0, 0, 0, 2, false, false, true); static ParVector interfaceVMDID ("VMDID", "The PDG code for the particle to be used for the VMD factor.", &PScalar4FermionsDecayer::_VMDid, 0, 0, 0, -10000, 10000, false, false, true); static ParVector interfaceVMDmass ("VMDmass", "The mass to use for the particle in the VMD factor", &PScalar4FermionsDecayer::_VMDmass, 1.*MeV, -1, ZERO, -10000*MeV, 10000*MeV, false, false, true); static ParVector interfaceVMDwidth ("VMDwidth", "The width to use for the particle in the VMD factor", &PScalar4FermionsDecayer::_VMDwidth, 1.*MeV, -1, ZERO, -10000*MeV, 10000*MeV, false, false, true); } void PScalar4FermionsDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // set up the spin information for the decay products ScalarWaveFunction::constructSpinInfo(const_ptr_cast(&part), incoming,true); // set up the spin information for the decay products for(unsigned int ix=0;ix<2;++ix) { SpinorBarWaveFunction:: constructSpinInfo(_wavebar[ix],decay[2*ix ],outgoing,true); SpinorWaveFunction:: constructSpinInfo(_wave[ix] ,decay[2*ix+1],outgoing,true); } } double PScalar4FermionsDecayer::me2(const int,const Particle & part, const tPDVector &, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half,PDT::Spin1Half))); bool identical((_outgoing1[imode()]==_outgoing2[imode()])); if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(_rho,const_ptr_cast(&part),incoming); } // calculate the spinors for(unsigned int ix=0;ix<2;++ix) { _wavebar[ix].resize(2); _wave [ix].resize(2); for(unsigned int ihel=0;ihel<2;++ihel) { _wavebar[ix][ihel] = HelicityFunctions::dimensionedSpinorBar(-momenta[2*ix ],ihel,Helicity::outgoing); _wave [ix][ihel] = HelicityFunctions::dimensionedSpinor (-momenta[2*ix+1],ihel,Helicity::outgoing); } } // // momenta of the outgoing photons Lorentz5Momentum poff[4]; poff[0]=momenta[0]+momenta[1]; poff[0].rescaleMass(); poff[1]=momenta[2]+momenta[3]; poff[1].rescaleMass(); if(identical) { poff[2]=momenta[2]+momenta[1]; poff[2].rescaleMass(); poff[3]=momenta[0]+momenta[3]; poff[3].rescaleMass(); } // compute the currents for the two leptonic decays LorentzPolarizationVectorE current[4][2][2]; unsigned int it,ix,iy,iz; for(iz=0;iz<2;++iz) { it = iz==0 ? 1 : 0; for(ix=0;ix<2;++ix) { for(iy=0;iy<2;++iy) { current[iz][iy][ix] = _wave[iz][ix].vectorCurrent(_wavebar[iz][iy]); // the second two currents if(identical) current[iz+2][iy][ix] = _wave[it][ix].vectorCurrent(_wavebar[iz][iy]); } } } // invariants Energy2 m12(sqr(poff[0].mass())); Energy2 m34(sqr(poff[1].mass())); Energy2 m14(ZERO), m23(ZERO); complex prop1(1./m12/m34),prop2(0./sqr(MeV2)); Complex ii(0.,1.); if(identical) { m14=poff[2].mass()*poff[2].mass(); m23=poff[3].mass()*poff[3].mass(); prop2=1./m14/m23; } // the VMD factor if needed if(_includeVMD[imode()]>0) { Energy2 mrho2(_VMDmass[imode()]*_VMDmass[imode()]); Energy2 mwrho(_VMDmass[imode()]*_VMDwidth[imode()]); prop1 = prop1*(-mrho2+ii*mwrho)/(m12-mrho2+ii*mwrho)* (-mrho2+ii*mwrho)/(m34-mrho2+ii*mwrho); if(identical) { prop2 = prop2*(-mrho2+ii*mwrho)/(m14-mrho2+ii*mwrho)* (-mrho2+ii*mwrho)/(m23-mrho2+ii*mwrho); } } // prefactor Complex pre(_coupling[imode()]*4.*Constants::pi *SM().alphaEM()*part.mass()); Complex diag; // now compute the matrix element LorentzVector > eps; vector ispin(5,0); for(ispin[1]=0; ispin[1]<2;++ispin[1]) { for(ispin[2]=0;ispin[2]<2;++ispin[2]) { for(ispin[3]=0;ispin[3]<2;++ispin[3]) { for(ispin[4]=0;ispin[4]<2;++ispin[4]) { // the first diagram eps = epsilon(current[0][ispin[1]][ispin[2]],poff[1], current[1][ispin[3]][ispin[4]]); - diag = prop1*(eps*poff[0]); + diag = Complex(prop1*(eps*poff[0])); // exchanged diagram if identical particles // (sign due normal ordering) if(identical) { eps = epsilon(current[2][ispin[1]][ispin[4]],poff[3], current[3][ispin[3]][ispin[2]]); - diag-= prop2*(eps*poff[2]); + diag-= Complex(prop2*(eps*poff[2])); } (*ME())(ispin)=pre*diag; } } } } double me=ME()->contract(_rho).real(); if(identical) me *= 0.25; return me; } // output the setup info for the particle database void PScalar4FermionsDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class DecayIntegrator::dataBaseOutput(output,false); for(unsigned int ix=0;ix<_incoming.size();++ix) { if(ix<_initsize) { output << "newdef " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "newdef " << name() << ":Outgoing1 " << ix << " " << _outgoing1[ix] << "\n"; output << "newdef " << name() << ":Outgoing2 " << ix << " " << _outgoing2[ix] << "\n"; output << "newdef " << name() << ":Coupling " << ix << " " << _coupling[ix]*MeV << "\n"; output << "newdef " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; output << "newdef " << name() << ":IncludeVMD " << ix << " " << _includeVMD[ix] << "\n"; output << "newdef " << name() << ":VMDID " << ix << " " << _VMDid[ix] << "\n"; output << "newdef " << name() << ":VMDmass " << ix << " " << _VMDmass[ix]/MeV << "\n"; output << "newdef " << name() << ":VMDwidth " << ix << " " << _VMDwidth[ix]/MeV << "\n"; } else { output << "insert " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "insert " << name() << ":Outgoing1 " << ix << " " << _outgoing1[ix] << "\n"; output << "insert " << name() << ":Outgoing2 " << ix << " " << _outgoing2[ix] << "\n"; output << "insert " << name() << ":Coupling " << ix << " " << _coupling[ix]*MeV << "\n"; output << "insert " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; output << "insert " << name() << ":IncludeVMD " << ix << " " << _includeVMD[ix] << "\n"; output << "insert " << name() << ":VMDID " << ix << " " << _VMDid[ix] << "\n"; output << "insert " << name() << ":VMDmass " << ix << " " << _VMDmass[ix]/MeV << "\n"; output << "insert " << name() << ":VMDwidth " << ix << " " << _VMDwidth[ix]/MeV << "\n"; } } if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Decay/ScalarMeson/PScalarPScalarVectorDecayer.cc b/Decay/ScalarMeson/PScalarPScalarVectorDecayer.cc --- a/Decay/ScalarMeson/PScalarPScalarVectorDecayer.cc +++ b/Decay/ScalarMeson/PScalarPScalarVectorDecayer.cc @@ -1,283 +1,283 @@ // -*- C++ -*- // // PScalarPScalarVectorDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the PScalarPScalarVectorDecayer class. // #include "PScalarPScalarVectorDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" #include "ThePEG/Helicity/HelicityFunctions.h" using namespace Herwig; using namespace ThePEG::Helicity; void PScalarPScalarVectorDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()) { for(unsigned int ix=0;ix<_incoming.size();++ix) if(mode(ix)) _maxweight[ix] = mode(ix)->maxWeight(); } } PScalarPScalarVectorDecayer::PScalarPScalarVectorDecayer() { _incoming = { 100111, 100211, 100211, 100311, 100321, 100311, 100321, 100311, 100321, 100311, 100321, 100331, 100331,9020221,9020221, 10221, 10221,9030221,9030221}; _outgoingP = { -211, 111, 211, 311, 321, 321, 311, 111, 111, -211, 211, -321, -311, -321, -311, -211, 111, -211, 111}; _outgoingV = { 213, 213, 113, 113, 113, -213, 213, 313, 323, 323, 313, 323, 313, 323, 313, 20213, 20113, 20213, 20113}; _coupling = { 3.57, 3.57, 3.57, 1., 1., 1.41, 1.41, 1.55, 1.55, 2.19, 2.19, 2.92, 2.92, 0.956, 0.956, 2.68, 2.68, 1.147, 1.147}; _maxweight = { 4.5, 4.5, 4.5, 4., 4., 4., 4., 2., 2., 2., 2., 3.5, 3.5, 4., 4., 4.5, 4.5, 3.2, 3.2}; // initial size of the arrays _initsize=_incoming.size(); // intermediates generateIntermediates(false); } void PScalarPScalarVectorDecayer::doinit() { DecayIntegrator::doinit(); // check the parameters arew consistent unsigned int isize=_coupling.size(); if(isize!=_incoming.size() || isize!=_outgoingP.size()|| isize!=_outgoingV.size() || isize!=_maxweight.size()) throw InitException() << "Inconsistent parameters in PScalarPScalarVectorDecayer" << Exception::abortnow; // set up the integration channels PhaseSpaceModePtr mode; for(unsigned int ix=0;ix<_incoming.size();++ix) { tPDPtr in = getParticleData( _incoming[ix]); tPDVector out = {getParticleData(_outgoingP[ix]), getParticleData(_outgoingV[ix])}; if(in&&out[0]&&out[1]) mode = new_ptr(PhaseSpaceMode(in,out,_maxweight[ix])); else mode=PhaseSpaceModePtr(); addMode(mode); } } int PScalarPScalarVectorDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { if(children.size()!=2) return -1; int id(parent->id()); int idbar = parent->CC() ? parent->CC()->id() : id; int id1(children[0]->id()); int id1bar = children[0]->CC() ? children[0]->CC()->id() : id1; int id2(children[1]->id()); int id2bar = children[1]->CC() ? children[1]->CC()->id() : id2; int imode(-1); unsigned int ix(0); cc=false; do { if(id ==_incoming[ix]) { if((id1 ==_outgoingP[ix]&&id2 ==_outgoingV[ix])|| (id2 ==_outgoingP[ix]&&id1 ==_outgoingV[ix])) imode=ix; } if(idbar==_incoming[ix]) { if((id1bar==_outgoingP[ix]&&id2bar==_outgoingV[ix])|| (id2bar==_outgoingP[ix]&&id1bar==_outgoingV[ix])) { imode=ix; cc=true; } } ++ix; } while(ix<_incoming.size()&&imode<0); return imode; } void PScalarPScalarVectorDecayer::persistentOutput(PersistentOStream & os) const { os << _coupling << _incoming << _outgoingP << _outgoingV << _maxweight; } void PScalarPScalarVectorDecayer::persistentInput(PersistentIStream & is, int) { is >> _coupling >> _incoming >> _outgoingP >> _outgoingV >> _maxweight; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigPScalarPScalarVectorDecayer("Herwig::PScalarPScalarVectorDecayer", "HwSMDecay.so"); void PScalarPScalarVectorDecayer::Init() { static ClassDocumentation documentation ("The PScalarPScalarVectorDecayer class is designed for" " the decay of a pseduoscalar meson to two spin-1 particles."); static ParVector interfaceIncoming ("Incoming", "The PDG code for the incoming particle", &PScalarPScalarVectorDecayer::_incoming, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutgoingScalar ("OutgoingPScalar", "The PDG code for the outgoing pseudoscalar meson", &PScalarPScalarVectorDecayer::_outgoingP, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutgoingVector ("OutgoingVector", "The PDG code for the outgoing vector meson", &PScalarPScalarVectorDecayer::_outgoingV, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceCoupling ("Coupling", "The coupling for the decay mode", &PScalarPScalarVectorDecayer::_coupling, 0, 0, 0, 0., 100., false, false, true); static ParVector interfaceMaxWeight ("MaxWeight", "The maximum weight for the decay mode", &PScalarPScalarVectorDecayer::_maxweight, 0, 0, 0, 0., 100., false, false, true); } void PScalarPScalarVectorDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // set up the spin information for the decay products ScalarWaveFunction::constructSpinInfo(const_ptr_cast(&part), incoming,true); ScalarWaveFunction::constructSpinInfo(decay[0],outgoing,true); VectorWaveFunction::constructSpinInfo(_vectors,decay[1], outgoing,true,false); } double PScalarPScalarVectorDecayer::me2(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin1))); if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(_rho,const_ptr_cast(&part),incoming); } _vectors.resize(3); bool massless = outgoing[1]->id()==ParticleID::gamma; for(unsigned int ix=0;ix<3;++ix) { if(massless && ix==1) continue; _vectors[ix] = HelicityFunctions::polarizationVector(-momenta[1],ix,Helicity::outgoing); } // calculate the matrix element Lorentz5Momentum psum(part.momentum()+momenta[0]); for(unsigned int ix=0;ix<3;++ix) { - (*ME())(0,0,ix)=_coupling[imode()]/part.mass()*(_vectors[ix]*psum); + (*ME())(0,0,ix) = Complex(_coupling[imode()]/part.mass()*(_vectors[ix]*psum)); } double me=ME()->contract(_rho).real(); // test of the matrix element // Energy pcm = Kinematics::pstarTwoBodyDecay(part.mass(),momenta[0].mass(), // momenta[1].mass()); // double test = 4.*sqr(_coupling[imode()]*pcm/momenta[1].mass()); // cerr << "testing matrix element for " << part.PDGName() << " -> " // << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " " // << me << " " << (me-test)/(me+test) << "\n"; // output the answer return me; } // specify the 1-2 matrix element to be used in the running width calculation bool PScalarPScalarVectorDecayer::twoBodyMEcode(const DecayMode & dm, int & mecode, double & coupling) const { int imode(-1); int id(dm.parent()->id()); int idbar = dm.parent()->CC() ? dm.parent()->CC()->id() : id; ParticleMSet::const_iterator pit(dm.products().begin()); int id1((**pit).id()); int id1bar = (**pit).CC() ? (**pit).CC()->id() : id1; ++pit; int id2((**pit).id()); int id2bar = (**pit).CC() ? (**pit).CC()->id() : id2; unsigned int ix(0); bool order(true); do { if(id ==_incoming[ix]) { if(id1==_outgoingP[ix]&&id2==_outgoingV[ix]) { imode=ix; order=true; } if(id2==_outgoingP[ix]&&id1==_outgoingV[ix]) { imode=ix; order=false; } } if(idbar==_incoming[ix]&&imode<0) { if(id1bar==_outgoingP[ix]&&id2bar==_outgoingV[ix]) { imode=ix; order=true; } if(id2bar==_outgoingP[ix]&&id1bar==_outgoingV[ix]) { imode=ix; order=false; } } ++ix; } while(ix<_incoming.size()&&imode<0); coupling=_coupling[imode]; mecode=10; return order; } // output the setup information for the particle database void PScalarPScalarVectorDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class DecayIntegrator::dataBaseOutput(output,false); // the rest of the parameters for(unsigned int ix=0;ix<_incoming.size();++ix) { if(ix<_initsize) { output << "newdef " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "newdef " << name() << ":OutgoingPScalar " << ix << " " << _outgoingP[ix] << "\n"; output << "newdef " << name() << ":OutgoingVector " << ix << " " << _outgoingV[ix] << "\n"; output << "newdef " << name() << ":Coupling " << ix << " " << _coupling[ix] << "\n"; output << "newdef " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; } else { output << "insert " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "insert " << name() << ":OutgoingPScalar " << ix << " " << _outgoingP[ix] << "\n"; output << "insert " << name() << ":OutgoingVector " << ix << " " << _outgoingV[ix] << "\n"; output << "insert " << name() << ":Coupling " << ix << " " << _coupling[ix] << "\n"; output << "insert " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; } } if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Decay/ScalarMeson/PScalarVectorVectorDecayer.cc b/Decay/ScalarMeson/PScalarVectorVectorDecayer.cc --- a/Decay/ScalarMeson/PScalarVectorVectorDecayer.cc +++ b/Decay/ScalarMeson/PScalarVectorVectorDecayer.cc @@ -1,275 +1,275 @@ // -*- C++ -*- // // PScalarVectorVectorDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the PScalarVectorVectorDecayer class. // #include "PScalarVectorVectorDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "ThePEG/Helicity/epsilon.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" using namespace Herwig; using namespace ThePEG::Helicity; void PScalarVectorVectorDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()) { for(unsigned int ix=0;ix<_incoming.size();++ix) _maxweight[ix] = mode(ix)->maxWeight(); } } PScalarVectorVectorDecayer::PScalarVectorVectorDecayer() : _incoming(10), _outgoing1(10), _outgoing2(10), _coupling(10), _maxweight(10) { // decay eta -> omega gamma _incoming[0] = 331; _outgoing1[0] = 223; _outgoing2[0] = 22; _coupling[0] = 0.1412/GeV; _maxweight[0] = 1.2; // decay pi -> gamma gamma _incoming[1] = 111; _outgoing1[1] = 22; _outgoing2[1] = 22; _coupling[1] = 0.0178/GeV; _maxweight[1] = 1.1; // decay eta -> gamma gamma _incoming[2] = 221; _outgoing1[2] = 22; _outgoing2[2] = 22; _coupling[2] = 0.0176/GeV; _maxweight[2] = 1.1; // decay eta' -> gamma gamma _incoming[3] = 331; _outgoing1[3] = 22; _outgoing2[3] = 22; _coupling[3] = 0.0221/GeV; _maxweight[3] = 1.1; // decay eta_c -> rho rho _incoming[4] = 441; _outgoing1[4] = 213; _outgoing2[4] = -213; _coupling[4] = 0.0525/GeV; _maxweight[4] = 2.7; _incoming[5] = 441; _outgoing1[5] = 113; _outgoing2[5] = 113; _coupling[5] = 0.0371/GeV; _maxweight[5] = 2.7; // decay eta-c -> phi phi _incoming[6] = 441; _outgoing1[6] = 333; _outgoing2[6] = 333; _coupling[6] = 0.0267/GeV; _maxweight[6] = 9.; // decay eta-c -> gamma gamma _incoming[7] = 441; _outgoing1[7] = 22; _outgoing2[7] = 22; _coupling[7] = 0.00521/GeV; _maxweight[7] = 1.2; // decay eta_c -> K* K* _incoming[8] = 441; _outgoing1[8] = 323; _outgoing2[8] = -323; _coupling[8] = 0.0308/GeV; _maxweight[8] = 5.3; _incoming[9] = 441; _outgoing1[9] = 313; _outgoing2[9] = -313; _coupling[9] = 0.0308/GeV; _maxweight[9] = 5.3; // initial size of the vectors _initsize = _incoming.size(); // intermediates generateIntermediates(false); } void PScalarVectorVectorDecayer::doinit() { DecayIntegrator::doinit(); // check the parameters arew consistent unsigned int isize(_coupling.size()); if(isize!=_incoming.size() || isize!=_outgoing1.size()|| isize!=_outgoing2.size() || isize!=_maxweight.size()) throw InitException() << "Inconsistent parameters in PScalarVectorVectorDecayer" << Exception::abortnow; // set up the integration channels vector wgt; PhaseSpaceModePtr mode; for(unsigned int ix=0;ix<_incoming.size();++ix) { tPDPtr in = getParticleData( _incoming[ix]); tPDVector out = {getParticleData(_outgoing1[ix]), getParticleData(_outgoing2[ix])}; if(in&&out[0]&&out[1]) mode = new_ptr(PhaseSpaceMode(in,out,_maxweight[ix])); else mode=PhaseSpaceModePtr(); addMode(mode); } } int PScalarVectorVectorDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { cc = false; if(children.size()!=2) return -1; int id(parent->id()); int id1(children[0]->id()); int id2(children[1]->id()); unsigned int ix(0); int imode(-1); do { if(_incoming[ix]==id) { if((id1==_outgoing1[ix]&&id2==_outgoing2[ix])|| (id2==_outgoing1[ix]&&id1==_outgoing2[ix])) imode=ix; } ++ix; } while(imode<0&&ix<_incoming.size()); return imode; } void PScalarVectorVectorDecayer::persistentOutput(PersistentOStream & os) const { os << ounit(_coupling,1/GeV) << _incoming << _outgoing1 << _outgoing2 << _maxweight; } void PScalarVectorVectorDecayer::persistentInput(PersistentIStream & is, int) { is >> iunit(_coupling,1/GeV) >> _incoming >> _outgoing1 >> _outgoing2 >> _maxweight; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigPScalarVectorVectorDecayer("Herwig::PScalarVectorVectorDecayer", "HwSMDecay.so"); void PScalarVectorVectorDecayer::Init() { static ClassDocumentation documentation ("The PScalarVectorVectorDecayer class is designed for" " the decay of a pseduoscalar meson to two spin-1 particles."); static ParVector interfaceIncoming ("Incoming", "The PDG code for the incoming particle", &PScalarVectorVectorDecayer::_incoming, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutcoming1 ("FirstOutgoing", "The PDG code for the first outgoing particle", &PScalarVectorVectorDecayer::_outgoing1, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutcoming2 ("SecondOutgoing", "The PDG code for the second outgoing particle", &PScalarVectorVectorDecayer::_outgoing2, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceCoupling ("Coupling", "The coupling for the decay mode", &PScalarVectorVectorDecayer::_coupling, 1/GeV, 0, ZERO, ZERO, 10000/GeV, false, false, true); static ParVector interfaceMaxWeight ("MaxWeight", "The maximum weight for the decay mode", &PScalarVectorVectorDecayer::_maxweight, 0, 0, 0, 0., 200., false, false, true); } void PScalarVectorVectorDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // set up the spin information for the decay products ScalarWaveFunction::constructSpinInfo(const_ptr_cast(&part), incoming,true); bool photon[2]={false,false}; for(unsigned int ix=0;ix<2;++ix) photon[ix] = decay[ix]->id()==ParticleID::gamma; for(unsigned int ix=0;ix<2;++ix) VectorWaveFunction::constructSpinInfo(_vectors[ix],decay[ix], outgoing,true,photon[ix]); } double PScalarVectorVectorDecayer::me2(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin0,PDT::Spin1,PDT::Spin1))); bool photon[2]={false,false}; for(unsigned int ix=0;ix<2;++ix) photon[ix] = outgoing[ix]->id()==ParticleID::gamma; if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(_rho,const_ptr_cast(&part),incoming); } for(unsigned int ix=0;ix<2;++ix) { _vectors[ix].resize(3); for(unsigned int ihel=0;ihel<3;++ihel) { if(photon[ix] && ihel==1) continue; _vectors[ix][ihel] = HelicityFunctions::polarizationVector(-momenta[ix],ihel,Helicity::outgoing); } } // now compute the matrix element InvEnergy2 fact(_coupling[imode()]/part.mass()); for(unsigned int ix=0;ix<3;++ix) { if(photon[0] && ix==1) continue; for(unsigned int iy=0;iy<3;++iy) { if(photon[1] && iy==1) continue; - (*ME())(0,ix,iy)=fact*epsilon(_vectors[0][ix],momenta[1], - _vectors[1][iy])*momenta[0]; + (*ME())(0,ix,iy)=Complex(fact*epsilon(_vectors[0][ix],momenta[1], + _vectors[1][iy])*momenta[0]); } } double output = ME()->contract(_rho).real(); // test of the matrix element // double test = 2.*sqr(fact*part.mass())* // sqr(Kinematics::pstarTwoBodyDecay(part.mass(),momenta[0].mass(),momenta[1].mass())); // cerr << "testing the matrix element for " << part.PDGName() << " -> " // << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " " // << output << " " << (output-test)/(output+test) << "\n"; return output; } // specify the 1-2 matrix element to be used in the running width calculation bool PScalarVectorVectorDecayer::twoBodyMEcode(const DecayMode & dm, int & itype, double & coupling) const { int imode(-1); int id(dm.parent()->id()); ParticleMSet::const_iterator pit(dm.products().begin()); int id1((**pit).id());++pit; int id2((**pit).id()); unsigned int ix(0); do { if(_incoming[ix]==id) { if((id1==_outgoing1[ix]&&id2==_outgoing2[ix])|| (id2==_outgoing1[ix]&&id1==_outgoing2[ix])) imode=ix; } ++ix; } while(imode<0&&ix<_incoming.size()); coupling=_coupling[imode]*dm.parent()->mass(); itype = 3; return id1==_outgoing1[imode]&&id2==_outgoing2[imode]; } // output the setup info for the particle database void PScalarVectorVectorDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class DecayIntegrator::dataBaseOutput(output,false); for(unsigned int ix=0;ix<_incoming.size();++ix) { if(ix<_initsize) { output << "newdef " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "newdef " << name() << ":FirstOutgoing " << ix << " " << _outgoing1[ix] << "\n"; output << "newdef " << name() << ":SecondOutgoing " << ix << " " << _outgoing2[ix] << "\n"; output << "newdef " << name() << ":Coupling " << ix << " " << _coupling[ix]*GeV << "\n"; output << "newdef " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; } else { output << "insert " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "insert " << name() << ":FirstOutgoing " << ix << " " << _outgoing1[ix] << "\n"; output << "insert " << name() << ":SecondOutgoing " << ix << " " << _outgoing2[ix] << "\n"; output << "insert " << name() << ":Coupling " << ix << " " << _coupling[ix]*GeV << "\n"; output << "insert " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; } } if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Decay/ScalarMeson/ScalarMesonTensorScalarDecayer.cc b/Decay/ScalarMeson/ScalarMesonTensorScalarDecayer.cc --- a/Decay/ScalarMeson/ScalarMesonTensorScalarDecayer.cc +++ b/Decay/ScalarMeson/ScalarMesonTensorScalarDecayer.cc @@ -1,281 +1,281 @@ // -*- C++ -*- // // ScalarMesonTensorScalarDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ScalarMesonTensorScalarDecayer class. // #include "ScalarMesonTensorScalarDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Helicity/WaveFunction/TensorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" using namespace Herwig; using namespace ThePEG::Helicity; void ScalarMesonTensorScalarDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()) { for(unsigned int ix=0;ix<_incoming.size();++ix) if(mode(ix)) _maxweight[ix] = mode(ix)->maxWeight(); } } ScalarMesonTensorScalarDecayer::ScalarMesonTensorScalarDecayer() : _incoming(3), _outgoingT(3), _outgoingS(3), _coupling(3), _maxweight(3) { // D+ -> f_2 pi _incoming[0] = 411; _outgoingT[0] = 225; _outgoingS[0] = 211; _coupling[0] = 8.23E-7/GeV; _maxweight[0] = 5; // chi_c0 -> K*_0 K*_2 _incoming[1] = 10441; _outgoingT[1] = 325; _outgoingS[1] = -10321; _coupling[1] = 0.0217/GeV; _maxweight[1] = 5; _incoming[2] = 10441; _outgoingT[2] = 315; _outgoingS[2] = -10311; _coupling[2] = 0.0217/GeV; _maxweight[2] = 5; // initial size of the arrays _initsize=_incoming.size(); // intermediates generateIntermediates(false); } void ScalarMesonTensorScalarDecayer::doinit() { DecayIntegrator::doinit(); // check the parameters arew consistent unsigned int isize=_coupling.size(); if(isize!=_incoming.size() || isize!=_outgoingT.size()|| isize!=_outgoingS.size() || isize!=_maxweight.size()) throw InitException() << "Inconsistent parameters in " << "ScalarMesonTensorScalarDecayer" << Exception::abortnow; // set up the integration channels for(unsigned int ix=0;ix<_incoming.size();++ix) { tPDPtr in = getParticleData(_incoming[ix]); tPDVector out = {getParticleData(_outgoingT[ix]), getParticleData(_outgoingS[ix])}; PhaseSpaceModePtr mode; if(in&&out[0]&&out[1]) mode = new_ptr(PhaseSpaceMode(in,out,_maxweight[ix])); else mode=PhaseSpaceModePtr(); addMode(mode); } } int ScalarMesonTensorScalarDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { if(children.size()!=2) return -1; int id0(parent->id()); int id0bar = parent->CC() ? parent->CC()->id() : id0; int id1(children[0]->id()); int id1bar = children[0]->CC() ? children[0]->CC()->id() : id1; int id2(children[1]->id()); int id2bar = children[1]->CC() ? children[1]->CC()->id() : id2; unsigned int ix(0); int imode(-1); do { if(id0 ==_incoming[ix]) { if((id1 ==_outgoingT[ix]&&id2 ==_outgoingS[ix])|| (id2 ==_outgoingT[ix]&&id1 ==_outgoingS[ix])) { imode=ix; cc=false; } } if(id0bar==_incoming[ix]&&imode<0) { if((id1bar==_outgoingT[ix]&&id2bar==_outgoingS[ix])|| (id2bar==_outgoingT[ix]&&id1bar==_outgoingS[ix])) { imode=ix; cc=true; } } ++ix; } while(ix<_incoming.size()&&imode<0); return imode; } void ScalarMesonTensorScalarDecayer::persistentOutput(PersistentOStream & os) const { os << ounit(_coupling,1/GeV) << _incoming << _outgoingT << _outgoingS << _maxweight; } void ScalarMesonTensorScalarDecayer::persistentInput(PersistentIStream & is, int) { is >> iunit(_coupling,1/GeV) >> _incoming >> _outgoingT >> _outgoingS >> _maxweight; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigScalarMesonTensorScalarDecayer("Herwig::ScalarMesonTensorScalarDecayer", "HwSMDecay.so"); void ScalarMesonTensorScalarDecayer::Init() { static ClassDocumentation documentation ("The ScalarMesonTensorScalarDecayer class is designed for" " the decay of a pseduoscalar meson to two spin-1 particles."); static ParVector interfaceIncoming ("Incoming", "The PDG code for the incoming particle", &ScalarMesonTensorScalarDecayer::_incoming, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutcomingT ("OutgoingTensor", "The PDG code for the outgoing tensor", &ScalarMesonTensorScalarDecayer::_outgoingT, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutcomingS ("OutgoingScalar", "The PDG code for the outgoing scalar", &ScalarMesonTensorScalarDecayer::_outgoingS, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceCoupling ("Coupling", "The coupling for the decay mode", &ScalarMesonTensorScalarDecayer::_coupling, 1/GeV, 0, ZERO, ZERO, 100./GeV, false, false, true); static ParVector interfaceMaxWeight ("MaxWeight", "The maximum weight for the decay mode", &ScalarMesonTensorScalarDecayer::_maxweight, 0, 0, 0, 0., 100., false, false, true); } void ScalarMesonTensorScalarDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // set up the spin information for the decay products ScalarWaveFunction::constructSpinInfo(const_ptr_cast(&part), incoming,true); TensorWaveFunction::constructSpinInfo(_tensors,decay[0], outgoing,true,false); ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true); } double ScalarMesonTensorScalarDecayer::me2(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin0,PDT::Spin2,PDT::Spin0))); if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(_rho,const_ptr_cast(&part),incoming); } TensorWaveFunction twave(momenta[0],outgoing[0],Helicity::outgoing); _tensors.resize(5); for(unsigned int ihel=0;ihel<5;++ihel) { twave.reset(ihel); _tensors[ihel] = twave.wave(); } // calculate the matrix element InvEnergy2 fact(_coupling[imode()]/part.mass()); LorentzPolarizationVectorE vtemp; for(unsigned int ix=0;ix<5;++ix) { vtemp = _tensors[ix]*part.momentum(); - (*ME())(0,ix,0) = fact * momenta[1].dot(vtemp); + (*ME())(0,ix,0) = Complex(fact * momenta[1].dot(vtemp)); } double me = ME()->contract(_rho).real(); // // test of the matrix element // Energy pcm = Kinematics::pstarTwoBodyDecay(part.mass(),momenta[0].mass(), // momenta[1].mass()); // double test = 2.*pow<4,1>(pcm)*sqr(_coupling[imode()]*part.mass())/ // 3./pow<4,1>(momenta[0].mass()); // cerr << "testing matrix element for " << part.PDGName() << " -> " // << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " " // << me << " " << (me-test)/(me+test) << "\n"; // output the answer return me; } // specify the 1-2 matrix element to be used in the running width calculation bool ScalarMesonTensorScalarDecayer::twoBodyMEcode(const DecayMode & dm, int & itype, double & coupling) const { int id(dm.parent()->id()); int idbar = dm.parent()->CC() ? dm.parent()->CC()->id() : id; ParticleMSet::const_iterator pit(dm.products().begin()); int id1((**pit).id()); int id1bar = (**pit).CC() ? (**pit).CC()->id() : id1; ++pit; int id2((**pit).id()); int id2bar = (**pit).CC() ? (**pit).CC()->id() : id2; unsigned int ix(0); bool order(false); int imode(-1); do { if(id ==_incoming[ix]) { if(id1==_outgoingT[ix]&&id2==_outgoingS[ix]) { imode=ix; order=true; } if(id2==_outgoingT[ix]&&id1==_outgoingS[ix]) { imode=ix; order=false; } } if(idbar==_incoming[ix]&&imode<0) { if(id1bar==_outgoingT[ix]&&id2bar==_outgoingS[ix]) { imode=ix; order=true; } if(id2bar==_outgoingT[ix]&&id1bar==_outgoingS[ix]) { imode=ix; order=false; } } ++ix; } while(ix<_incoming.size()&&imode<0); coupling=_coupling[imode]*dm.parent()->mass(); itype = 11; return order; } // output the setup information for the particle database void ScalarMesonTensorScalarDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class DecayIntegrator::dataBaseOutput(output,false); // the rest of the parameters for(unsigned int ix=0;ix<_incoming.size();++ix) { if(ix<_initsize) { output << "newdef " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "newdef " << name() << ":OutgoingTensor " << ix << " " << _outgoingT[ix] << "\n"; output << "newdef " << name() << ":OutgoingScalar " << ix << " " << _outgoingS[ix] << "\n"; output << "newdef " << name() << ":Coupling " << ix << " " << _coupling[ix]*GeV << "\n"; output << "newdef " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; } else { output << "insert " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "insert " << name() << ":OutgoingTensor " << ix << " " << _outgoingT[ix] << "\n"; output << "insert " << name() << ":OutgoingScalar " << ix << " " << _outgoingS[ix] << "\n"; output << "insert " << name() << ":Coupling " << ix << " " << _coupling[ix]*GeV << "\n"; output << "insert " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; } } if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Decay/Tau/TauDecayer.cc b/Decay/Tau/TauDecayer.cc --- a/Decay/Tau/TauDecayer.cc +++ b/Decay/Tau/TauDecayer.cc @@ -1,360 +1,360 @@ // -*- C++ -*- // // TauDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the TauDecayer class. // // Author: Peter Richardson // #include "TauDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Helicity/VectorSpinInfo.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "Herwig/Decay/DecayVertex.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include "ThePEG/Helicity/FermionSpinInfo.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Helicity/HelicityFunctions.h" using namespace Herwig; using namespace ThePEG::Helicity; void TauDecayer::doinit() { DecayIntegrator::doinit(); // make sure the current got initialised current_->init(); // set up the phase-space channels tPDPtr tau = getParticleData(ParticleID::tauminus); tPDPtr nu = getParticleData(ParticleID::nu_tau); Energy mtau(tau->mass()); vector channelwgts; modeMap_.clear(); for(unsigned int ix=0;ixnumberOfModes();++ix) { // get the external particles for this mode tPDVector out = {nu}; int iq(0),ia(0); tPDVector ptemp = current_->particles(-3,ix,iq,ia); out.insert(std::end(out), std::begin(ptemp), std::end(ptemp)); // the maximum weight double maxweight = wgtMax_.size()>numberModes() ? wgtMax_[numberModes()] : 1.; // create the mode PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(tau,out,maxweight)); // create the first piece of the channel PhaseSpaceChannel channel((PhaseSpaceChannel(mode),0,1)); if(!current_->createMode(-3,tcPDPtr(),FlavourInfo(),ix,mode,1,0,channel,mtau)) continue; // the channel weights // the weights for the channel if(wgtLoc_.size()>numberModes()&& wgtLoc_[numberModes()]+mode->channels().size()<=weights_.size()) { channelwgts=vector(weights_.begin()+wgtLoc_[numberModes()], weights_.begin()+wgtLoc_[numberModes()]+mode->channels().size()); } else { channelwgts.resize(mode->channels().size(),1./(mode->channels().size())); } modeMap_.push_back(ix); // special for the two body modes if(out.size()==2) { channelwgts.clear(); mode=new_ptr(PhaseSpaceMode(tau,out,maxweight)); } mode->setWeights(channelwgts); // need to do the weights addMode(mode); } current_->reset(); current_->touch(); current_->update(); } void TauDecayer::doinitrun() { current_->initrun(); DecayIntegrator::doinitrun(); if(initialize()) { weights_.clear(); wgtLoc_.clear(); wgtMax_.clear(); for(unsigned int ix=0;ixmaxWeight()); wgtLoc_.push_back(weights_.size()); for(unsigned int iy=0;iychannels().size();++iy) { weights_.push_back(mode(ix)->channels()[iy].weight()); } } } } bool TauDecayer::accept(tcPDPtr parent, const tPDVector & children) const { bool allowed(false); // find the neutrino int idnu(0),idtemp,idin(parent->id()); vector idother; tPDVector::const_iterator pit = children.begin(); tPDVector::const_iterator pend = children.end(); for( ; pit!=pend;++pit) { idtemp=(**pit).id(); if(abs(idtemp)==16) idnu=idtemp; else idother.push_back(idtemp); } if((idnu==ParticleID::nu_tau && idin==ParticleID::tauminus)|| (idnu==ParticleID::nu_taubar && idin==ParticleID::tauplus )) { allowed=current_->accept(idother); } return allowed; } int TauDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { int imode(-1); tPDVector::const_iterator pit = children.begin(); tPDVector::const_iterator pend = children.end(); int idtemp;vector idother; for( ; pit!=pend;++pit) { idtemp=(**pit).id(); if(abs(idtemp)!=16) idother.push_back(idtemp); } unsigned int itemp=current_->decayMode(idother); for(unsigned int ix=0;ixid()==ParticleID::tauplus; return imode; } void TauDecayer::persistentOutput(PersistentOStream & os) const { os << modeMap_ << current_ << wgtLoc_ << wgtMax_ << weights_ << polOpt_ << tauMpol_ << tauPpol_; } void TauDecayer::persistentInput(PersistentIStream & is, int) { is >> modeMap_ >> current_ >> wgtLoc_ >> wgtMax_ >> weights_ >> polOpt_ >> tauMpol_ >> tauPpol_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigTauDecayer("Herwig::TauDecayer", "HwTauDecay.so"); void TauDecayer::Init() { static ClassDocumentation documentation ("The TauDecayer class is designed to use a weak current" " to perform the decay of the tau."); static Reference interfaceWeakCurrent ("WeakCurrent", "The reference for the decay current to be used.", &TauDecayer::current_, false, false, true, false, false); static ParVector interfaceWeightLocation ("WeightLocation", "The locations of the weights for a given channel in the vector", &TauDecayer::wgtLoc_, 0, 0, 0, 0, 10000, false, false, true); static ParVector interfaceWeightMax ("MaximumWeight", "The maximum weight for a given channel.", &TauDecayer::wgtMax_, 0, 0, 0, 0., 100., false, false, true); static ParVector interfaceWeights ("Weights", "The weights for the integration.", &TauDecayer::weights_, 0, 0, 0, 0., 1., false, false, true); static Switch interfacePolarizationOption ("PolarizationOption", "Option of forcing the polarization of the tau leptons, N.B. you" " should only use this option for making distributions for" " comparision if you really know what you are doing.", &TauDecayer::polOpt_, false, false, false); static SwitchOption interfacePolarizationOptionDefault (interfacePolarizationOption, "Default", "Don't force the polarization use the full spin density matrices" " to get the right answer", false); static SwitchOption interfacePolarizationOptionForce (interfacePolarizationOption, "Force", "Force the polarizations", true); static Parameter interfaceTauMinusPolarization ("TauMinusPolarization", "The polarization of the tau-, left=-1, right=+1 if this is forced.", &TauDecayer::tauMpol_, 0.0, -1.0, 1.0, false, false, Interface::limited); static Parameter interfaceTauPlusPolarization ("TauPlusPolarization", "The polarization of the tau+, left=-1, right=+1 if this is forced.", &TauDecayer::tauPpol_, 0.0, -1.0, 1.0, false, false, Interface::limited); } void TauDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { if(part.id()==ParticleID::tauminus) { SpinorWaveFunction :: constructSpinInfo(inSpin_,const_ptr_cast(&part),incoming,true); SpinorBarWaveFunction:: constructSpinInfo(inBar_,decay[0],outgoing,true); } else { SpinorBarWaveFunction:: constructSpinInfo(inBar_ ,const_ptr_cast(&part),incoming,true); SpinorWaveFunction:: constructSpinInfo(inSpin_,decay[0],outgoing,true); } current_->constructSpinInfo(ParticleVector(decay.begin()+1,decay.end())); } // combine the currents to give the matrix element double TauDecayer::me2(const int ichan, const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { // map the mode to those in the current int mode(modeMap_[imode()]); // extract info on the decaying particle if(meopt==Initialize) { // spin density matrix for the decaying particle rho_ = RhoDMatrix(PDT::Spin1Half); if(part.id()==ParticleID::tauminus) SpinorWaveFunction ::calculateWaveFunctions(inSpin_,rho_, const_ptr_cast(&part), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(inBar_ ,rho_, const_ptr_cast(&part), incoming); // fix rho if no correlations fixRho(rho_); if(polOpt_) { rho_(0,1) = rho_(1,0) = 0.; if(part.id()==ParticleID::tauminus) { rho_(0,0) = 0.5*(1.-tauMpol_); rho_(1,1) = 0.5*(1.+tauMpol_); } else { rho_(0,0) = 0.5*(1.+tauPpol_); rho_(1,1) = 0.5*(1.-tauPpol_); } } // work out the mapping for the hadron vector constants_ = vector(outgoing.size()+1); iSpin_ = vector(outgoing.size()); int itemp(1); unsigned int ix(outgoing.size()); do { --ix; iSpin_[ix] = outgoing[ix]->iSpin(); itemp *= iSpin_[ix]; constants_[ix] = itemp; } while(ix>0); constants_[outgoing.size()] = 1; constants_[0 ] = constants_[1]; } if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,iSpin_))); // calculate the spinors for the decay products if(part.id()==ParticleID::tauminus) { inBar_.resize(2); for(unsigned int ihel=0;ihel<2;++ihel) inBar_[ihel] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ihel,Helicity::outgoing); } else { inSpin_.resize(2); for(unsigned int ihel=0;ihel<2;++ihel) inSpin_[ihel] = HelicityFunctions::dimensionedSpinor (-momenta[0],ihel,Helicity::outgoing); } // calculate the hadron current Energy q; vector hadron(current_->current(tcPDPtr(),FlavourInfo(),mode,ichan,q,tPDVector(outgoing.begin()+1,outgoing.end()), vector(momenta.begin()+1,momenta.end()),meopt)); // prefactor double pre = sqr(pow(part.mass()/q,int(outgoing.size()-3))); // calculate the lepton current LorentzPolarizationVectorE lepton[2][2]; for(unsigned ix=0;ix<2;++ix) { for(unsigned iy=0;iy<2;++iy) { if(part.id()==15) lepton[ix][iy]=2.*inSpin_[ix].leftCurrent(inBar_[iy]); else lepton[iy][ix]=2.*inSpin_[ix].leftCurrent(inBar_[iy]); } } // compute the matrix element vector ihel(outgoing.size()+1); for(unsigned int hhel=0;hhel1;--ix) { ihel[ix]=(hhel%constants_[ix-1])/constants_[ix]; } // loop over the helicities of the tau and neutrino and set up the matrix // element for(ihel[1]=0;ihel[1]<2;++ihel[1]){ for(ihel[0]=0;ihel[0]<2;++ihel[0]) { - (*ME())(ihel)= lepton[ihel[0]][ihel[1]].dot(hadron[hhel])* - SM().fermiConstant(); + (*ME())(ihel)= Complex(lepton[ihel[0]][ihel[1]].dot(hadron[hhel])* + SM().fermiConstant()); } } } // multiply by the CKM element int iq,ia; current_->decayModeInfo(mode,iq,ia); double ckm(1.); if(iq<=6) { if(iq%2==0) ckm = SM().CKM(iq/2-1,(abs(ia)-1)/2); else ckm = SM().CKM(abs(ia)/2-1,(iq-1)/2); } return 0.5*pre*ckm*(ME()->contract(rho_)).real(); } // output the setup information for the particle database void TauDecayer::dataBaseOutput(ofstream & output,bool header) const { unsigned int ix; if(header) output << "update decayers set parameters=\""; DecayIntegrator::dataBaseOutput(output,false); for(ix=0;ixdataBaseOutput(output,false,true); output << "newdef " << name() << ":WeakCurrent " << current_->name() << " \n"; output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";\n"; } diff --git a/Decay/TensorMeson/TensorMesonVectorPScalarDecayer.cc b/Decay/TensorMeson/TensorMesonVectorPScalarDecayer.cc --- a/Decay/TensorMeson/TensorMesonVectorPScalarDecayer.cc +++ b/Decay/TensorMeson/TensorMesonVectorPScalarDecayer.cc @@ -1,356 +1,356 @@ // -*- C++ -*- // // TensorMesonVectorPScalarDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the TensorMesonVectorPScalarDecayer class. // #include "TensorMesonVectorPScalarDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Helicity/WaveFunction/TensorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "ThePEG/Helicity/epsilon.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" #include "ThePEG/Helicity/HelicityFunctions.h" using namespace Herwig; using namespace ThePEG::Helicity; void TensorMesonVectorPScalarDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()){ for(unsigned int ix=0;ix<_incoming.size();++ix) if(mode(ix)) _maxweight[ix] = mode(ix)->maxWeight(); } } TensorMesonVectorPScalarDecayer::TensorMesonVectorPScalarDecayer() : _incoming(31), _outgoingV(31), _outgoingP(31), _coupling(31), _maxweight(31) { // a_2 -> rho pi _incoming[0] = 115; _outgoingV[0] = 213; _outgoingP[0] = -211; _coupling[0] = 21.1/GeV2; _maxweight[0] = 10.; _incoming[1] = 215; _outgoingV[1] = 113; _outgoingP[1] = 211; _coupling[1] = 21.1/GeV2; _maxweight[1] = 9.; _incoming[2] = 215; _outgoingV[2] = 213; _outgoingP[2] = 111; _coupling[2] = 21.1/GeV2; _maxweight[2] = 9.; // a_2+/- -> gamma pi+/- _incoming[3] = 215; _outgoingV[3] = 22; _outgoingP[3] = 211; _coupling[3] = 0.551/GeV2; _maxweight[3] = 2.; // k_2 -> K_2 omega _incoming[4] = 315; _outgoingV[4] = 223; _outgoingP[4] = 311; _coupling[4] = 11.66/GeV2; _maxweight[4] = 17.; _incoming[5] = 325; _outgoingV[5] = 223; _outgoingP[5] = 321; _coupling[5] = 11.66/GeV2; _maxweight[5] = 20.5; // k_2+/- -> K+/- gamma _incoming[6] = 325; _outgoingV[6] = 22; _outgoingP[6] = 321; _coupling[6] = 0.553/GeV2; _maxweight[6] = 2.2; // B_c2 -> B_c gamma _incoming[7] = 545; _outgoingV[7] = 22; _outgoingP[7] = 541; _coupling[7] = 0.651/GeV2; _maxweight[7] = 2.; // K_2 -> K rho _incoming[8] = 325; _outgoingV[8] = 113; _outgoingP[8] = 321; _coupling[8] = 10.14/GeV2; _maxweight[8] = 9.; _incoming[9] = 325; _outgoingV[9] = 213; _outgoingP[9] = 311; _coupling[9] = 14.33/GeV2; _maxweight[9] = 9.; _incoming[10] = 315; _outgoingV[10] = 113; _outgoingP[10] = 311; _coupling[10] = 10.14/GeV2; _maxweight[10] = 9.; _incoming[11] = 315; _outgoingV[11] = -213; _outgoingP[11] = 321; _coupling[11] = 14.33/GeV2; _maxweight[11] = 9.; // K_2 -> K* pi _incoming[12] = 325; _outgoingV[12] = 323; _outgoingP[12] = 111; _coupling[12] = 9.733/GeV2; _maxweight[12] = 13; _incoming[13] = 325; _outgoingV[13] = 313; _outgoingP[13] = 211; _coupling[13] = 13.77/GeV2; _maxweight[13] = 11; _incoming[14] = 315; _outgoingV[14] = 313; _outgoingP[14] = 111; _coupling[14] = 9.733/GeV2; _maxweight[14] = 8.; _incoming[15] = 315; _outgoingV[15] = 323; _outgoingP[15] = -211; _coupling[15] = 13.77/GeV2; _maxweight[15] = 8.; // D_2 -> D* pi _incoming[16] = 425; _outgoingV[16] = 423; _outgoingP[16] = 111; _coupling[16] = 8.035/GeV2; _maxweight[16] = 2.2; _incoming[17] = 425; _outgoingV[17] = 413; _outgoingP[17] = -211; _coupling[17] = 11.670/GeV2; _maxweight[17] = 2.4; _incoming[18] = 415; _outgoingV[18] = 413; _outgoingP[18] = 111; _coupling[18] = 6.801/GeV2; _maxweight[18] = 2.4; _incoming[19] = 415; _outgoingV[19] = 423; _outgoingP[19] = 211; _coupling[19] = 9.527/GeV2; _maxweight[19] = 2.; // D_s2 -> D* K _incoming[20] = 435; _outgoingV[20] = 423; _outgoingP[20] = 321; _coupling[20] = 13.10/GeV2; _maxweight[20] = 2.2; _incoming[21] = 435; _outgoingV[21] = 413; _outgoingP[21] = 311; _coupling[21] = 13.10/GeV2; _maxweight[21] = 2.5; // B_2 -> B* pi _incoming[22] = 525; _outgoingV[22] = 523; _outgoingP[22] = 111; _coupling[22] = 4.99/GeV2; _maxweight[22] = 2.1; _incoming[23] = 525; _outgoingV[23] = 513; _outgoingP[23] = 211; _coupling[23] = 7.059/GeV2; _maxweight[23] = 2.1; _incoming[24] = 515; _outgoingV[24] = 513; _outgoingP[24] = 111; _coupling[24] = 4.99/GeV2; _maxweight[24] = 2.1; _incoming[25] = 515; _outgoingV[25] = 523; _outgoingP[25] = -211; _coupling[25] = 7.059/GeV2; _maxweight[25] = 2.1; // D_s2 _incoming[26] = 435; _outgoingV[26] = 423; _outgoingP[26] = 321; _coupling[26] = 13.09/GeV2; _maxweight[26] = 2.2; _incoming[27] = 435; _outgoingV[27] = 413; _outgoingP[27] = 311; _coupling[27] = 13.09/GeV2; _maxweight[27] = 2.5; // B_s2 _incoming[28] = 535; _outgoingV[28] = 523; _outgoingP[28] = -321; _coupling[28] = 7.29/GeV2; _maxweight[28] = 2.4; _incoming[29] = 535; _outgoingV[29] = 513; _outgoingP[29] = -311; _coupling[29] = 9.43/GeV2; _maxweight[29] = 2.1; // upsilon_2(1d) to chi_b gamma _incoming[30] = 20555; _outgoingV[30] = 22; _outgoingP[30] = 10551; _coupling[30] = 1.11/GeV2; _maxweight[30] = 2.4; // initial size of the arrays _initsize=_incoming.size(); // intermediates generateIntermediates(false); } void TensorMesonVectorPScalarDecayer::doinit() { DecayIntegrator::doinit(); // check consistence of the parameters unsigned int isize=_incoming.size(); if(isize!=_outgoingV.size()||isize!=_outgoingP.size()|| isize!=_maxweight.size()||isize!=_coupling.size()) throw InitException() << "Inconsistent parameters TensorMesonVectorPScalarDecayer" << Exception::abortnow; // set up the integration channels PhaseSpaceModePtr mode; for(unsigned int ix=0;ix<_incoming.size();++ix) { tPDPtr in = getParticleData(_incoming[ix]); tPDVector out = {getParticleData(_outgoingV[ix]), getParticleData(_outgoingP[ix])}; if(in&&out[0]&&out[1]) mode=new_ptr(PhaseSpaceMode(in,out,_maxweight[ix])); else mode=PhaseSpaceModePtr(); addMode(mode); } } int TensorMesonVectorPScalarDecayer::modeNumber(bool & cc, tcPDPtr parent, const tPDVector & children) const { if(children.size()!=2) return -1; int id(parent->id()); int idbar = parent->CC() ? parent->CC()->id() : id; int id1(children[0]->id()); int id1bar = children[0]->CC() ? children[0]->CC()->id() : id1; int id2(children[1]->id()); int id2bar = children[1]->CC() ? children[1]->CC()->id() : id2; int imode(-1); unsigned int ix(0); cc=false; do { if(id ==_incoming[ix]) { if((id1 ==_outgoingP[ix]&&id2 ==_outgoingV[ix])|| (id2 ==_outgoingP[ix]&&id1 ==_outgoingV[ix])) imode=ix; } if(idbar==_incoming[ix]) { if((id1bar==_outgoingP[ix]&&id2bar==_outgoingV[ix])|| (id2bar==_outgoingP[ix]&&id1bar==_outgoingV[ix])) { imode=ix; cc=true; } } ++ix; } while(ix<_incoming.size()&&imode<0); return imode; } void TensorMesonVectorPScalarDecayer::persistentOutput(PersistentOStream & os) const { os << _incoming << _outgoingV << _outgoingP << _maxweight << ounit(_coupling,1/GeV2); } void TensorMesonVectorPScalarDecayer::persistentInput(PersistentIStream & is, int) { is >> _incoming >> _outgoingV >> _outgoingP >> _maxweight >> iunit(_coupling,1/GeV2); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigTensorMesonVectorPScalarDecayer("Herwig::TensorMesonVectorPScalarDecayer", "HwTMDecay.so"); void TensorMesonVectorPScalarDecayer::Init() { static ClassDocumentation documentation ("The TensorMesonVectorPScalarDecayer class implements the" " decay of a tensor meson to a spin-1 particle and a pseduoscalar meson"); static ParVector interfaceIncoming ("Incoming", "The PDG code for the incoming particle", &TensorMesonVectorPScalarDecayer::_incoming, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutcomingV ("OutgoingVector", "The PDG code for the outgoing spin-1particle", &TensorMesonVectorPScalarDecayer::_outgoingV, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutcomingP ("OutgoingScalar", "The PDG code for the outgoing pseudoscalar meson", &TensorMesonVectorPScalarDecayer::_outgoingP, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceCoupling ("Coupling", "The coupling for the decay mode", &TensorMesonVectorPScalarDecayer::_coupling, 1/GeV2, 0, ZERO, ZERO, 100./GeV2, false, false, true); static ParVector interfaceMaxWeight ("MaxWeight", "The maximum weight for the decay mode", &TensorMesonVectorPScalarDecayer::_maxweight, 0, 0, 0, 0., 1000., false, false, true); } void TensorMesonVectorPScalarDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { TensorWaveFunction::constructSpinInfo(_tensors,const_ptr_cast(&part), incoming,true,false); // set up the spin information for the decay products VectorWaveFunction::constructSpinInfo(_vectors,decay[0],outgoing,true, decay[0]->id()==ParticleID::gamma); ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true); } // matrix elememt for the process double TensorMesonVectorPScalarDecayer::me2(const int,const Particle & part, const tPDVector &, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin2,PDT::Spin1,PDT::Spin0))); // check for photons bool photon(_outgoingV[imode()]==ParticleID::gamma); // stuff for incoming particle if(meopt==Initialize) { _rho = RhoDMatrix(PDT::Spin2); TensorWaveFunction:: calculateWaveFunctions(_tensors,_rho,const_ptr_cast(&part), incoming,false); } _vectors.resize(3); for(unsigned int ix=0;ix<3;++ix) { if(photon && ix==1) continue; _vectors[ix] = HelicityFunctions::polarizationVector(-momenta[0],ix,Helicity::outgoing); } InvEnergy3 fact(_coupling[imode()]/part.mass()); // calculate the matrix element for(unsigned int inhel=0;inhel<5;++inhel) { for(unsigned int vhel=0;vhel<3;++vhel){ if(vhel==1&&photon) (*ME())(inhel,vhel,0)=0.; else { LorentzVector > vtemp= fact*epsilon(momenta[0],_vectors[vhel],momenta[1]); - (*ME())(inhel,vhel,0)= (momenta[1]*_tensors[inhel]).dot(vtemp); + (*ME())(inhel,vhel,0) = Complex((momenta[1]*_tensors[inhel]).dot(vtemp)); } } } double me = ME()->contract(_rho).real(); // test of the answer // Energy pcm = Kinematics::pstarTwoBodyDecay(part.mass(),momenta[0].mass(), // momenta[1].mass()); // double test = Energy4(pow<4,1>(2*pcm))*sqr( _coupling[imode()])/80.; // cout << "testing matrix element for " << part.PDGName() << " -> " // << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " " // << me << " " << test << " " << (me-test)/(me+test) << endl; // return the answer return me; } bool TensorMesonVectorPScalarDecayer::twoBodyMEcode(const DecayMode & dm,int & mecode, double & coupling) const { int imode(-1); int id(dm.parent()->id()); int idbar = dm.parent()->CC() ? dm.parent()->CC()->id() : id; ParticleMSet::const_iterator pit(dm.products().begin()); int id1((**pit).id()); int id1bar = (**pit).CC() ? (**pit).CC()->id() : id1; ++pit; int id2((**pit).id()); int id2bar = (**pit).CC() ? (**pit).CC()->id() : id2; unsigned int ix(0); bool order(false); do { if(id ==_incoming[ix]) { if(id1==_outgoingP[ix]&&id2==_outgoingV[ix]) { imode=ix; order=true; } if(id2==_outgoingP[ix]&&id1==_outgoingV[ix]) { imode=ix; order=false; } } if(idbar==_incoming[ix]&&imode<0) { if(id1bar==_outgoingP[ix]&&id2bar==_outgoingV[ix]) { imode=ix; order=true; } if(id2bar==_outgoingP[ix]&&id1bar==_outgoingV[ix]) { imode=ix; order=false; } } ++ix; } while(ix<_incoming.size()&&imode<0); coupling=_coupling[imode]*sqr(dm.parent()->mass()); mecode=8; return order; } void TensorMesonVectorPScalarDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class DecayIntegrator::dataBaseOutput(output,false); // the rest of the parameters for(unsigned int ix=0;ix<_incoming.size();++ix) { if(ix<_initsize) { output << "newdef " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "newdef " << name() << ":OutgoingVector " << ix << " " << _outgoingV[ix] << "\n"; output << "newdef " << name() << ":OutgoingScalar " << ix << " " << _outgoingP[ix] << "\n"; output << "newdef " << name() << ":Coupling " << ix << " " << _coupling[ix]*GeV2 << "\n"; output << "newdef " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; } else { output << "insert " << name() << ":Incoming " << ix << " " << _incoming[ix] << "\n"; output << "insert " << name() << ":OutgoingVector " << ix << " " << _outgoingV[ix] << "\n"; output << "insert " << name() << ":OutgoingScalar " << ix << " " << _outgoingP[ix] << "\n"; output << "insert " << name() << ":Coupling " << ix << " " << _coupling[ix]*GeV2 << "\n"; output << "insert " << name() << ":MaxWeight " << ix << " " << _maxweight[ix] << "\n"; } } if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Decay/VectorMeson/VectorMeson2BaryonsDecayer.h b/Decay/VectorMeson/VectorMeson2BaryonsDecayer.h --- a/Decay/VectorMeson/VectorMeson2BaryonsDecayer.h +++ b/Decay/VectorMeson/VectorMeson2BaryonsDecayer.h @@ -1,230 +1,230 @@ // -*- C++ -*- #ifndef Herwig_VectorMeson2BaryonsDecayer_H #define Herwig_VectorMeson2BaryonsDecayer_H // // This is the declaration of the VectorMeson2BaryonsDecayer class. // #include "Herwig/Decay/DecayIntegrator.h" #include "Herwig/Decay/PhaseSpaceMode.h" #include "ThePEG/Helicity/LorentzPolarizationVector.h" #include "ThePEG/Helicity/LorentzSpinorBar.h" #include "ThePEG/Helicity/LorentzRSSpinorBar.h" namespace Herwig { using namespace ThePEG; /** \ingroup Decay * * The VectorMeson2BaryonsDecayer class is designed for the decay * of a vector meson to a baryon-antibaryon pair. * * In this case the matrix element is taken to have the form * \f[\mathcal{M} = e_g \epsilon_\mu \bar{u}(p_f) \left[\gamma^\mu \right]u(p_{\bar{f}}).\f] * * The incoming vector mesons together with their decay products and the coupling * \f$e_g\f$, \f$G_E\f$ and the phase \f$\phi_E\f$ can be specified using the interfaces for the class. * The maximum weights * for the decays can be calculated using the Initialize interface of the * DecayIntegrator class or specified using the interface. * * The incoming and outgoing particles, couplings and maximum weights for * many of the common \f$V\to f\bar{f}\f$ decays are specified in the default * constructor. * * @see DecayIntegrator * @see \ref VectorMeson2BaryonsDecayerInterfaces "The interfaces" * defined for VectorMeson2BaryonsDecayer. * */ class VectorMeson2BaryonsDecayer: public DecayIntegrator { public: /** * The default constructor. */ VectorMeson2BaryonsDecayer(); /** * 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; /** * 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 outgoing The particles produced in the decay * @param momenta The momenta of 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. */ double me2(const int ichan,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const; /** * Construct the SpinInfos for the particles produced in the decay */ virtual void constructSpinInfo(const Particle & part, ParticleVector outgoing) const; /** * Output the setup information for the particle database * @param os The stream to output the information to * @param header Whether or not to output the information for MySQL */ virtual void dataBaseOutput(ofstream & os,bool header) 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 and * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object to the begining of the run phase. */ virtual void doinitrun(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - VectorMeson2BaryonsDecayer & operator=(const VectorMeson2BaryonsDecayer &); + VectorMeson2BaryonsDecayer & operator=(const VectorMeson2BaryonsDecayer &) = delete; private: /** * \f$G_M\f$ coupling */ vector gm_; /** * \f$G_E\f$ coupling */ vector ge_; /** * Relative phase of \f$G_E\f$/\f$G_M\f$ */ vector phi_; /** * the PDG codes for the incoming particles */ vector incoming_; /** * the PDG codes for the outgoing fermion */ vector outgoingf_; /** * the PDG codes for the outgoing antifermion. */ vector outgoinga_; /** * maximum weight for a decay */ vector maxweight_; /** * Initial size of the vectors */ unsigned int initsize_; /** * Spin density matrix */ mutable RhoDMatrix rho_; /** * Polarization vectors for the decaying particle */ mutable vector vectors_; /** * Spinors for the decay products */ mutable vector > wave_; /** * barred spinors for the decay products */ mutable vector > wavebar_; /** * Spinors for the decay products */ mutable vector > wave2_; /** * barred spinors for the decay products */ mutable vector > wave2bar_; }; } #endif /* Herwig_VectorMeson2BaryonsDecayer_H */ diff --git a/Decay/VectorMeson/f1FourPiDecayer.cc b/Decay/VectorMeson/f1FourPiDecayer.cc --- a/Decay/VectorMeson/f1FourPiDecayer.cc +++ b/Decay/VectorMeson/f1FourPiDecayer.cc @@ -1,313 +1,313 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the f1FourPiDecayer class. // #include "f1FourPiDecayer.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/epsilon.h" using namespace Herwig; f1FourPiDecayer::f1FourPiDecayer() : gRhoPiPi_(6.), ga1RhoPi_(4.8*GeV), gf1a1Pi_(9.77/GeV), maxWeight_({1.,1.}) { generateIntermediates(true); } // normally not implemented but do it here to get rid of the a_1 which // can't be on-shell ParticleVector f1FourPiDecayer::decay(const Particle & parent, const tPDVector & children) const { ParticleVector output = DecayIntegrator::decay(parent,children); ParticleVector::iterator it =output.begin(); for(;it!=output.end();++it) { long id = (**it).id(); if(id==20113 || id==20213 || id==-20213) break; } if(it!=output.end()) { PPtr a1 = *it; output.erase(it); for(PPtr child : a1->children()) { output.push_back(child); a1->abandonChild(child); } } return output; } IBPtr f1FourPiDecayer::clone() const { return new_ptr(*this); } IBPtr f1FourPiDecayer::fullclone() const { return new_ptr(*this); } void f1FourPiDecayer::persistentOutput(PersistentOStream & os) const { os << gRhoPiPi_ << ounit(ga1RhoPi_,GeV) << ounit(gf1a1Pi_,1./GeV) << ounit(ma1_,GeV) << ounit(ga1_,GeV) << ounit(mrho_,GeV) << ounit(grho_,GeV) << maxWeight_; } void f1FourPiDecayer::persistentInput(PersistentIStream & is, int) { is >> gRhoPiPi_ >> iunit(ga1RhoPi_,GeV) >> iunit(gf1a1Pi_,1./GeV) >> iunit(ma1_,GeV) >> iunit(ga1_,GeV) >> iunit(mrho_,GeV) >> iunit(grho_,GeV) >> maxWeight_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigf1FourPiDecayer("Herwig::f1FourPiDecayer", "HwVMDecay.so"); void f1FourPiDecayer::Init() { static ClassDocumentation documentation ("The f1FourPiDecayer class implements a simple model for " "f1 -> pipirho via an intermediate a_1"); static ParVector interfaceMaxWeight ("MaxWeight", "Maximum weights for the decays", &f1FourPiDecayer::maxWeight_, 2, 1.0, 0.0, 100.0, false, false, Interface::limited); static Parameter interfacegRhoPiPi ("gRhoPiPi", "The coupling of the rho to two pions", &f1FourPiDecayer::gRhoPiPi_, 6., 0.0, 10.0, false, false, Interface::limited); static Parameter interfacega1RhoPi ("ga1RhoPi", "Coupling of the a_1 to rho pi", &f1FourPiDecayer::ga1RhoPi_, GeV, 4.8*GeV, 0.*GeV, 20.*GeV, false, false, Interface::limited); static Parameter interfacegf1a1Pi ("gf1a1Pi", "Coupling of f_1 to a_1 pi", &f1FourPiDecayer::gf1a1Pi_, 1./GeV, 1.0/GeV, 0.0/GeV, 10.0/GeV, false, false, Interface::limited); } void f1FourPiDecayer::doinit() { DecayIntegrator::doinit(); // pointers to the particles we need as external particles tPDPtr f1 = getParticleData(ParticleID::f_1); tPDPtr a1p = getParticleData(ParticleID::a_1plus); tPDPtr a1m = getParticleData(ParticleID::a_1minus); tPDPtr a10 = getParticleData(ParticleID::a_10); tPDPtr pip = getParticleData(ParticleID::piplus); tPDPtr pim = getParticleData(ParticleID::piminus); tPDPtr pi0 = getParticleData(ParticleID::pi0); tPDPtr rhop = getParticleData(ParticleID::rhoplus); tPDPtr rhom = getParticleData(ParticleID::rhominus); tPDPtr rho0 = getParticleData(ParticleID::rho0); // decay mode f_1 -> pi+ pi- pi+ pi- PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(f1,{pip,pim,pip,pim},maxWeight_[0])); mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,rho0,2,3,2,4)); mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,1,1,rho0,2,3,2,4)); mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,3,1,2,1,rho0,2,1,2,4)); mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,3,1,rho0,2,1,2,4)); mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,4,1,rho0,2,3,2,2)); mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,4,1,1,1,rho0,2,3,2,2)); mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,3,1,4,1,rho0,2,1,2,2)); mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,4,1,3,1,rho0,2,1,2,2)); addMode(mode); // decay mode f_1 -> pi+ pi0 pi- pi0 mode = new_ptr(PhaseSpaceMode(f1,{pip,pi0,pim,pi0},maxWeight_[0])); mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,rhop,2,3,2,4)); mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,1,1,rhop,2,3,2,4)); mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,4,1,rhop,2,3,2,2)); mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,4,1,1,1,rhop,2,3,2,2)); mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,3,1,4,1,rhom,2,1,2,2)); mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,4,1,3,1,rhom,2,1,2,2)); mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,3,1,2,1,rhom,2,1,2,4)); mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,3,1,rhom,2,1,2,4)); addMode(mode); // masses of intermediates ma1_ = a1p->mass(); ga1_ = a1p->width(); mrho_ = rhop->mass(); grho_ = rhop->width(); } void f1FourPiDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()) { for(unsigned int ix=0;ixmaxWeight(); } } int f1FourPiDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { if(children.size()!=4 || parent->id()!=ParticleID::f_1) return -1; // check the pions int npi0(0),npiplus(0),npiminus(0); for(auto child : children) { int idtemp= child->id(); if(idtemp==ParticleID::piplus) ++npiplus; else if(idtemp==ParticleID::piminus) ++npiminus; else if(idtemp==ParticleID::pi0) ++npi0; } cc = false; // f_1 -> 2pi+ 2pi-mode if(npiplus==2 && npiminus==2) return 0; // f_1 -> pi+ pi- pi0 pi0 mode else if (npiplus ==1 && npiminus==1 && npi0==2) return 1; // not found return -1; } void f1FourPiDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { VectorWaveFunction::constructSpinInfo(vector_,const_ptr_cast(&part), incoming,true,false); // set up the spin information for the pions for(unsigned int ix=0;ix<4;++ix) ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true); } double f1FourPiDecayer::me2(const int ichan, const Particle & part, const tPDVector & , const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin0,PDT::Spin0,PDT::Spin0,PDT::Spin0))); useMe(); // polarization vectors for incoming if(meopt==Initialize) { VectorWaveFunction::calculateWaveFunctions(vector_,rho_, const_ptr_cast(&part), incoming,false); } // breit wigners Lorentz5Momentum pa1[4] = {part.momentum()-momenta[0],part.momentum()-momenta[1], part.momentum()-momenta[2],part.momentum()-momenta[3]}; complex bwa1[4]; for(unsigned int ix=0;ix<4;++ix) { pa1[ix].rescaleMass(); bwa1[ix] = Resonance::BreitWignera1(pa1[ix].mass2(),ma1_,ga1_)/sqr(ma1_); } // compute the matrix element (as a current and then contract) LorentzVector > current; // decay mode f_1 -> pi+ pi- pi+pi- double sym(0.5); if(imode()==0) { sym=0.25; complex bwrho[4] = {Resonance::BreitWignerPWave((momenta[0]+momenta[1]).m2(),mrho_,grho_,momenta[0].mass(),momenta[1].mass())/sqr(mrho_), Resonance::BreitWignerPWave((momenta[0]+momenta[3]).m2(),mrho_,grho_,momenta[0].mass(),momenta[3].mass())/sqr(mrho_), Resonance::BreitWignerPWave((momenta[2]+momenta[1]).m2(),mrho_,grho_,momenta[2].mass(),momenta[1].mass())/sqr(mrho_), Resonance::BreitWignerPWave((momenta[2]+momenta[3]).m2(),mrho_,grho_,momenta[2].mass(),momenta[3].mass())/sqr(mrho_)}; if(ichan<=0) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,rho0,2,3,2,4)); current += bwa1[0]*bwrho[3]*epsilon(part.momentum(),pa1[0],momenta[2]-momenta[3]); } if(ichan<0||ichan==1) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,1,1,rho0,2,3,2,4)); current -= bwa1[1]*bwrho[3]*epsilon(part.momentum(),pa1[1],momenta[2]-momenta[3]); } if(ichan<0||ichan==2) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,3,1,2,1,rho0,2,1,2,4)); current += bwa1[2]*bwrho[1]*epsilon(part.momentum(),pa1[2],momenta[0]-momenta[3]); } if(ichan<0||ichan==3) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,3,1,rho0,2,1,2,4)); current -= bwa1[1]*bwrho[1]*epsilon(part.momentum(),pa1[1],momenta[0]-momenta[3]); } if(ichan<0||ichan==4) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,4,1,rho0,2,3,2,2)); current += bwa1[0]*bwrho[2]*epsilon(part.momentum(),pa1[0],momenta[2]-momenta[1]); } if(ichan<0||ichan==5) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,4,1,1,1,rho0,2,3,2,2)); current -= bwa1[3]*bwrho[2]*epsilon(part.momentum(),pa1[3],momenta[2]-momenta[1]); } if(ichan<0||ichan==6) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,3,1,4,1,rho0,2,1,2,2)); current += bwa1[2]*bwrho[0]*epsilon(part.momentum(),pa1[2],momenta[0]-momenta[1]); } if(ichan<0||ichan==7) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,4,1,3,1,rho0,2,1,2,2)); current -= bwa1[3]*bwrho[0]*epsilon(part.momentum(),pa1[3],momenta[0]-momenta[1]); } } // decay mode f_1 -> pi+ pi0 pi- pi0 else { complex bwrho[4] = {Resonance::BreitWignerPWave((momenta[0]+momenta[1]).m2(),mrho_,grho_,momenta[0].mass(),momenta[1].mass())/sqr(mrho_), Resonance::BreitWignerPWave((momenta[0]+momenta[3]).m2(),mrho_,grho_,momenta[0].mass(),momenta[3].mass())/sqr(mrho_), Resonance::BreitWignerPWave((momenta[2]+momenta[1]).m2(),mrho_,grho_,momenta[2].mass(),momenta[1].mass())/sqr(mrho_), Resonance::BreitWignerPWave((momenta[2]+momenta[3]).m2(),mrho_,grho_,momenta[2].mass(),momenta[3].mass())/sqr(mrho_)}; double f1 = (momenta[2].mass2()-momenta[3].mass2())/sqr(mrho_); double f2 = 1+f1; f1 = 1.-f1; if(ichan<=0) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,rhop,2,3,2,4)); current += bwa1[0]*bwrho[3]*epsilon(part.momentum(),pa1[0],f1*momenta[2]-f2*momenta[3]); } if(ichan<0||ichan==1) { // mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,1,1,rhop,2,3,2,4)); current -= bwa1[1]*bwrho[3]*epsilon(part.momentum(),pa1[1],f1*momenta[2]-f2*momenta[3]); } if(ichan<0||ichan==2) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,4,1,rhop,2,3,2,2)); current += bwa1[0]*bwrho[2]*epsilon(part.momentum(),pa1[0],f1*momenta[2]-f2*momenta[1]); } if(ichan<0||ichan==3) { // mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,4,1,1,1,rhop,2,3,2,2)); current -= bwa1[3]*bwrho[2]*epsilon(part.momentum(),pa1[3],f1*momenta[2]-f2*momenta[1]); } if(ichan<0||ichan==4) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,3,1,4,1,rhom,2,1,2,2)); current += bwa1[2]*bwrho[0]*epsilon(part.momentum(),pa1[2],f1*momenta[0]-f2*momenta[1]); } if(ichan<0||ichan==5) { // mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,4,1,3,1,rhom,2,1,2,2)); current -= bwa1[3]*bwrho[0]*epsilon(part.momentum(),pa1[3],f1*momenta[0]-f2*momenta[1]); } if(ichan<0||ichan==6) { // mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,3,1,2,1,rhom,2,1,2,4)); current += bwa1[2]*bwrho[1]*epsilon(part.momentum(),pa1[2],f1*momenta[0]-f2*momenta[3]); } if(ichan<0||ichan==7) { // mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,3,1,rhom,2,1,2,4)); current -= bwa1[1]*bwrho[1]*epsilon(part.momentum(),pa1[1],f1*momenta[0]-f2*momenta[3]); } } // contract the current double pre = gRhoPiPi_*gf1a1Pi_*ga1RhoPi_; for(unsigned int ihel=0;ihel<3;++ihel) - (*ME())(ihel,0,0,0,0) = pre*current.dot(vector_[ihel])*part.mass(); + (*ME())(ihel,0,0,0,0) = Complex(pre*current.dot(vector_[ihel])*part.mass()); // matrix element return sym*ME()->contract(rho_).real(); } void f1FourPiDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class DecayIntegrator::dataBaseOutput(output,false); output << "newdef " << name() << ":gRhoPiPi " << gRhoPiPi_ << "\n"; output << "newdef " << name() << ":ga1RhoPi " << ga1RhoPi_/GeV << "\n"; output << "newdef " << name() << ":gf1a1Pi " << gf1a1Pi_*GeV << "\n"; for(unsigned int ix=0;ixchildren()) { output.push_back(child); a1->abandonChild(child); } } return output; } IBPtr f1RhoPiPiDecayer::clone() const { return new_ptr(*this); } IBPtr f1RhoPiPiDecayer::fullclone() const { return new_ptr(*this); } void f1RhoPiPiDecayer::persistentOutput(PersistentOStream & os) const { os << ounit(ga1RhoPi_,GeV) << ounit(gf1a1Pi_,1./GeV) << ounit(ma1_,GeV) << ounit(ga1_,GeV) << maxWeight_; } void f1RhoPiPiDecayer::persistentInput(PersistentIStream & is, int) { is >> iunit(ga1RhoPi_,GeV) >> iunit(gf1a1Pi_,1./GeV) >> iunit(ma1_,GeV) >> iunit(ga1_,GeV) >> maxWeight_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigf1RhoPiPiDecayer("Herwig::f1RhoPiPiDecayer", "HwVMDecay.so"); void f1RhoPiPiDecayer::Init() { static ClassDocumentation documentation ("The f1RhoPiPiDecayer class implements a simple model for " "f1 -> pipirho via an intermediate a_1"); static ParVector interfaceMaxWeight ("MaxWeight", "Maximum weights for the decays", &f1RhoPiPiDecayer::maxWeight_, 2, 1.0, 0.0, 100.0, false, false, Interface::limited); static Parameter interfacega1RhoPi ("ga1RhoPi", "Coupling of the a_1 to rho pi", &f1RhoPiPiDecayer::ga1RhoPi_, GeV, 4.8*GeV, 0.*GeV, 20.*GeV, false, false, Interface::limited); static Parameter interfacegf1a1Pi ("gf1a1Pi", "Coupling of f_1 to a_1 pi", &f1RhoPiPiDecayer::gf1a1Pi_, 1./GeV, 1.0/GeV, 0.0/GeV, 10.0/GeV, false, false, Interface::limited); } void f1RhoPiPiDecayer::doinit() { DecayIntegrator::doinit(); // pointers to the particles we need as external particles tPDPtr f1 = getParticleData(ParticleID::f_1); tPDPtr a1p = getParticleData(ParticleID::a_1plus); tPDPtr a1m = getParticleData(ParticleID::a_1minus); tPDPtr a10 = getParticleData(ParticleID::a_10); tPDPtr pip = getParticleData(ParticleID::piplus); tPDPtr pim = getParticleData(ParticleID::piminus); tPDPtr pi0 = getParticleData(ParticleID::pi0); tPDPtr rhop = getParticleData(ParticleID::rhoplus); tPDPtr rhom = getParticleData(ParticleID::rhominus); tPDPtr rho0 = getParticleData(ParticleID::rho0); // decay mode f_1 -> pi+ pi- rho0 PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(f1,{pip,pim,rho0},maxWeight_[0])); mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,3)); mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,1,1,3)); addMode(mode); // decay mode f_1 -> pi- pi0 rho+ mode = new_ptr(PhaseSpaceMode(f1,{pim,pi0,rhop},maxWeight_[1])); mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,1,1,2,1,3)); mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,1,1,3)); addMode(mode); ma1_ = a1p->mass(); ga1_ = a1p->width(); } void f1RhoPiPiDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()) { for(unsigned int ix=0;ixmaxWeight(); } } int f1RhoPiPiDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { if(children.size()!=3 || parent->id()!=ParticleID::f_1) return -1; // check the pions int npi0(0),npiplus(0),npiminus(0),idrho(0); for(auto child : children) { int idtemp= child->id(); if(idtemp==ParticleID::piplus) ++npiplus; else if(idtemp==ParticleID::piminus) ++npiminus; else if(idtemp==ParticleID::pi0) ++npi0; else idrho=idtemp; } cc = false; // f_1 -> pi+pi-rho0 mode if(idrho==113 && npiplus==1 && npiminus==1) return 0; // f_1 -> pi-pi0rho+ mode else if (idrho== 213 && npiminus==1 && npi0==1) return 1; else if (idrho==-213 && npiplus ==1 && npi0==1) { cc=true; return 1; } // not found return -1; } void f1RhoPiPiDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { VectorWaveFunction::constructSpinInfo(vectors_[0],const_ptr_cast(&part), incoming,true,false); // set up the spin information for the decay products // pions for(unsigned int ix=0;ix<2;++ix) ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true); // rho VectorWaveFunction::constructSpinInfo(vectors_[1],decay[2],outgoing,true,false); } double f1RhoPiPiDecayer::me2(const int ichan, const Particle & part, const tPDVector & , const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin0,PDT::Spin0,PDT::Spin1))); useMe(); // polarization vectors for incoming if(meopt==Initialize) { VectorWaveFunction::calculateWaveFunctions(vectors_[0],rho_, const_ptr_cast(&part), incoming,false); } // polarization vectors for rho vectors_[1].resize(3); for(unsigned int ix=0;ix<3;++ix) vectors_[1][ix] = HelicityFunctions::polarizationVector(-momenta[2],ix,Helicity::outgoing); // breit wigners Lorentz5Momentum pa1[2] = {part.momentum()-momenta[0], part.momentum()-momenta[1]}; for(unsigned int ix=0;ix<2;++ix) pa1[ix].rescaleMass(); complex bwa1[2] = {Resonance::BreitWignera1(pa1[0].mass2(),ma1_,ga1_)/sqr(ma1_), Resonance::BreitWignera1(pa1[1].mass2(),ma1_,ga1_)/sqr(ma1_)}; if(ichan>0) bwa1[ ichan == 0 ? 1 : 0 ] = ZERO; // compute the matrix element for(unsigned int ihel=0;ihel<3;++ihel) { LorentzVector > pol[2] = {epsilon(vectors_[0][ihel],part.momentum(),pa1[0]), epsilon(vectors_[0][ihel],part.momentum(),pa1[1])}; for(unsigned int ohel=0;ohel<3;++ohel) { - (*ME())(ihel,0,0,ohel) = gf1a1Pi_*ga1RhoPi_*(LorentzPolarizationVector(pol[0]*bwa1[0])- - LorentzPolarizationVector(pol[1]*bwa1[1])).dot(vectors_[1][ohel]); + (*ME())(ihel,0,0,ohel) = Complex(gf1a1Pi_*ga1RhoPi_*(LorentzPolarizationVector(pol[0]*bwa1[0])- + LorentzPolarizationVector(pol[1]*bwa1[1])).dot(vectors_[1][ohel])); } } // matrix element return ME()->contract(rho_).real(); } void f1RhoPiPiDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class DecayIntegrator::dataBaseOutput(output,false); output << "newdef " << name() << ":ga1RhoPi " << ga1RhoPi_/GeV << "\n"; output << "newdef " << name() << ":gf1a1Pi " << gf1a1Pi_*GeV << "\n"; for(unsigned int ix=0;ix describeHerwigThreePionCLEOCurrent("Herwig::ThreePionCLEOCurrent", "HwWeakCurrents.so"); HERWIG_INTERPOLATOR_CLASSDESC(ThreePionCLEOCurrent,Energy,Energy2) ThreePionCLEOCurrent::ThreePionCLEOCurrent() { addDecayMode(1,-1); addDecayMode(2,-2); addDecayMode(2,-1); addDecayMode(1,-1); addDecayMode(2,-2); addDecayMode(2,-1); setInitialModes(6); // rho masses and widths _rhomass = {0.7743*GeV,1.370*GeV}; _rhowidth = {0.1491*GeV,0.386*GeV}; // f_2 mass and width _f2mass = 1.275*GeV; _f2width = 0.185*GeV; // f_0(1370) mass and width _f0mass = 1.186*GeV; _f0width = 0.350*GeV; // sigma mass and width _sigmamass = 0.860*GeV; _sigmawidth = 0.880*GeV; // a1 mass and width _a1mass = 1.331*GeV; _a1width = 0.814*GeV; // parameters for the K K* contribution to the a_1 running width _mKstar = 894*MeV; _mK = 496*MeV; _gammk = 3.32; // pion decay constant _fpi = 130.7*MeV/sqrt(2.); // couplings and phases for the different channels // p-wave rho and rho prime using Constants::pi; _rhomagP = {1.,0.12}; _rhophaseP = {0.,0.99*pi}; // d-wave rho and rho prime _rhomagD = {0.37/GeV2,0.87/GeV2}; _rhophaseD = {-0.15*pi, 0.53*pi}; // f_2 _f2mag = 0.71/GeV2; _f2phase = 0.56*pi; _f2coup = ZERO; // sigma _sigmamag = 2.10; _sigmaphase = 0.23*pi; _sigmacoup = 0.; // f_0 _f0mag = 0.77; _f0phase = -0.54*pi; _f0coup = 0.; // initialize the a_1 width _initializea1=false; _a1opt=true; double a1q2in[200]={0 ,15788.6,31577.3,47365.9,63154.6,78943.2, 94731.9,110521 ,126309 ,142098 ,157886 ,173675 , 189464 ,205252 ,221041 ,236830 ,252618 ,268407 , 284196 ,299984 ,315773 ,331562 ,347350 ,363139 , 378927 ,394716 ,410505 ,426293 ,442082 ,457871 , 473659 ,489448 ,505237 ,521025 ,536814 ,552603 , 568391 ,584180 ,599969 ,615757 ,631546 ,647334 , 663123 ,678912 ,694700 ,710489 ,726278 ,742066 , 757855 ,773644 ,789432 ,805221 ,821010 ,836798 , 852587 ,868375 ,884164 ,899953 ,915741 ,931530 , 947319 ,963107 ,978896 ,994685 , 1.01047e+06,1.02626e+06,1.04205e+06,1.05784e+06, 1.07363e+06,1.08942e+06,1.10521e+06,1.12099e+06, 1.13678e+06,1.15257e+06,1.16836e+06,1.18415e+06, 1.19994e+06,1.21573e+06,1.23151e+06,1.24730e+06, 1.26309e+06,1.27888e+06,1.29467e+06,1.31046e+06, 1.32625e+06,1.34203e+06,1.35782e+06,1.37361e+06, 1.38940e+06,1.40519e+06,1.42098e+06,1.43677e+06, 1.45256e+06,1.46834e+06,1.48413e+06,1.49992e+06, 1.51571e+06,1.53150e+06,1.54729e+06,1.56308e+06, 1.57886e+06,1.59465e+06,1.61044e+06,1.62623e+06, 1.64202e+06,1.65781e+06,1.67360e+06,1.68939e+06, 1.70517e+06,1.72096e+06,1.73675e+06,1.75254e+06, 1.76833e+06,1.78412e+06,1.79991e+06,1.81569e+06, 1.83148e+06,1.84727e+06,1.86306e+06,1.87885e+06, 1.89464e+06,1.91043e+06,1.92621e+06,1.94200e+06, 1.95779e+06,1.97358e+06,1.98937e+06,2.00516e+06, 2.02095e+06,2.03674e+06,2.05252e+06,2.06831e+06, 2.08410e+06,2.09989e+06,2.11568e+06,2.13147e+06, 2.14726e+06,2.16304e+06,2.17883e+06,2.19462e+06, 2.21041e+06,2.22620e+06,2.24199e+06,2.25778e+06, 2.27356e+06,2.28935e+06,2.30514e+06,2.32093e+06, 2.33672e+06,2.35251e+06,2.36830e+06,2.38409e+06, 2.39987e+06,2.41566e+06,2.43145e+06,2.44724e+06, 2.46303e+06,2.47882e+06,2.49461e+06,2.51039e+06, 2.52618e+06,2.54197e+06,2.55776e+06,2.57355e+06, 2.58934e+06,2.60513e+06,2.62092e+06,2.63670e+06, 2.65249e+06,2.66828e+06,2.68407e+06,2.69986e+06, 2.71565e+06,2.73144e+06,2.74722e+06,2.76301e+06, 2.77880e+06,2.79459e+06,2.81038e+06,2.82617e+06, 2.84196e+06,2.85774e+06,2.87353e+06,2.88932e+06, 2.90511e+06,2.92090e+06,2.93669e+06,2.95248e+06, 2.96827e+06,2.98405e+06,2.99984e+06,3.01563e+06, 3.03142e+06,3.04721e+06,3.06300e+06,3.07879e+06, 3.09457e+06,3.11036e+06,3.12615e+06,3.14194e+06}; double a1widthin[200]={0,0,0,0,0,0,0,0, 0,0,0,0.00021256,0.0107225,0.0554708,0.150142,0.303848, 0.522655,0.81121,1.1736,1.61381,2.13606,2.74499,3.44583,4.24454, 5.14795,6.16391,7.3014,8.57079,9.98398,11.5547,13.2987,15.2344, 17.3827,19.7683,22.4195,25.3695,28.6568,32.3264,36.4311,41.0322, 46.201,52.0203,58.5847,66.0011,74.3871,83.8666,94.5615,106.578, 119.989,134.807,150.968,168.315,186.615,205.576,224.893,244.28, 263.499,282.364,300.748,318.569,335.781,352.367,368.327,383.677, 398.438,412.638,426.306,439.472,452.167,464.421,476.263,487.719, 498.815,509.576,520.024,530.179,540.063,549.693,559.621,568.26, 577.229,586.005,594.604,603.035,611.314,619.447,627.446,635.321, 643.082,650.736,658.288,665.75,673.127,680.427,687.659,694.82, 701.926,708.977,715.983,722.944,729.862,736.752,743.619,750.452, 757.271,764.076,770.874,777.658,784.444,791.233,798.027,804.838, 811.649,818.485,825.342,832.224,839.139,846.082,853.059,860.079, 867.143,874.248,881.409,919.527,945.28,965.514,983.228,999.471, 1014.69,1029.15,1043.05,1056.49,1069.57,1082.36,1094.88,1107.2, 1120.89,1131.4,1143.33,1155.15,1166.92,1178.61,1190.27,1201.92, 1213.55,1225.18,1236.81,1250.06,1260.16,1271.86,1283.64,1295.46, 1307.36,1319.3,1331.34,1343.45,1355.64,1367.93,1380.31,1392.77, 1405.35,1418.03,1430.83,1443.75,1457.17,1469.94,1483.22,1496.64, 1510.18,1523.86,1537.67,1551.64,1565.72,1579.99,1594.38,1608.92, 1623.63,1642.08,1653.51,1668.69,1684.03,1699.53,1715.21,1731.04, 1747.05,1763.23,1779.59,1796.12,1812.83,1829.72,1846.79,1864.04, 1881.49,1899.11,1916.93,1934.93,1953.13,1971.52,1990.12,2008.89}; if(_a1runwidth.empty()) { vector tmp1(a1widthin,a1widthin+200); std::transform(tmp1.begin(), tmp1.end(), back_inserter(_a1runwidth), [](double x){return x*GeV;}); vector tmp2(a1q2in,a1q2in+200); _a1runq2.clear(); std::transform(tmp2.begin(), tmp2.end(), back_inserter(_a1runq2), [](double x){return x*GeV2;}); } // zero parameters which will be calculated later to avoid problems _mpi0=ZERO; _mpic=ZERO; _fact=ZERO; _maxmass=ZERO; _maxcalc=ZERO; } void ThreePionCLEOCurrent::doinit() { WeakCurrent::doinit(); // parameters for the breit-wigners _mpic = getParticleData(ParticleID::piplus)->mass(); _mpi0 = getParticleData(ParticleID::pi0) ->mass(); // couplings for the different modes Complex ii(0.,1.); _rhocoupP.resize(_rhomagP.size()); for(unsigned int ix=0;ix<_rhomagP.size();++ix) _rhocoupP[ix]=_rhomagP[ix]*(cos(_rhophaseP[ix])+ii*sin(_rhophaseP[ix])); _rhocoupD.resize(_rhomagD.size()); for(unsigned int ix=0;ix<_rhomagD.size();++ix) _rhocoupD[ix]=_rhomagD[ix]*(cos(_rhophaseD[ix])+ii*sin(_rhophaseD[ix])); _f0coup=_f0mag*(cos(_f0phase)+ii*sin(_f0phase)); _f2coup=_f2mag*(cos(_f2phase)+ii*sin(_f2phase)); _sigmacoup=_sigmamag*(cos(_sigmaphase)+ii*sin(_sigmaphase)); // overall coupling _fact = 2.*sqrt(2.)/_fpi/3.; // initialise the a_1 running width calculation inita1Width(-1); } void ThreePionCLEOCurrent::persistentOutput(PersistentOStream & os) const { os << ounit(_rhomass,GeV) << ounit(_rhowidth,GeV) << ounit(_f2mass,GeV) << ounit(_f2width,GeV) << ounit(_f0mass,GeV) << ounit(_f0width,GeV) << ounit(_sigmamass,GeV) << ounit(_sigmawidth,GeV) << ounit(_mpi0,GeV) << ounit(_mpic,GeV) << ounit(_fpi,GeV) << ounit(_fact,1/GeV) << _rhomagP << _rhophaseP << _rhocoupP << ounit(_rhomagD,1/GeV2) << _rhophaseD << ounit(_rhocoupD,1/GeV2) <> iunit(_rhomass,GeV) >> iunit(_rhowidth,GeV) >> iunit(_f2mass,GeV) >> iunit(_f2width,GeV) >> iunit(_f0mass,GeV) >> iunit(_f0width,GeV) >> iunit(_sigmamass,GeV) >> iunit(_sigmawidth,GeV) >> iunit(_mpi0,GeV) >> iunit(_mpic,GeV) >> iunit(_fpi,GeV) >> iunit(_fact,1/GeV) >> _rhomagP >> _rhophaseP >> _rhocoupP >> iunit(_rhomagD,1/GeV2) >> _rhophaseD >> iunit(_rhocoupD,1/GeV2) >> iunit(_f2mag,1/GeV2) >> _f2phase >> iunit(_f2coup,1/GeV2) >> _f0mag >> _f0phase >> _f0coup >> _sigmamag >> _sigmaphase >> _sigmacoup >> iunit(_a1mass,GeV) >> iunit(_a1width,GeV) >> iunit(_a1runwidth,GeV) >> iunit(_a1runq2,GeV2) >> _initializea1 >> iunit(_mKstar,GeV) >> iunit(_mK,GeV) >> _gammk >> _a1opt >> iunit(_maxmass,GeV) >> iunit(_maxcalc,GeV) >> _a1runinter; } void ThreePionCLEOCurrent::Init() { static ClassDocumentation documentation ("The ThreePionCLEOCurrent class performs the decay of the" " tau to three pions using the currents from CLEO", "The decay of tau to three pions is modelled using the currents from " "\\cite{Asner:1999kj}.", " %\\cite{Asner:1999kj}\n" "\\bibitem{Asner:1999kj}\n" " D.~M.~Asner {\\it et al.} [CLEO Collaboration],\n" " ``Hadronic structure in the decay tau- --> nu/tau pi- pi0 pi0 and the sign\n" " %of the tau neutrino helicity,''\n" " Phys.\\ Rev.\\ D {\\bf 61}, 012002 (2000)\n" " [arXiv:hep-ex/9902022].\n" " %%CITATION = PHRVA,D61,012002;%%\n" ); static ParVector interfacerhomass ("RhoMasses", "The masses of the different rho resonnaces", &ThreePionCLEOCurrent::_rhomass, MeV, 0, ZERO, -10000*MeV, 10000*MeV, false, false, true); static ParVector interfacerhowidth ("RhoWidths", "The widths of the different rho resonnaces", &ThreePionCLEOCurrent::_rhowidth, MeV, 0, ZERO, -10000*MeV, 10000*MeV, false, false, true); static Parameter interfacef_2Mass ("f_2Mass", "The mass of the f_2 meson", &ThreePionCLEOCurrent::_f2mass, GeV, 1.275*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfacef_2Width ("f_2Width", "The width of the f_2 meson", &ThreePionCLEOCurrent::_f2width, GeV, 0.185*GeV, ZERO, 1.0*GeV, false, false, true); static Parameter interfacef_0Mass ("f_0Mass", "The mass of the f_0 meson", &ThreePionCLEOCurrent::_f0mass, GeV, 1.186*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfacef_0Width ("f_0Width", "The width of the f_0 meson", &ThreePionCLEOCurrent::_f0width, GeV, 0.350*GeV, ZERO, 1.0*GeV, false, false, true); static Parameter interfacesigmaMass ("sigmaMass", "The mass of the sigma meson", &ThreePionCLEOCurrent::_sigmamass, GeV, 0.860*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfacesigmaWidth ("sigmaWidth", "The width of the sigma meson", &ThreePionCLEOCurrent::_sigmawidth, GeV, 0.880*GeV, ZERO, 2.0*GeV, false, false, true); static Parameter interfacea1Mass ("a1Mass", "The mass of the a_1 meson", &ThreePionCLEOCurrent::_a1mass, GeV, 1.331*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfacea1Width ("a1Width", "The width of the a_1 meson", &ThreePionCLEOCurrent::_a1width, GeV, 0.814*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceKaonMass ("KaonMass", "The mass of the kaon", &ThreePionCLEOCurrent::_mK, GeV, 0.496*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceKStarMass ("KStarMass", "The mass of the k* meson", &ThreePionCLEOCurrent::_mKstar, GeV, 0.894*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceKaonCoupling ("KaonCoupling", "The relative coupling for the kaon in the a_1 running width", &ThreePionCLEOCurrent::_gammk, 3.32, 0.0, 10.0, false, false, true); static Parameter interfaceFpi ("Fpi", "The pion decay constant", &ThreePionCLEOCurrent::_fpi, MeV, 130.7*MeV/sqrt(2.), ZERO, 500.0*MeV, false, false, true); static ParVector interfacerhomagP ("RhoPWaveMagnitude", "The magnitude of the couplings for the p-wave rho currents", &ThreePionCLEOCurrent::_rhomagP, 0, 0, 0, 0, 10000, false, false, true); static ParVector interfacerhophaseP ("RhoPWavePhase", "The phase of the couplings for the p-wave rho currents", &ThreePionCLEOCurrent::_rhophaseP, 0, 0, 0, -Constants::twopi, Constants::twopi, false, false, true); static ParVector interfacerhomagD ("RhoDWaveMagnitude", "The magnitude of the couplings for the d-wave rho currents", &ThreePionCLEOCurrent::_rhomagD, 1/MeV2, 0, ZERO, ZERO, 10000/MeV2, false, false, true); static ParVector interfacerhophaseD ("RhoDWavePhase", "The phase of the couplings for the d-wave rho currents", &ThreePionCLEOCurrent::_rhophaseD, 0, 0, 0, -Constants::twopi, Constants::twopi, false, false, true); static Parameter interfacef0Phase ("f0Phase", "The phase of the f_0 scalar current", &ThreePionCLEOCurrent::_f0phase, 0.54*Constants::pi, -Constants::twopi, Constants::twopi, false, false, true); static Parameter interfacef2Phase ("f2Phase", "The phase of the f_2 tensor current", &ThreePionCLEOCurrent::_f2phase, 0.56*Constants::pi,-Constants::twopi, Constants::twopi, false, false, true); static Parameter interfacesigmaPhase ("sigmaPhase", "The phase of the sigma scalar current", &ThreePionCLEOCurrent::_sigmaphase, 0.23*Constants::pi, -Constants::twopi, Constants::twopi, false, false, true); static Parameter interfacef0Magnitude ("f0Magnitude", "The magnitude of the f_0 scalar current", &ThreePionCLEOCurrent::_f0mag, 0.77, 0.0, 10, false, false, true); static Parameter interfacef2Magnitude ("f2Magnitude", "The magnitude of the f_2 tensor current", &ThreePionCLEOCurrent::_f2mag, 1./GeV2, 0.71/GeV2, ZERO, 10./GeV2, false, false, true); static Parameter interfacesigmaMagnitude ("sigmaMagnitude", "The magnitude of the sigma scalar current", &ThreePionCLEOCurrent::_sigmamag, 2.1, 0.0, 10, false, false, true); static ParVector interfacea1RunningWidth ("a1RunningWidth", "The values of the a_1 width for interpolation to giving the running width.", &ThreePionCLEOCurrent::_a1runwidth, MeV, 0, ZERO, ZERO, 10000000*MeV, false, false, true); static ParVector interfacea1RunningQ2 ("a1RunningQ2", "The values of the q^2 for interpolation to giving the running width.", &ThreePionCLEOCurrent::_a1runq2, MeV2, 0, ZERO, ZERO, 10000000*MeV2, false, false, true); static Switch interfaceInitializea1 ("Initializea1", "Initialise the calculation of the a_1 running width", &ThreePionCLEOCurrent::_initializea1, false, false, false); static SwitchOption interfaceInitializea1Initialization (interfaceInitializea1, "Yes", "Initialize the calculation", true); static SwitchOption interfaceInitializea1NoInitialization (interfaceInitializea1, "No", "Use the default values", false); static Switch interfacea1WidthOption ("a1WidthOption", "Option for the treatment of the a1 width", &ThreePionCLEOCurrent::_a1opt, true, false, false); static SwitchOption interfacea1WidthOptionLocal (interfacea1WidthOption, "Local", "Use a calculation of the running width based on the parameters as" " interpolation table.", true); static SwitchOption interfacea1WidthOptionParam (interfacea1WidthOption, "Kuhn", "Use the parameterization of Kuhn and Santamaria for default parameters." " This should only be used for testing vs TAUOLA", false); } // initialisation of the a_1 width // (iopt=-1 initialises, iopt=0 starts the interpolation) void ThreePionCLEOCurrent::inita1Width(int iopt) { if(iopt==-1) { _maxcalc=_maxmass; if(!_initializea1||_maxmass==ZERO) return; // parameters for the table of values Energy2 step=sqr(_maxmass)/200.; // function to be integrated to give the matrix element // integrator to perform the integral vector inweights;inweights.push_back(0.5);inweights.push_back(0.5); vector intype;intype.push_back(2);intype.push_back(3); Energy mrho=getParticleData(ParticleID::rhoplus)->mass(); Energy wrho=getParticleData(ParticleID::rhoplus)->width(); vector inmass;inmass.push_back(mrho);inmass.push_back(mrho); vector inwidth;inwidth.push_back(wrho);inwidth.push_back(wrho); vector inpow(2,0.0); ThreeBodyAllOnCalculator widthgenN(inweights,intype,inmass,inwidth,inpow,*this,0,_mpi0,_mpi0,_mpic); ThreeBodyAllOnCalculator widthgenC(inweights,intype,inmass,inwidth,inpow,*this,1,_mpic,_mpic,_mpic); // normalisation constant to give physical width if on shell double a1const = _a1width/(widthgenN.partialWidth(sqr(_a1mass))+ widthgenC.partialWidth(sqr(_a1mass))); // loop to give the values _a1runq2.clear();_a1runwidth.clear(); for(Energy2 moff2=ZERO; moff2<=sqr(_maxmass); moff2+=step) { Energy moff=sqrt(moff2); _a1runq2.push_back(moff2); Energy charged=a1const*widthgenC.partialWidth(moff2); Energy neutral=a1const*widthgenN.partialWidth(moff2); Energy kaon = moff<=_mK+_mKstar ? ZERO : 2.870*_gammk*_gammk/8./Constants::pi* Kinematics::pstarTwoBodyDecay(moff,_mK,_mKstar)/moff2*GeV2; Energy total = charged + neutral + kaon; _a1runwidth.push_back(total); } } // set up the interpolator else if(iopt==0) { _a1runinter = make_InterpolatorPtr(_a1runwidth,_a1runq2,3); } } void ThreePionCLEOCurrent::CLEOFormFactor(int imode,int ichan, Energy2 q2,Energy2 s1, Energy2 s2, Energy2 s3, Complex & F1, Complex & F2, Complex & F3) const { useMe(); double fact=1.; if(imode<=1) { // identical particle factors fact = 1./sqrt(6.); // compute the breit wigners we need Complex sigbws1 = Resonance::BreitWignerSWave(s1,_sigmamass,_sigmawidth,_mpi0,_mpi0); Complex sigbws2 = Resonance::BreitWignerSWave(s2,_sigmamass,_sigmawidth,_mpi0,_mpi0); Complex sigbws3 = Resonance::BreitWignerSWave(s3,_sigmamass,_sigmawidth,_mpi0,_mpi0); Complex f0bws1 = Resonance::BreitWignerSWave(s1, _f0mass, _f0width,_mpi0,_mpi0); Complex f0bws2 = Resonance::BreitWignerSWave(s2, _f0mass, _f0width,_mpi0,_mpi0); Complex f0bws3 = Resonance::BreitWignerSWave(s3, _f0mass, _f0width,_mpi0,_mpi0); Complex f2bws1 = Resonance::BreitWignerDWave(s1, _f2mass, _f2width,_mpi0,_mpi0); Complex f2bws2 = Resonance::BreitWignerDWave(s2, _f2mass, _f2width,_mpi0,_mpi0); Complex f2bws3 = Resonance::BreitWignerDWave(s3, _f2mass, _f2width,_mpi0,_mpi0); if(ichan<0) { // the scalar terms F1=2./3.*(_sigmacoup*sigbws3+_f0coup*f0bws3) -2./3.*(_sigmacoup*sigbws2+_f0coup*f0bws2); F2=2./3.*(_sigmacoup*sigbws3+_f0coup*f0bws3) -2./3.*(_sigmacoup*sigbws1+_f0coup*f0bws1); F3=-2./3.*(_sigmacoup*sigbws1+_f0coup*f0bws1) +2./3.*(_sigmacoup*sigbws2+_f0coup*f0bws2); // the tensor terms complex Dfact1 = 1./18.*(4.*_mpi0*_mpi0-s1)*(q2+s1-_mpi0*_mpi0)/s1*f2bws1; complex Dfact2 = 1./18.*(4.*_mpi0*_mpi0-s2)*(q2+s2-_mpi0*_mpi0)/s2*f2bws2; complex Dfact3 = 1./18.*(4.*_mpi0*_mpi0-s3)*(q2-_mpi0*_mpi0+s3)/s3*f2bws3; - F1+=_f2coup*( 0.5*(s3-s2)*f2bws1-Dfact2+Dfact3); - F2+=_f2coup*( 0.5*(s3-s1)*f2bws2-Dfact1+Dfact3); - F3+=_f2coup*(-0.5*(s1-s2)*f2bws3-Dfact1+Dfact2); + F1+=Complex(_f2coup*( 0.5*(s3-s2)*f2bws1-Dfact2+Dfact3)); + F2+=Complex(_f2coup*( 0.5*(s3-s1)*f2bws2-Dfact1+Dfact3)); + F3+=Complex(_f2coup*(-0.5*(s1-s2)*f2bws3-Dfact1+Dfact2)); } else if(ichan==0) { F2=-2./3.*_sigmacoup*sigbws1; F3=-2./3.*_sigmacoup*sigbws1; } else if(ichan==1) { F1=-2./3.*_sigmacoup*sigbws2; F3=+2./3.*_sigmacoup*sigbws2; } else if(ichan==2) { F1= 2./3.*_sigmacoup*sigbws3; F2= 2./3.*_sigmacoup*sigbws3; } else if(ichan==3) { complex Dfact1 = 1./18.*(4.*_mpi0*_mpi0-s1)*(q2+s1-_mpi0*_mpi0)/s1*f2bws1; - F1+=_f2coup*0.5*(s3-s2)*f2bws1; - F2-=_f2coup*Dfact1; - F3-=_f2coup*Dfact1; + F1+=Complex(_f2coup*0.5*(s3-s2)*f2bws1); + F2-=Complex(_f2coup*Dfact1); + F3-=Complex(_f2coup*Dfact1); } else if(ichan==4) { complex Dfact2 = 1./18.*(4.*_mpi0*_mpi0-s2)*(q2+s2-_mpi0*_mpi0)/s2*f2bws2; - F2+=_f2coup*0.5*(s3-s1)*f2bws2; - F1-=_f2coup*Dfact2; - F3+=_f2coup*Dfact2; + F2+=Complex(_f2coup*0.5*(s3-s1)*f2bws2); + F1-=Complex(_f2coup*Dfact2); + F3+=Complex(_f2coup*Dfact2); } else if(ichan==5) { complex Dfact3 = 1./18.*(4.*_mpi0*_mpi0-s3)*(q2-_mpi0*_mpi0+s3)/s3*f2bws3; - F3+=-_f2coup*0.5*(s1-s2)*f2bws3; - F1+=_f2coup*Dfact3; - F2+=_f2coup*Dfact3; + F3+=Complex(-_f2coup*0.5*(s1-s2)*f2bws3); + F1+=Complex(_f2coup*Dfact3); + F2+=Complex(_f2coup*Dfact3); } else if(ichan==6) { F2=-2./3.*_f0coup*f0bws1; F3=-2./3.*_f0coup*f0bws1; } else if(ichan==7) { F1=-2./3.*_f0coup*f0bws2; F3=+2./3.*_f0coup*f0bws2; } else if(ichan==8) { F1= 2./3.*_f0coup*f0bws3; F2= 2./3.*_f0coup*f0bws3; } } // calculate the pi0 pi0 pi+ factor else if(imode==2) { // identical particle factors fact = 1./sqrt(2.); // compute the breit wigners we need Complex rhos1bw[3],rhos2bw[3]; for(unsigned int ix=0,N=max(_rhocoupP.size(),_rhocoupD.size());ixpi+pi-pi0 else if(imode==3||imode==4) { // compute the breit wigners we need Complex rhos1bw[3],rhos2bw[3]; for(unsigned int ix=0,N=max(_rhocoupP.size(),_rhocoupD.size());ix sfact1 = 1./18.*(4.*_mpic*_mpic-s1)*(q2+s1-_mpic*_mpic)/s1*f2bws1; complex sfact2 = 1./18.*(4.*_mpic*_mpic-s2)*(q2+s2-_mpic*_mpic)/s2*f2bws2; F1+=Complex(_f2coup*(0.5*(s3-s2)*f2bws1-sfact2)); F2+=Complex(_f2coup*(0.5*(s3-s1)*f2bws2-sfact1)); F3+=Complex(_f2coup*(-sfact1+sfact2)); } else if(ichan%2==0&&ichan<=4) { unsigned int ires=ichan/2; Energy2 Dfact2=1./3.*(s2-s3); if(ires<_rhocoupP.size()) F1-=_rhocoupP[ires]*rhos1bw[ires]; if(ires<_rhocoupD.size()) { F2-=Complex(Dfact2*_rhocoupD[ires]*rhos1bw[ires]); F3-=Complex(_rhocoupD[ires]*Dfact2*rhos1bw[ires]); } } else if(ichan%2==1&&ichan<=5) { unsigned int ires=(ichan-1)/2; Energy2 Dfact1=1./3.*(s1-s3); if(ires<_rhocoupP.size()) { F2-=_rhocoupP[ires]*rhos2bw[ires]; } if(ires<_rhocoupD.size()) { F1-=Complex(Dfact1*_rhocoupD[ires]*rhos2bw[ires]); F3+=Complex(_rhocoupD[ires]*Dfact1*rhos2bw[ires]); } } else if(ichan==6) { F2-=2./3.*_sigmacoup*sigbws1; F3-=2./3.*_sigmacoup*sigbws1; } else if(ichan==7) { F1-=2./3.*_sigmacoup*sigbws2; F3+=2./3.*_sigmacoup*sigbws2; } else if(ichan==8) { complex sfact1 = 1./18.*(4.*_mpic*_mpic-s1)*(q2+s1-_mpic*_mpic)/s1*f2bws1; F1+=Complex(_f2coup*0.5*(s3-s2)*f2bws1); F2-=Complex(_f2coup*sfact1); F3-=Complex(_f2coup*sfact1); } else if(ichan==9) { complex sfact2 = 1./18.*(4.*_mpic*_mpic-s2)*(q2+s2-_mpic*_mpic)/s2*f2bws2; F1-=Complex(_f2coup*sfact2); F2+=Complex(_f2coup*0.5*(s3-s1)*f2bws2); F3+=Complex(_f2coup*sfact2); } else if(ichan==10) { F2-=2./3.*_f0coup*f0bws1; F3-=2./3.*_f0coup*f0bws1; } else if(ichan==11) { F1-=2./3.*_f0coup*f0bws2; F3+=2./3.*_f0coup*f0bws2; } } else { throw Exception() << "ThreePionCLEOCurrent Unknown Decay" << imode << Exception::abortnow; } F1 *= fact; F2 *= fact; F3 *= fact; } // complete the construction of the decay mode for integration bool ThreePionCLEOCurrent::createMode(int icharge, tcPDPtr resonance, FlavourInfo flavour, unsigned int imode,PhaseSpaceModePtr mode, unsigned int iloc,int ires, PhaseSpaceChannel phase, Energy upp ) { // check the charge and resonance if(imode<=1||imode==3||imode==4) { if(icharge!=0) return false; if(resonance && resonance->id()!=ParticleID::a_10) return false; } else if(imode==2||imode==5) { if(abs(icharge)!=3) return false; if(resonance && abs(resonance->id())!=ParticleID::a_1plus) return false; } // check the total isospin if(flavour.I!=IsoSpin::IUnknown) { if(flavour.I!=IsoSpin::IOne) return false; } // check I_3 if(flavour.I3!=IsoSpin::I3Unknown) { switch(flavour.I3) { case IsoSpin::I3Zero: if(imode==2||imode==5) return false; break; case IsoSpin::I3One: if((imode!=2&&imode!=5) || icharge ==-3) return false; break; case IsoSpin::I3MinusOne: if((imode!=2&&imode!=5) || icharge == 3) return false; break; default: return false; } } if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return false; if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return false; if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return false; // get the particles and check the masses int iq(0),ia(0); tPDVector extpart=particles(1,imode,iq,ia); Energy min(ZERO); for(unsigned int ix=0;ixmassMin(); if(min>upp) return false; _maxmass=max(_maxmass,upp); // pointers to the particles we need tPDPtr a1m = getParticleData(ParticleID::a_1minus); tPDPtr a10 = getParticleData(ParticleID::a_10); // the different rho resonances tPDPtr rhom[3] = {getParticleData(-213),getParticleData(-100213),getParticleData(-30213)}; if(icharge==3) { for(unsigned int ix=0;ix<3;++ix) rhom[ix]=rhom[ix]->CC(); a1m = a1m->CC(); } tPDPtr rho0[3] = {getParticleData(113),getParticleData(100113),getParticleData(30113)}; // the sigma tPDPtr sigma = getParticleData(9000221); // the f_2 tPDPtr f2=getParticleData(225); // the f_0 tPDPtr f0=getParticleData(10221); assert(f2 && f0 && sigma); // a0 -> pi0 pi0 pi0 if(imode<=1) { for(unsigned int ix=0;ix<3;++ix) { tPDPtr temp; if(ix==0) temp = sigma; else if(ix==1) temp = f2; else if(ix==2) temp = f0; mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,temp,ires+1,iloc+1, ires+2,iloc+2,ires+2,iloc+3)); mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,temp,ires+1,iloc+2, ires+2,iloc+1,ires+2,iloc+3)); mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,temp,ires+1,iloc+3, ires+2,iloc+1,ires+2,iloc+2)); } } // decay mode a_1- -> pi0 pi0 pi- else if(imode==2) { for(unsigned int ix=0;ix<3;++ix) { // first rho+ channel mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,rhom[ix],ires+1,iloc+1, ires+2,iloc+2,ires+2,iloc+3)); // second rho+ channel mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,rhom[ix],ires+1,iloc+2, ires+2,iloc+1,ires+2,iloc+3)); } // I=0 channels for(unsigned int iy=0;iy<3;++iy) { tPDPtr temp; if(iy==0) temp = sigma; else if(iy==1) temp = f2; else if(iy==2) temp = f0; mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,temp,ires+1,iloc+3, ires+2,iloc+1,ires+2,iloc+2)); } } // decay mode a_10 -> pi+ pi- pi0 else if(imode==3||imode==4) { // rho modes for(unsigned int ix=0;ix<3;++ix) { // first rho channel mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,rhom[ix],ires+1,iloc+1, ires+2,iloc+2,ires+2,iloc+3)); // second channel mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,rhom[ix],ires+1,iloc+2, ires+2,iloc+1,ires+2,iloc+3)); } // I=0 channels for(unsigned int iy=0;iy<3;++iy) { tPDPtr temp; if(iy==0) temp = sigma; else if(iy==1) temp = f2; else if(iy==2) temp = f0; mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,temp,ires+1,iloc+3, ires+2,iloc+1,ires+2,iloc+2)); } } else if(imode==5) { for(unsigned int ix=0;ix<3;++ix) { // the neutral rho channels // first channel mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,rho0[ix],ires+1,iloc+1, ires+2,iloc+2,ires+2,iloc+3)); // interchanged channel mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,rho0[ix],ires+1,iloc+2, ires+2,iloc+1,ires+2,iloc+3)); } for(unsigned int iy=0;iy<3;++iy) { tPDPtr temp; if(iy==0) temp = sigma; else if(iy==1) temp = f2; else if(iy==2) temp = f0; // first channel mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,temp,ires+1,iloc+1, ires+2,iloc+2,ires+2,iloc+3)); // interchanged channel mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,temp,ires+1,iloc+2, ires+2,iloc+1,ires+2,iloc+3)); } } // reset the integration parameters for(unsigned int iy=0;iy<_rhomass.size();++iy) { mode->resetIntermediate(rho0[iy],_rhomass[iy],_rhowidth[iy]); mode->resetIntermediate(rhom[iy],_rhomass[iy],_rhowidth[iy]); } mode->resetIntermediate(sigma,_sigmamass,_sigmawidth); mode->resetIntermediate(f2,_f2mass,_f2width); mode->resetIntermediate(f0,_f0mass,_f0width); mode->resetIntermediate(a10,_a1mass,_a1width); mode->resetIntermediate(a10,_a1mass,_a1width); return true; } void ThreePionCLEOCurrent::dataBaseOutput(ofstream & output,bool header, bool create) const { if(header){output << "update decayers set parameters=\"";} if(create) { output << "create Herwig::ThreePionCLEOCurrent " << name() << " HwWeakCurrents.so\n"; } for(unsigned int ix=0;ix<_rhomass.size();++ix) { if(ix<2) { output << "newdef " << name() << ":RhoMasses " << ix << " " << _rhomass[ix]/MeV << "\n"; } else { output << "insert " << name() << ":RhoMasses " << ix << " " << _rhomass[ix]/MeV << "\n"; } } for(unsigned int ix=0;ix<_rhowidth.size();++ix) { if(ix<2) { output << "newdef " << name() << ":RhoWidths " << ix << " " << _rhowidth[ix]/MeV << "\n"; } else { output << "insert " << name() << ":RhoWidths " << ix << " " << _rhowidth[ix]/MeV << "\n"; } } output << "newdef " << name() << ":f_2Mass " << _f2mass/GeV << "\n"; output << "newdef " << name() << ":f_2Width " << _f2width/GeV << "\n"; output << "newdef " << name() << ":f_0Mass " << _f0mass/GeV << "\n"; output << "newdef " << name() << ":f_0Width " << _f0width/GeV << "\n"; output << "newdef " << name() << ":sigmaMass " << _sigmamass/GeV << "\n"; output << "newdef " << name() << ":sigmaWidth " << _sigmawidth/GeV << "\n"; output << "newdef " << name() << ":a1Mass " << _a1mass/GeV << "\n"; output << "newdef " << name() << ":a1Width " <<_a1width /GeV << "\n"; output << "newdef " << name() << ":KaonMass " << _mK/GeV << "\n"; output << "newdef " << name() << ":KStarMass " << _mKstar/GeV << "\n"; output << "newdef " << name() << ":KaonCoupling " << _gammk << "\n"; output << "newdef " << name() << ":Fpi " << _fpi/MeV << "\n"; output << "newdef " << name() << ":a1WidthOption " << _a1opt << "\n"; for(unsigned int ix=0;ix<_rhomagP.size();++ix) { if(ix<2) { output << "newdef " << name() << ":RhoPWaveMagnitude " << ix << " " << _rhomagP[ix] << "\n"; } else { output << "insert " << name() << ":RhoPWaveMagnitude " << ix << " " << _rhomagP[ix] << "\n"; } } for(unsigned int ix=0;ix<_rhophaseP.size();++ix) { if(ix<2) { output << "newdef " << name() << ":RhoPWavePhase " << ix << " " << _rhophaseP[ix] << "\n"; } else { output << "insert " << name() << ":RhoPWavePhase " << ix << " " << _rhophaseP[ix] << "\n"; } } for(unsigned int ix=0;ix<_rhomagD.size();++ix) { if(ix<2) { output << "newdef " << name() << ":RhoDWaveMagnitude " << ix << " " << _rhomagD[ix]*MeV2 << "\n"; } else { output << "insert " << name() << ":RhoDWaveMagnitude " << ix << " " << _rhomagD[ix]*MeV2 << "\n"; } } for(unsigned int ix=0;ix<_rhophaseD.size();++ix) { if(ix<2) { output << "newdef " << name() << ":RhoDWavePhase " << ix << " " << _rhophaseD[ix] << "\n"; } else { output << "insert " << name() << ":RhoDWavePhase " << ix << " " << _rhophaseD[ix] << "\n"; } } output << "newdef " << name() << ":f0Phase " << _f0phase << "\n"; output << "newdef " << name() << ":f2Phase " <<_f2phase << "\n"; output << "newdef " << name() << ":sigmaPhase " <<_sigmaphase << "\n"; output << "newdef " << name() << ":f0Magnitude " << _f0mag << "\n"; output << "newdef " << name() << ":f2Magnitude " << _f2mag*GeV2 << "\n"; output << "newdef " << name() << ":sigmaMagnitude " <<_sigmamag << "\n"; output << "newdef " << name() << ":Initializea1 " <<_initializea1 << "\n"; for(unsigned int ix=0;ix<_a1runwidth.size();++ix) { if(ix<200) { output << "newdef " << name() << ":a1RunningWidth " << ix << " " << _a1runwidth[ix]/MeV << "\n"; } else { output << "insert " << name() << ":a1RunningWidth " << ix << " " << _a1runwidth[ix]/MeV << "\n"; } } for(unsigned int ix=0;ix<_a1runq2.size();++ix) { if(ix<200) { output << "newdef " << name() << ":a1RunningQ2 " << ix << " " << _a1runq2[ix]/MeV2 << "\n"; } else { output << "insert " << name() << ":a1RunningQ2 " << ix << " " << _a1runq2[ix]/MeV2 << "\n"; } } WeakCurrent::dataBaseOutput(output,false,false); if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } void ThreePionCLEOCurrent::doinitrun() { // set up the running a_1 width inita1Width(0); WeakCurrent::doinitrun(); } void ThreePionCLEOCurrent::doupdate() { WeakCurrent::doupdate(); // update running width if needed if ( !touched() ) return; if(_maxmass!=_maxcalc) inita1Width(-1); } Energy ThreePionCLEOCurrent::a1width(Energy2 q2) const { Energy output; if(_a1opt) output=(*_a1runinter)(q2); else { double gam(0.); if(q2<0.1753*GeV2) { gam =0.; } else if(q2<0.823*GeV2) { double p=q2/GeV2-0.1753; gam = 5.80900*p*sqr(p)*(1.-3.00980*p+4.57920*sqr(p)); } else { double p=q2/GeV2; gam = -13.91400+27.67900*p-13.39300*sqr(p) +3.19240*sqr(p)*p-0.10487*sqr(sqr(p)); } if(q2<0.1676*GeV2) { gam+=0.; } else if(q2<0.823*GeV2) { double p=q2/GeV2-0.1676; gam+= 6.28450*p*sqr(p)*(1.-2.95950*p+4.33550*sqr(p)); } else { double p=q2/GeV2; gam+= -15.41100+32.08800*p-17.66600*sqr(p) +4.93550*sqr(p)*p-0.37498*sqr(sqr(p)); } Energy mkst=0.894*GeV,mk=0.496*GeV; Energy2 mk1sq=sqr(mkst+mk), mk2sq=sqr(mkst-mk); double c3pi=sqr(0.2384),ckst=sqr(4.7621)*c3pi; gam*=c3pi; if(q2>mk1sq) gam+=0.5*ckst*sqrt((q2-mk1sq)*(q2-mk2sq))/q2; gam = gam*_a1width*_a1mass/GeV2/1.331/0.814/1.0252088; output = gam*GeV2/sqrt(q2); } return output; } double ThreePionCLEOCurrent::threeBodyMatrixElement(const int iopt, const Energy2 q2, const Energy2 s3, const Energy2 s2, const Energy2 s1, const Energy, const Energy, const Energy) const { Energy p1[5],p2[5],p3[5]; Energy2 p1sq, p2sq, p3sq; Energy q=sqrt(q2); Energy2 mpi2c=_mpic*_mpic; Energy2 mpi20=_mpi0*_mpi0; // construct the momenta for the 2 neutral 1 charged mode Complex F1,F2,F3; if(iopt==0) { // construct the momenta of the decay products p1[0] = 0.5*(q2+mpi20-s1)/q; p1sq=p1[0]*p1[0]; p1[4]=sqrt(p1sq-mpi20); p2[0] = 0.5*(q2+mpi20-s2)/q; p2sq=p2[0]*p2[0]; p2[4]=sqrt(p2sq-mpi20); p3[0] = 0.5*(q2+mpi2c-s3)/q; p3sq=p3[0]*p3[0]; p3[4]=sqrt(p3sq-mpi2c); // take momentum of 1 parallel to z axis p1[1]=ZERO;p1[2]=ZERO;p1[3]=p1[4]; // construct 2 double cos2 = 0.5*(p1sq+p2sq-p3sq-2.*mpi20+mpi2c)/p1[4]/p2[4]; p2[1] = p2[4]*sqrt(1.-cos2*cos2); p2[2]=ZERO; p2[3]=-p2[4]*cos2; // construct 3 double cos3 = 0.5*(p1sq-p2sq+p3sq-mpi2c)/p1[4]/p3[4]; p3[1] =-p3[4]*sqrt(1.-cos3*cos3); p3[2]=ZERO; p3[3]=-p3[4]*cos3; // calculate the form factors CLEOFormFactor(1,-1,q2,s1,s2,s3,F1,F2,F3); } // construct the momenta for the 3 charged mode else { // construct the momenta of the decay products p1[0] = 0.5*(q2+mpi2c-s1)/q; p1sq=p1[0]*p1[0]; p1[4]=sqrt(p1sq-mpi2c); p2[0] = 0.5*(q2+mpi2c-s2)/q; p2sq=p2[0]*p2[0]; p2[4]=sqrt(p2sq-mpi2c); p3[0] = 0.5*(q2+mpi2c-s3)/q; p3sq=p3[0]*p3[0]; p3[4]=sqrt(p3sq-mpi2c); // take momentum of 1 parallel to z axis p1[1]=ZERO;p1[2]=ZERO;p1[3]=p1[4]; // construct 2 double cos2 = 0.5*(p1sq+p2sq-p3sq-mpi2c)/p1[4]/p2[4]; p2[1] = p2[4]*sqrt(1.-cos2*cos2); p2[2]=ZERO; p2[3]=-p2[4]*cos2; // construct 3 double cos3 = 0.5*(p1sq-p2sq+p3sq-mpi2c)/p1[4]/p3[4]; p3[1] =-p3[4]*sqrt(1.-cos3*cos3); p3[2]=ZERO; p3[3]=-p3[4]*cos3; // calculate the form factors CLEOFormFactor(0,-1,q2,s1,s2,s3,F1,F2,F3); } // construct a vector with the current complex current[4]; for(unsigned int ix=0;ix<4;++ix) current[ix] = F1*(p2[ix]-p3[ix])-F2*(p3[ix]-p1[ix])+F3*(p1[ix]-p2[ix]); complex dot1=current[0]*conj(current[0]); for(unsigned int ix=1;ix<4;++ix) dot1-=current[ix]*conj(current[ix]); complex dot2=current[0]*q; return(-dot1+dot2*conj(dot2)/q2).real() / sqr(_rhomass[0]); } // the hadronic currents vector ThreePionCLEOCurrent::current(tcPDPtr resonance, FlavourInfo flavour, const int imode, const int ichan, Energy & scale, const tPDVector & , const vector & momenta, DecayIntegrator::MEOption) const { useMe(); // check the isospin if(flavour.I!=IsoSpin::IUnknown && flavour.I!=IsoSpin::IOne) return vector(); // check I_3 if(flavour.I3!=IsoSpin::I3Unknown) { switch(flavour.I3) { case IsoSpin::I3Zero: if(imode==2||imode==5) return vector(); break; case IsoSpin::I3One: if(imode!=2&&imode!=5) return vector(); break; case IsoSpin::I3MinusOne: if(imode!=2&&imode!=5) return vector(); break; default: return vector(); } } if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return vector(); if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return vector(); if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return vector(); // calculate q2,s1,s2,s3 Lorentz5Momentum q; for(unsigned int ix=0;ix a1fact = _fact; if(!resonance) a1fact *= a1BreitWigner(q2); // current LorentzPolarizationVectorE vect = q.mass()*a1fact* ((F2-F1)*momenta[2] + (F1-F3)*momenta[1] + (F3-F2)*momenta[0]); // scalar piece Complex dot=(vect*q)/q2; vect -= dot*q; // return the answer return vector(1,vect); } bool ThreePionCLEOCurrent::accept(vector id) { if(id.size()!=3) return false; for(unsigned int ix=0;ix id) { if(id.size()!=3) return -1; int npip(0),npim(0),npi0(0); for(unsigned int ix=0;ixCC()) extpart[ix]=extpart[ix]->CC(); } } // return the answer return extpart; } diff --git a/Hadronization/SpinHadronizer.h b/Hadronization/SpinHadronizer.h --- a/Hadronization/SpinHadronizer.h +++ b/Hadronization/SpinHadronizer.h @@ -1,191 +1,191 @@ // -*- C++ -*- #ifndef Herwig_SpinHadronizer_H #define Herwig_SpinHadronizer_H // // This is the declaration of the SpinHadronizer class. // #include "ThePEG/Handlers/StepHandler.h" namespace Herwig { using namespace ThePEG; /** * The SpinHadronizer class is designed to be used as a post-hadronization handler to * give a simple model of spin transfer between the perturbative and non-perturbative * stages. * * @see \ref SpinHadronizerInterfaces "The interfaces" * defined for SpinHadronizer. */ class SpinHadronizer: public StepHandler { public: /** * The default constructor. */ SpinHadronizer() : omegaHalf_(2./3.), omegaThreeHalf_(0.2), minFlav_(3), maxFlav_(5), debug_(false), qPol_(6,make_pair(0.,0.)) {} public: /** @name Virtual functions required by the StepHandler class. */ //@{ /** * The main function called by the EventHandler class to * perform a step. Given the current state of an Event, this function * performs the event generation step and includes the result in a new * Step object int the Event record. * @param eh the EventHandler in charge of the Event generation. * @param tagged if not empty these are the only particles which should * be considered by the StepHandler. * @param hint a Hint object with possible information from previously * performed steps. * @throws Veto if the StepHandler requires the current step to be discarded. * @throws Stop if the generation of the current Event should be stopped * after this call. * @throws Exception if something goes wrong. */ virtual void handle(EventHandler & eh, const tPVector & tagged, const Hint & hint); //@} 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: /** * Functions to calculate the spins */ //@{ /** * Calculate the spin of a baryon */ void baryonSpin(tPPtr baryon); /** * Calculate the spin of a meson */ void mesonSpin(tPPtr meson); //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); //@} protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - SpinHadronizer & operator=(const SpinHadronizer &); + SpinHadronizer & operator=(const SpinHadronizer &) = delete; private: /** * Parameters */ //@{ /** * Falk-Peskin \f$\omega_\frac12\f$ parameter */ double omegaHalf_; /** * Falk-Peskin \f$\omega_\frac32\f$ parameter */ double omegaThreeHalf_; /** * Minimum quark flavour */ unsigned int minFlav_; /** * Maximum quark flavour */ unsigned int maxFlav_; /** * Print out debugging info */ bool debug_; /** * Polarization of the quarks */ vector > qPol_; //@} }; } #endif /* Herwig_SpinHadronizer_H */ diff --git a/MatrixElement/Matchbox/Utility/DensityOperator.h b/MatrixElement/Matchbox/Utility/DensityOperator.h --- a/MatrixElement/Matchbox/Utility/DensityOperator.h +++ b/MatrixElement/Matchbox/Utility/DensityOperator.h @@ -1,248 +1,248 @@ // -*- C++ -*- // // DensityOperator.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef Herwig_DensityOperator_H #define Herwig_DensityOperator_H // // This is the declaration of the DensityOperator class. // #include "ThePEG/Handlers/HandlerBase.h" #include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h" #include #include #include namespace Herwig { using namespace ThePEG; typedef boost::numeric::ublas::vector CVector; /** * Here is the documentation of the DensityOperator class. * * @see \ref DensityOperatorInterfaces "The interfaces" * defined for DensityOperator. */ class DensityOperator: public HandlerBase { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ DensityOperator(); /** * The destructor. */ virtual ~DensityOperator(); //@} public: /** * Clears theDensityOperatorMap. */ void clear(); /** * Prepare for the given sub process. */ void prepare(const cPDVector&); /** * Fill the density operator for the given hard subprocess, summing over all * helicity configurations. */ void fill(const Ptr::ptr, const cPDVector&, const vector& momenta); /** * Evolve the density operator, by * M_{n+1} = -\sum_{i,k}{-4*pi*alpha_s/Ti2*V_{ij,k} T_{i,n}M_nT_{k,n}^\dag}, * see arXiv:1206.0180 eq. (5), note that the pi*pj factor is assumed to be * included in V_{ij,k}. */ void evolve(const map,Complex>& Vijk, const cPDVector& before, const cPDVector& after, const map,map >& emissionsMap, const bool splitAGluon, const bool initialGluonSplitting); /** * Calculate the colour matrix element correction. * -(1+delta(i,gluon))/Ti^2 Tr(Sn+1 Ti Mn Tk^dagger)/Tr(Sn Mn) * where the bracket in front compensates for the gluon symmetry factor, * Ti^2 is C_f or C_a, Sn+1 is the matrix of scalar products, and * Ti is the radiation matrix. * The first arg contains (emitter index, spectator index, emission pid) * */ double colourMatrixElementCorrection(const std::tuple& ikemission, const cPDVector& particles); /** * Checking colour conservation for the colour matrix element corrections. */ void colourConservation(const cPDVector& particles); /** * Get the colour basis. */ Ptr::tptr colourBasis() { return theColourBasis; } /** * Get the colour basis. */ const Ptr::tptr colourBasis() const { return theColourBasis; } /** * Set the colour basis. */ void colourBasis(Ptr::ptr ptr) { theColourBasis = ptr; } /** * Get the correlator map. */ const map,pair >,double>& correlatorMap() const { return theCorrelatorMap; } /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * Number of colours used in colourNorm. */ double Nc; /** * QCD vertex normalization. */ double TR; /** * Normalization of colour charges \mathbf{T}_{ij}^2. */ double colourNorm(const cPDPtr particle); /** * Fast evaluation of Tij*Mn, where a Tij is the matrix from ColourBasis::charge, * which is a sparse matrix, and Mn is the density operator, a dense matrix. * */ matrix prodSparseDense(const compressed_matrix&, const matrix&); /** * Fast evaluation of TijMn*Tkdagger, where a TijMn is the result from the method * prodSparseDense, a dense matrix, and Tkdagger is the transponse conjugate of * the matrix from ColourBasis::charge, a sparse matrix. * */ matrix prodDenseSparse(const matrix&, const compressed_matrix&); /** * Boosts a vector of momenta to the rest frame of the initial pair * of particles (the first 2 elements of the argument vector). Returns * the boosted vectors */ vector boostToRestFrame(const vector& momenta); /** * Boosts a vector of momenta to the rest frame of the initial pair */ bool compareMomentum(const Lorentz5Momentum& p, const Lorentz5Momentum& q); /** * Mapping of colour structures to density operator matrices. * */ map,matrix > theDensityOperatorMap; /** * Mapping of colour structures and legs to colour correlators. */ map,pair >,double> theCorrelatorMap; /** * A map from the hard subprocess particles to a map of amplitude colour * basis order to the normal ordered colour basis. */ map > theColourBasisToColourBasisMap; /** * Colour basis used. */ Ptr::ptr theColourBasis; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - DensityOperator & operator=(const DensityOperator &); + DensityOperator & operator=(const DensityOperator &) = delete; }; } #endif /* Herwig_DensityOperator_H */ diff --git a/Models/DarkMatter/DMDMMediatorVertex.h b/Models/DarkMatter/DMDMMediatorVertex.h --- a/Models/DarkMatter/DMDMMediatorVertex.h +++ b/Models/DarkMatter/DMDMMediatorVertex.h @@ -1,113 +1,113 @@ // -*- C++ -*- #ifndef Herwig_DMDMMediatorVertex_H #define Herwig_DMDMMediatorVertex_H // // This is the declaration of the DMDMMediatorVertex class. // #include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h" namespace Herwig { using namespace ThePEG; /** * The DMDMMediatorVertex class implements the coupling of the dark matter to the mediator * * @see \ref DMDMMediatorVertexInterfaces "The interfaces" * defined for DMDMMediatorVertex. */ class DMDMMediatorVertex: public FFVVertex { public: /** * The default constructor. */ DMDMMediatorVertex(); /** * Calculate the couplings. * @param q2 The scale \f$q^2\f$ for the coupling at the vertex. * @param part1 The ParticleData pointer for the first particle. * @param part2 The ParticleData pointer for the second particle. * @param part3 The ParticleData pointer for the third particle. */ virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - DMDMMediatorVertex & operator=(const DMDMMediatorVertex &); + DMDMMediatorVertex & operator=(const DMDMMediatorVertex &) = delete; private: /** * DM coupling to the dark mediator */ double cDMmed_; }; } #endif /* Herwig_DMDMMediatorVertex_H */ diff --git a/Models/DarkMatter/MEDM2Mesons.cc b/Models/DarkMatter/MEDM2Mesons.cc --- a/Models/DarkMatter/MEDM2Mesons.cc +++ b/Models/DarkMatter/MEDM2Mesons.cc @@ -1,266 +1,266 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEDM2Mesons class. // #include "MEDM2Mesons.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "Herwig/Models/General/BSMModel.h" using namespace Herwig; typedef LorentzVector > LorentzPolarizationVectorInvE; MEDM2Mesons::MEDM2Mesons() : cDMmed_(0.), cSMmed_({0.0,0.0,0.0}) {} Energy2 MEDM2Mesons::scale() const { return sHat(); } unsigned int MEDM2Mesons::orderInAlphaS() const { return 0; } unsigned int MEDM2Mesons::orderInAlphaEW() const { return 0; } IBPtr MEDM2Mesons::clone() const { return new_ptr(*this); } IBPtr MEDM2Mesons::fullclone() const { return new_ptr(*this); } void MEDM2Mesons::doinit() { // make sure the current got initialised current_->init(); // max energy Energy Emax = generator()->maximumCMEnergy(); // loop over the modes int nmode=0; for(unsigned int imode=0;imodenumberOfModes();++imode) { // get the external particles for this mode int iq(0),ia(0); tPDVector out = current_->particles(0,imode,iq,ia); current_->decayModeInfo(imode,iq,ia); if(iq==2&&ia==-2) continue; PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(incomingA_,out,1., incomingB_,Emax)); PhaseSpaceChannel channel(mode); if(!current_->createMode(0,tcPDPtr(), FlavourInfo(), imode,mode,0,-1,channel,Emax)) continue; modeMap_[imode] = nmode; addMode(mode); ++nmode; } // cast the model Ptr::ptr model = dynamic_ptr_cast::ptr>(generator()->standardModel()); bool foundDM(false),foundU(false),foundD(false),foundS(false); // find the vertices we need and extract the couplings for(unsigned int ix = 0; ix < model->numberOfVertices(); ++ix ) { VertexBasePtr vertex = model->vertex(ix); if(vertex->getNpoint()!=3) continue; for(unsigned int iloc = 0;iloc < 3; ++iloc) { vector ext = vertex->search(iloc, Mediator_->id()); if(ext.empty()) continue; for(unsigned int ioff=0;ioffid() && ext[ioff+1]==incomingB_->id()) || (ext[ioff]==incomingB_->id() && ext[ioff+1]==incomingA_->id())) { foundDM = true; vertex->setCoupling(sqr(Emax),incomingA_,incomingB_,Mediator_); cDMmed_ = vertex->norm(); } else if(abs(ext[ioff])==1 && abs(ext[ioff+1])==1 && ext[ioff]==-ext[ioff+1]) { foundD = true; vertex->setCoupling(sqr(Emax),getParticleData(1),getParticleData(-1),Mediator_); cSMmed_[0] = vertex->norm(); } else if(abs(ext[ioff])==2 && abs(ext[ioff+1])==2 && ext[ioff]==-ext[ioff+1]) { foundU = true; vertex->setCoupling(sqr(Emax),getParticleData(2),getParticleData(-2),Mediator_); cSMmed_[1] = vertex->norm(); } else if(abs(ext[ioff])==3 && abs(ext[ioff+1])==3 && ext[ioff]==-ext[ioff+1]) { foundS = true; vertex->setCoupling(sqr(Emax),getParticleData(3),getParticleData(-3),Mediator_); cSMmed_[2] = vertex->norm(); } } } } if(!foundDM) { throw InitException() << "Cannot find DM coupling in MEDM2Mesons::doinit()"; } if(!foundD) { throw InitException() << "Cannot find down quark coupling in MEDM2Mesons::doinit()"; } if(!foundU) { throw InitException() << "Cannot find up quark coupling in MEDM2Mesons::doinit()"; } if(!foundS) { throw InitException() << "Cannot find strange quark coupling in MEDM2Mesons::doinit()"; } MEMultiChannel::doinit(); } void MEDM2Mesons::persistentOutput(PersistentOStream & os) const { os << current_ << modeMap_ << incomingA_ << incomingB_ << Mediator_ << cDMmed_ << cSMmed_; } void MEDM2Mesons::persistentInput(PersistentIStream & is, int) { is >> current_ >> modeMap_ >> incomingA_ >> incomingB_ >> Mediator_ >> cDMmed_ >> cSMmed_; } //The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigMEDM2Mesons("Herwig::MEDM2Mesons", "Herwig.so"); void MEDM2Mesons::Init() { static ClassDocumentation documentation ("The MEDM2Mesons class simulates the annhilation of" " DM particles to mesons at low energy"); static Reference interfaceWeakCurrent ("WeakCurrent", "The reference for the decay current to be used.", &MEDM2Mesons::current_, false, false, true, false, false); static Reference interfaceIncomingA ("IncomingA", "First incoming particle", &MEDM2Mesons::incomingA_, false, false, true, false, false); static Reference interfaceIncomingB ("IncomingB", "Second incoming particle", &MEDM2Mesons::incomingB_, false, false, true, false, false); static Reference interfaceMediator_ ("Mediator", "DM mediator", &MEDM2Mesons::Mediator_, false, false, true, false, false); } double MEDM2Mesons::me2(const int ichan) const { // compute the incoming current LorentzPolarizationVectorInvE lepton[2][2]; if(incomingA_->iSpin()==PDT::Spin1Half && incomingB_->iSpin()==PDT::Spin1Half) { SpinorWaveFunction em_in( meMomenta()[0],mePartonData()[0],incoming); SpinorBarWaveFunction ep_in( meMomenta()[1],mePartonData()[1],incoming); vector f1; vector a1; for(unsigned int ix=0;ix<2;++ix) { em_in.reset(ix); f1.push_back(em_in); ep_in.reset(ix); a1.push_back(ep_in); } // this should be coupling of DM to mediator/ mediator propagator complex mmed = Mediator_->mass(); complex mmed2 = sqr(mmed); complex mwid = Mediator_->width(); complex prop = sHat()-mmed2+Complex(0.,1.)*mmed*mwid; complex pre = cDMmed_/prop; for(unsigned ix=0;ix<2;++ix) { for(unsigned iy=0;iy<2;++iy) { lepton[ix][iy]= pre*f1[ix].dimensionedWave().vectorCurrent(a1[iy].dimensionedWave()); } } } // TODO think about other spins for the DM else assert(false); // work out the mapping for the hadron vector int nOut = int(meMomenta().size())-2; vector constants(nOut+1); vector iSpin(nOut); vector hadrons(nOut); int itemp(1); int ix(nOut); do { --ix; iSpin[ix] = mePartonData()[ix+2]->iSpin(); itemp *= iSpin[ix]; constants[ix] = itemp; hadrons[ix] = mePartonData()[ix+2]->id(); } while(ix>0); constants[nOut] = 1; // calculate the matrix element me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,iSpin)); // calculate the hadron current unsigned int imode = current_->decayMode(hadrons); Energy q = sqrt(sHat()); vector momenta(meMomenta() .begin()+2, meMomenta().end()); tPDVector out = mode(modeMap_.at(imode))->outgoing(); if(ichan<0) iMode(modeMap_.at(imode)); // get the hadronic currents for the I=1 and I=0 components vector hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero), imode,ichan,q,out,momenta,DecayIntegrator::Calculate)); vector hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero), imode,ichan,q,out,momenta,DecayIntegrator::Calculate)); vector hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::ssbar), imode,ichan,q,out,momenta,DecayIntegrator::Calculate)); // compute the matrix element vector ihel(meMomenta().size()); double output(0.); unsigned int hI0_size = hadronI0.size(); unsigned int hI1_size = hadronI1.size(); unsigned int hss_size = hadronssbar.size(); unsigned int maxsize = max(max(hadronI0.size(),hadronI1.size()),hss_size); for(unsigned int hhel=0;hhel0;--ix) { ihel[ix+1]=(hhel%constants[ix-1])/constants[ix]; } // loop over the helicities of the incoming particles for(ihel[1]=0;ihel[1]<2;++ihel[1]){ for(ihel[0]=0;ihel[0]<2;++ihel[0]) { Complex amp = 0.; // work on coefficients for the I1 and I0 bits if(hI0_size != 0 ) - amp += (cSMmed_[0]+cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI0[hhel])); + amp += Complex((cSMmed_[0]+cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI0[hhel]))); if(hI1_size !=0) - amp += (cSMmed_[0]-cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI1[hhel])); + amp += Complex((cSMmed_[0]-cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI1[hhel]))); if(hss_size !=0) - amp += cSMmed_[2]* (lepton[ihel[0]][ihel[1]].dot(hadronssbar[hhel])); + amp += Complex(cSMmed_[2]* (lepton[ihel[0]][ihel[1]].dot(hadronssbar[hhel]))); me_(ihel)= amp; output += std::norm(amp); } } } // symmetry factors map ncount; double symmetry(1.); for(tPDPtr o : out) ncount[o->id()]+=1; for(map::const_iterator it=ncount.begin();it!=ncount.end();++it) { symmetry *= it->second; } // prefactors output *= 0.25*sqr(pow(sqrt(sHat())/q,int(momenta.size()-2))); return output/symmetry; } void MEDM2Mesons::constructVertex(tSubProPtr) { } diff --git a/Models/DarkMatter/MEDM2Mesons.h b/Models/DarkMatter/MEDM2Mesons.h --- a/Models/DarkMatter/MEDM2Mesons.h +++ b/Models/DarkMatter/MEDM2Mesons.h @@ -1,177 +1,177 @@ // -*- C++ -*- #ifndef Herwig_MEDM2Mesons_H #define Herwig_MEDM2Mesons_H // // This is the declaration of the MEDM2Mesons class. // #include "Herwig/MatrixElement/MEMultiChannel.h" #include "Herwig/Decay/WeakCurrents/WeakCurrent.h" #include "Herwig/MatrixElement/ProductionMatrixElement.h" namespace Herwig { using namespace ThePEG; /** * The MEDM2Mesons class implements dark matter annhilation via a vector current to * mesons at low energies * * @see \ref MEDM2MesonsInterfaces "The interfaces" * defined for MEDM2Mesons. */ class MEDM2Mesons: public MEMultiChannel { public: /** * The default constructor. */ MEDM2Mesons(); public: /** @name Virtual functions required by the MEBase class. */ //@{ /** * Return the order in \f$\alpha_S\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaS() const; /** * Return the order in \f$\alpha_{EW}\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaEW() const; /** * Return the scale associated with the last set phase space point. */ virtual Energy2 scale() const; //@} /** * Construct the vertex of spin correlations. */ virtual void constructVertex(tSubProPtr); 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: /** * Return the matrix element squared for a given mode and phase-space channel. * @param ichan The channel we are calculating the matrix element for. */ virtual double me2(const int ichan) const; protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - MEDM2Mesons & operator=(const MEDM2Mesons &); + MEDM2Mesons & operator=(const MEDM2Mesons &) = delete; private : /** * Hadronic current etc */ //@{ /** * the hadronic current */ WeakCurrentPtr current_; /** * The matrix element */ mutable ProductionMatrixElement me_; /** * Map for the modes */ map modeMap_; //@} /** * DM */ //@{ /** * Incoming Particles */ PDPtr incomingA_, incomingB_; //@} /** * DM coupling to the dark mediator */ Complex cDMmed_; /** * SM couplings to the dark mediator */ vector cSMmed_; /** * DM vector mediator */ PDPtr Mediator_; }; } #endif /* Herwig_MEDM2Mesons_H */ diff --git a/PDF/HwRemDecayer.h b/PDF/HwRemDecayer.h --- a/PDF/HwRemDecayer.h +++ b/PDF/HwRemDecayer.h @@ -1,747 +1,747 @@ // -*- C++ -*- // // HwRemDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_HwRemDecayer_H #define HERWIG_HwRemDecayer_H // // This is the declaration of the HwRemDecayer class. // #include "ThePEG/PDT/RemnantDecayer.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/EventRecord/SubProcess.h" #include "ThePEG/PDF/BeamParticleData.h" #include "Herwig/Shower/ShowerAlpha.h" #include "Herwig/PDT/StandardMatchers.h" #include "ThePEG/PDT/StandardMatchers.h" #include "HwRemDecayer.fh" namespace Herwig { using namespace ThePEG; /** * The HwRemDecayer class is responsible for the decay of the remnants. Additional * secondary scatters have to be evolved backwards to a gluon, the * first/hard interaction has to be evolved back to a valence quark. * This is all generated inside this class, * which main methods are then called by the ShowerHandler. * * A simple forced splitting algorithm is used. * This takes the Remnant object produced from the PDF and backward * evolution (hadron - parton) and produce partons with the remaining * flavours and with the correct colour connections. * * The algorithim operates by starting with the parton which enters the hard process. * If this is from the sea there is a forced branching to produce the antiparticle * from a gluon branching. If the parton entering the hard process was a gluon, or * a gluon was produced from the first step of the algorithm, there is then a further * branching back to a valence parton. After these partons have been produced a quark or * diquark is produced to give the remaining valence content of the incoming hadron. * * The forced branching are generated using a scale between QSpac and EmissionRange times * the minimum scale. The energy fractions are then distributed using * \f[\frac{\alpha_S}{2\pi}\frac{P(z)}{z}f(x/z,\tilde{q})\f] * with the massless splitting functions. * * \author Manuel B\"ahr * * @see \ref HwRemDecayerInterfaces "The interfaces" * defined for HwRemDecayer. */ class HwRemDecayer: public RemnantDecayer { public: /** Typedef to store information about colour partners */ typedef vector > PartnerMap; public: /** * The default constructor. */ HwRemDecayer() : allowTop_(false), multiPeriph_(true), quarkPair_(false), ptmin_(-1.*GeV), beta_(ZERO), maxtrySoft_(10), colourDisrupt_(1.0), ladderbFactor_(0.0), ladderPower_(-0.08), ladderNorm_(1.0), ladderMult_(1.0), gaussWidth_(0.1), valOfN_(0), initTotRap_(0), _kinCutoff(0.75*GeV), _forcedSplitScale(2.5*GeV), _range(1.1), _zbin(0.05),_ybin(0.), _nbinmax(100), DISRemnantOpt_(0), PtDistribution_(0), pomeronStructure_(0), mg_(ZERO) {} /** @name Virtual functions required by the Decayer class. */ //@{ /** * Check if this decayer can perfom the decay specified by the * given decay mode. * @return true if this decayer can handle the given mode, otherwise false. */ virtual bool accept(const DecayMode &) const { return true; } /** * Return true if this decayer can handle the extraction of the \a * extracted parton from the given \a particle. */ virtual bool canHandle(tcPDPtr particle, tcPDPtr parton) const; /** * Return true if this decayed can extract more than one parton from * a particle. */ virtual bool multiCapable() const { return true; } /** * Perform a decay for a given DecayMode and a given Particle instance. * @param dm the DecayMode describing the decay. * @param p the Particle instance to be decayed. * @param step the step we are working on. * @return a ParticleVector containing the decay products. */ virtual ParticleVector decay(const DecayMode & dm, const Particle & p, Step & step) const; //@} public: /** * struct that is used to catch exceptions which are thrown * due to energy conservation issues of additional soft scatters */ struct ExtraSoftScatterVeto {}; /** @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(); /** * Do several checks and initialization, for remnantdecay inside ShowerHandler. */ void initialize(pair rems, tPPair beam, Step & step, Energy forcedSplitScale); /** * Initialize the soft scattering machinery. * @param ptmin = the pt cutoff used in the UE model * @param beta = slope of the soft pt-spectrum */ void initSoftInteractions(Energy ptmin, InvEnergy2 beta); /** * Perform the acual forced splitting. * @param partons is a pair of ThePEG::Particle pointers which store the final * partons on which the shower ends. * @param pdfs are pointers to the pdf objects for both beams * @param first is a flage wether or not this is the first or a secondary interation */ void doSplit(pair partons, pair pdfs, bool first); /** * Perform the final creation of the diquarks. Set the remnant masses and do * all colour connections. * @param colourDisrupt = variable to control how many "hard" scatters * are colour isolated * @param softInt = parameter for the number of soft scatters */ void finalize(double colourDisrupt=0.0, unsigned int softInt=0); /** * Find the children */ void findChildren(tPPtr,vector &) const; protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** 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 {return new_ptr(*this);} //@} 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() { Interfaced::doinit(); _ybin=0.25/_zbin; mg_ = getParticleData(ParticleID::g)->constituentMass(); } //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - HwRemDecayer & operator=(const HwRemDecayer &); + HwRemDecayer & operator=(const HwRemDecayer &) = delete; public: /** * Simple struct to store info about baryon quark and di-quark * constituents. */ struct HadronContent { /** * manually extract the valence flavour \a id. */ inline void extract(int id) { for(unsigned int i=0; iid() == ParticleID::gamma || (hadron->id() == ParticleID::pomeron && pomeronStructure==1) || hadron->id() == ParticleID::reggeon) { flav[0] = id; flav[1] = -id; extracted = 0; flav.resize(2); } else if (hadron->id() == ParticleID::pomeron && pomeronStructure==0) { extracted = 0; } else { extracted = i; } break; } } } /** * Return a proper particle ID assuming that \a id has been removed * from the hadron. */ long RemID() const; /** * Method to determine whether \a parton is a quark from the sea. * @return TRUE if \a parton is neither a valence quark nor a gluon. */ bool isSeaQuark(tcPPtr parton) const { return ((parton->id() != ParticleID::g) && ( !isValenceQuark(parton) ) ); } /** * Method to determine whether \a parton is a valence quark. */ bool isValenceQuark(tcPPtr parton) const { return isValenceQuark(parton->id()); } /** * Method to determine whether \a parton is a quark from the sea. * @return TRUE if \a parton is neither a valence quark nor a gluon. */ bool isSeaQuarkData(tcPDPtr partonData) const { return ((partonData->id() != ParticleID::g) && ( !isValenceQuarkData(partonData) ) ); } /** * Method to determine whether \a parton is a valence quark. */ bool isValenceQuarkData(tcPDPtr partonData) const { int id(sign*partonData->id()); return find(flav.begin(),flav.end(),id) != flav.end(); } /** * Method to determine whether \a parton is a valence quark. */ bool isValenceQuark(int id) const { return find(flav.begin(),flav.end(),sign*id) != flav.end(); } /** The valence flavours of the corresponding baryon. */ vector flav; /** The array index of the extracted particle. */ int extracted; /** -1 if the particle is an anti-particle. +1 otherwise. */ int sign; /** The ParticleData objects of the hadron */ tcPDPtr hadron; /** Pomeron treatment */ unsigned int pomeronStructure; }; /** * Return the hadron content objects for the incoming particles. */ const pair& content() const { return theContent; } /** * Return a HadronContent struct from a PPtr to a hadron. */ HadronContent getHadronContent(tcPPtr hadron) const; /** * Set the hadron contents. */ void setHadronContent(tPPair beam) { theContent.first = getHadronContent(beam.first); theContent.second = getHadronContent(beam.second); } private: /** * Do the forced Splitting of the Remnant with respect to the * extracted parton \a parton. * @param parton = PPtr to the parton going into the subprocess. * @param content = HadronContent struct to keep track of flavours. * @param rem = Pointer to the ThePEG::RemnantParticle. * @param used = Momentum vector to keep track of remaining momenta. * @param partners = Vector of pairs filled with tPPtr to the particles * which should be colour connected. * @param pdf pointer to the PDF Object which is used for this particle * @param first = Flag for the first interaction. */ void split(tPPtr parton, HadronContent & content, tRemPPtr rem, Lorentz5Momentum & used, PartnerMap & partners, tcPDFPtr pdf, bool first); /** * Merge the colour lines of two particles * @param p1 = Pointer to particle 1 * @param p2 = Pointer to particle 2 * @param anti = flag to indicate, if (anti)colour was extracted as first parton. */ void mergeColour(tPPtr p1, tPPtr p2, bool anti) const; /** * Set the colour connections. * @param partners = Object that holds the information which particles to connect. * @param anti = flag to indicate, if (anti)colour was extracted as first parton. * @param disrupt parameter for disruption of the colour structure */ void fixColours(PartnerMap partners, bool anti, double disrupt) const; /** * Set the momenta of the Remnants properly and boost the decay particles. */ void setRemMasses() const; /** * This creates a parton from the remaining flavours of the hadron. The * last parton used was a valance parton, so only 2 (or 1, if meson) flavours * remain to be used. */ PPtr finalSplit(const tRemPPtr rem, long remID, Lorentz5Momentum usedMomentum) const { // Create the remnant and set its momentum, also reset all of the decay // products from the hadron PPtr remnant = new_ptr(Particle(getParticleData(remID))); Lorentz5Momentum prem(rem->momentum()-usedMomentum); prem.setMass(getParticleData(remID)->constituentMass()); prem.rescaleEnergy(); remnant->set5Momentum(prem); // Add the remnant to the step, but don't do colour connections thestep->addDecayProduct(rem,remnant,false); return remnant; } /** * This takes the particle and find a splitting for np -> p + child and * creates the correct kinematics and connects for such a split. This * Splitting has an upper bound on qtilde given by the energy argument * @param rem The Remnant * @param child The PDG code for the outgoing particle * @param oldQ The maximum scale for the evolution * @param oldx The fraction of the hadron's momentum carried by the last parton * @param pf The momentum of the last parton at input and after branching at output * @param p The total emitted momentum * @param content The content of the hadron */ PPtr forceSplit(const tRemPPtr rem, long child, Energy &oldQ, double &oldx, Lorentz5Momentum &pf, Lorentz5Momentum &p, HadronContent & content) const; /** * Check if a particle is a parton from a hadron or not * @param parton The parton to be tested */ bool isPartonic(tPPtr parton) const; /** @name Soft interaction methods. */ //@{ /** * Produce pt values according to dN/dp_T = N p_T exp(-beta_*p_T^2) */ Energy softPt() const; /** * Get the 2 pairs of 5Momenta for the scattering. Needs calling of * initSoftInteractions. */ void softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2, Lorentz5Momentum &g1, Lorentz5Momentum &g2) const; /** * Create N soft gluon interactions */ void doSoftInteractions(unsigned int N){ if(!multiPeriph_){ doSoftInteractions_old(N);} //outdated model for soft interactions else{ doSoftInteractions_multiPeriph(N); // Multiperipheral model } } /** * Create N soft gluon interactions (old version) */ void doSoftInteractions_old(unsigned int N); /** * Create N soft gluon interactions with multiperhpheral kinematics */ void doSoftInteractions_multiPeriph(unsigned int N); /** * Phase space generation for the ladder partons */ bool doPhaseSpaceGenerationGluons(vector &softGluons, Energy energy, unsigned int &its) const; /** * This returns the rotation matrix needed to rotate p into the z axis */ LorentzRotation rotate(const LorentzMomentum &p) const; /** * Methods to generate random distributions also all stolen form UA5Handler **/ template inline T gaussDistribution(T mean, T stdev) const{ double x = rnd(); x = sqrt(-2.*log(x)); double y; randAzm(x,x,y); return mean + stdev*x; } /** * This returns a random number with a flat distribution * [-A,A] plus gaussian tail with stdev B * TODO: Should move this to Utilities * @param A The width of the flat part * @param B The standard deviation of the gaussian tail * @return the randomly generated value */ inline double randUng(double A, double B) const{ double prun; if(A == 0.) prun = 0.; else prun = 1./(1.+B*1.2533/A); if(rnd() < prun) return 2.*(rnd()-0.5)*A; else { double temp = gaussDistribution(0.,B); if(temp < 0) return temp - abs(A); else return temp + abs(A); } } template inline void randAzm(T pt, T &px, T &py) const{ double c,s,cs; while(true) { c = 2.*rnd()-1.; s = 2.*rnd()-1.; cs = c*c+s*s; if(cs <= 1.&&cs!=0.) break; } T qt = pt/cs; px = (c*c-s*s)*qt; py = 2.*c*s*qt; } inline Energy randExt(Energy AM0,InvEnergy B) const{ double r = rnd(); // Starting value Energy am = AM0-log(r)/B; for(int i = 1; i<20; ++i) { double a = exp(-B*(am-AM0))/(1.+B*AM0); double f = (1.+B*am)*a-r; InvEnergy df = -B*B*am*a; Energy dam = -f/df; am += dam; if(am theanti; /** * variable to sum up the x values of the extracted particles */ pair theX; /**Pair of HadronContent structs to know about the quark content of the beams*/ pair theContent; /**Pair of Lorentz5Momentum to keep track of the forced splitting product momenta*/ pair theUsed; /** * Pair of PartnerMap's to store the particles, which will be colour * connected in the end. */ pair theMaps; /** * Variable to hold a pointer to the current step. The variable is used to * determine, wether decay(const DecayMode & dm, const Particle & p, Step & step) * has been called in this event or not. */ StepPtr thestep; /** * Pair of Remnant pointers. This is needed to boost * in the Remnant-Remnant CMF after all have been decayed. */ pair theRems; /** * The beam particle data for the current incoming hadron */ mutable tcPPtr theBeam; /** * the beam data */ mutable Ptr::const_pointer theBeamData; /** * The PDF for the current initial-state shower */ mutable tcPDFPtr _pdf; private: /** * Switch to control handling of top quarks in proton */ bool allowTop_; /** * Switch to control using multiperipheral kinemaics */ bool multiPeriph_; /** * True if kinematics is to be calculated for quarks */ bool quarkPair_; /** @name Soft interaction variables. */ //@{ /** * Pair of soft Remnant pointers, i.e. Diquarks. */ tPPair softRems_; /** * ptcut of the UE model */ Energy ptmin_; /** * slope of the soft pt-spectrum: dN/dp_T = N p_T exp(-beta*p_T^2) */ InvEnergy2 beta_; /** * Maximum number of attempts for the regeneration of an additional * soft scattering, before the number of scatters is reduced. */ unsigned int maxtrySoft_; /** * Variable to store the relative number of colour disrupted * connections to additional soft subprocesses. */ double colourDisrupt_; /** * Variable to store the additive factor of the multiperipheral ladder multiplicity. */ double ladderbFactor_; /** * Variable of the parameterization of the ladder multiplicity. */ double ladderPower_; /** * Variable of the parameterization of the ladder multiplicity. */ double ladderNorm_; double ladderMult_; /** * Variable to store the gaussian width of the * fluctuation of the longitudinal momentum * fraction. */ double gaussWidth_; /** * Variable to store the current total multiplicity of a ladder. */ double valOfN_; /** * Variable to store the initial total rapidity between of the remnants. */ double initTotRap_; //@} /** @name Forced splitting variables. */ //@{ /** * The kinematic cut-off */ Energy _kinCutoff; /** * The PDF freezing scale as set in ShowerHandler */ Energy _forcedSplitScale; /** * Range for emission */ double _range; /** * Size of the bins in z for the interpolation */ double _zbin; /** * Size of the bins in y for the interpolation */ double _ybin; /** * Maximum number of bins for the z interpolation */ int _nbinmax; /** * Pointer to the object calculating the QCD coupling */ ShowerAlphaPtr _alphaS; /** * Pointer to the object calculating the QED coupling */ ShowerAlphaPtr _alphaEM; /** * Option for the DIS remnant */ unsigned int DISRemnantOpt_; /** * Option for the pT generation */ unsigned int PtDistribution_; /** * Option for the treatment of the pomeron structure */ unsigned int pomeronStructure_; //@} /** * The gluon constituent mass. */ Energy mg_; }; } #endif /* HERWIG_HwRemDecayer_H */ diff --git a/PDF/MultiPartonExtractor.h b/PDF/MultiPartonExtractor.h --- a/PDF/MultiPartonExtractor.h +++ b/PDF/MultiPartonExtractor.h @@ -1,114 +1,114 @@ // -*- C++ -*- #ifndef Herwig_MultiPartonExtractor_H #define Herwig_MultiPartonExtractor_H // // This is the declaration of the MultiPartonExtractor class. // #include "ThePEG/PDF/PartonExtractor.h" #include namespace Herwig { using namespace ThePEG; /** * The MultiPartonExtractor class inherits from the PartonExtractor of ThePEG * but allows more control over the PDFs used in the case that there are multiple * stages of parton extraction * * @see \ref MultiPartonExtractorInterfaces "The interfaces" * defined for MultiPartonExtractor. */ class MultiPartonExtractor: public PartonExtractor { public: /** * The default constructor. */ MultiPartonExtractor() {}; /** * Return a vector of possible pairs of parton bins which can be * produced within a given maximum total particle-particle * invariant mass squared, \a maxEnergy sBin. */ virtual PartonPairVec getPartons(Energy maxEnergy, const cPDPair &, const Cuts &) 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: /** * Add parton bins to pbins for the given incoming particle and the * specified cuts. */ virtual void addPartons(tPBPtr incoming ,const PDFCuts & cuts, std::deque pdf ,PartonVector & pbins) const; protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - MultiPartonExtractor & operator=(const MultiPartonExtractor &); + MultiPartonExtractor & operator=(const MultiPartonExtractor &) = delete; /** * PDFBase object to override first PDF */ vector firstPDF_; /** * PDFBase object to override second PDF */ vector secondPDF_; }; } #endif /* Herwig_MultiPartonExtractor_H */ diff --git a/Shower/Dipole/Kernels/ColourMatrixElementCorrection.h b/Shower/Dipole/Kernels/ColourMatrixElementCorrection.h --- a/Shower/Dipole/Kernels/ColourMatrixElementCorrection.h +++ b/Shower/Dipole/Kernels/ColourMatrixElementCorrection.h @@ -1,171 +1,171 @@ // -*- C++ -*- // // ColourMatrixElementCorrection.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef Herwig_ColourMatrixElementCorrection_H #define Herwig_ColourMatrixElementCorrection_H // // This is the declaration of the ColourMatrixElementCorrection class. // #include "Herwig/Shower/Dipole/Base/DipoleSplittingReweight.h" #include namespace Herwig { using namespace ThePEG; /** * \ingroup DipoleShower * \author Johan Thoren, Simon Platzer * * \brief ColourMatrixElementCorrection is implementing colour matrix element * corrections through the weighted Sudakov algorithm * * @see \ref ColourMatrixElementCorrectionInterfaces "The interfaces" * defined for ColourMatrixElementCorrection. */ class ColourMatrixElementCorrection: public DipoleSplittingReweight { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ ColourMatrixElementCorrection(); /** * The destructor. */ virtual ~ColourMatrixElementCorrection(); //@} public: /** * Calculate and cache the colour matrix element correction factor * for the given splitting type. */ double cmec(const DipoleSplittingInfo&) const; /** * Return the reweighting factor for the given splitting type. */ virtual double evaluate(const DipoleSplittingInfo& s) const { return cmec(s); } /** * Return the absolute value of the colour matrix element correction * as an enhancement hint for the sampling of the un-reweighted * splitting kernel. */ virtual double hint(const DipoleSplittingInfo& s) const { if ( hintOnly(s) ) return cmec(s); return abs(cmec(s))*lambda; } /** * Return true, if the reweight can be entirely absorbed into the hint. A * possible detuning will be switched off. */ virtual bool hintOnly(const DipoleSplittingInfo& s) const { return cmec(s) > 0.; } /** * Set the factor in front of enhance used by the veto algorithm. */ virtual void reweightFactor(const double c) { assert(c > 0.0); lambda = c; } /** * Set the factor in front of enhance used by the veto algorithm. */ virtual void negativeScaling(const double c) { assert(c >= 0.0); negCMECScaling = c; } public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * Factor to shuffle the magnitude of the CMEC between the splitting kernel * and the weight in the reweighted veto algorithm. */ double lambda; /** * Scaling factor multiplying all of the negative colour matrix element * corrections. The physically sensible value is 1.0, but this factor can * be used to examine the effects of the negative contributions. */ double negCMECScaling; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - ColourMatrixElementCorrection & operator=(const ColourMatrixElementCorrection &); + ColourMatrixElementCorrection & operator=(const ColourMatrixElementCorrection &) = delete; }; } #endif /* Herwig_ColourMatrixElementCorrection_H */ diff --git a/Shower/Dipole/Kinematics/IFMassiveKinematics.h b/Shower/Dipole/Kinematics/IFMassiveKinematics.h --- a/Shower/Dipole/Kinematics/IFMassiveKinematics.h +++ b/Shower/Dipole/Kinematics/IFMassiveKinematics.h @@ -1,279 +1,279 @@ // -*- C++ -*- // // IFLightKinematics.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_IFLightKinematics_H #define HERWIG_IFLightKinematics_H // // This is the declaration of the IFLightKinematics class. // #include "DipoleSplittingKinematics.h" namespace Herwig { using namespace ThePEG; /** * \ingroup DipoleShower * \author Simon Platzer, Martin Stoll * * \brief IFMassiveKinematics implements massless splittings * off an initial-final dipole. * * @see \ref IFMassiveKinematicsInterfaces "The interfaces" * defined for IFMassiveKinematics. */ class IFMassiveKinematics: public DipoleSplittingKinematics { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ IFMassiveKinematics(); /** * The destructor. */ virtual ~IFMassiveKinematics(); //@} public: /** * Return the boundaries in between the evolution * variable random number is to be sampled; the lower * cuoff is assumed to correspond to the infrared cutoff. */ virtual pair kappaSupport(const DipoleSplittingInfo& dIndex) const; /** * Return the boundaries in between the momentum * fraction random number is to be sampled. */ virtual pair xiSupport(const DipoleSplittingInfo& dIndex) const; /** * Return the boundaries on the momentum fraction */ virtual pair zBoundaries(Energy, const DipoleSplittingInfo&, const DipoleSplittingKernel&) const { return {0.0,1.0}; } /** * Return the dipole scale associated to the * given pair of emitter and spectator. This * should be the invariant mass or absolute value * final/final or initial/initial and the absolute * value of the momentum transfer for intial/final or * final/initial dipoles. */ virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter, const Lorentz5Momentum& pSpectator) const; /** * Return the maximum pt for the given dipole scale. */ virtual Energy ptMax(Energy dScale, double emX, double specX, const DipoleSplittingInfo& dInfo, const DipoleSplittingKernel& split) const; /** * Return the maximum pt for the given dipole scale. */ virtual Energy ptMax(Energy dScale, double, double, const DipoleIndex& dIndex, const DipoleSplittingKernel& split, tPPtr emitter, tPPtr) const; /** * Return the maximum pt for the given dipole scale. */ virtual Energy ptMax(Energy, double, double, const DipoleIndex&, const DipoleSplittingKernel&) const { // Only the DipoleSplittingInfo version should be used for massive // dipoles, for now anyway. assert(false); return ZERO; } /** * Return the maximum virtuality for the given dipole scale. */ virtual Energy QMax(Energy dScale, double emX, double specX, const DipoleSplittingInfo& dInfo, const DipoleSplittingKernel& split) const; /** * Return the maximum virtuality for the given dipole scale. */ virtual Energy QMax(Energy, double, double, const DipoleIndex&, const DipoleSplittingKernel&) const { // Only the DipoleSplittingInfo version should be used for massive // dipoles, for now anyway. assert(false); return ZERO; } /** * Return the pt given a virtuality. */ virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const; /** * Return the virtuality given a pt. */ virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const; /** * Return the random number associated to * the given pt. */ virtual double ptToRandom(Energy pt, Energy dScale, double emX, double specX, const DipoleIndex& dIndex, const DipoleSplittingKernel&) const; /** * Generate splitting variables given three random numbers * and the momentum fractions of the emitter and spectator. * Return true on success. */ virtual bool generateSplitting(double kappa, double xi, double phi, DipoleSplittingInfo& dIndex, const DipoleSplittingKernel&); /** * Generate the full kinematics given emitter and * spectator momentum and a previously completeted * DipoleSplittingInfo object. */ virtual void generateKinematics(const Lorentz5Momentum& pEmitter, const Lorentz5Momentum& pSpectator, const DipoleSplittingInfo& dInfo); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initIFMassiveKinematics; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - IFMassiveKinematics & operator=(const IFMassiveKinematics &); + IFMassiveKinematics & operator=(const IFMassiveKinematics &) = delete; private: /** * Wether or not to choose the `collinear' scheme */ bool theCollinearScheme; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of IFMassiveKinematics. */ template <> struct BaseClassTrait { /** Typedef of the first base class of IFMassiveKinematics. */ typedef Herwig::DipoleSplittingKinematics NthBase; }; /** This template specialization informs ThePEG about the name of * the IFMassiveKinematics class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::IFMassiveKinematics"; } /** * The name of a file containing the dynamic library where the class * IFMassiveKinematics is implemented. It may also include several, space-separated, * libraries if the class IFMassiveKinematics depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_IFMassiveKinematics_H */ diff --git a/Shower/Dipole/SpinCorrelations/DipoleShowerParticle.h b/Shower/Dipole/SpinCorrelations/DipoleShowerParticle.h --- a/Shower/Dipole/SpinCorrelations/DipoleShowerParticle.h +++ b/Shower/Dipole/SpinCorrelations/DipoleShowerParticle.h @@ -1,182 +1,182 @@ // -*- C++ -*- // // DipoleShowerParticle.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_DipoleShowerParticle_H #define HERWIG_DipoleShowerParticle_H // // This is the declaration of the DipoleShowerParticle class. // #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Helicity/WaveFunction/WaveFunctionBase.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "DipoleShowerVertex.h" namespace Herwig { using namespace ThePEG; /** \ingroup DipoleShower * * \author Stephen Webster * */ class DipoleShowerParticle : public Base { public: /** * Default constructor */ DipoleShowerParticle() {} /** * Default destructor **/ ~DipoleShowerParticle() {} public: /** * Reset the member variables of the object. */ void clear(); /** * Set up the decay vertex for the emitter. */ void prepare( PPtr& part, const Helicity::Direction emmDir, const Helicity::Direction specDir, const Lorentz5Momentum& pVector, const Lorentz5Momentum& nVector ); /** * Return the associated decay vertex, * a DipoleShowerVertex. **/ DSVertexPtr decayVertex() { return theDecayVertex; } /** * Create fermion decay basis states. * It returns the decay basis states in the * decay frame as required for the mapping. */ vector > createFermionDecayStates(); /** * Create vector decay basis states. */ void createVectorDecayStates(); /** * Create fermion production basis states * for the given particle produced in the splitting. */ void createNewFermionSpinInfo( PPtr& outgoing, Helicity::Direction dir); /** * Create vector production basis states * for the given particle produced in the splitting. */ void createNewVectorSpinInfo( PPtr& outgoing, Helicity::Direction dir); /** * Create the mappings between the production * and decay states for the fermion and * store them in the associated decay vertex. * (No longer applicable) reason for passing the * decay states as an argument: * Previously used a check on zero values for computing * the mapping, rather than a 1e-5, this would only * work when using the original decay state as calculated * in the decay frame (i.e. without transforming to the * lab frame and back). Now it simply avoids doing an * unnecessary rotation of the decay basis */ void setFermionMapping( const vector>& decayBasis ); /** * Create the mappings between the production * and decay states the boson and * store them in the associated decay vertex. */ void setVectorMapping(); public: /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); private: /** * The pptr to this particle. */ PPtr theParticle; /** * The dipole shower vertex associated * with this particle. */ DSVertexPtr theDecayVertex; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - DipoleShowerParticle & operator=(const DipoleShowerParticle &); + DipoleShowerParticle & operator=(const DipoleShowerParticle &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of DipoleShowerParticle. */ template <> struct BaseClassTrait { /** Typedef of the first base class of DipoleShowerParticle. */ typedef Base NthBase; }; /** This template specialization informs ThePEG about the name of * the DipoleShowerParticle class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::DipoleShowerParticle"; } /** * The name of a file containing the dynamic library where the class * DipoleShowerParticle is implemented. It may also include several, space-separated, * libraries if the class DipoleShowerParticle depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_DipoleShowerParticle_H */ diff --git a/Shower/Dipole/SpinCorrelations/DipoleShowerVertex.h b/Shower/Dipole/SpinCorrelations/DipoleShowerVertex.h --- a/Shower/Dipole/SpinCorrelations/DipoleShowerVertex.h +++ b/Shower/Dipole/SpinCorrelations/DipoleShowerVertex.h @@ -1,285 +1,285 @@ // -*- C++ -*- // // DipoleShowerVertex.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_DipoleShowerVertex_H #define HERWIG_DipoleShowerVertex_H // // This is the declaration of the DipoleShowerVertex class. // #include "ThePEG/EventRecord/HelicityVertex.h" #include "Herwig/Decay/DecayMatrixElement.h" #include "DipoleShowerVertex.fh" namespace Herwig { using namespace ThePEG; /** \ingroup DipoleShower * * This class represents the vertex for a given splitting * in the dipole shower. * * \author Stephen Webster * */ class DipoleShowerVertex: public HelicityVertex { public: /** * Default constructor */ DipoleShowerVertex(); /** * Default destructor **/ ~DipoleShowerVertex() {} public: /** * Return the matrix element for this vertex. */ inline const DecayMEPtr ME() const { return theMatrixElement; } /** * Set the matrix element */ inline void ME(DecayMEPtr in) { theMatrixElement = in; } public: /** * Method to calculate the \f$\rho\f$ matrix for one of the decay products * in the frame of this splitting vertex. * * @param iprod The splitting product to compute the \f$\rho\f$ matrix for. */ RhoDMatrix getRhoMatrix(int iprod, bool ) const; /** * Method to calculate the \f$D\f$ matrix for the decaying * particle / the incoming to the vertex, in the frame of * the vertex. The argument is a dummy argument. */ RhoDMatrix getDMatrix(int) const; /** * Get the lorentz rotation from the working frame * to the frame of the splitting. */ LorentzRotation boostToSplitting(); /** * Set the p vector for this splitting */ void pVector( const Lorentz5Momentum& emitterMom ) { thePVector = emitterMom; } /** * Set the n vector for this splitting */ void nVector( const Lorentz5Momentum& nMom ) { theNVector = nMom; } /** * Set the emitter,Spectator Config (II,IF,FF,FI - F=true, I=false) */ void dipoleConfig(const pair& newConfig) { theDipoleConfig = newConfig; } /** * Return the p vector for this splitting */ Lorentz5Momentum pVector() const { return thePVector; } /** * Return the n/spectator vector for this splitting */ Lorentz5Momentum nVector() const { return theNVector; } /** * Return the emitter,Spectator Config (II,IF,FF,FI - F=true, I=false) */ const pair& dipoleConfig() const { return theDipoleConfig; } /** * Set the decay state to production state * mapping for this vertex. */ void mappingD2P( RhoDMatrix& mapping ) { theMappingDecay2Prod = mapping; } /** * Return the mapping from the decay * states to the production states. */ RhoDMatrix mappingD2P() { return theMappingDecay2Prod; } /** * Set the production state to decay state * mapping for this vertex. */ void mappingP2D( RhoDMatrix& mapping ) { theMappingProd2Decay = mapping; } /** * Return the mapping from the production * states to the decay states. */ RhoDMatrix mappingP2D() { return theMappingProd2Decay; } /** * Set the new to old spectator mapping * for this vertex. */ void mappingSpecNewToOld( RhoDMatrix& mapping ) { theMappingSpectatorNewToOld = mapping; } /** * Return the new to old spectator mapping * for this vertex. */ RhoDMatrix mappingSpecNewToOld() { return theMappingSpectatorNewToOld; } /** * Set the new to old spectator mapping * for this vertex. */ void mappingSpecOldToNew( RhoDMatrix& mapping ) { theMappingSpectatorOldToNew = mapping; } /** * Return the new to old spectator mapping * for this vertex. */ RhoDMatrix mappingSpecOldToNew() { return theMappingSpectatorOldToNew; } public: /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); private: /** * Storage of the decay matrix element. */ DecayMEPtr theMatrixElement; /** * The p vector of the 'splitting basis' * associated with this vertex. **/ Lorentz5Momentum thePVector; /** * The n vector of the 'splitting basis' * associated with this vertex. **/ Lorentz5Momentum theNVector; /** * Initial/final config {emitter, spectator} */ pair theDipoleConfig; /** * An indicator flag to record if the * boost to shower for this vertex has been done */ bool theBoostCalculated; /** * The lorentz transformation from the * working frame to this splitting. */ LorentzRotation theBoostToSplitting; /** * The mapping from the decay basis states * to the production basis states. */ RhoDMatrix theMappingDecay2Prod; /** * The mapping from the production basis states * to the decay basis states. */ RhoDMatrix theMappingProd2Decay; /** * The mapping from the new spectator basis * states to the old spectator basis states. */ RhoDMatrix theMappingSpectatorNewToOld; /** * The mapping from the old spectator basis * states to the new spectator basis states. */ RhoDMatrix theMappingSpectatorOldToNew; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - DipoleShowerVertex & operator=(const DipoleShowerVertex &); + DipoleShowerVertex & operator=(const DipoleShowerVertex &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of DipoleShowerVertex. */ template <> struct BaseClassTrait { /** Typedef of the first base class of DipoleShowerVertex. */ typedef HelicityVertex NthBase; }; /** This template specialization informs ThePEG about the name of * the DipoleShowerVertex class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::DipoleShowerVertex"; } /** * The name of a file containing the dynamic library where the class * DipoleShowerVertex is implemented. It may also include several, space-separated, * libraries if the class DipoleShowerVertex depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_DipoleShowerVertex_H */ diff --git a/Shower/Dipole/SpinCorrelations/DipoleVertexRecord.h b/Shower/Dipole/SpinCorrelations/DipoleVertexRecord.h --- a/Shower/Dipole/SpinCorrelations/DipoleVertexRecord.h +++ b/Shower/Dipole/SpinCorrelations/DipoleVertexRecord.h @@ -1,207 +1,207 @@ // -*- C++ -*- #ifndef Herwig_DipoleVertexRecord_H #define Herwig_DipoleVertexRecord_H // // This is the declaration of the DipoleVertexRecord class. // #include "ThePEG/Config/ThePEG.h" #include "Herwig/Shower/Dipole/Base/DipoleSplittingInfo.h" #include "Herwig/Shower/Dipole/Base/Dipole.h" #include "ThePEG/EventRecord/RhoDMatrix.h" #include "DipoleShowerParticle.h" namespace Herwig { using namespace ThePEG; /** * Here is the documentation of the DipoleVertexRecord class. */ class DipoleVertexRecord: public Base { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ DipoleVertexRecord() {} /** * The destructor. */ virtual ~DipoleVertexRecord() { clear(); } //@} public: /** * Prepare the emitter and spectator * for the spin correlations computations. **/ void prepareSplitting( const DipoleSplittingInfo& dInfo, const Dipole& dip); /** * Correctly initialise the decay matrix * to a delta matrix for an external particle. */ void initDecayMatrix(PPtr& particle, Helicity::Direction dir); /** * Compute the spin density matrix for the given emitter. * This tracks the path between the given emitter and * the previous emitter, calculating a rho/decay matrix * at each vertex as appropriate. */ RhoDMatrix emitterDensityMatrix(PPtr emitter); /** * Generate the spin-correlated azimuthal angle for a splitting. */ void generatePhi(DipoleSplittingInfo& dInfo, Dipole& dip); /** * Identify the type of particle and use the appropriate function * to set up the spin info. * Required for e.g. MPI */ void createSpinInfo(PPtr& part, const Helicity::Direction& dir); /** * Create and set up fermion spin info. * Required for e.g. MPI */ void createFermionSpinInfo(PPtr& part, const Helicity::Direction& dir); /** * Create and set up vector spin info. * Required for e.g. MPI */ void createVectorSpinInfo(PPtr& part, const Helicity::Direction& dir); /** * Update the vertex record following a splitting. */ void update(const DipoleSplittingInfo& dInfo); /** * For spectators. Set new particle spin info the that of the * old particle. Update the spin info to include any momentum changes. */ void updateSpinInfo( PPtr& oldPart, PPtr& newPart ); /** * Set the stopUpdate flag in the spin info of a particle * incoming to the current decay. */ void prepareParticleDecay( const PPtr& parent ); /** * Update the spin info of the incoming to the decay * following showering of the decay. */ void updateParticleDecay(); /** * SW 06/02/2019: Required for NearestNeighbourDipoleAnalysis tests. * Access the emitter info record. */ //map emitterInfoRecord() const { //return theEmitterInfoRecord; //} /** * SW 06/02/2019: Required for NearestNeighbourDipoleAnalysis tests. * Add a splitting to the emitter info record. */ // void addToRecord(const DipoleSplittingInfo& dInfo) { // assert(dInfo.emitter()); // theEmitterInfoRecord[dInfo.emitter()] = dInfo; // } /** * Clear the vertex record: Give up ownership * on any object involved in the evolution. */ virtual void clear(); /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); private: /** * The current emitter. */ DipoleShowerParticle theCurrentEmitter; /** * SW 06/02/2019: Required for NearestNeighbourDipoleAnalysis tests. * Record of the splittings as * required for the testing analysis. */ //map theEmitterInfoRecord; /** * The spin info of a particle incoming to the decay * under consideration. */ tcSpinPtr theDecayParentSpinInfo; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - DipoleVertexRecord & operator=(const DipoleVertexRecord &); + DipoleVertexRecord & operator=(const DipoleVertexRecord &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of DipoleVertexRecord. */ template <> struct BaseClassTrait { /** Typedef of the first base class of DipoleVertexRecord. */ typedef Base NthBase; }; /** This template specialization informs ThePEG about the name of * the DipoleVertexRecord class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::DipoleVertexRecord"; } /** * The name of a file containing the dynamic library where the class * DipoleVertexRecord is implemented. It may also include several, space-separated, * libraries if the class DipoleVertexRecord depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* Herwig_DipoleVertexRecord_H */