diff --git a/Decay/Baryon/BaryonFactorizedDecayer.cc b/Decay/Baryon/BaryonFactorizedDecayer.cc --- a/Decay/Baryon/BaryonFactorizedDecayer.cc +++ b/Decay/Baryon/BaryonFactorizedDecayer.cc @@ -1,809 +1,809 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the BaryonFactorizedDecayer class. // #include "BaryonFactorizedDecayer.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/PDT/DecayMode.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorBarWaveFunction.h" #include "Herwig/Decay/DecayVertex.h" #include "ThePEG/Helicity/FermionSpinInfo.h" #include "ThePEG/Helicity/RSFermionSpinInfo.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" using namespace Herwig; using namespace ThePEG::Helicity; BaryonFactorizedDecayer::BaryonFactorizedDecayer() { // default values taken from PRD56, 2799 _a1c= 1.1; _a2c=-0.5; _a1b= 1.0; _a2b= 0.28; // intermediates generateIntermediates(true); } void BaryonFactorizedDecayer::doinitrun() { _current->initrun(); _form->initrun(); DecayIntegrator::doinitrun(); _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()); } } void BaryonFactorizedDecayer::doinit() { DecayIntegrator::doinit(); // get the CKM matrix (unsquared for interference) Complex ckmmat[3][3]; vector< vector > CKM(_theCKM->getUnsquaredMatrix(SM().families())); for(unsigned int ix=0;ix<3;++ix) { for(unsigned int iy=0;iy<3;++iy) { ckmmat[ix][iy]=CKM[ix][iy]; } } // make sure the current and form factor got initialised _current->init(); _form->init(); // find all the possible modes vector tformmap,tcurrmap; vector inquark,outquark,currq,curra; tPDVector incoming; vector outgoing; for(unsigned int iform=0;iform<_form->numberOfFactors();++iform) { // particles from the form factor int id0,id1; _form->particleID (iform,id0,id1); int spect1,spect2,inq,outq,ispin,ospin; _form->formFactorInfo(iform,ispin,ospin,spect1,spect2,inq,outq); // particles from the form factor tPDPtr in = getParticleData(id0); tPDPtr out = getParticleData(id1); // the charge of the decay products int Wcharge = in->iCharge()-out->iCharge(); // max mass for the particles in the current Energy min = in->massMax()-out->massMin(); for(unsigned int icurr=0;icurr<_current->numberOfModes();++icurr) { // get the particles from the current int iq,ia; _current->decayModeInfo(icurr,iq,ia); tPDVector ptemp=_current->particles(Wcharge,icurr,iq,ia); tPDVector outV = {out}; outV.insert(std::end(outV), std::begin(ptemp), std::end(ptemp)); Energy minb=ZERO; for(unsigned int iz=0;izmassMin(); // valid mode if(outV.size()>1&&minb0&&inq%2!=iq%2)|| (inq<0&&abs(inq)%2!=abs(ia)%2))))) { tformmap.push_back(iform);tcurrmap.push_back(icurr); incoming.push_back(in); outgoing.push_back(outV); inquark.push_back(inq);outquark.push_back(outq); currq.push_back(iq);curra.push_back(ia); } // if the meson is neutral try the CC mode if(Wcharge==0&&iq!=-ia&&((inq>0&&inq%2!=iq%2)|| (inq<0&&abs(inq)%2!=abs(ia)%2))) { ptemp=_current->particles(Wcharge,icurr,-ia,-iq); tPDVector outV = {out}; outV.insert(std::end(outV), std::begin(ptemp), std::end(ptemp)); Energy minb=ZERO; for(unsigned int iz=0;izmassMin(); if(outV.size()>1&&minb modeloc; vector modecc; findModes(ix,incoming,outgoing,modeloc,modecc); // if more than two outgoing particles only allow one diagram if ( outgoing[ix].size() > 2 && !modeloc.empty() ) {break;} // create the mode and set the particles as for the first instance PhaseSpaceModePtr mode=new_ptr(PhaseSpaceMode(incoming[ix],outgoing[ix],1.)); PhaseSpaceChannel channel((PhaseSpaceChannel(mode),0,1)); Energy min = incoming[ix]->massMax()-outgoing[ix][0]->massMin(); int Wcharge = incoming[ix]->iCharge()-outgoing[ix][0]->iCharge(); bool done = _current->createMode(Wcharge,tcPDPtr(),FlavourInfo(), tcurrmap[ix],mode,1,0,channel,min); if(!done){throw InitException() << "Failed to construct mode in " << "BaryonFactorizedDecayer::doinit()." << Exception::abortnow;} // set the parameters for the additional modes vectorttform,ttcurr; ttform.push_back(tformmap[ix]);ttcurr.push_back(tcurrmap[ix]); for(unsigned int iy=0;iy tCKM; Complex ckm; for(unsigned int iy=0;iyparticleID(ttform[iy],id0,id1); int Wcharge = getParticleData(id0)->iCharge()-getParticleData(id1)->iCharge(); Complex ckm=1.; if(Wcharge!=0) { if(abs(iq)%2==0){ckm *= conj(ckmmat[abs(iq)/2-1][(abs(ia)-1)/2]);} else{ckm *= conj(ckmmat[abs(ia)/2-1][(abs(iq)-1)/2]);} if(abs(inq)%2==0){ckm *= ckmmat[abs(inq)/2-1][(abs(outq)-1)/2];} else{ckm *= ckmmat[abs(outq)/2-1][(abs(inq)-1)/2];} if(abs(inq)==5){ckm*=_a1b;} else{ckm*=_a1c;} } else { if(inq>0) { if(abs(inq)%2==0){ckm *= ckmmat[abs(inq)/2-1][(abs(iq)-1)/2];} else{ckm *= ckmmat[abs(iq)/2-1][(abs(inq)-1)/2];} if(abs(outq)%2==0) {ckm *= conj(ckmmat[abs(outq)/2-1][(abs(ia)-1)/2]);} else{ckm *= conj(ckmmat[abs(ia)/2-1][(abs(outq)-1)/2]);} } else { if(abs(inq)%2==0){ckm *= ckmmat[abs(inq)/2-1][(abs(ia)-1)/2];} else{ckm *= ckmmat[abs(ia)/2-1][(abs(inq)-1)/2];} if(abs(outq)%2==0) {ckm *= conj(ckmmat[abs(outq)/2-1][(abs(iq)-1)/2]);} else{ckm *= conj(ckmmat[abs(iq)/2-1][(abs(outq)-1)/2]);} } if(abs(inq)==5){ckm*=_a2b;} else{ckm*=_a2c;} } if((abs(inq)%2==0&&inq<0)||(abs(inq)%2!=0&&inq>0)){ckm=conj(ckm);} tCKM.push_back(ckm); } // add the parameters for the mode to the list _currentmap.push_back(ttcurr); _formmap.push_back(ttform); _factCKM.push_back(tCKM); double maxweight(0.); // add the mode to the list if(_wgtmax.size()>numberModes()) maxweight=_wgtmax[numberModes()]; // the weights for the channel vector channelwgts; if(_wgtloc.size()>numberModes()&& _wgtloc[numberModes()]+mode->channels().size()<=_weights.size()) { vector::iterator start=_weights.begin()+_wgtloc[numberModes()]; vector::iterator end = start+mode->channels().size(); channelwgts=vector(start,end); } else { channelwgts.resize(mode->channels().size(),1./(mode->channels().size())); } // don't need channels for two body decays if(outgoing[ix].size()==2) { channelwgts.clear(); mode = new_ptr(PhaseSpaceMode(incoming[ix],outgoing[ix],maxweight)); } else { mode->maxWeight(maxweight); mode->setWeights(channelwgts); } addMode(mode); // resize the duplicate modes to remove them for(unsigned int iy=0;iyid()),ibaryon,foundb,id0,id1; vector idall,idother; tPDVector::const_iterator pit = children.begin(); tPDVector::const_iterator pend = children.end(); for( ; pit!=pend;++pit){idall.push_back((**pit).id());} // loop over the particles in the form factor do { _form->particleID(iform,id0,id1); ibaryon=0; if(id0==idin){ibaryon=id1;} else if(id0==-idin){ibaryon=-id1;} if(ibaryon!=0) { foundb=false; idother.clear(); for(ix=0;ixaccept(idother);} } ++iform; } while(!allowed&&iform<_form->numberOfFactors()); return allowed; } int BaryonFactorizedDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { unsigned int ix,iy; int idin(parent->id()),ibaryon,foundb,id0,id1,icurr(-1),iform(0); vector idall,idother; tPDVector::const_iterator pit = children.begin(); tPDVector::const_iterator pend = children.end(); for( ; pit!=pend;++pit){idall.push_back((**pit).id());} // loop over the particles in the form factor do { _form->particleID(iform,id0,id1); ibaryon=0; if(id0==idin){ibaryon=id1;} else if(id0==-idin){ibaryon=-id1;} ++iform; foundb=false; idother.clear(); for(ix=0;ixdecayMode(idother);} } while(icurr<0&&iformnumberOfFactors())); // now find the mode int imode=-1; ix=0; --iform; do { for(iy=0;iy<_currentmap[ix].size();++iy) {if(int(_currentmap[ix][iy])==icurr&&int(_formmap[ix][iy])==iform){imode=ix;}} ++ix; } while(imode<0&&ix> _current >> _form >> _a1b >> _a2b >>_a1c >>_a2c >> _currentmap >> _formmap >> _factCKM >> _wgtloc >> _wgtmax >> _weights >> _theCKM; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigBaryonFactorizedDecayer("Herwig::BaryonFactorizedDecayer", "HwBaryonDecay.so"); void BaryonFactorizedDecayer::Init() { static ClassDocumentation documentation ("The BaryonFactorizedDecayer class combines the baryon form factor and a" " weak current to perform a decay in the naive factorization approximation."); static Reference interfaceWeakCurrent ("Current", "The reference for the decay current to be used.", &BaryonFactorizedDecayer::_current, false, false, true, false, false); static ParVector interfaceWeightLocation ("WeightLocation", "The locations of the weights for a given channel in the vector", &BaryonFactorizedDecayer::_wgtloc, 0, 0, 0, 0, 10000, false, false, true); static ParVector interfaceWeightMax ("MaximumWeight", "The maximum weight for a given channel.", &BaryonFactorizedDecayer::_wgtmax, 0, 0, 0, 0., 100., false, false, true); static ParVector interfaceWeights ("Weights", "The weights for the integration.", &BaryonFactorizedDecayer::_weights, 0, 0, 0, 0., 1., false, false, true); static Reference interfaceFormFactor ("FormFactor", "The form-factor", &BaryonFactorizedDecayer::_form, true, true, true, false, false); static Parameter interfacea1Bottom ("a1Bottom", "The factorization paramter a_1 for decays of bottom baryons", &BaryonFactorizedDecayer::_a1b, 1., -10.0, 10.0, false, false, true); static Parameter interfacea2Bottom ("a2Bottom", "The factorization paramter a_2 for decays of bottom baryons", &BaryonFactorizedDecayer::_a2b, 0.28, -10.0, 10.0, false, false, true); static Parameter interfacea1Charm ("a1Charm", "The factorization paramter a_1 for decays of charm baryons", &BaryonFactorizedDecayer::_a1c, 1.1, -10.0, 10.0, false, false, true); static Parameter interfacea2Charm ("a2Charm", "The factorization paramter a_2 for decays of charm baryons", &BaryonFactorizedDecayer::_a2c, -0.5, -10.0, 10.0, false, false, true); static Reference interfaceCKM ("CKM", "Reference to the Standard Model object", &BaryonFactorizedDecayer::_theCKM, false, false, true, false); } void BaryonFactorizedDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // for the decaying particle if(part.id()>0) { SpinorWaveFunction:: constructSpinInfo(_inHalf,const_ptr_cast(&part),incoming,true); } else { SpinorBarWaveFunction:: constructSpinInfo(_inHalfBar,const_ptr_cast(&part),incoming,true); } // decay product // spin 1/2 if(decay[0]->dataPtr()->iSpin()==PDT::Spin1Half) { if(part.id()>0) { SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true); } else { SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true); } } // spin 3/2 else if(decay[0]->dataPtr()->iSpin()==PDT::Spin3Half) { if(part.id()>0) { RSSpinorBarWaveFunction::constructSpinInfo(_inThreeHalfBar, decay[0],outgoing,true); } else { RSSpinorWaveFunction::constructSpinInfo(_inThreeHalf, decay[0],outgoing,true); } } else assert(false); // and the stuff from the current _current->constructSpinInfo(ParticleVector(decay.begin()+1,decay.end())); } double BaryonFactorizedDecayer::me2(const int ichan, const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { double me(0.); assert(part.dataPtr()->iSpin()==PDT::Spin1Half); if(outgoing[0]->iSpin()==PDT::Spin1Half) me=halfHalf(ichan,part,outgoing,momenta,meopt); else if(outgoing[0]->iSpin()==PDT::Spin3Half) me=halfThreeHalf(ichan,part,outgoing,momenta,meopt); else assert(false); return me; } // matrix element for a 1/2 -> 1/2 decay double BaryonFactorizedDecayer::halfHalf(const int ichan, const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { Energy scale; // extract the spins of the particles vector spin; for(unsigned ix=0;ixiSpin()); if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,spin))); if(meopt==Initialize) { // spinors and rho if(part.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&part), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&part), incoming); } ME()->zero(); // spinors for the decay product if(part.id()>0) { _inHalfBar.resize(2); for(unsigned int ihel=0;ihel<2;++ihel) _inHalfBar[ihel] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ihel,Helicity::outgoing); } else { _inHalf.resize(2); for(unsigned int ihel=0;ihel<2;++ihel) _inHalf[ihel] = HelicityFunctions::dimensionedSpinor (-momenta[0],ihel,Helicity::outgoing); } // get the information on the form-factor int id0(part.id()),id1(outgoing[0]->id()); // work out the value of q and calculate the form factors Lorentz5Momentum q(part.momentum()-momenta[0]); q.rescaleMass(); Energy m0(part.mass()),m1(momenta[0].mass()); Energy2 q2(q.mass2()); Lorentz5Momentum sum(part.momentum()+momenta[0]); // calculate the baryonic part of the current for the decay double pre(0.); for(unsigned int mode=0;mode<_formmap[imode()].size();++mode) { Complex f1v,f2v,f3v,f1a,f2a,f3a; // calculate the form factor piece _form->SpinHalfSpinHalfFormFactor(q2,_formmap[imode()][mode],id0,id1,m0,m1, - f1v,f2v,f3v,f1a,f2a,f3a); + f1v,f2v,f3v,f1a,f2a,f3a,FlavourInfo()); Complex left = f1v-f1a-f2v-double((m0-m1)/(m0+m1))*f2a; Complex right = f1v+f1a-f2v+double((m0-m1)/(m0+m1))*f2a; vector baryon(4); for(unsigned int ix=0;ix<2;++ix) { for(unsigned int iy=0;iy<2;++iy) { LorentzPolarizationVectorE vtemp = _inHalf[ix].generalCurrent(_inHalfBar[iy],left,right); complex vspin=_inHalf[ix].scalar(_inHalfBar[iy]); complex aspin=_inHalf[ix].pseudoScalar(_inHalfBar[iy]); // the momentum like pieces if(part.id()>0) { vtemp+= (f2v*vspin+f2a*aspin)/(m0+m1)*sum; vtemp+= (f3v*vspin+f3a*aspin)/(m0+m1)*q; } else { vtemp+= (f2v*vspin-f2a*aspin)/(m0+m1)*sum; vtemp+= (f3v*vspin-f3a*aspin)/(m0+m1)*q; } if(part.id()>0) baryon[2*ix+iy]=vtemp; else baryon[2*iy+ix]=vtemp; } } // construct the weak current vector hadron = _current->current(tcPDPtr(),FlavourInfo(), _currentmap[imode()][mode],ichan,scale, tPDVector(outgoing.begin()+1,outgoing.end()), vector(momenta.begin()+1,momenta.end()),meopt); pre=pow(part.mass()/scale,int(outgoing.size()-3));pre*=pre; vector constants(outgoing.size()+1),ihel(outgoing.size()+1); int itemp=1; unsigned int ibar=0; for(int iz=int(outgoing.size()-1);iz>=0;--iz) { if(abs(outgoing[iz]->id())!=id1) { itemp *= outgoing[iz]->iSpin(); constants[iz]=itemp; } else ibar=iz; constants[outgoing.size()]=1; constants[ibar]=constants[ibar+1]; } for(unsigned int mhel=0;mhel0;--ix) { if(ix-1!=ibar){ihel[ix]=(lhel%constants[ix-1])/constants[ix];}} (*ME())(ihel) += Complex(hadron[lhel].dot(baryon[mhel])* _factCKM[imode()][mode]*SM().fermiConstant()); } } } // return the answer return 0.5*pre*(ME()->contract(_rho)).real(); } // matrix element for a 1/2 -> 3/2 decay double BaryonFactorizedDecayer::halfThreeHalf(const int ichan, const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { // spins Energy scale; vector spin(outgoing.size()); for(unsigned int ix=0;ixiSpin(); if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,spin))); // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&part), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&part), incoming); } ME()->zero(); // spinors for the decay product LorentzPolarizationVector in=UnitRemoval::InvE*part.momentum(); if(part.id()>0) { RSSpinorBarWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalfBar.resize(4); _inHalfBar.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalfBar[ihel]=swave.dimensionedWf(); _inHalfBar[ihel] = _inThreeHalfBar[ihel].dot(in); } } else { RSSpinorWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalf.resize(4); _inHalf.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalf[ihel]=swave.dimensionedWf(); _inHalf[ihel] = _inThreeHalf[ihel].dot(in); } } // get the information on the form-factor int id0(part.id()),id1(outgoing[0]->id()); // work out the value of q and calculate the form factors Lorentz5Momentum q(part.momentum()-momenta[0]); q.rescaleMass(); Energy m0(part.mass()),m1(momenta[0].mass()); Energy2 q2(q.mass2()); Lorentz5Momentum sum(part.momentum()+momenta[0]); InvEnergy ms(1./(m0+m1)); InvEnergy2 ms2(ms*ms); double pre(0.); for(unsigned int mode=0;mode<_formmap[imode()].size();++mode) { // calculate the form factors Complex f1v,f2v,f3v,f4v,f1a,f2a,f3a,f4a; _form->SpinHalfSpinThreeHalfFormFactor(q2,_formmap[imode()][mode],id0,id1,m0,m1, - f1v,f2v,f3v,f4v,f1a,f2a,f3a,f4a); + f1v,f2v,f3v,f4v,f1a,f2a,f3a,f4a,FlavourInfo()); complex lS1,lS2,rS1,rS2; Complex left,right; complex lV,rV; if(part.id()>0) { left = f1a-f1v; right = f1a+f1v; lS1 = ms2*(f3a-f4a-f3v+f4v); rS1 = ms2*(f3a-f4a+f3v-f4v); lS2 = ms2*(f4a-f4v); rS2 = ms2*(f4a+f4v); lV = ms*(f2a-f2v); rV = ms*(f2a+f2v); } else { left = conj(f1a+f1v); right = conj(f1a-f1v); lS1 = ms2*conj(f3a-f4a+f3v-f4v); rS1 = ms2*conj(f3a-f4a-f3v+f4v); lS2 = ms2*conj(f4a-f4v); rS2 = ms2*conj(f4a+f4v); lV = ms *conj(f2a-f2v); rV = ms *conj(f2a+f2v); } // construct the vectors for the decay LorentzPolarizationVectorE baryon[4][2]; for(unsigned int iya=0;iya<4;++iya) { for(unsigned int ixa=0;ixa<2;++ixa) { unsigned int ix,iy; if(outgoing[0]->id()>0) { ix=iya; iy=ixa; } else { ix=ixa; iy=iya; } // scalar like terms complex lfact = _inHalf[iy].leftScalar( _inHalfBar[ix]); complex rfact = _inHalf[iy].rightScalar(_inHalfBar[ix]); Complex scalar1 = (lS1*lfact+rS1*rfact)*UnitRemoval::E; Complex scalar2 = (lS2*lfact+rS2*rfact)*UnitRemoval::E; LorentzPolarizationVector svec = _inHalf[iy].generalCurrent(_inHalfBar[ix],lV/ms,rV/ms)*ms; LorentzPolarizationVectorE tvec; if(part.id()>0) { tvec=_inThreeHalfBar[ix].generalCurrent(_inHalf[iy],left,right); } else { tvec=_inThreeHalf[iy].generalCurrent(_inHalfBar[ix],left,right); } baryon[iya][ixa] = tvec+svec*UnitRemoval::E +scalar1*momenta[0]+scalar2*part.momentum(); } } vector hadron = _current->current(tcPDPtr(),FlavourInfo(),_currentmap[imode()][mode],ichan,scale, tPDVector(outgoing.begin()+1,outgoing.end()), vector(momenta.begin()+1,momenta.end()),meopt); // prefactor pre = pow(part.mass()/scale,int(outgoing.size()-3)); pre *= pre; // work out the mapping for the hadron vector vector constants(outgoing.size()+1),ihel(outgoing.size()+1); int itemp = 1; int ibar = 0; for(int ix=int(outgoing.size()-1);ix>=0;--ix) { if(abs(outgoing[ix]->id())!=id1) { itemp*=outgoing[ix]->iSpin(); constants[ix]=itemp; } else{ibar=ix;} } constants[outgoing.size()]=1; constants[ibar]=constants[ibar+1]; for(unsigned int iya=0;iya<4;++iya) { ihel[1]=iya; for(unsigned int ixa=0;ixa<2;++ixa) { ihel[0]=ixa; for(unsigned int lhel=0;lhel0;--ix) {if(ix-1!=ibar){ihel[ix]=(lhel%constants[ix-1])/constants[ix];}} (*ME())(ihel) += Complex(hadron[lhel].dot(baryon[iya][ixa])* _factCKM[imode()][mode]*SM().fermiConstant()); } } } } // return the answer return 0.5*pre*(ME()->contract(_rho)).real(); } void BaryonFactorizedDecayer::findModes(unsigned int imode, tPDVector & incoming, vector & outgoing, vector & loc, vector & cc) { // get the id's for the mode // incoming int id_in = incoming[imode]->id(); int idbar_in = incoming[imode]->CC() ? incoming[imode]->CC()->id() : incoming[imode]->id(); // outgoing vector id_out,idbar_out; for(unsigned int ix=0;ixid()); if(outgoing[imode][ix]->CC()) idbar_out.push_back(outgoing[imode][ix]->CC()->id()); else idbar_out.push_back(id_out[ix]); } // loop over the modes for(unsigned int ix=0;ixid()==id_in&&outgoing[ix].size()==id_out.size()) { vector done(id_out.size(),false); unsigned int nfound = 0; for(unsigned int iy=0;iyid(); unsigned int iz = 0; bool found = false; do { if(idtemp==id_out[iz]&&!done[iz]) { done[iz]=true; found=true; } ++iz; } while(izid()==idbar_in&&outgoing[ix].size()==idbar_out.size()) { vector done(id_out.size(),false); unsigned int nfound=0; for(unsigned int iy=0;iyid(); unsigned int iz=0; bool found = false; do { if(idtemp==idbar_out[iz]&&!done[iz]) { done[iz]=true; found=true; } ++iz; } while(izname() << " \n"; for(ix=0;ix<_wgtloc.size();++ix) {output << "insert " << name() << ":WeightLocation " << ix << " " << _wgtloc[ix] << "\n";} for(ix=0;ix<_wgtmax.size();++ix) {output << "insert " << name() << ":MaximumWeight " << ix << " " << _wgtmax[ix] << "\n";} for(ix=0;ix<_weights.size();++ix) {output << "insert " << name() << ":Weights " << ix << " " << _weights[ix] << "\n";} _current->dataBaseOutput(output,false,true); output << "newdef " << name() << ":Current " << _current->name() << " \n"; _form->dataBaseOutput(output,false,true); output << "newdef " << name() << ":FormFactor " << _form->name() << " \n"; if(header){output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;} } diff --git a/Decay/Baryon/SemiLeptonicBaryonDecayer.cc b/Decay/Baryon/SemiLeptonicBaryonDecayer.cc --- a/Decay/Baryon/SemiLeptonicBaryonDecayer.cc +++ b/Decay/Baryon/SemiLeptonicBaryonDecayer.cc @@ -1,499 +1,499 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the SemiLeptonicBaryonDecayer class. // #include "SemiLeptonicBaryonDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorBarWaveFunction.h" #include "ThePEG/Helicity/LorentzPolarizationVector.h" #include "ThePEG/Helicity/FermionSpinInfo.h" #include "ThePEG/Helicity/RSFermionSpinInfo.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include "ThePEG/Helicity/HelicityFunctions.h" using namespace Herwig; using namespace ThePEG::Helicity; SemiLeptonicBaryonDecayer::SemiLeptonicBaryonDecayer() { // intermediates generateIntermediates(true); } void SemiLeptonicBaryonDecayer::doinitrun() { _current->initrun(); _form->initrun(); DecayIntegrator::doinitrun(); if(initialize()) { _maxwgt.clear(); for(unsigned int ix=0;ixmaxWeight()); } } void SemiLeptonicBaryonDecayer::doinit() { DecayIntegrator::doinit(); // make sure the current got initialised _current->init(); // and the form factors _form->init(); // the channels _modemap.clear(); for(unsigned int ix=0;ix<_form->numberOfFactors();++ix) { int id0(0),id1(0); // get the external particles for this mode _form->particleID(ix,id0,id1); int inspin,spect1,spect2,inquark,outquark,outspin; _form->formFactorInfo(ix,inspin,outspin,spect1,spect2,inquark,outquark); // incoming and outgoing particles tPDPtr in = getParticleData(id0); tPDPtr out = getParticleData(id1); // charge of the W int Wcharge =(in->iCharge()-out->iCharge()); Energy min = in->mass()+in->widthUpCut() -out->mass()+out->widthLoCut(); _modemap.push_back(numberModes()); for(unsigned int iy=0;iy<_current->numberOfModes();++iy) { int iq(0),ia(0); tPDVector outV = {out}; tPDVector ptemp=_current->particles(Wcharge,iy,iq,ia); outV.insert(std::end(outV), std::begin(ptemp), std::end(ptemp)); // create the mode PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,outV,1.)); // create the first piece of the channel PhaseSpaceChannel channel((PhaseSpaceChannel(mode),0,1)); // and the rest bool done = _current->createMode(Wcharge,tcPDPtr(),FlavourInfo(), iy,mode,1,0,channel,min); // check the result if(done&&abs(Wcharge)==3&&inspin==2&&(outspin==2||outspin==4)) { // the maximum weight double maxweight = _maxwgt.size()>numberModes() ? _maxwgt[numberModes()] : 2.; mode->maxWeight(maxweight); addMode(mode); } } } } bool SemiLeptonicBaryonDecayer::accept(tcPDPtr parent, const tPDVector & children) const { // find the non-lepton int ibar(0),idtemp,idin(parent->id()); vector idother; bool dummy; tPDVector::const_iterator pit = children.begin(); tPDVector::const_iterator pend = children.end(); for( ; pit!=pend;++pit) { idtemp=(**pit).id(); if(abs(idtemp)>16) ibar=idtemp; else idother.push_back(idtemp); } // check that the form factor exists if(_form->formFactorNumber(idin,ibar,dummy)<0) return false; // and the current return _current->accept(idother); } int SemiLeptonicBaryonDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { // find the ids of the particles for the decay current tPDVector::const_iterator pit = children.begin(); tPDVector::const_iterator pend = children.end(); int idtemp,ibar(0),idin(parent->id()); vector idother; cc=false; for( ; pit!=pend;++pit) { idtemp=(**pit).id(); if(abs(idtemp)>16) ibar=idtemp; else idother.push_back(idtemp); } return _modemap[_form->formFactorNumber(idin,ibar,cc)] +_current->decayMode(idother); } void SemiLeptonicBaryonDecayer::persistentOutput(PersistentOStream & os) const { os << _current << _form << _maxwgt << _modemap; } void SemiLeptonicBaryonDecayer::persistentInput(PersistentIStream & is, int) { is >> _current >> _form >> _maxwgt >> _modemap; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigSemiLeptonicBaryonDecayer("Herwig::SemiLeptonicBaryonDecayer", "HwBaryonDecay.so"); void SemiLeptonicBaryonDecayer::Init() { static ClassDocumentation documentation ("The SemiLeptonicBaryonDecayer class is designed for" " the semi-leptonic decay of the baryons."); static Reference interfaceCurrent ("Current", "The current for the leptons produced in the decay.", &SemiLeptonicBaryonDecayer::_current, true, true, true, false, false); static Reference interfaceFormFactor ("FormFactor", "The form factor", &SemiLeptonicBaryonDecayer::_form, true, true, true, false, false); static ParVector interfaceMaximumWeight ("MaximumWeight", "The maximum weights for the decays", &SemiLeptonicBaryonDecayer::_maxwgt, 0, 0, 0, 0, 10000, false, false, true); } void SemiLeptonicBaryonDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // for the decaying particle if(part.id()>0) { SpinorWaveFunction:: constructSpinInfo(_inHalf,const_ptr_cast(&part),incoming,true); } else { SpinorBarWaveFunction:: constructSpinInfo(_inHalfBar,const_ptr_cast(&part),incoming,true); } // decay product // spin 1/2 if(decay[0]->dataPtr()->iSpin()==PDT::Spin1Half) { if(part.id()>0) { SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true); } else { SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true); } } // spin 3/2 else if(decay[0]->dataPtr()->iSpin()==PDT::Spin3Half) { if(part.id()>0) { RSSpinorBarWaveFunction::constructSpinInfo(_inThreeHalfBar, decay[0],outgoing,true); } else { RSSpinorWaveFunction::constructSpinInfo(_inThreeHalf, decay[0],outgoing,true); } } else assert(false); // and the stuff from the current _current->constructSpinInfo(ParticleVector(decay.begin()+1,decay.end())); } // combine the currents and form-factors to give the matrix element double SemiLeptonicBaryonDecayer::me2(const int , const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { assert(part.dataPtr()->iSpin()==PDT::Spin1Half); double me(0.); if(outgoing[0]->iSpin()==PDT::Spin1Half) me = halfHalf(part,outgoing,momenta,meopt); else if(outgoing[0]->iSpin()==PDT::Spin3Half) me = halfThreeHalf(part,outgoing,momenta,meopt); else assert(false); return me; } // matrix element for a 1/2 -> 1/2 semi-leptonic decay double SemiLeptonicBaryonDecayer::halfHalf(const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&part), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&part), incoming); // work out the mapping for the lepton vector _constants.resize(outgoing.size()+1); _ispin.resize(outgoing.size()); int itemp(1); _ibar=0; for(int ix=int(outgoing.size()-1);ix>=0;--ix) { _ispin[ix]=outgoing[ix]->iSpin(); if(abs(outgoing[ix]->id())<=16) { itemp*=_ispin[ix]; _constants[ix]=itemp; } else _ibar=ix; } _constants[outgoing.size()]=1; _constants[_ibar]=_constants[_ibar+1]; } if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,_ispin))); // spinors for the decay product if(part.id()>0) { _inHalfBar.resize(2); for(unsigned int ihel=0;ihel<2;++ihel) _inHalfBar[ihel] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ihel,Helicity::outgoing); } else { _inHalf.resize(2); for(unsigned int ihel=0;ihel<2;++ihel) _inHalf[ihel] = HelicityFunctions::dimensionedSpinor (-momenta[0],ihel,Helicity::outgoing); } // get the information on the form-factor int spinin(0),spinout(0),spect1,spect2,inquark,outquark; int id0(part.id()),id1(outgoing[0]->id()); bool cc; int iloc(_form->formFactorNumber(id0,id1,cc)); _form->formFactorInfo(iloc,spinin,spinout,spect1,spect2,inquark,outquark); // work out the value of q and calculate the form factors Lorentz5Momentum q(part.momentum()-momenta[0]); q.rescaleMass(); Energy m0(part.mass()),m1(momenta[0].mass()); Energy2 q2(q.mass2()); Lorentz5Momentum sum(part.momentum()+momenta[0]); // calculate the form factors Complex f1v,f2v,f3v,f1a,f2a,f3a; _form->SpinHalfSpinHalfFormFactor(q2,iloc,id0,id1,m0,m1, - f1v,f2v,f3v,f1a,f2a,f3a); + f1v,f2v,f3v,f1a,f2a,f3a,FlavourInfo()); // calculate the hadronic current for the decay vector hadron(4); Complex left =f1v-f1a-f2v-double((m0-m1)/(m0+m1))*f2a; Complex right =f1v+f1a-f2v+double((m0-m1)/(m0+m1))*f2a; LorentzPolarizationVectorE vtemp; for(unsigned int ix=0;ix<2;++ix) { for(unsigned int iy=0;iy<2;++iy) { vtemp = _inHalf[ix].generalCurrent(_inHalfBar[iy],left,right); complex vspin = _inHalf[ix].scalar(_inHalfBar[iy]); complex aspin = _inHalf[ix].pseudoScalar(_inHalfBar[iy]); // the momentum like pieces if(part.id()>0) { vtemp+= (f2v*vspin+f2a*aspin)/(m0+m1)*sum; vtemp+= (f3v*vspin+f3a*aspin)/(m0+m1)*q; } else { vtemp-= (f2v*vspin-f2a*aspin)/(m0+m1)*sum; vtemp+= (f3v*vspin-f3a*aspin)/(m0+m1)*q; } if(part.id()>0) hadron[2*ix+iy]=vtemp; else hadron[2*iy+ix]=vtemp; } } // construct the lepton current Energy scale; int mode((abs(outgoing[1]->id())-11)/2); vector lepton(_current->current(tcPDPtr(),FlavourInfo(), mode,-1,scale,tPDVector(outgoing.begin()+1,outgoing.end()), vector(momenta.begin()+1,momenta.end()), meopt)); // matrix element vector ihel(outgoing.size()+1); unsigned int mhel,ix,lhel; for(mhel=0;mhel0;--ix) { if(ix-1!=_ibar) ihel[ix]=(lhel%_constants[ix-1])/_constants[ix]; } (*ME())(ihel) = Complex(lepton[lhel].dot(hadron[mhel])*SM().fermiConstant()); } } // ckm factor double ckm(1.); if(inquark<=6) { if(inquark%2==0) ckm = SM().CKM(inquark/2-1,(abs(outquark)-1)/2); else ckm = SM().CKM(abs(outquark)/2-1,(inquark-1)/2); } // return the answer return 0.5*(ME()->contract(_rho)).real()*ckm; } // matrix element for a 1/2 -> 3/2 semi-leptonic decay double SemiLeptonicBaryonDecayer::halfThreeHalf(const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&part), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&part), incoming); // work out the mapping for the lepton vector _constants.resize(outgoing.size()+1); _ispin.resize(outgoing.size()); int itemp(1); _ibar=0; for(int ix=int(outgoing.size()-1);ix>=0;--ix) { _ispin[ix]=outgoing[ix]->iSpin(); if(abs(outgoing[ix]->id())<=16) { itemp*=_ispin[ix]; _constants[ix]=itemp; } else _ibar=ix; } _constants[outgoing.size()]=1; _constants[_ibar]=_constants[_ibar+1]; } if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,_ispin))); // spinors for the decay product LorentzPolarizationVector in=UnitRemoval::InvE*part.momentum(); if(part.id()>0) { RSSpinorBarWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalfBar.resize(4); _inHalfBar.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalfBar[ihel]=swave.dimensionedWf(); _inHalfBar[ihel] = _inThreeHalfBar[ihel].dot(in); } } else { RSSpinorWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalf.resize(4); _inHalf.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalf[ihel]=swave.dimensionedWf(); _inHalf[ihel] = _inThreeHalf[ihel].dot(in); } } // get the information on the form-factor int spinin(0),spinout(0),inquark,outquark,spect1,spect2; int id0(part.id()),id1(outgoing[0]->id()); bool cc; int iloc(_form->formFactorNumber(id0,id1,cc)); _form->formFactorInfo(iloc,spinin,spinout,spect1,spect2,inquark,outquark); // work out the value of q and calculate the form factors Lorentz5Momentum q(part.momentum()-momenta[0]); q.rescaleMass(); Energy m0(part.mass()),m1(momenta[0].mass()); Energy2 q2(q.mass2()); Lorentz5Momentum sum(part.momentum()+momenta[0]); // calculate the form factors Complex f1v,f2v,f3v,f4v,f1a,f2a,f3a,f4a; _form->SpinHalfSpinThreeHalfFormFactor(q2,iloc,id0,id1,m0,m1, - f1v,f2v,f3v,f4v,f1a,f2a,f3a,f4a); + f1v,f2v,f3v,f4v,f1a,f2a,f3a,f4a,FlavourInfo()); LorentzPolarizationVector vtemp; complex lS1,lS2,rS1,rS2; complex lV,rV; Complex left,right; InvEnergy ms(1./(m0+m1)); InvEnergy2 ms2(ms*ms); if(part.id()>0) { left = f1a-f1v; right = f1a+f1v; lS1 = ms2*(f3a-f4a-f3v+f4v); rS1 = ms2*(f3a-f4a+f3v-f4v); lS2 = ms2*(f4a-f4v); rS2 = ms2*(f4a+f4v); lV = ms*(f2a-f2v); rV = ms*(f2a+f2v); } else { left = conj(f1a+f1v); right = conj(f1a-f1v); lS1 = ms2*conj(f3a-f4a+f3v-f4v); rS1 = ms2*conj(f3a-f4a-f3v+f4v); lS2 = ms2*conj(f4a-f4v); rS2 = ms2*conj(f4a+f4v); lV = ms *conj(f2a-f2v); rV = ms *conj(f2a+f2v); } // calculate the hadronic current for the decay LorentzPolarizationVectorE hadron[4][2]; // construct the vectors for the decay Complex scalar1,scalar2; complex lfact,rfact; LorentzPolarizationVectorE tvec; LorentzPolarizationVector svec; for(unsigned int iya=0;iya<4;++iya) { for(unsigned int ixa=0;ixa<2;++ixa) { unsigned int ix=iya,iy=ixa; if(outgoing[0]->id()<0) swap(ix,iy); // scalar like terms lfact = _inHalf[iy].leftScalar(_inHalfBar[ix]); rfact = _inHalf[iy].rightScalar(_inHalfBar[ix]); scalar1 = Complex((lS1*lfact+rS1*rfact)*UnitRemoval::E); scalar2 = Complex((lS2*lfact+rS2*rfact)*UnitRemoval::E); svec = _inHalf[iy].generalCurrent(_inHalfBar[ix],lV/ms,rV/ms)*ms; if(part.id()>0) tvec=_inThreeHalfBar[ix].generalCurrent(_inHalf[iy],left,right); else tvec=_inThreeHalf[iy].generalCurrent(_inHalfBar[ix],left,right); hadron[iya][ixa] = tvec+svec*UnitRemoval::E+scalar1*momenta[0] +scalar2*part.momentum(); } } // construct the lepton current Energy scale; int mode((abs(outgoing[1]->id())-11)/12); vector lepton(_current->current(tcPDPtr(),FlavourInfo(),mode,-1,scale,tPDVector(outgoing.begin()+1,outgoing.end()), vector(momenta.begin()+1,momenta.end()),meopt)); vector ihel(outgoing.size()+1); for(unsigned int iya=0;iya<4;++iya) { ihel[1]=iya; for(unsigned int ixa=0;ixa<2;++ixa) { ihel[0]=ixa; for(unsigned int lhel=0;lhelcontract(_rho)).real()*ckm; } // output the setup information for the particle database void SemiLeptonicBaryonDecayer::dataBaseOutput(ofstream & output,bool header) const { if(header) output << "update decayers set parameters=\""; DecayIntegrator::dataBaseOutput(output,false); for(unsigned int ix=0;ix<_maxwgt.size();++ix) { output << "insert " << name() << ":MaximumWeight " << ix << " " << _maxwgt[ix] << " \n"; } _current->dataBaseOutput(output,false,true); output << "newdef " << name() << ":Current " << _current->name() << " \n"; _form->dataBaseOutput(output,false,true); output << "newdef " << name() << ":FormFactor " << _form->name() << " \n"; if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Decay/FormFactors/BaryonFormFactor.cc b/Decay/FormFactors/BaryonFormFactor.cc --- a/Decay/FormFactors/BaryonFormFactor.cc +++ b/Decay/FormFactors/BaryonFormFactor.cc @@ -1,163 +1,163 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the BaryonFormFactor class. // #include "BaryonFormFactor.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" namespace Herwig { using namespace ThePEG; void BaryonFormFactor::doinit() { Interfaced::doinit(); // check the consistency of the parameters unsigned int isize(_incomingid.size()); if(isize!=_outgoingid.size() || isize!=_incomingJ.size()||isize!=_outgoingJ.size()|| isize!=_spectator1.size()|| isize!=_spectator2.size()|| isize!=_inquark.size()|| isize!=_outquark.size()) {throw InitException() << "Inconsistent parameters in BaryonFormFactor::doinit() " << Exception::abortnow;} } void BaryonFormFactor::persistentOutput(PersistentOStream & os) const { os << _incomingid << _outgoingid << _incomingJ << _outgoingJ << _spectator1 << _spectator2 << _inquark << _outquark << _numbermodes; } void BaryonFormFactor::persistentInput(PersistentIStream & is, int) { is >> _incomingid >> _outgoingid >> _incomingJ >> _outgoingJ >> _spectator1 >> _spectator2 >> _inquark >> _outquark >> _numbermodes; } // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractClass describeHerwigBaryonFormFactor("Herwig::BaryonFormFactor", "Herwig.so"); void BaryonFormFactor::Init() { static ClassDocumentation documentation ("The BaryonFormFactor class is the base class for" " the implementation of the form factors for weak decays of baryon."); static ParVector interfaceIncoming ("Incoming", "The id of the incoming baryons.", &BaryonFormFactor::_incomingid, 0, 0, 0, 0, 1000000, false, false, true); static ParVector interfaceOutgoing ("Outgoing", "The id of the outgoing baryons.", &BaryonFormFactor::_outgoingid, 0, 0, 0, 0, 1000000, false, false, true); static ParVector interfaceInSpin ("InSpin", "The spin of the incoming baryons to decide which form factors to use.", &BaryonFormFactor::_incomingJ, 0, 0, 0, 0, 2, false, false, true); static ParVector interfaceOutSpin ("OutSpin", "The spin of the outgoing baryons to decide which form factors to use.", &BaryonFormFactor::_outgoingJ, 0, 0, 0, 0, 4, false, false, true); static ParVector interfaceSpectator1 ("Spectator1", "The first specator quark.", &BaryonFormFactor::_spectator1, 0, 0, 0, 0, 5, false, false, true); static ParVector interfaceSpectator2 ("Spectator2", "The second specator quark.", &BaryonFormFactor::_spectator2, 0, 0, 0, 0, 5, false, false, true); static ParVector interfaceInQuark ("InQuark", "The PDG code for the decaying quark.", &BaryonFormFactor::_inquark, 0, 0, 0, 0, 5, false, false, true); static ParVector interfaceOutQuark ("OutQuark", "The PDG code for the quark produced in the decay.", &BaryonFormFactor::_outquark, 0, 0, 0, 0, 5, false, false, true); } // form factor for spin-1/2 to spin-1/2 void BaryonFormFactor::SpinHalfSpinHalfFormFactor(Energy2,int,int,int,Energy,Energy, Complex &,Complex &,Complex &, Complex &,Complex &,Complex &, - Virtuality) { + FlavourInfo,Virtuality) { throw Exception() << "Error in BaryonFormFactor::SpinHalfSpinHalfFormFactor" - << " not implemented" + << " not implemented" << fullName() << Exception::abortnow; } void BaryonFormFactor::SpinHalfSpinThreeHalfFormFactor(Energy2,int,int,int,Energy, Energy,Complex &,Complex &, Complex &,Complex &,Complex &, Complex &,Complex &,Complex &, - Virtuality) { + FlavourInfo,Virtuality) { throw Exception() << "Error in BaryonFormFactor::SpinHalfSpinThreeHalfFormFactor" - << " not implemented" + << " not implemented" << fullName() << Exception::abortnow; } // output the information for the database void BaryonFormFactor::dataBaseOutput(ofstream & output,bool header,bool create) const { if(header) output << "update decayers set parameters=\""; if(create) output << "create Herwig::BaryonFormFactor " << name() << " \n"; for(unsigned int ix=0;ix<_incomingid.size();++ix) { if(ix<_numbermodes) { output << "newdef " << name() << ":Incoming " << ix << " " << _incomingid[ix] << endl; output << "newdef " << name() << ":Outgoing " << ix << " " << _outgoingid[ix] << endl; output << "newdef " << name() << ":InSpin " << ix << " " << _incomingJ[ix] << endl; output << "newdef " << name() << ":OutSpin " << ix << " " << _outgoingJ[ix] << endl; output << "newdef " << name() << ":Spectator1 " << ix << " " << _spectator1[ix] << endl; output << "newdef " << name() << ":Spectator2 " << ix << " " << _spectator2[ix] << endl; output << "newdef " << name() << ":InQuark " << ix << " " << _inquark[ix] << endl; output << "newdef " << name() << ":OutQuark " << ix << " " << _outquark[ix]<< endl; } else { output << "insert " << name() << ":Incoming " << ix << " " << _incomingid[ix] << endl; output << "insert " << name() << ":Outgoing " << ix << " " << _outgoingid[ix] << endl; output << "insert " << name() << ":InSpin " << ix << " " << _incomingJ[ix] << endl; output << "insert " << name() << ":OutSpin " << ix << " " << _outgoingJ[ix] << endl; output << "insert " << name() << ":Spectator1 " << ix << " " << _spectator1[ix] << endl; output << "insert " << name() << ":Spectator2 " << ix << " " << _spectator2[ix] << endl; output << "insert " << name() << ":InQuark " << ix << " " << _inquark[ix] << endl; output << "insert " << name() << ":OutQuark " << ix << " " << _outquark[ix]<< endl; } } if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } } diff --git a/Decay/FormFactors/BaryonFormFactor.h b/Decay/FormFactors/BaryonFormFactor.h --- a/Decay/FormFactors/BaryonFormFactor.h +++ b/Decay/FormFactors/BaryonFormFactor.h @@ -1,406 +1,408 @@ // -*- C++ -*- #ifndef HERWIG_BaryonFormFactor_H #define HERWIG_BaryonFormFactor_H // // This is the declaration of the BaryonFormFactor class. // #include "ThePEG/Interface/Interfaced.h" #include "BaryonFormFactor.fh" #include "ThePEG/Config/Complex.h" #include "Herwig/Decay/IsoSpin.h" namespace Herwig { using namespace ThePEG; /** \ingroup Decay * The BaryonFormFactor class is the base class for the implementation * of the form-factors for the weak decay of a baryon. * * This is designed so that the form factors can be used for both the semi-leptonic * decays and using factorization for hadronic weak decays. * * The form factors are given below for the decay \f$X(p_0)\to Y(p_1)\f$ with * \f$q_\mu=(p_0-p_1)_\mu\f$. * * * The form-factors are defined to be * * \f[\bar{u}(p_1) \left[ \gamma^\mu \left(F^V_1+F^A_1 \gamma_5\right) * +\frac{i}{(m_0+m_1)}\sigma_{\mu\nu}q^\nu\left(F^V_2+F^A_2\gamma_5\right) * +\frac1{(m_0+m_1)}q^\mu\left(F^V_3+F^A_3\gamma_5\right)\right] u(p_0) \f] * * for the \f$\frac12\to\frac12\f$ transition and * * \f[\bar{u}^\alpha(p_1) \left[ g_{\alpha\mu}\left(G^V_1+G^A_1 \gamma_5\right) * +\frac1{(m_0+m_1)}p_{0\alpha}\gamma_\mu\left(G^V_2+G^A_2\gamma_5\right) * +\frac1{(m_0+m_1)^2}p_{0\alpha}p_{1\mu}\left(G^V_3+G^A_3\gamma_5\right)\right.\f] * \f[ \left. * +\frac1{(m_0+m_1)^2}p_{0\alpha}q_\mu\left(G^V_4+G^A_4\gamma_5\right)\right] * \gamma_5 u(p_0) * \f] * for the \f$\frac12\to\frac32\f$ transition. * * These definitions differ from those in the liturature because we have divided some * terms by the sum of the baryon masses to ensure that the form-factors are all * dimensionless. * * In many applications, particularly for the decay of baryons containing a heavy * quark, an alternative version of the form factors in terms of the velocities * of the baryons is used. This form is * * \f[\bar{u}(v') \left[ \gamma^\mu \left(F_1-G_1 \gamma_5\right) * v^\mu \left(F_2-G_2 \gamma_5\right) * {v'}^\mu \left(F_3-G_3\gamma_5\right)\right] u(v) \f] * * for the \f$\frac12\to\frac12\f$ transition and * * \f[\bar{u}^\alpha(v') \left[ v_\alpha\gamma_\mu\left(N_1-K_1 \gamma_5\right) * +v_\alpha v^\mu \left(N_2-K_2 \gamma_5\right) * +v_\alpha {v'}^\mu \left(N_3-K_3 \gamma_5\right) * +g_\alpha^\mu \left(N_4-K_4 \gamma_5\right)\right] * \gamma_5 u(v) * \f] * for the \f$\frac12\to\frac32\f$ transition. * * In terms of these form factors the form factors we use are * * \f[F^V_1= F_1+\frac12(m_0+m_1)\left(\frac{F_2}{m_0}+\frac{F_3}{m_1}\right)\f] * \f[F^V_2=\frac12(m_0+m_1)\left(\frac{F_2}{m_0}+\frac{F_3}{m_1}\right)\f] * \f[F^V_3=\frac12(m_0+m_1)\left(\frac{F_2}{m_0}-\frac{F_3}{m_1}\right)\f] * \f[F^A_1=-G_1+\frac12(m_0-m_1)\left(\frac{G_2}{m_0}+\frac{G_3}{m_1}\right)\f] * \f[F^A_2=-\frac12(m_0+m_1)\left(\frac{G_2}{m_0}+\frac{G_3}{m_1}\right)\f] * \f[F^A_3=\frac12(m_0+m_1)\left(-\frac{G_2}{m_0}+\frac{G_3}{m_1}\right)\f] * * for the \f$\frac12\to\frac12\f$ transition and * * \f[G_1^V = N_4\f] * \f[G_2^V = N_1\frac{(m_0+m_1)}{m_0}\f] * \f[G_3^V = \frac{(m_0+m_1)^2}{m_0}\left(\frac{N_2}{m_0}+\frac{N_3}{m_1}\right)\f] * \f[G_4^V = N_2\frac{(m_0+m_1)^2}{m^2_0}\f] * \f[G_1^A =-K_4\f] * \f[G_2^A =-K_1\frac{(m_0+m_1)}{m_0}\f] * \f[G_3^A =-\frac{(m_0+m_1)^2}{m_0}\left(\frac{K_2}{m_0}+\frac{K_3}{m_1}\right)\f] * \f[G_4^A =-K_2\frac{(m_0+m_1)^2}{m^2_0}\f] * * for the \f$\frac12\to\frac32\f$ transition. * * * @see ScalarFormFactor * */ class BaryonFormFactor: public Interfaced { public: /** * Whether the form factor is space- or time-like */ enum Virtuality {TimeLike, SpaceLike}; public: /** * Default constructor */ BaryonFormFactor() : _numbermodes() {} /** @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); //@} /** * Standard Init function used to initialize the interfaces. */ 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 {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);} //@} public: /** @name Functions to give information about the form factors available. */ //@{ /** * Find the location for a given pair of particle. * \param in PDG code for the incoming baryon. * \param out PDG code for the outgoing baryon. * \param cc particles or charge conjugates stored in form factor. * @return The location in the vectors storing the data. */ int formFactorNumber(int in,int out,bool & cc) const { int output(-1); unsigned int ix(0); if(_incomingid.size()==0){return output;} do { if(_incomingid[ix]== in && _outgoingid[ix]== out) { cc=false; output=ix; } else if (_incomingid[ix]==-in && _outgoingid[ix]==-out) { cc=true; output=ix; } ++ix; } while(ix<_incomingid.size()&&output<0); return output; } /** * Get the particle ids for an entry. * @param iloc The location in the list. * @param id0 The PDG code for the incoming baryon. * @param id1 The PDG code for the outgoing baryon. */ void particleID(int iloc,int& id0,int& id1) { id0=_incomingid[iloc]; id1=_outgoingid[iloc]; } /** * Information on the form factor. * @param iloc The location in the list. * @param ispin The spin of the incoming baryon. * @param ospin The spin of the outgoing baryon. * @param spect1 The PDG code of the first spectator quark. * @param spect2 The PDG code of the second spectator quark. * @param inquark The PDG code for decaying incoming quark. * @param outquark The PDG code for the outgoing quark produced in the decay. */ void formFactorInfo(int iloc,int & ispin,int & ospin,int & spect1, int & spect2, int & inquark,int & outquark) { ispin = _incomingJ[iloc]; ospin = _outgoingJ[iloc]; spect1 = _spectator1[iloc]; spect2 = _spectator2[iloc]; inquark = _inquark[iloc]; outquark = _outquark[iloc]; } /** * Information on the form factor. * @param in The PDG code of the incoming baryon. * @param out The PDG code of the outgoing baryon. * @param ispin The spin of the incoming baryon. * @param ospin The spin of the outgoing baryon. * @param spect1 The PDG code of the first spectator quark. * @param spect2 The PDG code of the second spectator quark. * @param inquark The PDG code for decaying incoming quark. * @param outquark The PDG code for the outgoing quark produced in the decay. */ void formFactorInfo(int in,int out,int & ispin,int & ospin,int & spect1, int & spect2, int & inquark,int & outquark) { bool dummy; unsigned int ix=formFactorNumber(in,out,dummy); formFactorInfo(ix,ispin,ospin,spect1,spect2,inquark,outquark); } /** * number of form factors */ unsigned int numberOfFactors() const {return _incomingid.size();} //@} public: /** @name Form Factors */ //@{ /** * The form factor for the weak decay of a spin \f$\frac12\f$ baryon to a * spin \f$\frac12\f$ baryon. * This method is virtual and must be implementented in classes * inheriting from this which include spin\f$\frac12\f$ to spin \f$\frac12\f$ * form factors. * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form factor list. * @param id0 The PDG code of the incoming baryon. * @param id1 The PDG code of the outgoing baryon. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param f1v The form factor \f$F^V_1\f$. * @param f2v The form factor \f$F^V_2\f$. * @param f3v The form factor \f$F^V_3\f$. * @param f1a The form factor \f$F^A_1\f$. * @param f2a The form factor \f$F^A_2\f$. * @param f3a The form factor \f$F^A_3\f$. * @param virt Whether the q2 is space or timelike */ virtual void SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo flavour, Virtuality virt=SpaceLike); /** * The form factor for the weak decay of a spin \f$\frac12\f$ baryon to a * spin \f$\frac32\f$ baryon. * This method is virtual and must be implementented in classes * inheriting from this which include spin\f$\frac12\f$ to spin \f$\frac32\f$ * form factors. * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form factor list. * @param id0 The PDG code of the incoming baryon. * @param id1 The PDG code of the outgoing baryon. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param g1v The form factor \f$G^V_1\f$. * @param g2v The form factor \f$G^V_2\f$. * @param g3v The form factor \f$G^V_3\f$. * @param g4v The form factor \f$G^V_4\f$. * @param g1a The form factor \f$G^A_1\f$. * @param g2a The form factor \f$G^A_2\f$. * @param g3a The form factor \f$G^A_3\f$. * @param g4a The form factor \f$G^A_4\f$. * @param virt Whether the q2 is space or timelike */ virtual void SpinHalfSpinThreeHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & g1v,Complex & g2v,Complex & g3v, Complex & g4v,Complex & g1a,Complex & g2a, Complex & g3a,Complex & g4a, + FlavourInfo flavour, Virtuality virt=SpaceLike); //@} /** * 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 * @param create Whether or not to add a statement creating the object */ virtual void dataBaseOutput(ofstream & os,bool header,bool create) const; protected: /** * Add a form factor to the list. * @param in The PDG code of the incoming baryon. * @param out The PDG code of the outgoing baryon. * @param inspin The spin of the incoming baryon. * @param outspin The spin of the outgoing baryon. * @param spect1 The PDG code of the first spectator quark. * @param spect2 The PDG code of the second spectator quark. * @param inquark The PDG code for decaying incoming quark. * @param outquark The PDG code for the outgoing quark produced in the decay. */ void addFormFactor(int in,int out,int inspin, int outspin, int spect1, int spect2, int inquark,int outquark) { _incomingid.push_back(in); _outgoingid.push_back(out); _incomingJ.push_back(inspin); _outgoingJ.push_back(outspin); _spectator1.push_back(spect1); _spectator2.push_back(spect2); _inquark.push_back(inquark); _outquark.push_back(outquark); } /** * Set initial number of modes * @param nmodes The number of modes. */ void initialModes(unsigned int nmodes) {_numbermodes=nmodes;} /** * Get the initial number of modes */ unsigned int initialModes() const {return _numbermodes;} 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(); //@} private: /** * Private and non-existent assignment operator. */ BaryonFormFactor & operator=(const BaryonFormFactor &) = delete; private: /** * the id's of the incoming particles */ vector _incomingid; /** * the id's of the outgoing particles */ vector _outgoingid; /** * spin of the incoming particle */ vector _incomingJ; /** * spin of the outgoing particle */ vector _outgoingJ; /** * the id of the first spectator quark */ vector _spectator1; /** * the id of the second spectator quark */ vector _spectator2; /** * the id of the decaying quark */ vector _inquark; /** * the id of the outgoing quark */ vector _outquark; /** * The initial number of modes */ unsigned int _numbermodes; }; } #endif /* HERWIG_BaryonFormFactor_H */ diff --git a/Decay/FormFactors/BaryonSimpleFormFactor.cc b/Decay/FormFactors/BaryonSimpleFormFactor.cc --- a/Decay/FormFactors/BaryonSimpleFormFactor.cc +++ b/Decay/FormFactors/BaryonSimpleFormFactor.cc @@ -1,240 +1,241 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the BaryonSimpleFormFactor class. // #include "BaryonSimpleFormFactor.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; using namespace ThePEG; BaryonSimpleFormFactor::BaryonSimpleFormFactor() { // axial vector coupling in beta decay _gA=1.25; // SU(3) breaking paramters // D/(D+F) ratio _alphaD=0.6; // eta_V parameter _etaV=0.970; // eta_A parameter _etaA=1.080; // SU(3) breaking factors for the electric dipole moment _rhoE=0.094; // SU(3) breaking factors for the magnetic dipole moment _rhoM=0.860; // the various decay modes handled by the model addFormFactor(2112,2212,2,2,2,1,1,2); addFormFactor(3222,3122,2,2,3,2,1,2); addFormFactor(3112,3122,2,2,3,1,1,2); addFormFactor(3112,3212,2,2,3,1,1,2); addFormFactor(3212,3222,2,2,3,2,1,2); addFormFactor(3312,3322,2,2,3,3,1,2); addFormFactor(3122,2212,2,2,2,1,3,2); addFormFactor(3212,2212,2,2,2,1,3,2); addFormFactor(3112,2112,2,2,1,1,3,2); addFormFactor(3312,3122,2,2,3,1,3,2); addFormFactor(3312,3212,2,2,3,1,3,2); addFormFactor(3322,3222,2,2,3,2,3,2); // set the inital number of form factors initialModes(numberOfFactors()); } void BaryonSimpleFormFactor::doinit() { BaryonFormFactor::doinit(); _f1.clear(); _f2.clear(); _g1.clear(); _g2.clear(); _f1.resize(numberOfFactors()); _f2.resize(numberOfFactors()); _g1.resize(numberOfFactors()); _g2.resize(numberOfFactors()); // calculate the couplings for the different modes int id0,id1; double root23(sqrt(2./3.)),root2(sqrt(2.)),root32(sqrt(3./2.)); for(unsigned int ix=0;ix> _gA >> _alphaD >> _etaV >> _etaA >> _rhoE >> _rhoM >> _f1 >> _f2 >> _g1 >> _g2;} // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigBaryonSimpleFormFactor("Herwig::BaryonSimpleFormFactor", "HwFormFactors.so"); void BaryonSimpleFormFactor::Init() { static ClassDocumentation documentation ("The BaryonSimpleFormFactor class implements the" " quark model calculation of the form-factors from PRD25, 206", "The BaryonSimpleFormFactor class which implements the results" " of \\cite{Donoghue:1981uk}was used for the weak decays of the light baryons.", "\\bibitem{Donoghue:1981uk}\n" "J.~F.~Donoghue and B.~R.~Holstein,\n" "Phys.\\ Rev.\\ D {\\bf 25} (1982) 206.\n" "%%CITATION = PHRVA,D25,206;%%\n"); static Parameter interfaceg_A ("g_A", "The axial-vector coupling in neutron beta decay.", &BaryonSimpleFormFactor::_gA, 1.25, 0.0, 10.0, false, false, true); static Parameter interfacealpha_D ("alpha_D", "SU(3) breaking parameter which is the ratio D/(D+F). ", &BaryonSimpleFormFactor::_alphaD, 0.6, 0.0, 1.0, false, false, true); static Parameter interfaceeta_V ("eta_V", "The eta_V SU(3) breaking parameter", &BaryonSimpleFormFactor::_etaV, .97, 0.0, 10.0, false, false, true); static Parameter interfaceeta_A ("eta_A", "The eta_A SU(3) breaking parameter", &BaryonSimpleFormFactor::_etaA, 1.08, 0.0, 10.0, false, false, true); static Parameter interfacerho_E ("rho_E", "The SU(3) breaking parameter for the electric dipole moment.", &BaryonSimpleFormFactor::_rhoE, 0.094, 0.0, 10.0, false, false, true); static Parameter interfacerho_M ("rho_M", "The SU(3) breaking parameter for the magentic dipole moment.", &BaryonSimpleFormFactor::_rhoM, 0.860, 0.0, 10.0, false, false, true); } void BaryonSimpleFormFactor:: SpinHalfSpinHalfFormFactor(Energy2,int iloc, int,int,Energy,Energy, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo, Virtuality virt) { assert(virt==SpaceLike); useMe(); f1v = _f1[iloc]; f1a = -_g1[iloc]; f2v = -_f2[iloc]; f2a = _g2[iloc]; f3v = 0.; f3a = 0.; } void BaryonSimpleFormFactor::dataBaseOutput(ofstream& output,bool header, bool create) const { if(header) output << "update decayers set parameters=\""; if(create) output << "create Herwig::BaryonSimpleFormFactor " << name() << " \n"; output << "newdef " << name() << ":g_A " << _gA << " \n"; output << "newdef " << name() << ":alpha_D " << _alphaD << " \n"; output << "newdef " << name() << ":eta_V " << _etaV << " \n"; output << "newdef " << name() << ":eta_A " << _etaA << " \n"; output << "newdef " << name() << ":rho_E " << _rhoE << " \n"; output << "newdef " << name() << ":rho_M " << _rhoM << " \n"; BaryonFormFactor::dataBaseOutput(output,false,false); if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Decay/FormFactors/BaryonSimpleFormFactor.h b/Decay/FormFactors/BaryonSimpleFormFactor.h --- a/Decay/FormFactors/BaryonSimpleFormFactor.h +++ b/Decay/FormFactors/BaryonSimpleFormFactor.h @@ -1,180 +1,181 @@ // -*- C++ -*- #ifndef HERWIG_BaryonSimpleFormFactor_H #define HERWIG_BaryonSimpleFormFactor_H // // This is the declaration of the BaryonSimpleFormFactor class. // #include "BaryonFormFactor.h" namespace Herwig { using namespace ThePEG; /** \ingroup Decay * * The BaryonSimpleFormFactor class is a simple model for the form-factors * for the semi-leptonic decay of the light (i.e. uds) baryons. The form-factors * are assumed to be constant and are taken from the quark model results * of PRD25, 206 (1982). * * @ see BaryonFormFactor */ class BaryonSimpleFormFactor: public BaryonFormFactor { public: /** * Default constructor */ BaryonSimpleFormFactor(); /** @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); //@} /** * Standard Init function used to initialize the interfaces. */ static void Init(); public: /** @name Form Factors */ //@{ /** * The form factor for the weak decay of a spin \f$\frac12\f$ baryon to a * spin \f$\frac12\f$ baryon. * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form factor list. * @param id0 The PDG code of the incoming baryon. * @param id1 The PDG code of the outgoing baryon. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param f1v The form factor \f$F^V_1\f$. * @param f2v The form factor \f$F^V_2\f$. * @param f3v The form factor \f$F^V_3\f$. * @param f1a The form factor \f$F^A_1\f$. * @param f2a The form factor \f$F^A_2\f$. * @param f3a The form factor \f$F^A_3\f$. * @param virt Whether the q2 is space or timelike */ virtual void SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo flavour, Virtuality virt=SpaceLike); //@} /** * 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 * @param create Whether or not to add a statement creating the object */ virtual void dataBaseOutput(ofstream & os,bool header,bool create) 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 and * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * Private and non-existent assignment operator. */ BaryonSimpleFormFactor & operator=(const BaryonSimpleFormFactor &) = delete; private: /** * The axial vector coupling in \f$\beta\f$ decay. */ double _gA; /** * The \f$\frac{D}{D+F}\f$ ratio. */ double _alphaD; /** * The \f$\eta_V\f$ \f$SU(3)\f$ breaking parameter. */ double _etaV; /** * The \f$\eta_V\f$ \f$SU(3)\f$ breaking parameter. */ double _etaA; /** * The \f$\rho_E\f$ \f$SU(3)\f$ breaking parameter for the electic dipole moment. */ double _rhoE; /** * The \f$\rho_M\f$ \f$SU(3)\f$ breaking parameter for the electic dipole moment. */ double _rhoM; /** * The calculated constants for the \f$f_1\f$ form factor. */ vector _f1; /** * The calculated constants for the \f$f_2\f$ form factor. */ vector _f2; /** * The calculated constants for the \f$g_1\f$ form factor. */ vector _g1; /** * The calculated constants for the \f$g_2\f$ form factor. */ vector _g2; }; } #endif /* HERWIG_BaryonSimpleFormFactor_H */ diff --git a/Decay/FormFactors/BaryonThreeQuarkModelFormFactor.cc b/Decay/FormFactors/BaryonThreeQuarkModelFormFactor.cc --- a/Decay/FormFactors/BaryonThreeQuarkModelFormFactor.cc +++ b/Decay/FormFactors/BaryonThreeQuarkModelFormFactor.cc @@ -1,652 +1,654 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the BaryonThreeQuarkModelFormFactor class. // #include "BaryonThreeQuarkModelFormFactor.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Utilities/GaussianIntegrator.h" #include "ThePEG/StandardModel/StandardModelBase.h" using namespace Herwig; // function for the integral of the partial width namespace { using namespace Herwig; struct BaryonCMatrixElement { typedef ThePEG::Ptr::pointer BaryonThreeQuarkModelFormFactorPtr; BaryonCMatrixElement(BaryonThreeQuarkModelFormFactorPtr in, Energy m0, Energy m1, int type, int mass, int id0,int id1) { _formFactor=in; _m0=m0; _m1=m1; _type=type; _mass=mass; _id0=id0; _id1=id1; } // calculate the integrand Energy operator() (double x) const { return _formFactor->widthIntegrand(x,_m0,_m1,_type,_mass,_id0,_id1); } /** Argument type for GaussianIntegrator */ typedef double ArgType; /** Return type for GaussianIntegrator */ typedef Energy ValType; private: // private variables BaryonThreeQuarkModelFormFactorPtr _formFactor; Energy _m0,_m1; int _type,_mass; int _id0,_id1; }; } BaryonThreeQuarkModelFormFactor::BaryonThreeQuarkModelFormFactor() : _initialize(false), _order(50),_mlight(420*MeV),_mstrange(570*MeV), _LambdaQ(2.5*GeV),_Lambdaqq(0.71*GeV),_Lambdasq(850*MeV),_Lambdass(1.0*GeV) { // modes handled by this form factor // lambda_b addFormFactor(5122,4122,2,2,1,2,5,4); // xi_b addFormFactor(5232,4232,2,2,2,3,5,4); addFormFactor(5132,4132,2,2,1,3,5,4); // sigma_b addFormFactor(5112,4112,2,2,1,1,5,4); addFormFactor(5212,4212,2,2,2,1,5,4); addFormFactor(5222,4222,2,2,2,2,5,4); // omega_b addFormFactor(5332,4332,2,2,3,3,5,4); // sigma_b-> sigma_c* addFormFactor(5112,4114,2,4,1,1,5,4); addFormFactor(5212,4214,2,4,2,1,5,4); addFormFactor(5222,4224,2,4,2,2,5,4); // omega_b -> omega_c* addFormFactor(5332,4334,2,4,3,3,5,4); // set the inital number of form factors initialModes(numberOfFactors()); } void BaryonThreeQuarkModelFormFactor::doinit() { BaryonFormFactor::doinit(); // initialization in needed if(_initialize) { GaussianIntegrator integrator; _C0.clear();_C1.clear();_C2.clear(); double pre(0.),root(2.*sqrt(6.)); double gamma1(1),gamma2(1),gamma3(sqrt(acos(-1.))); unsigned int ix,iy; for(iy=0;iy<2;++iy) { _mu2 = iy==0 ? sqr(_mlight /_LambdaQ) : sqr(_mstrange/_LambdaQ); for(ix=0;ix<=_order;++ix) { if(ix>0) gamma1*=ix; if(ix%2==1) { gamma2*=(ix+1)/2.0; gamma3*=ix/2.0; } if(ix%2==0) pre=pow(root,double(ix))/12.*gamma2/gamma1; else pre=pow(root,double(ix))/12.*gamma3/gamma1; // for the xi_0 function _a=_mu2;_b=2.;_N=ix; _C0.push_back(pre*integrator.value(*this,0.,1.)); // for the xi_1 function _a=_mu2;_b=1.; _C1.push_back(pre*integrator.value(*this,0.,1.)); // for the xi_2 function _a=0.; _b=0.; _C2.push_back(pre*integrator.value(*this,0.,1.)); } } generator()->log() << "Checking results of BaryonThreeQuarkModelFormFactor" << "vs results from the original paper\n"; // first matrix element Energy m0=getParticleData(5122)->mass(); Energy m1=getParticleData(4122)->mass(); double omegamax=0.5*(m0*m0+m1*m1)/m0/m1; BaryonCMatrixElement int1(this,m0,m1,1,0,5122,4122); Energy width = integrator.value(int1,1.,omegamax); generator()->log() << "Lambda_b0->Lambda_c+ decay" << width/6.582119E-22/MeV << "\n"; // second matrix element m0=getParticleData(5222)->mass(); m1=getParticleData(4222)->mass(); omegamax=0.5*(m0*m0+m1*m1)/m0/m1; int1 = BaryonCMatrixElement(this,m0,m1,2,0,5222,4222); width = integrator.value(int1,1.,omegamax); generator()->log() << "Sigma_b+->Sigma_c++ decay" << width/6.582119E-22/MeV << "\n"; // third matrix element m0=getParticleData(5232)->mass(); m1=getParticleData(4232)->mass(); omegamax=0.5*(m0*m0+m1*m1)/m0/m1; int1 = BaryonCMatrixElement(this,m0,m1,1,1,5232,4232); width = integrator.value(int1,1.,omegamax); generator()->log() << "Xi_b0->Xi_c+ decay" << width/6.582119E-22/MeV << "\n"; // fourth matrix element m0=getParticleData(5332)->mass(); m1=getParticleData(4332)->mass(); omegamax=0.5*(m0*m0+m1*m1)/m0/m1; int1 = BaryonCMatrixElement(this,m0,m1,2,2,5332,4332); width = integrator.value(int1,1.,omegamax); generator()->log() << "Omega_b-->Omega_c0 decay" << width/6.582119E-22/MeV << "\n"; // fifth matrix element m0=getParticleData(5222)->mass(); m1=getParticleData(4224)->mass(); omegamax=0.5*(m0*m0+m1*m1)/m0/m1; int1 = BaryonCMatrixElement(this,m0,m1,3,0,5222,4224); width = integrator.value(int1,1.,omegamax); generator()->log() << "Sigma_vb+->Sigma_c*++ decay" << width/6.582119E-22/MeV << "\n"; // fourth matrix element m0=getParticleData(5332)->mass(); m1=getParticleData(4334)->mass(); omegamax=0.5*(m0*m0+m1*m1)/m0/m1; int1 = BaryonCMatrixElement(this,m0,m1,3,2,5332,4334); width = integrator.value(int1,1.,omegamax); generator()->log() << "Omega_b-->Omega_c*0 decay" << width/6.582119E-22/MeV << "\n"; // output some plots for testing double lambdabar = -999.999; ofstream output("ThreeQuark.top"); output << "newdef font duplex" << endl; output << "title top \"Figure 3 from paper \"" << endl; output << "newdef limits x 1 1.44 y 0.5 1" << endl; for(unsigned int ix=0;ix<5;++ix) { double omegamin(1.),omegamax(1.44), step((omegamax-omegamin)/100.),omega(1.),xi; unsigned int ioff(0); if(ix==0){lambdabar=600*MeV/_LambdaQ;} else if(ix==1){lambdabar=650*MeV/_LambdaQ;} else if(ix==2){lambdabar=710*MeV/_LambdaQ;} else if(ix==3){lambdabar=750*MeV/_LambdaQ;} else if(ix==4){lambdabar=800*MeV/_LambdaQ;} for(;omega phi(phiFunction(omega)); double power(1.),numer(0.),denom(0.); for(unsigned int iy=0;iy<=_order;++iy) { numer+=phi[iy]*power*_C0[iy+ioff]; denom+=power*_C0[iy+ioff]; power*=lambdabar; } xi=numer/denom; output << omega << " " << xi << endl; } if(ix==0){output << "join red" << endl;} else if(ix==1){output << "join green" << endl;} else if(ix==2){output << "join blue" << endl;} else if(ix==3){output << "join cyan" << endl;} else if(ix==4){output << "join magenta" << endl;} } output << "new frame " << endl; output << "newdef font duplex" << endl; output << "title top \"Figure 6 from paper \"" << endl; output << "newdef limits x 1 1.4 y 0.5 1" << endl; for(unsigned int ix=0;ix<5;++ix) { double omegamin(1.),omegamax(1.45),step((omegamax-omegamin)/100.),omega(1.); unsigned int ioff(0); if(ix==0){lambdabar=600*MeV/_LambdaQ;} else if(ix==1){lambdabar=650*MeV/_LambdaQ;} else if(ix==2){lambdabar=710*MeV/_LambdaQ;} else if(ix==3){lambdabar=750*MeV/_LambdaQ;} else if(ix==4){lambdabar=800*MeV/_LambdaQ;} for(;omega phi(phiFunction(omega)); double power(1.),numer[2]={0.,0.},denom(0.); for(unsigned int iy=0;iy<=_order;++iy) { numer[0]+=phi[iy]*power*_C1[iy+ioff]; denom+=power*_C1[iy+ioff]; numer[1]+=power*_C2[iy+ioff]*(phi[iy]-phi[iy+2]); power*=lambdabar; } numer[1]/=(omega-1.); double xi1(numer[0]/denom); output << omega << " " << xi1 << endl; } if(ix==0){output << "join red" << endl;} else if(ix==1){output << "join green" << endl;} else if(ix==2){output << "join blue" << endl;} else if(ix==3){output << "join cyan" << endl;} else if(ix==4){output << "join magenta" << endl;} } output << "new frame " << endl; output << "newdef font duplex" << endl; output << "title top \"Figure 7 from paper \"" << endl; output << "newdef limits x 1 1.33 y 0.4 1" << endl; for(unsigned int ix=0;ix<5;++ix) { double omegamin(1.),omegamax(1.45),step((omegamax-omegamin)/100.),omega(1.); unsigned int ioff(_order+1); if(ix==0){lambdabar=800*MeV/_LambdaQ;} else if(ix==1){lambdabar=900*MeV/_LambdaQ;} else if(ix==2){lambdabar=1000*MeV/_LambdaQ;} else if(ix==3){lambdabar=1050*MeV/_LambdaQ;} else if(ix==4){lambdabar=1100*MeV/_LambdaQ;} for(;omega phi(phiFunction(omega)); double power(1.),numer[2]={0.,0.},denom(0.); for(unsigned int iy=0;iy<=_order;++iy) { numer[0]+=phi[iy]*power*_C1[iy+ioff]; denom+=power*_C1[iy+ioff]; numer[1]+=power*_C2[iy+ioff]*(phi[iy]-phi[iy+2]); power*=lambdabar; } numer[1]/=(omega-1.); double xi1(numer[0]/denom); output << omega << " " << xi1 << endl; } if(ix==0){output << "join red" << endl;} else if(ix==1){output << "join green" << endl;} else if(ix==2){output << "join blue" << endl;} else if(ix==3){output << "join cyan" << endl;} else if(ix==4){output << "join magenta" << endl;} } } } void BaryonThreeQuarkModelFormFactor::persistentOutput(PersistentOStream & os) const { os << _initialize << _order << ounit(_mlight,MeV) << ounit(_mstrange,MeV) << ounit(_LambdaQ,MeV) << ounit(_Lambdaqq,MeV) << ounit(_Lambdasq,MeV) << ounit(_Lambdass,MeV) << _C0 << _C1 << _C2;} void BaryonThreeQuarkModelFormFactor::persistentInput(PersistentIStream & is, int) { is >> _initialize >> _order >> iunit(_mlight,MeV) >> iunit(_mstrange,MeV) >> iunit(_LambdaQ,MeV) >> iunit(_Lambdaqq,MeV) >> iunit(_Lambdasq,MeV) >> iunit(_Lambdass,MeV) >> _C0 >> _C1 >> _C2;} // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigBaryonThreeQuarkModelFormFactor("Herwig::BaryonThreeQuarkModelFormFactor", "HwFormFactors.so"); void BaryonThreeQuarkModelFormFactor::Init() { static ClassDocumentation documentation ("The BaryonThreeQuarkModelFormFactor class implements" " the form-factors for semi-leptonic decay of baryon containing a" " heavy quark from PRD56, 348.", "The form factors from \\cite{Ivanov:1996fj} were used.", "%\\cite{Ivanov:1996fj}\n" "\\bibitem{Ivanov:1996fj}\n" " M.~A.~Ivanov, V.~E.~Lyubovitskij, J.~G.~Korner and P.~Kroll,\n" " ``Heavy baryon transitions in a relativistic three-quark model,''\n" " Phys.\\ Rev.\\ D {\\bf 56} (1997) 348\n" " [arXiv:hep-ph/9612463].\n" " %%CITATION = PHRVA,D56,348;%%\n" ); static Parameter interfaceOrder ("Order", "The order of terms to include in the series expansion of the form-factor.", &BaryonThreeQuarkModelFormFactor::_order, 10, 0, 1000, false, false, true); static Parameter interfaceLightMass ("LightMass", "The mass of the light quark", &BaryonThreeQuarkModelFormFactor::_mlight, GeV, .42*GeV, ZERO, 2.0*GeV, false, false, true); static Parameter interfaceStrangeMass ("StrangeMass", "The mass of the strange quark", &BaryonThreeQuarkModelFormFactor::_mstrange, GeV, .57*GeV, ZERO, 2.0*GeV, false, false, true); static Parameter interfaceLambdaQ ("LambdaQ", "Heavy Baryon Size Parameter", &BaryonThreeQuarkModelFormFactor::_LambdaQ, GeV, 2.5*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceLambdaqq ("Lambdaqq", "The size parameter for light quarks", &BaryonThreeQuarkModelFormFactor::_Lambdaqq, GeV, 0.71*GeV, ZERO, 2.0*GeV, false, false, true); static Parameter interfaceLambdasq ("Lambdasq", "The size parameter for one strange quark", &BaryonThreeQuarkModelFormFactor::_Lambdasq, GeV, 0.85*GeV, ZERO, 2.0*GeV, false, false, true); static Parameter interfaceLambdass ("Lambdass", "The size parameter with two strange quarks.", &BaryonThreeQuarkModelFormFactor::_Lambdass, GeV, 1.0*GeV, ZERO, 2.0*GeV, false, false, true); static ParVector interfaceC0 ("C0", "The coefficient of the expansion for xi_0.", &BaryonThreeQuarkModelFormFactor::_C0, 0, 0, 0, -1E20, 1E20, false, false, true); static ParVector interfaceC1 ("C1", "The coefficient of the expansion for xi_1.", &BaryonThreeQuarkModelFormFactor::_C1, 0, 0, 0, -1E20, 1E20, false, false, true); static ParVector interfaceC2 ("C2", "The coefficient of the expansion for xi_2.", &BaryonThreeQuarkModelFormFactor::_C2, 0, 0, 0, -1E20, 1E20, false, false, true); static Switch interfaceInitialize ("Initialize", "Initialize the coefficient for the expansion of the form-factor", &BaryonThreeQuarkModelFormFactor::_initialize, false, false, false); static SwitchOption interfaceInitializeInitialize (interfaceInitialize, "Yes", "Perform the initialize", true); static SwitchOption interfaceInitializeNoInitialize (interfaceInitialize, "No", "No initialization", false); } void BaryonThreeQuarkModelFormFactor:: SpinHalfSpinHalfFormFactor(Energy2 q2,int,int id0,int id1,Energy m0,Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo , Virtuality virt) { assert(virt==SpaceLike); useMe(); // this model is based on heavy quark theory // therefore most of the factors are zero Complex g1v(0.),g1a(0.),g2v(0.),g2a(0.),g3a(0.),g3v(0.); // work out which light quark constant to use double lambdabar;unsigned int ioff(0); if(abs(id1)==4332) { lambdabar=_Lambdass/_LambdaQ; ioff=_order+1; } else if(abs(id1)==4232||abs(id1)==4132) { lambdabar=_Lambdasq/_LambdaQ; ioff=_order+1; } else { lambdabar=_Lambdaqq/_LambdaQ; } // the omega value double omega = 0.5/m0/m1*(m0*m0+m1*m1-q2); // calculate the form factor // the phi functions vector phi(phiFunction(omega)); // use the xi0 functions if(abs(id0)==5122||abs(id0)==5232||abs(id0)==5132) { double power(1.),numer(0.),denom(0.); for(unsigned int ix=0;ix<=_order;++ix) { numer+=phi[ix]*power*_C0[ix+ioff]; denom+=power*_C0[ix+ioff]; power*=lambdabar; } g1v=numer/denom; g1a=g1v; } else { double power(1.),numer[2]={0.,0.},denom(0.); for(unsigned int ix=0;ix<=_order;++ix) { numer[0]+=phi[ix]*power*_C1[ix+ioff]; denom+=power*_C1[ix+ioff]; numer[1]+=power*_C2[ix+ioff]*(phi[ix]-phi[ix+2]); power*=lambdabar; } numer[1]/=(omega-1.); double xi1(numer[0]/denom),xi2(numer[1]/denom); // the couplings in the velocity form g1v = -(omega*xi1-(sqr(omega)-1.)*xi2)/3.; g1a = g1v; g2v = 2./3.*(xi1-(omega-1.)*xi2); g3v = g2v; g2a = 2./3.*(xi1-(omega+1.)*xi2); g3a =-g2a; } // convert to our form f1v = g1v+Complex(0.5*(m0+m1)*(g2v/m0+g3v/m1)); f1a =-g1a+Complex(0.5*(m0-m1)*(g2a/m0+g3a/m1)); f2v = Complex(0.5*(m0+m1)*( g2v/m0+g3v/m1)); f3v = Complex(0.5*(m0+m1)*( g2v/m0-g3v/m1)); f2a =-Complex(0.5*(m0+m1)*( g2a/m0+g3a/m1)); f3a = Complex(0.5*(m0+m1)*(-g2a/m0+g3a/m1)); } void BaryonThreeQuarkModelFormFactor:: SpinHalfSpinThreeHalfFormFactor(Energy2 q2,int,int,int id1,Energy m0, Energy m1, Complex & f1v,Complex & f2v, Complex & f3v,Complex & f4v,Complex & f1a, Complex & f2a,Complex & f3a,Complex & f4a, + FlavourInfo , Virtuality virt) { assert(virt==SpaceLike); useMe(); // work out which light quark constant to use double lambdabar;unsigned int ioff(0); if(abs(id1)==4334) { lambdabar=_Lambdass/_LambdaQ; ioff=_order+1; } else if(abs(id1)==4234||abs(id1)==4134||abs(id1)==3324) { lambdabar=_Lambdasq/_LambdaQ; ioff=_order+1; } else { lambdabar=_Lambdaqq/_LambdaQ; } // the omega value double omega=0.5/m0/m1*(m0*m0+m1*m1-q2); // calculate the form factor // the phi functions vector phi=phiFunction(omega); double power(1.),numer[2]={0.,0.},denom(0.); for(unsigned int ix=0;ix<=_order;++ix) { numer[0]+=phi[ix]*power*_C1[ix+ioff]; denom+=power*_C1[ix+ioff]; numer[1]+=power*_C2[ix+ioff]*(phi[ix]-phi[ix+2]); power*=lambdabar; } numer[1]/=(omega-1.); double xi1=numer[0]/denom; double xi2=numer[1]/denom; // couplings in the velocity form Complex g1v,g2v,g3v,g4v,g1a,g2a,g3a,g4a; double orr(1./sqrt(3.)); Energy msum(m0+m1); Energy2 msum2(msum*msum); g1v = 2.*orr*xi1; g1a = -g1v; g2v = orr*(xi1-(omega-1)*xi2); g2a = orr*(xi1-(omega+1)*xi2); g3v = 0.; g3a = 0.; g4v = -2.*orr*xi2; g4a = -g4v; // convert to our form f1v = g1v; f1a =-g1a; f2v = g2v*msum/m0; f2a =-g2a*msum/m0; f3v = Complex(msum2/m0*(g3v/m0+g4v/m1)); f3a =-Complex(msum2/m0*(g3a/m0+g4a/m1)); f4v = Complex(msum2/m0/m0*g3v); f4a =-Complex(msum2/m0/m0*g3a); } void BaryonThreeQuarkModelFormFactor::dataBaseOutput(ofstream & output,bool header, bool create) const { if(header) output << "update decayers set parameters=\""; if(create) output << "create Herwig::BaryonThreeQuarkModelFormFactor " << name() << " \n"; output << "newdef " << name() << ":Order " << _order << " \n"; output << "newdef " << name() << ":LightMass " << _mlight/GeV << " \n"; output << "newdef " << name() << ":StrangeMass " << _mstrange/GeV << " \n"; output << "newdef " << name() << ":LambdaQ " << _LambdaQ/GeV << " \n"; output << "newdef " << name() << ":Lambdaqq " << _Lambdaqq/GeV << " \n"; output << "newdef " << name() << ":Lambdasq " << _Lambdasq/GeV << " \n"; output << "newdef " << name() << ":Lambdass " << _Lambdass/GeV << " \n"; // the number of terms to include in the sum for the form-factors for(unsigned int ix=0;ix<_C0.size();++ix) output << "insert " << name() << ":C0 " << ix << " " << _C0[ix] << " \n"; for(unsigned int ix=0;ix<_C1.size();++ix) output << "insert " << name() << ":C1 " << ix << " " << _C1[ix] << " \n"; for(unsigned int ix=0;ix<_C2.size();++ix) output << "insert " << name() << ":C2 " << ix << " " << _C2[ix] << " \n"; BaryonFormFactor::dataBaseOutput(output,false,false); if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } double BaryonThreeQuarkModelFormFactor::operator ()(double x) const { // convert the integration variable double y=(1.-x)/x; double output =exp(-24.*_mu2*y)*y; // the integrals double I,Im2; SN(y,_N,Im2,I); double Nfact=0.5*_N+1.0; output *= (_a+Nfact/24./(1.+y))*I+1./12./(1.+y)*(_b-Nfact)*Im2; return output; } void BaryonThreeQuarkModelFormFactor::SN(double y, int N, double & SNm2, double & SN) const { // special cases for the low lying values if(N==0) { double root=sqrt((1.+y)*(3.+y)); SN = 0.5/y*sqrt((1.+y)/(3.+y))*log((root+y)/(root-y)); SNm2 = 0.5/(y+3.)*(SN+(1.+y)/(3.+4.*y)); } else if(N==1) { SN = sqrt(1.+y)/y*asin(y/sqrt(1.+y)/sqrt(3.+y)); SNm2 = 1./(3.+y)*sqrt((1.+y)/(3.+4.*y)); } else if(N==2) { double root=sqrt((1.+y)*(3.+y)); SN = 1.; SNm2 = 0.5/y*sqrt((1.+y)/(3.+y))*log((root+y)/(root-y)); } // the general case else { int ix; double root; if(N%2==0) { SN=1.; ix=2; root=1.; } else { SN = sqrt(1.+y)/y*asin(y/sqrt(1.+y)/sqrt(3.+y)); ix=1; root=sqrt((1.+y)/(3.+4.*y)); } do { SNm2=SN; ix+=2; root*=(3.+4.*y)/(1.+y); SN=1.0/(ix-1.0)*(root+(ix-2.0)*(y+3.)*SNm2); } while(ix BaryonThreeQuarkModelFormFactor::phiFunction(double omega) { vector output; double root(sqrt(omega*omega-1.)); output.push_back(1./root*log(omega+root)); if(omega<1.00001) output.back()=1.; if(_order>0) output.push_back(2./(omega+1.)); if(_order<2) return output; for(unsigned int ix=2;ix<=_order+2;++ix) { output.push_back(2./ix/(omega+1.)*(1.+(ix-1)*output[ix-2])); } return output; } Energy BaryonThreeQuarkModelFormFactor::widthIntegrand(double omega,Energy m0, Energy m1, int type, int ,int id0, int id1) { // prefactors double omegamax=0.5*(m0*m0+m1*m1)/m0/m1; double pi=acos(-1.); InvEnergy kw=sqr(generator()->standardModel()->fermiConstant()) /8./pi/pi/pi*m1*m1*m1/6.*(omegamax-omega)*sqrt(omega*omega-1.); Energy2 q2 = sqr(m0)+sqr(m1)-2.*m0*m1*omega; if(type<=2) { Complex f1v,f2v,f3v,f1a,f2a,f3a; - SpinHalfSpinHalfFormFactor(q2,0,id0,id1,m0,m1,f1v,f2v,f3v,f1a,f2a,f3a); + SpinHalfSpinHalfFormFactor(q2,0,id0,id1,m0,m1,f1v,f2v,f3v,f1a,f2a,f3a,FlavourInfo()); Complex left =f1v-f1a-f2v-double((m0-m1)/(m0+m1))*f2a; Complex right =f1v+f1a-f2v+double((m0-m1)/(m0+m1))*f2a; double g1v = 0.5*( left+right).real(); double g2v = m0*(f2v+f3v).real()/(m0+m1); double g3v = m1*(f2v-f3v).real()/(m0+m1); double g1a = -0.5*(+right-left).real(); double g2a = -m0*(f2a+f3a).real()/(m0+m1); double g3a = -m1*(f2a-f3a).real()/(m0+m1); Energy Hpp = -2.*sqrt(m0*m1*(omega-1.))*g1v+2.*sqrt(m0*m1*(omega+1))*g1a; Energy Hmm = -2.*sqrt(m0*m1*(omega-1.))*g1v-2.*sqrt(m0*m1*(omega+1))*g1a; Energy Hp0 = (sqrt(2.*m0*m1*(omega-1))*((m0+m1)*g1v+m1*(omega+1)*g2v+m0*(omega+1)*g3v)- sqrt(2.*m0*m1*(omega+1))*((m0-m1)*g1a-m1*(omega-1)*g2a-m0*(omega-1)*g3a))/sqrt(q2); Energy Hm0 = (sqrt(2.*m0*m1*(omega-1))*((m0+m1)*g1v+m1*(omega+1)*g2v+m0*(omega+1)*g3v)+ sqrt(2.*m0*m1*(omega+1))*((m0-m1)*g1a-m1*(omega-1)*g2a-m0*(omega-1)*g3a))/sqrt(q2); return kw*sqr(0.04)*(sqr(Hpp)+sqr(Hmm)+sqr(Hp0)+sqr(Hm0)); } else { Complex f1v,f2v,f3v,f4v,f1a,f2a,f3a,f4a; double g1v,g2v,g3v,g4v,g1a,g2a,g3a,g4a; SpinHalfSpinThreeHalfFormFactor(q2,0,id0,id1,m0,m1, f1v,f2v,f3v,f4v, - f1a,f2a,f3a,f4a); + f1a,f2a,f3a,f4a,FlavourInfo()); g1v = f1v.real(); g1a = -f1a.real(); g2v = m0/(m0+m1)*f2v.real(); g2a =-m0/(m0+m1)*f2a.real(); g3v = sqr(m0/(m0+m1))*f4v.real(); g3a =-sqr(m0/(m0+m1))*f4a.real(); g4v = m0*m1/sqr(m0+m1)*(f3v.real()-f4v.real()); g4a =-m0*m1/sqr(m0+m1)*(f3a.real()-f4a.real()); Energy HppC = +sqrt(2./3.)*sqrt(m0*m1*(omega-1))*(g1v-2.*(omega+1)*g2v) -sqrt(2./3.)*sqrt(m0*m1*(omega+1))*(g1a-2.*(omega-1)*g2a); Energy HmmC = -sqrt(2./3.)*sqrt(m0*m1*(omega-1))*(g1v-2.*(omega+1)*g2v) -sqrt(2./3.)*sqrt(m0*m1*(omega+1))*(g1a-2.*(omega-1)*g2a); Energy HppbC = -sqrt(2.*m0*m1*(omega-1))*g1v -sqrt(2.*m0*m1*(omega+1))*g1a; Energy HmmbC = sqrt(2.*m0*m1*(omega-1))*g1v -sqrt(2.*m0*m1*(omega+1))*g1a; Energy Hp0C = ( -2./sqrt(3.)*sqrt(m0*m1*(omega-1))* ((m0*omega-m1)*g1v-(m0-m1)*(omega+1)*g2v+m1*(sqr(omega)-1)*g3v+m0*(sqr(omega)-1)*g4v) -2./sqrt(3.)*sqrt(m0*m1*(omega+1))* ((m0*omega-m1)*g1a+(m0+m1)*(omega-1)*g2a+m1*(sqr(omega)-1)*g3a+m0*(sqr(omega)-1)*g4a))/sqrt(q2); Energy Hm0C = ( +2./sqrt(3.)*sqrt(m0*m1*(omega-1))* ((m0*omega-m1)*g1v-(m0-m1)*(omega+1)*g2v+m1*(sqr(omega)-1)*g3v+m0*(sqr(omega)-1)*g4v) -2./sqrt(3.)*sqrt(m0*m1*(omega+1))* ((m0*omega-m1)*g1a+(m0+m1)*(omega-1)*g2a+m1*(sqr(omega)-1)*g3a+m0*(sqr(omega)-1)*g4a))/sqrt(q2); return kw*sqr(0.04)*(sqr(HppC)+sqr(HmmC)+sqr(HppbC)+sqr(HmmbC)+sqr(Hp0C)+sqr(Hm0C)); } } diff --git a/Decay/FormFactors/BaryonThreeQuarkModelFormFactor.h b/Decay/FormFactors/BaryonThreeQuarkModelFormFactor.h --- a/Decay/FormFactors/BaryonThreeQuarkModelFormFactor.h +++ b/Decay/FormFactors/BaryonThreeQuarkModelFormFactor.h @@ -1,292 +1,294 @@ // -*- C++ -*- #ifndef HERWIG_BaryonThreeQuarkModelFormFactor_H #define HERWIG_BaryonThreeQuarkModelFormFactor_H // // This is the declaration of the BaryonThreeQuarkModelFormFactor class. // #include "BaryonFormFactor.h" #include "ThePEG/PDT/ParticleData.h" namespace Herwig { using namespace ThePEG; /** \ingroup Decay * * The BaryonThreeQuarkModelFormFactor class implements the * form factors for the semi-leptonic decay of baryons containing a heavy quark * from the relativistic three-quark model calculation of PRD56, 348. * * As the only formulae in the paper are for the heavy-to-heavy i.e. bottom * to charm decay this there are the only modes included, although the paper * also includes charm decays and bottom decays to light quarks. * * The form factors are calculated by numerical computing the integrals from * PRD56, 348 to obtain the coefficients for the expansion of the form factors. * * @see BaryonFormFactor */ class BaryonThreeQuarkModelFormFactor: public BaryonFormFactor { public: /** * Default constructor */ BaryonThreeQuarkModelFormFactor(); /** @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); //@} /** * Standard Init function used to initialize the interfaces. */ 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 {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);} //@} public: /** @name Form Factors */ //@{ /** * The form factor for the weak decay of a spin \f$\frac12\f$ baryon to a * spin \f$\frac12\f$ baryon. * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form factor list. * @param id0 The PDG code of the incoming baryon. * @param id1 The PDG code of the outgoing baryon. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param f1v The form factor \f$F^V_1\f$. * @param f2v The form factor \f$F^V_2\f$. * @param f3v The form factor \f$F^V_3\f$. * @param f1a The form factor \f$F^A_1\f$. * @param f2a The form factor \f$F^A_2\f$. * @param f3a The form factor \f$F^A_3\f$. * @param virt Whether the q2 is space or timelike */ virtual void SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo flavour, Virtuality virt=SpaceLike); /** * The form factor for the weak decay of a spin \f$\frac12\f$ baryon to a * spin \f$\frac32\f$ baryon. * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form factor list. * @param id0 The PDG code of the incoming baryon. * @param id1 The PDG code of the outgoing baryon. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param g1v The form factor \f$G^V_1\f$. * @param g2v The form factor \f$G^V_2\f$. * @param g3v The form factor \f$G^V_3\f$. * @param g4v The form factor \f$G^V_4\f$. * @param g1a The form factor \f$G^A_1\f$. * @param g2a The form factor \f$G^A_2\f$. * @param g3a The form factor \f$G^A_3\f$. * @param g4a The form factor \f$G^A_4\f$. * @param virt Whether the q2 is space or timelike */ virtual void SpinHalfSpinThreeHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & g1v,Complex & g2v,Complex & g3v, Complex & g4v,Complex & g1a,Complex & g2a, Complex & g3a,Complex & g4a, + FlavourInfo flavour, Virtuality virt=SpaceLike); //@} /** * 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 * @param create Whether or not to add a statement creating the object */ virtual void dataBaseOutput(ofstream & os,bool header,bool create) const; /** * The integrand for the coefficients of the expansion. This is a function of the * integration variable \f$x\f$ which is chosen to transform the integrand over \f$y\f$ * which is from \f$0\f$ to \f$\infty\f$ to an integral between 0 and 1. This means * that \f$y=\frac{1-x}{x}\f$. * @param x The integration variable. */ double operator ()(double x) const; /** Argument type for GaussianIntegrator */ typedef double ArgType; /** Return type for GaussianIntegrator */ typedef double ValType; /** * The integrand for the semi-analytic calculation of the semi-leptonic width. * This is included for testing purposes. * @param omega The \f$\omega\f$ parameter of the heavy quark form-factors. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param type The type of the decay * @param imass The baryonic mass parameter to use. * @param id0 PDG code for the decaying particle * @param id1 PDG code for the decay product */ Energy widthIntegrand(double omega,Energy m0, Energy m1, int type, int imass, int id0,int id1); protected: /** @name Function needed to calculate the form factors */ //@{ /** * Returns the function \f$\Phi_N\f$ function of PRD56, 348 as a function * of \f$\omega\f$. */ vector phiFunction(double); /** * The integral of a power of the the \f$S\f$ function of PRD56, 348 with respect * to \f$\theta\f$. This is used to * calculate the coefficients of the expansion to compute the form factors. * The integral is calculated by computing a low power of the integrand and then * using recursion relations to calculate the pwoer requested. * @param y The \f$y\f$ variable of the function which is integrate over numerically. * @param N The power to which the integrand is raised. * @param SNm2 The integral with the function raised to the power \f$N/2-1\f$. * @param SN The integral with the function raised to the power \f$N\f$. */ void SN(double y,int N,double & SNm2,double & SN) 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(); //@} private: /** * Private and non-existent assignment operator. */ BaryonThreeQuarkModelFormFactor & operator=(const BaryonThreeQuarkModelFormFactor &) = delete; private: /** @name Parameters for the form factors */ //@{ /** * Initialization of the expansion coefficients */ bool _initialize; /** * Order of the expansion for the form factors. */ unsigned int _order; /** * Mass of the light quarks used in the calculation of the form factors. */ Energy _mlight; /** * Mass of the strange quark used in the calculation of the form factors. */ Energy _mstrange; /** * The heavy quark \f$\Lambda_Q\f$ parameter for the calculation of the form factors. */ Energy _LambdaQ; /** * The \f$\Lambda_{qq}\f$ parameter for the calculation of the form factors. */ Energy _Lambdaqq; /** * The \f$\Lambda_{sq}\f$ parameter for the calculation of the form factors. */ Energy _Lambdasq; /** * The \f$\Lambda_{ss}\f$ parameter for the calculation of the form factors. */ Energy _Lambdass; /** * Coefficients for the expansion of the \f$\xi(\omega)\f$ form factor. */ vector _C0; /** * Coefficients for the expansion of the \f$\xi_1(\omega)\f$ form factor. */ vector _C1; /** * Coefficients for the expansion of the \f$\xi_2(\omega)\f$ form factor. */ vector _C2; /** * Coefficient of the first term in the integrand for the coefficient calculation. */ double _a; /** * Coefficient of the second term in the integrand for the coefficient calculation. */ double _b; /** * The \f$\mu^2\f$ parameter for the coefficient calculation. */ double _mu2; /** * the order of the coefficient being calculated. */ int _N; //@} }; } #endif /* HERWIG_BaryonThreeQuarkModelFormFactor_H */ diff --git a/Decay/FormFactors/ChengHeavyBaryonFormFactor.cc b/Decay/FormFactors/ChengHeavyBaryonFormFactor.cc --- a/Decay/FormFactors/ChengHeavyBaryonFormFactor.cc +++ b/Decay/FormFactors/ChengHeavyBaryonFormFactor.cc @@ -1,441 +1,443 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the ChengHeavyBaryonFormFactor class. // #include "ChengHeavyBaryonFormFactor.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; using namespace ThePEG; ChengHeavyBaryonFormFactor::ChengHeavyBaryonFormFactor() { // consituent quark masses _mu = 338*MeV; _md = 322*MeV; _ms = 510*MeV; _mc = 1.6*GeV; _mb = 5.0*GeV; // masses for the q^2 dependence _mVbc = 6.34*GeV; _mVbs = 5.42*GeV; _mVcs = 2.11*GeV; _mVbd = 5.32*GeV; _mVcu = 2.01*GeV; _mAbc = 6.73*GeV; _mAbs = 5.86*GeV; _mAcs = 2.54*GeV; _mAbd = 5.71*GeV; _mAcu = 2.42*GeV; double one3(1./sqrt(3.)),one2(1./sqrt(2.)); // lambda_b to lambda_c addFormFactor(5122,4122,2,2,1,2,5,4);_Nfi.push_back(1. );_eta.push_back(1.); // lambda_b to lambda addFormFactor(5122,3122,2,2,1,2,5,3);_Nfi.push_back(one3 );_eta.push_back(1.); // lambda_b to n addFormFactor(5122,2112,2,2,1,2,5,1);_Nfi.push_back(one2 );_eta.push_back(1.); // xi_b to xi_c addFormFactor(5232,4232,2,2,2,3,5,4);_Nfi.push_back(1. );_eta.push_back(1.); addFormFactor(5132,4132,2,2,1,3,5,4);_Nfi.push_back(1. );_eta.push_back(1.); // xi_b to xi addFormFactor(5232,3322,2,2,2,3,5,3);_Nfi.push_back(one2 );_eta.push_back(1.); addFormFactor(5132,3312,2,2,1,3,5,3);_Nfi.push_back(one2 );_eta.push_back(1.); // xi_b to sigma addFormFactor(5232,3212,2,2,2,3,5,1);_Nfi.push_back(0.5 );_eta.push_back(1.); addFormFactor(5132,3112,2,2,1,3,5,1);_Nfi.push_back(0.5 );_eta.push_back(1.); // xi_b to lambda addFormFactor(5232,3122,2,2,2,3,5,1);_Nfi.push_back(one3/2.);_eta.push_back(1.); // omega_b to omega_c addFormFactor(5332,4332,2,2,3,3,5,4);_Nfi.push_back(1. );_eta.push_back(-1./3.); // omega_b to xi addFormFactor(5332,3312,2,2,3,3,5,1);_Nfi.push_back(one3 );_eta.push_back(-1./3.); // omega_b to omega_c* addFormFactor(5332,4334,2,4,3,3,5,4);_Nfi.push_back(1. );_eta.push_back(0.); // omega_b to omega addFormFactor(5332,3334,2,4,3,3,5,3);_Nfi.push_back(1. );_eta.push_back(0.); // omega_b to xi* addFormFactor(5332,3314,2,4,3,3,5,1);_Nfi.push_back(one3 );_eta.push_back(0.); // omega_c to omega addFormFactor(4332,3334,2,4,3,3,4,3);_Nfi.push_back(1. );_eta.push_back(0.); // omega_c to xi* addFormFactor(4332,3324,2,4,3,3,4,2);_Nfi.push_back(one3 );_eta.push_back(0.); // lambda_c to lambda_0 addFormFactor(4122,3122,2,2,1,2,4,3);_Nfi.push_back(1./sqrt(3.));_eta.push_back(1.); // xi_c to xi addFormFactor(4232,3322,2,2,2,3,4,3);_Nfi.push_back(1./sqrt(3.));_eta.push_back(1.); addFormFactor(4132,3312,2,2,1,3,4,3);_Nfi.push_back(1./sqrt(3.));_eta.push_back(1.); // initial number of form factors initialModes(numberOfFactors()); } void ChengHeavyBaryonFormFactor::doinit() { BaryonFormFactor::doinit(); // check the parameters are consistent unsigned int isize(numberOfFactors()); if(isize!=_eta.size()||isize!=_Nfi.size()) {throw InitException() << "Inconsistent paramters in ChengHeavyBaryon" << "FormFactor::doinit() " << Exception::abortnow;} Energy mi,mf,mq,mQ,lambda,delta,msum; int id0,id1,inspin,outspin,isp1,isp2,inq,outq; for(unsigned int ix=0;ixmass(); mf=getParticleData(id1)->mass(); msum=mi+mf; // masses of the incoming and outgoing quarks if((id0==4122&&id1==3122)||(id0==4232&&id1==3322)||(id0==4132&&id1==3312)|| (id0==4332&&id1==3334)) {mq=_ms;mQ=_mc;} else if((id0==4332&&id1==3322)||(id0==4332&&id1==3324)) {mq=_mu;mQ=_mc;} else if((id0==5122&&id1==4122)||(id0==5232&&id1==4232)||(id0==5132&&id1==4132)|| (id0==5332&&id1==4332)||(id0==5332&&id1==4334)) {mq=_mc;mQ=_mb;} else if((id0==5122&&id1==3122)||(id0==5132&&id1==3312)||(id0==5232&&id1==3322)|| (id0==5332&&id1==3334)) {mq=_ms;mQ=_mb;} else if((id0==5122&&id1==2112)||(id0==5132&&id1==3112)||(id0==5232&&id1==3212)|| (id0==5332&&id1==3312)||(id0==5232&&id1==3122)||(id0==5332&&id1==3314)) {mq=_md;mQ=_mb;} else {throw InitException() << "Unknown decay in ChengHeavyBaryon" << "FormFactor::doinit() " << Exception::abortnow;} // parameters lambda = mf-mq; delta = mi-mf; // compute the form-factors if(inspin==2&&outspin==2) { _f1.push_back(_Nfi[ix]*(1.-0.5*delta/mi +0.25*delta/mi/mq*(1.-0.5*lambda/mf)* (mi+mf-_eta[ix]*delta) -0.125*delta/mi/mf/mQ*lambda*(mi+mf+_eta[ix]*delta))); _f2.push_back(_Nfi[ix]*msum*(0.5/mi+0.25/mi/mq*(1.-0.5*lambda/mf)* (delta-(mi+mf)*_eta[ix]) -0.125*lambda/mi/mf/mQ*(delta+(mi+mf)*_eta[ix]))); _f3.push_back(_Nfi[ix]*msum*(0.5/mi-0.25/mi/mq*(1.-0.5*lambda/mf)* (mi+mf-_eta[ix]*delta) +0.125*lambda/mi/mf/mQ*(mi+mf+_eta[ix]*delta))); _g1.push_back(_Nfi[ix]*_eta[ix]*(1.+0.25*delta*lambda*(1./mi/mq-1./mf/mQ))); _g2.push_back(-0.25*msum*_Nfi[ix]*_eta[ix]*lambda*(1./mi/mq-1./mf/mQ)); _g3.push_back(-0.25*msum*_Nfi[ix]*_eta[ix]*lambda*(1./mi/mq+1./mf/mQ)); } else if(inspin==2&&outspin==4) { _f1.push_back(2.*_Nfi[ix]/sqrt(3.)*(1.+0.5*lambda*(1./mq+1./mQ))); _f2.push_back(_Nfi[ix]*msum/sqrt(3.)/mi*(1.+0.5*lambda*(1./mq+1./mQ))); _f3.push_back(-_Nfi[ix]*msum*msum/mi/mf/sqrt(3.)* (1.+0.5*lambda*(1./mq+1./mQ))); _g1.push_back(-2./sqrt(3.)*_Nfi[ix]); _g2.push_back(-_Nfi[ix]*msum/sqrt(3.)*lambda/mq/mi); _g3.push_back(-_f3.back()); } else {throw InitException() << "Unknown spin combination in ChengHeavyBaryon" << "FormFactor::doinit() " << Exception::abortnow;} } for(unsigned int ix=0;ixmass(); tcPDPtr part1=getParticleData(id1); Energy m1=part1->mass(); if ( part1->iSpin() == 2 ) { Complex f1v,f2v,f3v,f4v,f1a,f2a,f3a,f4a; // dummy variables - SpinHalfSpinHalfFormFactor(ZERO,ix,id0,id1,m0,m1,f1v,f2v,f3v,f1a,f2a,f3a); + SpinHalfSpinHalfFormFactor(ZERO,ix,id0,id1,m0,m1,f1v,f2v,f3v,f1a,f2a,f3a,FlavourInfo()); } else { Complex f1v,f2v,f3v,f4v,f1a,f2a,f3a,f4a; // dummy variables SpinHalfSpinThreeHalfFormFactor(ZERO,ix,id0,id1,m0,m1,f1v,f2v,f3v, - f4v,f1a,f2a,f3a,f4a); + f4v,f1a,f2a,f3a,f4a,FlavourInfo()); } } } void ChengHeavyBaryonFormFactor::persistentOutput(PersistentOStream & os) const { os << ounit(_mu,MeV) << ounit(_md,MeV) << ounit(_ms,MeV) << ounit(_mc,MeV) << ounit(_mb,MeV) << _Nfi << _eta << _f1 << _f2 << _f3 << _g1 << _g2 << _g3 << ounit(_mVbc,MeV) << ounit(_mVbs,MeV) << ounit(_mVcs,MeV) << ounit(_mVbd,MeV) << ounit(_mVcu,MeV) << ounit(_mAbc,MeV) << ounit(_mAbs,MeV) << ounit(_mAcs,MeV) << ounit(_mAbd,MeV) << ounit(_mAcu,MeV); } void ChengHeavyBaryonFormFactor::persistentInput(PersistentIStream & is, int) { is >> iunit(_mu,MeV) >> iunit(_md,MeV) >> iunit(_ms,MeV) >> iunit(_mc,MeV) >> iunit(_mb,MeV) >> _Nfi >> _eta >> _f1 >> _f2 >> _f3 >> _g1 >> _g2 >> _g3 >> iunit(_mVbc,MeV) >> iunit(_mVbs,MeV) >> iunit(_mVcs,MeV) >> iunit(_mVbd,MeV) >> iunit(_mVcu,MeV) >> iunit(_mAbc,MeV) >> iunit(_mAbs,MeV) >> iunit(_mAcs,MeV) >> iunit(_mAbd,MeV) >> iunit(_mAcu,MeV); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigChengHeavyBaryonFormFactor("Herwig::ChengHeavyBaryonFormFactor", "HwFormFactors.so"); void ChengHeavyBaryonFormFactor::Init() { static ClassDocumentation documentation ("The ChengHeavyBaryonFormFactor class is the implementation" " of the form-factors of PRD53, 1457 and PRD56, 2799 for the weak decay of" "baryons containing a heavy quark. This model can be used for either" "semi-leptonic decays, or with the factorization approximation for" " non-leptonic weak decays", "The weak decay of baryons containing a heavy quark used form factors from " "\\cite{Cheng:1995fe,Cheng:1996cs}.", "%\\cite{Cheng:1995fe}\n" "\\bibitem{Cheng:1995fe}\n" " H.~Y.~Cheng and B.~Tseng,\n" " %``1/M corrections to baryonic form-factors in the quark model,''\n" " Phys.\\ Rev.\\ D {\\bf 53} (1996) 1457\n" " [Erratum-ibid.\\ D {\\bf 55} (1997) 1697]\n" " [arXiv:hep-ph/9502391].\n" " %%CITATION = PHRVA,D53,1457;%%\n" "%\\cite{Cheng:1996cs}\n" "\\bibitem{Cheng:1996cs}\n" " H.~Y.~Cheng,\n" " %``Nonleptonic weak decays of bottom baryons,''\n" " Phys.\\ Rev.\\ D {\\bf 56} (1997) 2799\n" " [arXiv:hep-ph/9612223].\n" " %%CITATION = PHRVA,D56,2799;%%\n" ); static Parameter interfaceUpMass ("DownMass", "The consituent mass of the down quark", &ChengHeavyBaryonFormFactor::_md, GeV, 0.322*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceDownMass ("UpMass", "The consituent mass of the up quark", &ChengHeavyBaryonFormFactor::_mu, GeV, 0.338*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfacStrangeMass ("StrangeMass", "The consituent mass of the strange quark", &ChengHeavyBaryonFormFactor::_ms, GeV, 0.510*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceCharmMass ("CharmMass", "The consituent mass of the charm quark", &ChengHeavyBaryonFormFactor::_mc, GeV, 1.6*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceBottomMass ("BottomMass", "The consituent mass of the bottom quark", &ChengHeavyBaryonFormFactor::_mb, GeV, 5.0*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceVectorMassbc ("VectorMassbc", "The vector mass for the b->c transitions.", &ChengHeavyBaryonFormFactor::_mVbc, GeV, 6.34*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceAxialMassbc ("AxialMassbc", "The axial-vector mass for the b->c transitions.", &ChengHeavyBaryonFormFactor::_mAbc, GeV, 6.73*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceVectorMassbs ("VectorMassbs", "The vector mass for the b->s transitions.", &ChengHeavyBaryonFormFactor::_mVbs, GeV, 5.42*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceAxialMassbs ("AxialMassbs", "The axial-vector mass for the b->s transitions.", &ChengHeavyBaryonFormFactor::_mAbs, GeV, 5.86*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceVectorMassbd ("VectorMassbd", "The vector mass for the b->d transitions.", &ChengHeavyBaryonFormFactor::_mVbd, GeV, 5.32*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceAxialMassbd ("AxialMassbd", "The axial-vector mass for the b->d transitions.", &ChengHeavyBaryonFormFactor::_mAbd, GeV, 5.71*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceVectorMasscs ("VectorMasscs", "The vector mass for the c->s transitions.", &ChengHeavyBaryonFormFactor::_mVcs, GeV, 2.11*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceAxialMasscs ("AxialMasscs", "The axial-vector mass for the c->s transitions.", &ChengHeavyBaryonFormFactor::_mAcs, GeV, 2.54*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceVectorMasscu ("VectorMasscu", "The vector mass for the c->u transitions.", &ChengHeavyBaryonFormFactor::_mVcu, GeV, 2.01*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceAxialMasscu ("AxialMasscu", "The axial-vector mass for the c->u transitions.", &ChengHeavyBaryonFormFactor::_mAcu, GeV, 2.42*GeV, ZERO, 10.0*GeV, false, false, true); static ParVector interfaceNfi ("Nfi", "The prefactor for a given form factor", &ChengHeavyBaryonFormFactor::_Nfi, -1, 1.0, -10.0, 10.0, false, false, true); static ParVector interfaceEta ("Eta", "The eta parameter for the form factor", &ChengHeavyBaryonFormFactor::_eta, -1, 0.0, -10.0, 10.0, false, false, true); } // form factor for spin-1/2 to spin-1/2 void ChengHeavyBaryonFormFactor:: SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc,int id0,int id1, Energy m0,Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo , Virtuality virt) { assert(virt==SpaceLike); useMe(); id0=abs(id0); id1=abs(id1); // masses for the energy dependence of the form-factors Energy mV(ZERO),mA(ZERO); if((id0==4122&&id1==3122)||(id0==4232&&id1==3322)||(id0==4132&&id1==3312)|| (id0==4332&&id1==3334)) {mA=_mAcs;mV=_mVcs;} else if((id0==4332&&id1==3322)||(id0==4332&&id1==3324)) {mA=_mAcu;mV=_mVcu;} else if((id0==5122&&id1==4122)||(id0==5232&&id1==4232)||(id0==5132&&id1==4132)|| (id0==5332&&id1==4332)||(id0==5332&&id1==4334)) {mA=_mAbc;mV=_mVbc;} else if((id0==5122&&id1==3122)||(id0==5132&&id1==3312)||(id0==5232&&id1==3322)|| (id0==5332&&id1==3334)) {mA=_mAbs;mV=_mVbs;} else if((id0==5122&&id1==2112)||(id0==5132&&id1==3112)||(id0==5232&&id1==3212)|| (id0==5332&&id1==3312)||(id0==5232&&id1==3122)||(id0==5332&&id1==3314)) {mA=_mAbd;mV=_mVbd;} Energy delta=m0-m1; double Vfact = (1.-delta*delta/mV/mV)/(1.-q2/mV/mV); Vfact *=Vfact; double Afact = (1.-delta*delta/mA/mA)/(1.-q2/mA/mA); Afact *=Afact; f1v = _f1[iloc]*Vfact; f2v = _f2[iloc]*Vfact; f3v = _f3[iloc]*Vfact; f1a =-_g1[iloc]*Afact; f2a =-_g2[iloc]*Afact; f3a =-_g3[iloc]*Afact; } // form factor for spin-1/2 to spin-3/2 void ChengHeavyBaryonFormFactor:: SpinHalfSpinThreeHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & g1v,Complex & g2v,Complex & g3v, Complex & g4v,Complex & g1a,Complex & g2a, - Complex & g3a,Complex & g4a, Virtuality virt) { + Complex & g3a,Complex & g4a, + FlavourInfo , Virtuality virt) { assert(virt==SpaceLike); useMe(); id0=abs(id0); id1=abs(id1); // masses for the energy dependence of the form-factors Energy mV(ZERO),mA(ZERO); if((id0==4122&&id1==3122)||(id0==4232&&id1==3322)||(id0==4132&&id1==3312)|| (id0==4332&&id1==3334)) {mA=_mAcs;mV=_mVcs;} else if((id0==4332&&id1==3322)||(id0==4332&&id1==3324)) {mA=_mAcu;mV=_mVcu;} else if((id0==5122&&id1==4122)||(id0==5232&&id1==4232)||(id0==5132&&id1==4132)|| (id0==5332&&id1==4332)||(id0==5332&&id1==4334)) {mA=_mAbc;mV=_mVbc;} else if((id0==5122&&id1==3122)||(id0==5132&&id1==3312)||(id0==5232&&id1==3322)|| (id0==5332&&id1==3334)) {mA=_mAbs;mV=_mVbs;} else if((id0==5122&&id1==2112)||(id0==5132&&id1==3112)||(id0==5232&&id1==3212)|| (id0==5332&&id1==3312)||(id0==5232&&id1==3122)||(id0==5332&&id1==3314)) {mA=_mAbd;mV=_mVbd;} Energy delta=m0-m1; double Vfact = (1.-delta*delta/mV/mV)/(1.-q2/mV/mV); Vfact *=Vfact; double Afact = (1.-delta*delta/mA/mA)/(1.-q2/mA/mA); Afact *=Afact; g1v = _f1[iloc]*Vfact; g2v = _f2[iloc]*Vfact; g3v = _f3[iloc]*Vfact; g4v = 0.; g1a =-_g1[iloc]*Afact; g2a =-_g2[iloc]*Afact; g3a =-_g3[iloc]*Afact; g4a = 0.; } // output the information for the database void ChengHeavyBaryonFormFactor::dataBaseOutput(ofstream& output,bool header, bool create) const { if(header){output << "update decayers set parameters=\"";} if(create) {output << "create Herwig::ChengHeavyBaryonFormFactor " << name() << " \n";} output << "newdef " << name() << ":DownMass " << _md/GeV << " \n"; output << "newdef " << name() << ":UpMass " << _mu/GeV << " \n"; output << "newdef " << name() << ":StrangeMass " << _ms/GeV << " \n"; output << "newdef " << name() << ":CharmMass " << _mc/GeV << " \n"; output << "newdef " << name() << ":BottomMass " << _mb/GeV << " \n"; output << "newdef " << name() << ":VectorMassbc " << _mVbc/GeV << " \n"; output << "newdef " << name() << ":AxialMassbc " << _mAbc/GeV << " \n"; output << "newdef " << name() << ":VectorMassbs " << _mVbs/GeV << " \n"; output << "newdef " << name() << ":AxialMassbs " << _mAbs/GeV << " \n"; output << "newdef " << name() << ":VectorMassbd " << _mVbd/GeV << " \n"; output << "newdef " << name() << ":AxialMassbd " << _mAbd/GeV << " \n"; output << "newdef " << name() << ":VectorMasscs " << _mVcs/GeV << " \n"; output << "newdef " << name() << ":AxialMasscs " << _mAcs/GeV << " \n"; output << "newdef " << name() << ":VectorMasscu " << _mVcu/GeV << " \n"; output << "newdef " << name() << ":AxialMasscu " << _mAcu/GeV << " \n"; for(unsigned int ix=0;ix _Nfi; /** * The ratio of the flavour and spin overlap for the outgoing over incoming * baryon, \f$\eta\f$ from PRD56, 2799. */ vector _eta; /** * the form factor \f$f_1\f$ evaluated at the maximum value of \f$q^2\f$. */ vector _f1; /** * the form factor \f$f_2\f$ evaluated at the maximum value of \f$q^2\f$. */ vector _f2; /** * the form factor \f$f_3\f$ evaluated at the maximum value of \f$q^2\f$. */ vector _f3; /** * the form factor \f$g_1\f$ evaluated at the maximum value of \f$q^2\f$. */ vector _g1; /** * the form factor \f$g_2\f$ evaluated at the maximum value of \f$q^2\f$. */ vector _g2; /** * the form factor \f$g_3\f$ evaluated at the maximum value of \f$q^2\f$. */ vector _g3; /** * The pole mass for the vector form factors for the \f$b\to c\f$ transitiion. */ Energy _mVbc; /** * The pole mass for the vector form factors for the \f$b\to s\f$ transitiion. */ Energy _mVbs; /** * The pole mass for the vector form factors for the \f$c\to s\f$ transitiion. */ Energy _mVcs; /** * The pole mass for the vector form factors for the \f$b\to d\f$ transitiion. */ Energy _mVbd; /** * The pole mass for the vector form factors for the \f$c\to u\f$ transitiion. */ Energy _mVcu; /** * The pole mass for the axial-vector form factors for the \f$b\to c\f$ transitiion. */ Energy _mAbc; /** * The pole mass for the axial-vector form factors for the \f$b\to s\f$ transitiion. */ Energy _mAbs; /** * The pole mass for the axial-vector form factors for the \f$c\to s\f$ transitiion. */ Energy _mAcs; /** * The pole mass for the axial-vector form factors for the \f$b\to d\f$ transitiion. */ Energy _mAbd; /** * The pole mass for the axial-vector form factors for the \f$c\to u\f$ transitiion. */ Energy _mAcu; }; } #endif /* HERWIG_ChengHeavyBaryonFormFactor_H */ diff --git a/Decay/FormFactors/CzyzNucleonFormFactor.cc b/Decay/FormFactors/CzyzNucleonFormFactor.cc --- a/Decay/FormFactors/CzyzNucleonFormFactor.cc +++ b/Decay/FormFactors/CzyzNucleonFormFactor.cc @@ -1,319 +1,339 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the CzyzNucleonFormFactor class. // #include "CzyzNucleonFormFactor.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Decay/ResonanceHelpers.h" using namespace Herwig; CzyzNucleonFormFactor::CzyzNucleonFormFactor() { // masses and widths of rho resonances rhoMasses_ = {775.49*MeV, 1465*MeV, 1720*MeV, 2.12*GeV,2.32647*GeV}; rhoWidths_ = {149.10*MeV, 400*MeV, 250*MeV, 0.3 *GeV,0.4473 *GeV}; // masses and width of omega resonances omegaMasses_ = {782.65*MeV, 1425*MeV, 1670*MeV, 2.0707 *GeV, 2.34795 *GeV}; omegaWidths_ = {8.49 *MeV, 215*MeV, 315*MeV, 1.03535*GeV, 1.173975*GeV}; // c_1 couplings c1Re_ = {1.,-0.467,-0.177, 0.301}; c1Im_ = {0.,-0.385, 0.149, 0.264}; // c_2 couplings c2Re_ = {1., 0.0521, -0.00308,-0.348}; c2Im_ = {0.,-3.04, 2.38,-0.104}; // c_3 couplings c3Re_ = {1.,-7.88,10.2}; c3Im_ = {0., 5.67, -1.94}; // c_4 couplings c4Re_ = {1.,-0.832, 0.405}; c4Im_ = {0., 0.308,-0.25}; // Magnetic moments mup_ = 2.792847356; mun_ = -1.9130427; // set up the form factors addFormFactor(2212,2212,2,2,2,1,1,1); addFormFactor(2112,2112,2,2,2,1,1,1); initialModes(numberOfFactors()); } IBPtr CzyzNucleonFormFactor::clone() const { return new_ptr(*this); } IBPtr CzyzNucleonFormFactor::fullclone() const { return new_ptr(*this); } void CzyzNucleonFormFactor::doinit() { BaryonFormFactor::doinit(); static const Complex ii(0.,1.); // calculate c_1 c1_.clear(); assert(c1Re_.size()==4 && c1Im_.size()==4); // c_1 1 -> 4 complex fact(ZERO); for(unsigned int ix=0;ix<4;++ix) { c1_.push_back(c1Re_[ix]+ii*c1Im_[ix]); fact += c1_[ix]*sqr(omegaMasses_[ix]); } c1_.push_back(-fact/sqr(omegaMasses_[4])); // calculate c_2 c2_.clear(); assert(c2Re_.size()==4 && c2Im_.size()==4); // c_2 1 -> 4 fact = ZERO; for(unsigned int ix=0;ix<4;++ix) { c2_.push_back(c2Re_[ix]+ii*c2Im_[ix]); fact += c2_[ix]*sqr(rhoMasses_[ix]); } c2_.push_back(-fact/sqr(rhoMasses_[4])); // calculate c_3 c3_.clear(); assert(c3Re_.size()==3 && c3Im_.size()==3); // c_3 1 -> 4 fact = ZERO; complex fact2(ZERO); for(unsigned int ix=0;ix<3;++ix) { c3_.push_back(c3Re_[ix]+ii*c3Im_[ix]); fact += c3_[ix]*sqr(omegaMasses_[ix]); fact2 += c3_[ix]*sqr(omegaMasses_[ix])* (sqr(omegaMasses_[ix])-sqr(omegaMasses_[4]) +ii*(omegaMasses_[4]*omegaWidths_[4]-omegaMasses_[ix]*omegaWidths_[ix])); } c3_.push_back(fact2/sqr(omegaMasses_[3])/ (sqr(omegaMasses_[4])-sqr(omegaMasses_[3]) +ii*(omegaMasses_[3]*omegaWidths_[3]-omegaMasses_[4]*omegaWidths_[4])) ); fact += c3_[3]*sqr(omegaMasses_[3]); c3_.push_back(-fact/sqr(omegaMasses_[4])); // c_4 1 -> 4 fact = ZERO; fact2 = ZERO; for(unsigned int ix=0;ix<3;++ix) { c4_.push_back(c4Re_[ix]+ii*c4Im_[ix]); fact += c4_[ix]*sqr(rhoMasses_[ix]); fact2 += c4_[ix]*sqr(rhoMasses_[ix])* (sqr(rhoMasses_[ix])-sqr(rhoMasses_[4]) +ii*(rhoMasses_[4]*rhoWidths_[4]-rhoMasses_[ix]*rhoWidths_[ix])); } c4_.push_back(fact2/sqr(rhoMasses_[3])/ (sqr(rhoMasses_[4])-sqr(rhoMasses_[3]) +ii*(rhoMasses_[3]*rhoWidths_[3]-rhoMasses_[4]*rhoWidths_[4])) ); fact += c4_[3]*sqr(rhoMasses_[3]); c4_.push_back(-fact/sqr(rhoMasses_[4])); // a and b parameters a_ = mup_-mun_-1; b_ = -mup_-mun_+1.; } void CzyzNucleonFormFactor::persistentOutput(PersistentOStream & os) const { os << ounit(rhoMasses_,GeV) << ounit(rhoWidths_,GeV) << ounit(omegaMasses_,GeV) << ounit(omegaWidths_,GeV) << c1Re_ << c1Im_ << c2Re_ << c2Im_ << c3Re_ << c3Im_ << c4Re_ << c4Im_ << c1_ << c2_ << c3_ << c4_ << mup_ << mun_ << a_ << b_; } void CzyzNucleonFormFactor::persistentInput(PersistentIStream & is, int) { is >> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV) >> iunit(omegaMasses_,GeV) >> iunit(omegaWidths_,GeV) >> c1Re_ >> c1Im_ >> c2Re_ >> c2Im_ >> c3Re_ >> c3Im_ >> c4Re_ >> c4Im_ >> c1_ >> c2_ >> c3_ >> c4_ >> mup_ >> mun_ >> a_ >> b_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigCzyzNucleonFormFactor("Herwig::CzyzNucleonFormFactor", "HwFormFactors.so"); void CzyzNucleonFormFactor::Init() { static ClassDocumentation documentation ("The CzyzNucleonFormFactor class implements the model of " "Phys.Rev. D90 (2014) no.11, 114021 for the nucleon form factor", "The nucleon form factor model of \\cite{Czyz:2014sha} was used", "\\bibitem{Czyz:2014sha}\n" "H.~Czyż, J.~H.~Kühn and S.~Tracz,\n" "%``Nucleon form factors and final state radiative corrections to $e^+e^- \\to p\\bar{p}γ$,''\n" "Phys.\\ Rev.\\ D {\\bf 90} (2014) no.11, 114021\n" "doi:10.1103/PhysRevD.90.114021\n" "[arXiv:1407.7995 [hep-ph]].\n" "%%CITATION = doi:10.1103/PhysRevD.90.114021;%%\n" "%5 citations counted in INSPIRE as of 25 Aug 2018\n"); static ParVector interfaceRhoMasses ("RhoMasses", "The masses of the rho mesons for the form factor", &CzyzNucleonFormFactor::rhoMasses_, GeV, 5, 1.0*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); static ParVector interfaceRhoWidths ("RhoWidths", "The widths of the rho mesons for the form factor", &CzyzNucleonFormFactor::rhoWidths_, GeV, 5, 1.0*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); static ParVector interfaceOmegaMasses ("OmegaMasses", "The masses of the omega mesons for the form factor", &CzyzNucleonFormFactor::omegaMasses_, GeV, 5, 1.0*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); static ParVector interfaceOmegaWidths ("OmegaWidths", "The widths of the omega mesons for the form factor", &CzyzNucleonFormFactor::omegaWidths_, GeV, 5, 1.0*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); static ParVector interfacec1Real ("c1Real", "The real part of the c_1 coupling", &CzyzNucleonFormFactor::c1Re_, 4, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec1Imag ("c1Imag", "The imaginary part of the c_1 coupling", &CzyzNucleonFormFactor::c1Im_, 4, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec2Real ("c2Real", "The real part of the c_2 coupling", &CzyzNucleonFormFactor::c2Re_, 4, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec2Imag ("c2Imag", "The imaginary part of the c_2 coupling", &CzyzNucleonFormFactor::c2Im_, 4, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec3Real ("c3Real", "The real part of the c_3 coupling", &CzyzNucleonFormFactor::c3Re_, 3, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec3Imag ("c3Imag", "The imaginary part of the c_3 coupling", &CzyzNucleonFormFactor::c3Im_, 3, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec4Real ("c4Real", "The real part of the c_4 coupling", &CzyzNucleonFormFactor::c4Re_, 3, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec4Imag ("c4Imag", "The imaginary part of the c_4 coupling", &CzyzNucleonFormFactor::c4Im_, 3, 1., -100., 100.0, false, false, Interface::limited); static Parameter interfaceMuProton ("MuProton", "The proton magnetic moment", &CzyzNucleonFormFactor::mup_, 2.792, 0.0, 10.0, false, false, Interface::limited); static Parameter interfaceMuNeutron ("MuNeutron", "The proton magnetic moment", &CzyzNucleonFormFactor::mun_, -1.913, 0.0, 10.0, false, false, Interface::limited); } void CzyzNucleonFormFactor:: SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc, int id0,int id1,Energy,Energy, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo flavour, Virtuality virt) { + f1a = f2a = f3a = f1v = f2v = f3v = 0.; assert(abs(id0)==abs(id1)); if(iloc==0) assert(abs(id0)==2212); else assert(abs(id0)==2112); assert(virt==TimeLike); + if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return; + if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return; + if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return; + bool I0 = true, I1 = true; + if(flavour.I!=IsoSpin::IUnknown) { + if(flavour.I==IsoSpin::IZero) { + I1=false; + if(flavour.I3!=IsoSpin::I3Unknown and flavour.I3!=IsoSpin::I3Zero) return; + } + else if(flavour.I==IsoSpin::IOne) { + I0=false; + if(flavour.I3!=IsoSpin::I3Unknown and flavour.I3!=IsoSpin::I3Zero) return; + } + } + if(flavour.I3!=IsoSpin::I3Unknown and flavour.I3!=IsoSpin::I3Zero) return; // calculate the form factors Complex F1S(0.),F1V(0.),F2S(0.),F2V(0.); Complex n1(0.),n2(0.),n3(0.),n4(0.); for(unsigned int ix=0;ix<5;++ix) { - F1S += c1_[ix]*Resonance::BreitWignerFW(q2,omegaMasses_[ix],omegaWidths_[ix]); - F1V += c2_[ix]*Resonance::BreitWignerFW(q2, rhoMasses_[ix], rhoWidths_[ix]); - F2S += c3_[ix]*Resonance::BreitWignerFW(q2,omegaMasses_[ix],omegaWidths_[ix]); - F2V += c4_[ix]*Resonance::BreitWignerFW(q2, rhoMasses_[ix], rhoWidths_[ix]); + if(I0) { + F1S += c1_[ix]*Resonance::BreitWignerFW(q2,omegaMasses_[ix],omegaWidths_[ix]); + F2S += c3_[ix]*Resonance::BreitWignerFW(q2,omegaMasses_[ix],omegaWidths_[ix]); + } + if(I1) { + F1V += c2_[ix]*Resonance::BreitWignerFW(q2, rhoMasses_[ix], rhoWidths_[ix]); + F2V += c4_[ix]*Resonance::BreitWignerFW(q2, rhoMasses_[ix], rhoWidths_[ix]); + } n1 += c1_[ix]; n2 += c2_[ix]; n3 += c3_[ix]; n4 += c4_[ix]; } F1S *= 0.5 /n1; F1V *= 0.5 /n2; F2S *= -0.5*b_/n3; F2V *= 0.5*a_/n4; - f1a = f2a = f3v = f3a = 0.; if(iloc==0) { f1v = F1S + F1V; f2v = F2S + F2V; } else { f1v = F1S - F1V; f2v = F2S - F2V; } } void CzyzNucleonFormFactor:: dataBaseOutput(ofstream& output,bool header, bool create) const { if(header) output << "update decayers set parameters=\""; if(create) output << "create Herwig::CzyzNucleonFormFactor " << name() << " \n"; for(unsigned int ix=0;ix<5;++ix) { output << "newdef " << name() << ":RhoMasses " << ix << " " << rhoMasses_[ix]/GeV << "\n"; output << "newdef " << name() << ":RhoWidths " << ix << " " << rhoWidths_[ix]/GeV << "\n"; output << "newdef " << name() << ":OmegaMasses " << ix << " " << omegaMasses_[ix]/GeV << "\n"; output << "newdef " << name() << ":OmegaWidths " << ix << " " << omegaWidths_[ix]/GeV << "\n"; } for(unsigned int ix=0;ix<4;++ix) { output << "newdef " << name() << ":c1Real " << ix << " " << c1Re_[ix] << "\n"; output << "newdef " << name() << ":c1Imag " << ix << " " << c1Im_[ix] << "\n"; output << "newdef " << name() << ":c2Real " << ix << " " << c2Re_[ix] << "\n"; output << "newdef " << name() << ":c2Imag " << ix << " " << c2Im_[ix] << "\n"; } for(unsigned int ix=0;ix<3;++ix) { output << "newdef " << name() << ":c3Real " << ix << " " << c3Re_[ix] << "\n"; output << "newdef " << name() << ":c3Imag " << ix << " " << c3Im_[ix] << "\n"; output << "newdef " << name() << ":c4Real " << ix << " " << c4Re_[ix] << "\n"; output << "newdef " << name() << ":c4Imag " << ix << " " << c4Im_[ix] << "\n"; } output << "newdef " << name() << ":MuProton " << mup_ << "\n"; output << "newdef " << name() << ":MuNeutron " << mun_ << "\n"; BaryonFormFactor::dataBaseOutput(output,false,false); if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Decay/FormFactors/CzyzNucleonFormFactor.h b/Decay/FormFactors/CzyzNucleonFormFactor.h --- a/Decay/FormFactors/CzyzNucleonFormFactor.h +++ b/Decay/FormFactors/CzyzNucleonFormFactor.h @@ -1,226 +1,227 @@ // -*- C++ -*- #ifndef Herwig_CzyzNucleonFormFactor_H #define Herwig_CzyzNucleonFormFactor_H // // This is the declaration of the CzyzNucleonFormFactor class. // #include "BaryonFormFactor.h" namespace Herwig { using namespace ThePEG; /** * The CzyzNucleonFormFactor class implements the model of * Phys.Rev. D90 (2014) no.11, 114021 for the proton and neutron form factors * * @see \ref CzyzNucleonFormFactorInterfaces "The interfaces" * defined for CzyzNucleonFormFactor. */ class CzyzNucleonFormFactor: public BaryonFormFactor { public: /** * The default constructor. */ CzyzNucleonFormFactor(); /** * The form factor for the weak decay of a spin \f$\frac12\f$ baryon to a * spin \f$\frac12\f$ baryon. * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form factor list. * @param id0 The PDG code of the incoming baryon. * @param id1 The PDG code of the outgoing baryon. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param f1v The form factor \f$F^V_1\f$. * @param f2v The form factor \f$F^V_2\f$. * @param f3v The form factor \f$F^V_3\f$. * @param f1a The form factor \f$F^A_1\f$. * @param f2a The form factor \f$F^A_2\f$. * @param f3a The form factor \f$F^A_3\f$. * @param virt Whether the q2 is space or timelike */ virtual void SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo flavour, Virtuality virt=SpaceLike); /** * 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 * @param create Whether or not to add a statement creating the object */ virtual void dataBaseOutput(ofstream & os,bool header,bool create) 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: /** * 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. */ CzyzNucleonFormFactor & operator=(const CzyzNucleonFormFactor &) = delete; private: /** * The masses of the \f$\rho\f$ resonances */ vector rhoMasses_; /** * The widths of the \f$\rho\f$ resonances */ vector rhoWidths_; /** * The masses of the \f$\omega\f$ resonances */ vector omegaMasses_; /** * The widths of the \f$\omega\f$ resonances */ vector omegaWidths_; /** * The real part of the \f$c_1\f$ couplings */ vector c1Re_; /** * The imaginary part of the \f$c_1\f$ couplings */ vector c1Im_; /** * The real part of the \f$c_2\f$ couplings */ vector c2Re_; /** * The imaginary part of the \f$c_2\f$ couplings */ vector c2Im_; /** * The real part of the \f$c_3\f$ couplings */ vector c3Re_; /** * The imaginary part of the \f$c_3\f$ couplings */ vector c3Im_; /** * The real part of the \f$c_4\f$ couplings */ vector c4Re_; /** * The imaginary part of the \f$c_4\f$ couplings */ vector c4Im_; /** * The complex \f$c_1\f$ coupling */ vector c1_; /** * The complex \f$c_2\f$ coupling */ vector c2_; /** * The complex \f$c_3\f$ coupling */ vector c3_; /** * The complex \f$c_4\f$ coupling */ vector c4_; /** * Proton magentic moment */ double mup_; /** * Neutron magentic moment */ double mun_; /** * \f$a\f$ parameter */ double a_; /** * \f$b\f$ parameter */ double b_; }; } #endif /* Herwig_CzyzNucleonFormFactor_H */ diff --git a/Decay/FormFactors/LambdabExcitedLambdacSumRuleFormFactor.cc b/Decay/FormFactors/LambdabExcitedLambdacSumRuleFormFactor.cc --- a/Decay/FormFactors/LambdabExcitedLambdacSumRuleFormFactor.cc +++ b/Decay/FormFactors/LambdabExcitedLambdacSumRuleFormFactor.cc @@ -1,146 +1,148 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the LambdabExcitedLambdacSumRuleFormFactor class. // #include "LambdabExcitedLambdacSumRuleFormFactor.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; using namespace ThePEG; LambdabExcitedLambdacSumRuleFormFactor::LambdabExcitedLambdacSumRuleFormFactor() { _xi1=0.29; _rho2=2.01; // modes handled by this form-factor // lambda_b to lambda_c1 addFormFactor(5122,14122,2,2,1,2,5,4); // lambda_b to lambda_c1* addFormFactor(5122,4124 ,2,4,1,2,5,4); } void LambdabExcitedLambdacSumRuleFormFactor:: persistentOutput(PersistentOStream & os) const { os << _xi1 << _rho2; } void LambdabExcitedLambdacSumRuleFormFactor:: persistentInput(PersistentIStream & is, int) { is >> _xi1 >> _rho2; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigLambdabExcitedLambdacSumRuleFormFactor("Herwig::LambdabExcitedLambdacSumRuleFormFactor", "HwFormFactors.so"); void LambdabExcitedLambdacSumRuleFormFactor::Init() { static ClassDocumentation documentation ("The LambdabExcitedLambdacSumRuleFormFactor class implements the" " form-factors for Lambda_b to Lambda_c1(*) from hep-ph/0012114.", "Lambda_b to Lambda_c1(*) used the formfactors from \\cite{Huang:2000xw}.", "%\\cite{Huang:2000xw}\n" "\\bibitem{Huang:2000xw}\n" " M.~Q.~Huang, J.~P.~Lee, C.~Liu and H.~S.~Song,\n" " %``Leading Isgur-Wise form factor of Lambda/b to Lambda/c1 transition using\n" " %QCD sum rules,''\n" " Phys.\\ Lett.\\ B {\\bf 502}, 133 (2001)\n" " [arXiv:hep-ph/0012114].\n" " %%CITATION = PHLTA,B502,133;%%\n" ); static Parameter interfaceXi ("Xi", "The intercept for the Isgur-Wise form-factor", &LambdabExcitedLambdacSumRuleFormFactor::_xi1, 0.29, 0.0, 10.0, false, false, true); static Parameter interfaceRho2 ("Rho2", "The slope parameter for the form-factor.", &LambdabExcitedLambdacSumRuleFormFactor::_rho2, 2.01, -10.0, 10.0, false, false, true); } void LambdabExcitedLambdacSumRuleFormFactor:: SpinHalfSpinHalfFormFactor(Energy2 q2,int,int,int,Energy m0,Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo , Virtuality virt) { assert(virt==SpaceLike); useMe(); double omega(.5/m0/m1*(m0*m0+m1*m1-q2)),orr(1./sqrt(3.)); // the universal form-factor double xi=_xi1*(1.-_rho2*(omega-1.)); // the couplings in the velocity form Complex g1v,g1a,g2v,g2a,g3a,g3v; g1v = orr*(omega-1.)*xi; g1a = orr*(omega+1.)*xi; g2v =-2.*orr*xi; g3v = 0.; g2a =-2.*orr*xi; g3a = 0.; // convert to our form f1a = g1v-0.5*(m0-m1)*(g2v/m0+g3v/m1); f1v =-g1a-0.5*(m0+m1)*(g2a/m0+g3a/m1); f2a = Complex(0.5*(m0+m1)*( g2v/m0+g3v/m1)); f2v =-Complex(0.5*(m0+m1)*( g2a/m0+g3a/m1)); f3a = Complex(0.5*(m0+m1)*( g2v/m0-g3v/m1)); f3v = Complex(0.5*(m0+m1)*(-g2a/m0+g3a/m1)); } void LambdabExcitedLambdacSumRuleFormFactor:: SpinHalfSpinThreeHalfFormFactor(Energy2 q2,int,int,int,Energy m0,Energy m1, Complex & f1v,Complex & f2v, Complex & f3v,Complex & f4v, Complex & f1a,Complex & f2a, Complex & f3a,Complex & f4a, + FlavourInfo , Virtuality virt) { assert(virt==SpaceLike); useMe(); // the omega value double omega(.5/m0/m1*(m0*m0+m1*m1-q2)); // the universal form-factor double xi(_xi1*(1.-_rho2*(omega-1.))); // calculate the form factor // the couplings in the velocity form Complex N1,N2,N3,N4,K1,K2,K3,K4; Energy msum(m0+m1);Energy2 msum2(msum*msum); // in the form of the heavy quark papers N1 = xi; K1 = xi; N2 = 0.; K2 = 0.; N3 = 0.; K3 = 0.; N4 = 0.; K4 = 0.; // convert to our form f1v =-N4; f1a = K4; f2v =-N1*msum/m0; f2a = K1*msum/m0; f3v =-Complex(msum2/m0*(N2/m0+N3/m1)); f3a = Complex(msum2/m0*(K2/m0+K3/m1)); f4v =-Complex(msum2/m0/m0*N2); f4a = Complex(msum2/m0/m0*K2); } void LambdabExcitedLambdacSumRuleFormFactor::dataBaseOutput(ofstream & output, bool header, bool create) const { if(header) output << "update decayers set parameters=\""; if(create) output << "create Herwig::LambdabExcitedLambdacSumRuleFormFactor " << name() << " \n"; output << "newdef " << name() << ":Xi " << _xi1 << " \n"; output << "newdef " << name() << ":Rho2 " << _rho2 << " \n"; BaryonFormFactor::dataBaseOutput(output,false,false); if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Decay/FormFactors/LambdabExcitedLambdacSumRuleFormFactor.h b/Decay/FormFactors/LambdabExcitedLambdacSumRuleFormFactor.h --- a/Decay/FormFactors/LambdabExcitedLambdacSumRuleFormFactor.h +++ b/Decay/FormFactors/LambdabExcitedLambdacSumRuleFormFactor.h @@ -1,151 +1,153 @@ // -*- C++ -*- #ifndef HERWIG_LambdabExcitedLambdacSumRuleFormFactor_H #define HERWIG_LambdabExcitedLambdacSumRuleFormFactor_H // This is the declaration of the LambdabExcitedLambdacSumRuleFormFactor class. #include "BaryonFormFactor.h" namespace Herwig { using namespace ThePEG; /** \ingroup Decay * * The LambdabExcitedLambdacSumRuleFormFactor class implements the * form-factors of hep-ph/0012114 for the decay of the \f$Lambda_b^0\f$ * to excited \f$\Lambda^+_c\f$ states. * * @see BaryonFormFactor. * */ class LambdabExcitedLambdacSumRuleFormFactor: public BaryonFormFactor { public: /** * Default constructor */ LambdabExcitedLambdacSumRuleFormFactor(); /** @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); //@} /** * Standard Init function used to initialize the interfaces. */ 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 {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);} //@} public: /** @name Form Factors */ //@{ /** * The form factor for the weak decay of a spin \f$\frac12\f$ baryon to a * spin \f$\frac12\f$ baryon. * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form factor list. * @param id0 The PDG code of the incoming baryon. * @param id1 The PDG code of the outgoing baryon. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param f1v The form factor \f$F^V_1\f$. * @param f2v The form factor \f$F^V_2\f$. * @param f3v The form factor \f$F^V_3\f$. * @param f1a The form factor \f$F^A_1\f$. * @param f2a The form factor \f$F^A_2\f$. * @param f3a The form factor \f$F^A_3\f$. */ virtual void SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo flavour, Virtuality virt=SpaceLike); /** * The form factor for the weak decay of a spin \f$\frac12\f$ baryon to a * spin \f$\frac32\f$ baryon. * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form factor list. * @param id0 The PDG code of the incoming baryon. * @param id1 The PDG code of the outgoing baryon. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param g1v The form factor \f$G^V_1\f$. * @param g2v The form factor \f$G^V_2\f$. * @param g3v The form factor \f$G^V_3\f$. * @param g4v The form factor \f$G^V_4\f$. * @param g1a The form factor \f$G^A_1\f$. * @param g2a The form factor \f$G^A_2\f$. * @param g3a The form factor \f$G^A_3\f$. * @param g4a The form factor \f$G^A_4\f$. * @param virt Whether the q2 is space or timelike */ virtual void SpinHalfSpinThreeHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & g1v,Complex & g2v,Complex & g3v, Complex & g4v,Complex & g1a,Complex & g2a, Complex & g3a,Complex & g4a, + FlavourInfo flavour, Virtuality virt=SpaceLike); //@} /** * 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 * @param create Whether or not to add a statement creating the object */ virtual void dataBaseOutput(ofstream & os,bool header,bool create) const; private: /** * Private and non-existent assignment operator. */ LambdabExcitedLambdacSumRuleFormFactor & operator=(const LambdabExcitedLambdacSumRuleFormFactor &) = delete; private: /** * The intercept for the form factor. */ double _xi1; /** * The slope for the form factor. */ double _rho2; }; } #endif /* HERWIG_LambdabExcitedLambdacSumRuleFormFactor_H */ diff --git a/Decay/FormFactors/LightBaryonQuarkModelFormFactor.cc b/Decay/FormFactors/LightBaryonQuarkModelFormFactor.cc --- a/Decay/FormFactors/LightBaryonQuarkModelFormFactor.cc +++ b/Decay/FormFactors/LightBaryonQuarkModelFormFactor.cc @@ -1,252 +1,253 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the LightBaryonQuarkModelFormFactor class. // #include "LightBaryonQuarkModelFormFactor.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" using namespace Herwig; void LightBaryonQuarkModelFormFactor::doinit() { BaryonFormFactor::doinit(); // check that the parameters are consistent unsigned int isize=numberOfFactors(); if(isize!=_f1.size()||isize!=_f2.size()||isize!=_g1.size()||isize!=_g2.size()|| isize!=_Lambdaf1.size()||isize!=_Lambdaf2.size()|| isize!=_Lambdag1.size()||isize!=_Lambdag2.size()) throw InitException() << "Inconsistent parameters in " << "LightBaryonQuarkModelFormFactor::doinit()" << Exception::abortnow; } LightBaryonQuarkModelFormFactor::LightBaryonQuarkModelFormFactor() { // the various decay modes handled by the model and the parameters // neutron to proton addFormFactor(2112,2212,2,2,2,1,1,2); _f1.push_back(1.00);_f2.push_back(1.81/GeV); _g1.push_back(1.25);_g2.push_back(ZERO); _Lambdaf1.push_back(0.69*GeV);_Lambdaf2.push_back(0.96*GeV); _Lambdag1.push_back(0.76*GeV);_Lambdag2.push_back(1.04*GeV); // sigma+ to lambda addFormFactor(3222,3122,2,2,3,2,2,1); _f1.push_back(0.00);_f2.push_back(1.04/GeV); _g1.push_back(0.60);_g2.push_back(ZERO); _Lambdaf1.push_back(-0.32*GeV);_Lambdaf2.push_back(-1.72*GeV); _Lambdag1.push_back( 0.77*GeV);_Lambdag2.push_back( 1.05*GeV); // sigma- to lambda addFormFactor(3112,3122,2,2,3,1,1,2); _f1.push_back(0.00);_f2.push_back(1.04/GeV); _g1.push_back(0.60);_g2.push_back(ZERO); _Lambdaf1.push_back(-0.32*GeV);_Lambdaf2.push_back(-1.72*GeV); _Lambdag1.push_back( 0.77*GeV);_Lambdag2.push_back( 1.05*GeV); // sigma- to sigma0 addFormFactor(3112,3212,2,2,3,1,1,2); _f1.push_back(1.41);_f2.push_back(0.76/GeV); _g1.push_back(0.69);_g2.push_back(ZERO); _Lambdaf1.push_back(0.60*GeV);_Lambdaf2.push_back(0.81*GeV); _Lambdag1.push_back(0.77*GeV);_Lambdag2.push_back(1.04*GeV); // sigma0 to sigma+ addFormFactor(3212,3222,2,2,3,2,1,2); _f1.push_back(-1.41);_f2.push_back(-0.76/GeV); _g1.push_back(-0.69);_g2.push_back( ZERO); _Lambdaf1.push_back(0.60*GeV);_Lambdaf2.push_back(0.81*GeV); _Lambdag1.push_back(0.77*GeV);_Lambdag2.push_back(1.04*GeV); // Xi- to Xi0 addFormFactor(3312,3322,2,2,3,3,1,2); _f1.push_back(-1.00);_f2.push_back(0.73/GeV); _g1.push_back( 0.24);_g2.push_back(ZERO); _Lambdaf1.push_back(0.56*GeV);_Lambdaf2.push_back(0.71*GeV); _Lambdag1.push_back(0.76*GeV);_Lambdag2.push_back(1.04*GeV); // lambda to proton addFormFactor(3122,2212,2,2,1,2,3,2); _f1.push_back(-1.19);_f2.push_back(-0.850/GeV); _g1.push_back(-0.99);_g2.push_back(-0.025/GeV); _Lambdaf1.push_back(0.71*GeV);_Lambdaf2.push_back(0.98*GeV); _Lambdag1.push_back(0.81*GeV);_Lambdag2.push_back(1.12*GeV); // sigma0 to proton addFormFactor(3212,2212,2,2,2,1,3,2); _f1.push_back(-0.69);_f2.push_back(0.44/GeV); _g1.push_back( 0.19);_g2.push_back(0.0043/GeV); _Lambdaf1.push_back(0.64*GeV);_Lambdaf2.push_back(0.84*GeV); _Lambdag1.push_back(0.83*GeV);_Lambdag2.push_back(1.16*GeV); // sigma- to neutron addFormFactor(3112,2112,2,2,1,1,3,2); _f1.push_back(-0.97);_f2.push_back(0.62/GeV); _g1.push_back(0.27);_g2.push_back(0.0061/GeV); _Lambdaf1.push_back(0.64*GeV);_Lambdaf2.push_back(0.90*GeV); _Lambdag1.push_back(0.83*GeV);_Lambdag2.push_back(1.16*GeV); // xi- to lambda addFormFactor(3312,3122,2,2,3,1,3,2); _f1.push_back(1.19);_f2.push_back(0.07/GeV); _g1.push_back(0.33);_g2.push_back(0.0076/GeV); _Lambdaf1.push_back(0.68*GeV);_Lambdaf2.push_back(0.89*GeV); _Lambdag1.push_back(0.81*GeV);_Lambdag2.push_back(1.10*GeV); // xi- to sigma0 addFormFactor(3312,3212,2,2,3,1,3,2); _f1.push_back(0.69);_f2.push_back(0.98/GeV); _g1.push_back(0.94);_g2.push_back(0.022/GeV); _Lambdaf1.push_back(0.75*GeV);_Lambdaf2.push_back(1.05*GeV); _Lambdag1.push_back(0.81*GeV);_Lambdag2.push_back(1.12*GeV); // xi0 to sigma+ addFormFactor(3322,3222,2,2,3,2,3,2); _f1.push_back(0.98);_f2.push_back(1.38/GeV); _g1.push_back(1.33);_g2.push_back(0.031/GeV); _Lambdaf1.push_back(0.75*GeV);_Lambdaf2.push_back(1.05*GeV); _Lambdag1.push_back(0.81*GeV);_Lambdag2.push_back(1.12*GeV); // set the inital number of form factors initialModes(numberOfFactors()); } void LightBaryonQuarkModelFormFactor::persistentOutput(PersistentOStream & os) const { os << _f1 << ounit(_f2,1/GeV) << _g1 << ounit(_g2,1/GeV) << ounit(_Lambdaf1,GeV) << ounit(_Lambdaf2,GeV) << ounit(_Lambdag1,GeV) << ounit(_Lambdag2,GeV); } void LightBaryonQuarkModelFormFactor::persistentInput(PersistentIStream & is, int) { is >> _f1 >> iunit(_f2,1/GeV) >> _g1 >> iunit(_g2,1/GeV) >> iunit(_Lambdaf1,GeV) >> iunit(_Lambdaf2,GeV) >> iunit(_Lambdag1,GeV) >> iunit(_Lambdag2,GeV); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigLightBaryonQuarkModelFormFactor("Herwig::LightBaryonQuarkModelFormFactor", "HwFormFactors.so"); void LightBaryonQuarkModelFormFactor::Init() { static ClassDocumentation documentation ("The LightBaryonQuarkModelFormFactor class implements" " the quark model calculation of hep-ph/9409272 for the form-factors" " for the light quarks", "The quark model calculation of \\cite{Schlumpf:1994fb} was used" "for the weak decay of the light baryons", "\\bibitem{Schlumpf:1994fb}\n" "F.~Schlumpf,\n" "Phys.\\ Rev.\\ D {\\bf 51} (1995) 2262 [arXiv:hep-ph/9409272].\n" "%%CITATION = PHRVA,D51,2262;%%\n"); static ParVector interfacef1 ("f1", "The form-factor f1 at zero q^2", &LightBaryonQuarkModelFormFactor::_f1, 0, 0, 0, -10., 10., false, false, true); static ParVector interfaceg1 ("g1", "The form-factor g1 at zero q^2", &LightBaryonQuarkModelFormFactor::_g1, 0, 0, 0, -10., 10., false, false, true); static ParVector interfacef2 ("f2", "The form-factor f2 at zero q^2", &LightBaryonQuarkModelFormFactor::_f2, 1./GeV, 0, ZERO, -10./GeV, 10./GeV, false, false, true); static ParVector interfaceg2 ("g2", "The form-factor g2 at zero q^2", &LightBaryonQuarkModelFormFactor::_g2, 1./GeV, 0, ZERO, -10./GeV, 10./GeV, false, false, true); static ParVector interfaceLambdaf1 ("Lambdaf1", "The first mass for the energy dependence of the f1 form-factor.", &LightBaryonQuarkModelFormFactor::_Lambdaf1, 1.*GeV, 0, ZERO, -10.*GeV, 10.*GeV, false, false, true); static ParVector interfaceLambdaf2 ("Lambdaf2", "The second mass for the energy dependence of the f1 form-factor.", &LightBaryonQuarkModelFormFactor::_Lambdaf2, 1.*GeV, 0, ZERO, -10.*GeV, 10.*GeV, false, false, true); static ParVector interfaceLambdag1 ("Lambdag1", "The first mass for the energy dependence of the g1 form-factor.", &LightBaryonQuarkModelFormFactor::_Lambdag1, 1.*GeV, 0, ZERO, -10.*GeV, 10.*GeV, false, false, true); static ParVector interfaceLambdag2 ("Lambdag2", "The second mass for the energy dependence of the g1 form-factor.", &LightBaryonQuarkModelFormFactor::_Lambdag2, 1.*GeV, 0, ZERO, -10.*GeV, 10.*GeV, false, false, true); } // form factor for spin-1/2 to spin-1/2 void LightBaryonQuarkModelFormFactor:: SpinHalfSpinHalfFormFactor(Energy2 q2,int mode,int, int, Energy m0, Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo , Virtuality virt) { assert(virt==SpaceLike); useMe(); // f_3 is zero f3v = 0.; f3a = 0.; // energy dependence of the f1 and g1 factors InvEnergy2 lam1,lam2; if(_Lambdaf1[mode]>ZERO) { lam1 = 1./sqr(_Lambdaf1[mode]); lam2 = 1./sqr(_Lambdaf2[mode]); f1v= _f1[mode]/(1.-q2*(lam1-q2*sqr(lam2))); } else { f1v = _f1[mode]*(1.+_Lambdaf1[mode]/GeV*q2/GeV2+_Lambdaf2[mode]/GeV*sqr(q2/GeV2)); } lam1 = 1./sqr(_Lambdag1[mode]); lam2 = 1./sqr(_Lambdag2[mode]); f1a=-_g1[mode]/(1.-q2*(lam1-q2*sqr(lam2))); // the f2 and g2 factors f2v =-(m0+m1)*_f2[mode]; f2a = (m0+m1)*_g2[mode]; } void LightBaryonQuarkModelFormFactor::dataBaseOutput(ofstream& output,bool header, bool create) const { if(header) output << "update decayers set parameters=\""; if(create) output << "create Herwig::LightBaryonQuarkModelFormFactor " << name() << " \n"; for(unsigned int ix=0;ix<_f1.size();++ix) { if(ixLightBaryonQuarkModelFormFactor class implements the * quark model calculation of hep-ph/9409272 for the form-factors for * the decay of baryons containing the light quarks. This can be adapted for * other models. * * @see BaryonFormFactor * */ class LightBaryonQuarkModelFormFactor: public BaryonFormFactor { public: /** * Default constructor */ LightBaryonQuarkModelFormFactor(); /** @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); //@} /** * Standard Init function used to initialize the interfaces. */ 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 {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);} //@} public: /** @name Form Factors */ //@{ /** * The form factor for the weak decay of a spin \f$\frac12\f$ baryon to a * spin \f$\frac12\f$ baryon. * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form factor list. * @param id0 The PDG code of the incoming baryon. * @param id1 The PDG code of the outgoing baryon. * @param m0 The mass of the incoming baryon. * @param m1 The mass of the outgoing baryon. * @param f1v The form factor \f$F^V_1\f$. * @param f2v The form factor \f$F^V_2\f$. * @param f3v The form factor \f$F^V_3\f$. * @param f1a The form factor \f$F^A_1\f$. * @param f2a The form factor \f$F^A_2\f$. * @param f3a The form factor \f$F^A_3\f$. * @param virt Whether the q2 is space or timelike */ virtual void SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc, int id0, int id1, Energy m0, Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, + FlavourInfo flavour, Virtuality virt=SpaceLike); //@} /** * 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 * @param create Whether or not to add a statement creating the object */ virtual void dataBaseOutput(ofstream & os,bool header,bool create) 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(); //@} private: /** * Private and non-existent assignment operator. */ LightBaryonQuarkModelFormFactor & operator=(const LightBaryonQuarkModelFormFactor &) = delete; private: /** * The form-factor \f$f_1\f$ at \f$q^2=0\f$. */ vector _f1; /** * The form-factor \f$f_2\f$ at \f$q^2=0\f$. */ vector _f2; /** * The form-factor \f$g_1\f$ at \f$q^2=0\f$. */ vector _g1; /** * The form-factor \f$g_2\f$ at \f$q^2=0\f$. */ vector _g2; /** * The mass for the energy dependence of the \f$f_1\f$ form factor. */ vector _Lambdaf1; /** * The mass for the energy dependence of the \f$f_2\f$ form factor. */ vector _Lambdaf2; /** * The mass for the energy dependence of the \f$g_1\f$ form factor. */ vector _Lambdag1; /** * The mass for the energy dependence of the \f$g_2\f$ form factor. */ vector _Lambdag2; }; } #endif /* HERWIG_LightBaryonQuarkModelFormFactor_H */ diff --git a/Decay/FormFactors/SingletonFormFactor.cc b/Decay/FormFactors/SingletonFormFactor.cc --- a/Decay/FormFactors/SingletonFormFactor.cc +++ b/Decay/FormFactors/SingletonFormFactor.cc @@ -1,235 +1,235 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the SingletonFormFactor class. // #include "SingletonFormFactor.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Repository/EventGenerator.h" using namespace Herwig; SingletonFormFactor::SingletonFormFactor() { // the charm and strange quark masses _mcharm = 1.80*GeV; _mstrange = 0.51*GeV; // mixing angles using Constants::pi; _thetalambda = 0.25*pi; _thetasigma = 0.25*pi; _thetaxi = 0.25*pi; _thetaxip = 0.25*pi; // the particles handled and the masses for them // lambda_b addFormFactor(5122,4122,2,2,1,2,5,4);_polemass.push_back(6.0*GeV); // sigma_b addFormFactor(5112,4112,2,2,1,1,5,4);_polemass.push_back(6.0*GeV); addFormFactor(5212,4212,2,2,2,1,5,4);_polemass.push_back(6.0*GeV); addFormFactor(5222,4222,2,2,2,2,5,4);_polemass.push_back(6.0*GeV); // omega_b addFormFactor(5332,4332,2,2,3,3,5,4);_polemass.push_back(6.4*GeV); // xi_b addFormFactor(5232,4232,2,2,2,3,5,4);_polemass.push_back(6.0*GeV); addFormFactor(5132,4132,2,2,1,3,5,4);_polemass.push_back(6.0*GeV); // lambda_c addFormFactor(4122,3122,2,2,1,2,4,3);_polemass.push_back(2.5*GeV); // sigma_c addFormFactor(4112,3112,2,2,1,1,4,3);_polemass.push_back(2.8*GeV); addFormFactor(4212,3212,2,2,2,1,4,3);_polemass.push_back(2.8*GeV); addFormFactor(4222,3222,2,2,2,2,4,3);_polemass.push_back(2.8*GeV); // xi_c addFormFactor(4232,3322,2,2,2,3,4,3);_polemass.push_back(2.8*GeV); addFormFactor(4132,3312,2,2,1,3,4,3);_polemass.push_back(2.8*GeV); // set the inital number of form factors initialModes(numberOfFactors()); } void SingletonFormFactor::doinit() { BaryonFormFactor::doinit(); if(numberOfFactors()!=_polemass.size()) throw InitException() << "Inconsistent parameters in SingletonFormFactor::doinit()" << Exception::abortnow; // calculate the constants for the form-factors int id0,id1; _xi.clear(); _nmM.clear(); _mquark.clear(); for(unsigned int ix=0;ix> iunit(_mcharm,GeV) >> iunit(_mstrange,GeV) >> _thetalambda >> _thetasigma >> _thetaxi >> _thetaxip >> iunit(_polemass,GeV) >> _xi >> _nmM >> iunit(_mquark,GeV); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigSingletonFormFactor("Herwig::SingletonFormFactor", "HwFormFactors.so"); void SingletonFormFactor::Init() { static ClassDocumentation documentation ("The SingletonFormFactor class implements the" " form-factors of PRD43, 2939 for the decay of spin-1/2 baryons" " containing one heavy quark.", "Spin-1/2 baryons with one heavy quark were decayed using " "the form factors in \\cite{Singleton:1990ye}.", "%\\cite{Singleton:1990ye}\n" "\\bibitem{Singleton:1990ye}\n" " R.~L.~Singleton,\n" " %``Semileptonic baryon decays with a heavy quark,''\n" " Phys.\\ Rev.\\ D {\\bf 43} (1991) 2939.\n" " %%CITATION = PHRVA,D43,2939;%%\n" ); static Parameter interfaceCharmMass ("CharmMass", "The mass of the charm quark", &SingletonFormFactor::_mcharm, GeV, 1.8*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceStrangeMass ("StrangeMass", "The mass of the strange quark", &SingletonFormFactor::_mstrange, GeV, 0.51*GeV, ZERO, 10.0*GeV, false, false, true); static Parameter interfaceThetaLambda ("ThetaLambda", "The mixing angle for the Lambda", &SingletonFormFactor::_thetalambda, 0.785398163, 0.0, 6.28318507, false, false, true); static Parameter interfaceThetaSigma ("ThetaSigma", "The mixing angle for the Sigma", &SingletonFormFactor::_thetasigma, 0.785398163, 0.0, 6.28318507, false, false, true); static Parameter interfaceThetaXi ("ThetaXi", "The mixing angle for the Xi", &SingletonFormFactor::_thetaxi, 0.785398163, 0.0, 6.28318507, false, false, true); static Parameter interfaceThetaXiPrime ("ThetaXiPrime", "The mixing angle for the Xi'", &SingletonFormFactor::_thetaxip, 0.785398163, 0.0, 6.28318507, false, false, true); static ParVector interfacePoleMass ("PoleMass", "The mass for the energy dependence of the form-factors.", &SingletonFormFactor::_polemass, 1.*GeV, 0, ZERO, -10.*GeV, 10.*GeV, false, false, true); } // form factor for spin-1/2 to spin-1/2 void SingletonFormFactor:: SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc,int, int, Energy m0, Energy m1, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, - Virtuality virt) { + FlavourInfo, Virtuality virt) { assert(virt==SpaceLike); useMe(); InvEnergy ratio(0.5/m0); // all factors divided by sqrt(4.*m0*m1) normalisation double gbar(m1*_xi[iloc]*_nmM[iloc]/_mquark[iloc]); double abar(_xi[iloc]*_nmM[iloc]); InvEnergy apbar(-ratio*(m1/_mquark[iloc]-1.)*_xi[iloc]*_nmM[iloc]); InvEnergy ambar(apbar); InvEnergy gpbar(-ratio*_nmM[iloc]*((_xi[iloc]*m1/_mquark[iloc]-1.) +0.5*(m1-m0)/_mquark[iloc]*(1.-_xi[iloc]))); InvEnergy gmbar(-ratio*_nmM[iloc]*((_xi[iloc]*m1/_mquark[iloc]-1.) +0.5*(m1+m0)/_mquark[iloc]*(1.-_xi[iloc]))); // energy dependence double ymax(sqr(1.-m1/m0)); double yres(sqr(_polemass[iloc]/m0)); double y(q2/sqr(m0)),efact((ymax-yres)/(y-yres)); f1v = efact*(gbar+(m0+m1)*gpbar); f2v = efact*gpbar*(m0+m1); f3v = efact*gmbar*(m0+m1); f1a = efact*(-abar+(m0-m1)*apbar); f2a =-efact*apbar*(m0+m1); f3a =-efact*ambar*(m0+m1); } void SingletonFormFactor::dataBaseOutput(ofstream & output,bool header, bool create) const { if(header) output << "update decayers set parameters=\""; if(create) output << "create Herwig::SingletonFormFactor " << name() << " \n"; output << "newdef " << name() << ":CharmMass " << _mcharm/GeV << " \n"; output << "newdef " << name() << ":StrangeMass " << _mstrange/GeV << " \n"; output << "newdef " << name() << ":ThetaLambda " << _thetalambda << " \n"; output << "newdef " << name() << ":ThetaSigma " << _thetasigma << " \n"; output << "newdef " << name() << ":ThetaXi " << _thetaxi << " \n"; output << "newdef " << name() << ":ThetaXiPrime " << _thetaxip << " \n"; for(unsigned int ix=0;ix<_polemass.size();++ix) { if(ix _polemass; /** * The \f$\xi\f$ parameter for the form factors. */ vector _xi; /** * The normalisation factor, \f$N_{mM}\f$, for the form factors. */ vector _nmM; /** * The mass of the quark for the form factor. */ vector _mquark; }; } #endif /* HERWIG_SingletonFormFactor_H */ diff --git a/Decay/WeakCurrents/WeakBaryonCurrent.cc b/Decay/WeakCurrents/WeakBaryonCurrent.cc --- a/Decay/WeakCurrents/WeakBaryonCurrent.cc +++ b/Decay/WeakCurrents/WeakBaryonCurrent.cc @@ -1,209 +1,202 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the WeakBaryonCurrent class. // #include "WeakBaryonCurrent.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/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "ThePEG/Helicity/HelicityFunctions.h" using namespace Herwig; WeakBaryonCurrent::WeakBaryonCurrent() {} IBPtr WeakBaryonCurrent::clone() const { return new_ptr(*this); } IBPtr WeakBaryonCurrent::fullclone() const { return new_ptr(*this); } void WeakBaryonCurrent::persistentOutput(PersistentOStream & os) const { os << formFactor_; } void WeakBaryonCurrent::persistentInput(PersistentIStream & is, int) { is >> formFactor_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigWeakBaryonCurrent("Herwig::WeakBaryonCurrent", "Herwig.so"); void WeakBaryonCurrent::Init() { static ClassDocumentation documentation ("The WeakBaryonCurrent class is a wrapper for the BaryonFormFactor" " so it can be used as a WeakCurrent"); static Reference interfaceFormFactor ("FormFactor", "The baryon form factor", &WeakBaryonCurrent::formFactor_, false, false, true, false, false); } void WeakBaryonCurrent::doinit() { // initialize the form factor formFactor_->init(); // set the modes for(unsigned int iloc=0;ilocnumberOfFactors();++iloc) { int ispin(0), ospin(0), spect1(0), spect2(0), inquark(0), outquark(0); formFactor_->formFactorInfo(iloc,ispin,ospin,spect1,spect2,inquark,outquark); addDecayMode(outquark,-inquark); } setInitialModes(formFactor_->numberOfFactors()); WeakCurrent::doinit(); } // complete the construction of the decay mode for integration bool WeakBaryonCurrent::createMode(int icharge, tcPDPtr , FlavourInfo flavour, unsigned int imode,PhaseSpaceModePtr mode, unsigned int iloc,int ires, PhaseSpaceChannel phase, Energy upp ) { - // todo isospin in the form factors - assert(false); - // // no isospin here - // if(Itotal!=IsoSpin::IUnknown || i3 !=IsoSpin::I3Unknown) return false; unsigned int iq(0),ia(0); tPDVector out = particles(icharge,imode,iq,ia); // make sure the the decays are kinematically allowed Energy min =out[0]->massMin()+out[1]->massMin(); if(min>=upp) return false; // set the resonances and check charge tPDPtr res; if(icharge==3) res=getParticleData(ParticleID::Wplus ); else if(icharge==-3) res=getParticleData(ParticleID::Wminus); else res=getParticleData(ParticleID::gamma ); // create the channel mode->addChannel((PhaseSpaceChannel(phase),ires,res,ires+1,iloc+1,ires+1,iloc+2)); // return if successful return true; } // the particles produced by the current tPDVector WeakBaryonCurrent::particles(int icharge, unsigned int imode, int , int ) { tPDVector extpart(2); int id0(0),id1(0); formFactor_->particleID(imode,id0,id1); extpart[0] = getParticleData(id0); if(extpart[0]->CC()) extpart[0]=extpart[0]->CC(); extpart[1] = getParticleData(id1); int charge = extpart[0]->iCharge()+extpart[1]->iCharge(); if(charge==icharge) return extpart; else if(charge==-icharge) { for(unsigned int ix=0;ix<2;++ix) if(extpart[ix]->CC()) extpart[ix]=extpart[ix]->CC(); return extpart; } else return tPDVector(); } void WeakBaryonCurrent::constructSpinInfo(ParticleVector decay) const { if(decay[0]->id()>0) { SpinorWaveFunction ::constructSpinInfo(wave_ ,decay[1],outgoing,true); SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true); } else { SpinorWaveFunction ::constructSpinInfo( wave_,decay[0],outgoing,true); SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[1],outgoing,true); } } // hadronic current vector WeakBaryonCurrent::current(tcPDPtr , - FlavourInfo flavour, - const int, const int, Energy & scale, - const tPDVector & outgoing, - const vector & momenta, - DecayIntegrator::MEOption) const { - // no isospin here - assert(false); - // if(Itotal!=IsoSpin::IUnknown || i3 !=IsoSpin::I3Unknown) return vector(); + FlavourInfo flavour, + const int, const int, Energy & scale, + const tPDVector & outgoing, + const vector & momenta, + DecayIntegrator::MEOption) const { useMe(); Lorentz5Momentum q = momenta[0]+momenta[1]; q.rescaleMass(); scale=q.mass(); int in = abs(outgoing[0]->id()); int out = abs(outgoing[1]->id()); Energy m1 = outgoing[0]->mass(); Energy m2 = outgoing[1]->mass(); bool cc = false; unsigned int imode = formFactor_->formFactorNumber(in,out,cc); // todo generalize to spin != 1/2 assert(outgoing[0]->iSpin()==PDT::Spin1Half && outgoing[1]->iSpin()==PDT::Spin1Half ); wave_.resize(2); wavebar_.resize(2); for(unsigned int ix=0;ix<2;++ix) { wavebar_[ix] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ix,Helicity::outgoing); wave_[ix] = HelicityFunctions::dimensionedSpinor (-momenta[1],ix,Helicity::outgoing); } // get the form factors Complex f1v(0.),f2v(0.),f3v(0.),f1a(0.),f2a(0.),f3a(0.); formFactor_->SpinHalfSpinHalfFormFactor(sqr(scale),imode,in,out,m1,m2, - f1v,f2v,f3v,f1a,f2a,f3a, + f1v,f2v,f3v,f1a,f2a,f3a,flavour, BaryonFormFactor::TimeLike); Complex left = f1v - f1a + f2v -double((m1-m2)/(m1+m2))*f2a; Complex right = f1v + f1a + f2v +double((m1-m2)/(m1+m2))*f2a; vector baryon; Lorentz5Momentum diff = momenta[0]-momenta[1]; for(unsigned int ohel1=0;ohel1<2;++ohel1) { for(unsigned int ohel2=0;ohel2<2;++ohel2) { LorentzPolarizationVectorE vtemp = wave_[ohel2].generalCurrent(wavebar_[ohel1],left,right); complex vspin=wave_[ohel2].scalar (wavebar_[ohel1]); complex aspin=wave_[ohel2].pseudoScalar(wavebar_[ohel1]); vtemp-= (f2v*vspin+f2a*aspin)/(m1+m2)*diff; vtemp+= (f3v*vspin+f3a*aspin)/(m1+m2)*q; baryon.push_back(vtemp); } } // return the answer return baryon; } bool WeakBaryonCurrent::accept(vector id) { assert(id.size()==2); int itemp[2] = {id[0],id[1]}; for(unsigned int ix=0;ix<2;++ix) if(itemp[ix]<0) itemp[ix]=-itemp[ix]; bool cc = false; return formFactor_->formFactorNumber(itemp[0],itemp[1],cc)>=0; } // the decay mode unsigned int WeakBaryonCurrent::decayMode(vector idout) { assert(idout.size()==2); int itemp[2] = {idout[0],idout[1]}; for(unsigned int ix=0;ix<2;++ix) if(itemp[ix]<0) itemp[ix]=-itemp[ix]; bool cc = false; return formFactor_->formFactorNumber(itemp[0],itemp[1],cc); } // output the information for the database void WeakBaryonCurrent::dataBaseOutput(ofstream & output,bool header, bool create) const { if(header) output << "update decayers set parameters=\""; if(create) output << "create Herwig::WeakBaryonCurrent " << name() << "\n"; output << "newdef " << name() << ":FormFactor " << formFactor_ << "\n"; WeakCurrent::dataBaseOutput(output,false,false); if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; }