diff --git a/Decay/WeakCurrents/TwoKaonCzyzCurrent.cc b/Decay/WeakCurrents/TwoKaonCzyzCurrent.cc --- a/Decay/WeakCurrents/TwoKaonCzyzCurrent.cc +++ b/Decay/WeakCurrents/TwoKaonCzyzCurrent.cc @@ -1,764 +1,765 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the TwoKaonCzyzCurrent class. // #include "TwoKaonCzyzCurrent.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 "Herwig/Decay/ResonanceHelpers.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include using namespace Herwig; HERWIG_INTERPOLATOR_CLASSDESC(TwoKaonCzyzCurrent,double,Energy2) TwoKaonCzyzCurrent::TwoKaonCzyzCurrent() // changed parameters from 1002.0279, Fit 2 // substituted by own fit : betaRho_(2.19680665014), betaOmega_(2.69362046884), betaPhi_(1.94518176513), - nMax_(200), etaPhi_(1.055), gammaOmega_(0.5), gammaPhi_(0.2), mpi_(140.*MeV), eMax_(10.*GeV) { + nMax_(200), etaPhi_(1.055), gammaOmega_(0.5), gammaPhi_(0.2), mpi_(140.*MeV), eMax_(-GeV) { using Constants::pi; // rho parameter rhoMag_ = {1.1148916618504967,0.050374779737077324, 0.014908906283692132,0.03902475997619905,0.038341465215871416}; rhoPhase_ = {0 , pi, pi, pi, pi}; rhoMasses_ = {775.49*MeV,1520.6995754050117*MeV,1740.9719246639341*MeV,1992.2811314327789*MeV}; rhoWidths_ = {149.4 *MeV,213.41728317817743*MeV, 84.12224414791908*MeV,289.9733272437917*MeV}; // omega parameters omegaMag_ = {1.3653229680598022, 0.02775156567495144, 0.32497165559032715,1.3993153161869765}; omegaPhase_ = {0 , pi, pi, 0,pi}; omegaMasses_ = {782.65*MeV,1414.4344268685891*MeV,1655.375231284883*MeV}; omegaWidths_ = {8.490000000000001*MeV, 85.4413887755723*MeV, 160.31760444832305*MeV}; // phi parameters phiMag_ = {0.965842498579515,0.002379766320723148,0.1956211640216197,0.16527771485190898}; phiPhase_ = {0. ,pi ,pi ,0. ,0.}; phiMasses_ = {1019.4209171596993*MeV,1594.759278457624*MeV,2156.971341201067*MeV}; phiWidths_ = {4.252653332329334*MeV, 28.741821847408196*MeV,673.7556174184005*MeV}; // set up for the modes in the base class addDecayMode(2,-1); addDecayMode(1,-1); addDecayMode(2,-2); addDecayMode(1,-1); addDecayMode(2,-2); setInitialModes(5); } IBPtr TwoKaonCzyzCurrent::clone() const { return new_ptr(*this); } IBPtr TwoKaonCzyzCurrent::fullclone() const { return new_ptr(*this); } void TwoKaonCzyzCurrent::persistentOutput(PersistentOStream & os) const { os << betaRho_ << betaOmega_ << betaPhi_ << rhoWgt_ << rhoMag_ << rhoPhase_ << ounit(rhoMasses_,GeV) << ounit(rhoWidths_,GeV) << phiWgt_ << phiMag_ << phiPhase_ << ounit(phiMasses_,GeV) << ounit(phiWidths_,GeV) << ounit(mass_,GeV) << ounit(width_,GeV) << coup_ << dh_ << ounit(hres_,GeV2) << ounit(h0_,GeV2) << nMax_ << etaPhi_ << gammaOmega_ << gammaPhi_ << ounit(mpi_,GeV) << ounit(eMax_,GeV) << fKI0Re_ << fKI0Im_ << fKI1Re_ << fKI1Im_; } void TwoKaonCzyzCurrent::persistentInput(PersistentIStream & is, int) { is >> betaRho_ >> betaOmega_ >> betaPhi_ >> rhoWgt_ >> rhoMag_ >> rhoPhase_ >> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV) >> phiWgt_ >> phiMag_ >> phiPhase_ >> iunit(phiMasses_,GeV) >> iunit(phiWidths_,GeV) >> iunit(mass_,GeV) >> iunit(width_,GeV) >> coup_ >> dh_ >> iunit(hres_,GeV2) >> iunit(h0_,GeV2) >> nMax_ >> etaPhi_ >> gammaOmega_ >> gammaPhi_ >> iunit(mpi_,GeV) >> iunit(eMax_,GeV) >> fKI0Re_ >> fKI0Im_ >> fKI1Re_ >> fKI1Im_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigTwoKaonCzyzCurrent("Herwig::TwoKaonCzyzCurrent", "HwWeakCurrents.so"); void TwoKaonCzyzCurrent::Init() { static ClassDocumentation documentation ("The TwoKaonCzyzCurrent class uses the currents from " "PRD 81 094014 for the weak current with two kaons", "The current for two kaons from \\cite{Czyz:2010hj} was used.", "%\\cite{Czyz:2010hj}\n" "\\bibitem{Czyz:2010hj}\n" "H.~Czyz, A.~Grzelinska and J.~H.~Kuhn,\n" "%``Narrow resonances studies with the radiative return method,''\n" "Phys.\\ Rev.\\ D {\\bf 81} (2010) 094014\n" "doi:10.1103/PhysRevD.81.094014\n" "[arXiv:1002.0279 [hep-ph]].\n" "%%CITATION = doi:10.1103/PhysRevD.81.094014;%%\n" "%28 citations counted in INSPIRE as of 30 Jul 2018\n"); static ParVector interfaceRhoMasses ("RhoMasses", "The masses of the different rho resonances for the pi pi channel", &TwoKaonCzyzCurrent::rhoMasses_, MeV, -1, 775.8*MeV, ZERO, 10000.*MeV, false, false, true); static ParVector interfaceRhoWidths ("RhoWidths", "The widths of the different rho resonances for the pi pi channel", &TwoKaonCzyzCurrent::rhoWidths_, MeV, -1, 150.3*MeV, ZERO, 1000.*MeV, false, false, true); static ParVector interfaceRhoMagnitude ("RhoMagnitude", "Magnitude of the weight of the different resonances for the pi pi channel", &TwoKaonCzyzCurrent::rhoMag_, -1, 0., 0, 0, false, false, Interface::nolimits); static ParVector interfaceRhoPhase ("RhoPhase", "Phase of the weight of the different resonances for the pi pi channel", &TwoKaonCzyzCurrent::rhoPhase_, -1, 0., 0, 0, false, false, Interface::nolimits); static ParVector interfaceOmegaMasses ("OmegaMasses", "The masses of the different omega resonances for the pi pi channel", &TwoKaonCzyzCurrent::omegaMasses_, MeV, -1, 775.8*MeV, ZERO, 10000.*MeV, false, false, true); static ParVector interfaceOmegaWidths ("OmegaWidths", "The widths of the different omega resonances for the pi pi channel", &TwoKaonCzyzCurrent::omegaWidths_, MeV, -1, 150.3*MeV, ZERO, 1000.*MeV, false, false, true); static ParVector interfaceOmegaMagnitude ("OmegaMagnitude", "Magnitude of the weight of the different resonances for the pi pi channel", &TwoKaonCzyzCurrent::omegaMag_, -1, 0., 0, 0, false, false, Interface::nolimits); static ParVector interfaceOmegaPhase ("OmegaPhase", "Phase of the weight of the different resonances for the pi pi channel", &TwoKaonCzyzCurrent::omegaPhase_, -1, 0., 0, 0, false, false, Interface::nolimits); static ParVector interfacePhiMasses ("PhiMasses", "The masses of the different phi resonances for the pi pi channel", &TwoKaonCzyzCurrent::phiMasses_, MeV, -1, 775.8*MeV, ZERO, 10000.*MeV, false, false, true); static ParVector interfacePhiWidths ("PhiWidths", "The widths of the different phi resonances for the pi pi channel", &TwoKaonCzyzCurrent::phiWidths_, MeV, -1, 150.3*MeV, ZERO, 1000.*MeV, false, false, true); static ParVector interfacePhiMagnitude ("PhiMagnitude", "Magnitude of the weight of the different resonances for the pi pi channel", &TwoKaonCzyzCurrent::phiMag_, -1, 0., 0, 0, false, false, Interface::nolimits); static ParVector interfacePhiPhase ("PhiPhase", "Phase of the weight of the different resonances for the pi pi channel", &TwoKaonCzyzCurrent::phiPhase_, -1, 0., 0, 0, false, false, Interface::nolimits); static Parameter interfacenMax ("nMax", "The maximum number of resonances to include in the sum," " should be approx infinity", &TwoKaonCzyzCurrent::nMax_, 200, 10, 10000, false, false, Interface::limited); static Parameter interfacebetaRho ("betaRho", "The beta parameter for the rho couplings", &TwoKaonCzyzCurrent::betaRho_, 2.23, 0.0, 100., false, false, Interface::limited); static Parameter interfacebetaOmega ("betaOmega", "The beta parameter for the rho couplings", &TwoKaonCzyzCurrent::betaOmega_, 2.23, 0.0, 100., false, false, Interface::limited); static Parameter interfacebetaPhi ("betaPhi", "The beta parameter for the phi couplings", &TwoKaonCzyzCurrent::betaPhi_, 1.97, 0.0, 100., false, false, Interface::limited); static Parameter interfaceEtaPhi ("EtaPhi", "The eta_phi mixing parameter", &TwoKaonCzyzCurrent::etaPhi_, 1.04, 0.0, 10.0, false, false, Interface::limited); static Parameter interfacegammaOmega ("gammaOmega", "The gamma parameter for the widths of omega resonances", &TwoKaonCzyzCurrent::gammaOmega_, 0.5, 0.0, 1.0, false, false, Interface::limited); static Parameter interfacegammaPhi ("gammaPhi", "The gamma parameter for the widths of phi resonances", &TwoKaonCzyzCurrent::gammaPhi_, 0.2, 0.0, 1.0, false, false, Interface::limited); } void TwoKaonCzyzCurrent::doinit() { WeakCurrent::doinit(); // check consistency of parametrers if(rhoMasses_.size() != rhoWidths_.size() || omegaMasses_.size() != omegaWidths_.size() || phiMasses_.size() != phiWidths_.size() ) throw InitException() << "Inconsistent parameters in TwoKaonCzyzCurrent" << "::doinit()" << Exception::abortnow; // weights for the rho channels if(rhoMag_.size()!=rhoPhase_.size()) throw InitException() << "The vectors containing the weights and phase for the" << " rho channel must be the same size in " << "TwoKaonCzyzCurrent::doinit()" << Exception::runerror; // combine mags and phase for(unsigned int ix=0;ixmass(); // rho masses and couplings double gamB(std::tgamma(2.-betaRho_)); mass_.push_back(vector()); width_.push_back(vector()); coup_.push_back(vector()); Complex total(0.); for(unsigned int ix=0;ix0) { gamB *= ((1.-betaRho_+double(ix)))/double(ix); } Complex c_n = std::tgamma(betaRho_-0.5) /(0.5+double(ix)) / sqrt(Constants::pi) * sin(Constants::pi*(betaRho_-1.-double(ix)))/Constants::pi*gamB; if(ix%2!=0) c_n *= -1.; // couplings coup_[0].push_back(c_n); total+=c_n; // set the masses and widths // calc for higher resonances if(ix>=rhoMasses_.size()) { mass_ [0].push_back(rhoMasses_[0]*sqrt(1.+2.*double(ix))); width_[0].push_back(rhoWidths_[0]/rhoMasses_[0]*mass_[0].back()); } // input for lower ones else { mass_ [0].push_back(rhoMasses_[ix]); width_[0].push_back(rhoWidths_[ix]); } // parameters for the gs propagators hres_.push_back(Resonance::Hhat(sqr(mass_[0].back()), mass_[0].back(),width_[0].back(),mpi_,mpi_)); dh_ .push_back(Resonance::dHhatds(mass_[0].back(),width_[0].back(),mpi_,mpi_)); h0_ .push_back(Resonance::H(ZERO,mass_[0].back(),width_[0].back(), mpi_,mpi_,dh_.back(),hres_.back())); } for(unsigned int ix=0;ix()); width_.push_back(vector()); coup_.push_back(vector()); total=0.; for(unsigned int ix=0;ix0) { gamB *= ((1.-betaOmega_+double(ix)))/double(ix); } Complex c_n = std::tgamma(betaOmega_-0.5) /(0.5+double(ix)) / sqrt(Constants::pi) * sin(Constants::pi*(betaOmega_-1.-double(ix)))/Constants::pi*gamB; if(ix%2!=0) c_n *= -1.; // couplings coup_[1].push_back(c_n); total+=c_n; // set the masses and widths // calc for higher resonances if(ix>=omegaMasses_.size()) { mass_ [1].push_back(omegaMasses_[0]*sqrt(1.+2.*double(ix))); width_[1].push_back(gammaOmega_*mass_[1].back()); } // input for lower ones else { mass_ [1].push_back(omegaMasses_[ix]); width_[1].push_back(omegaWidths_[ix]); } } for(unsigned int ix=0;ix()); width_.push_back(vector()); coup_.push_back(vector()); total=0.; for(unsigned int ix=0;ix0) { gamB *= ((1.-betaPhi_+double(ix)))/double(ix); } Complex c_n = std::tgamma(betaPhi_-0.5) /(0.5+double(ix)) / sqrt(Constants::pi) * sin(Constants::pi*(betaPhi_-1.-double(ix)))/Constants::pi*gamB; if(ix%2!=0) c_n *= -1.; // couplings coup_[2].push_back(c_n); total +=c_n; // set the masses and widths // calc for higher resonances if(ix>=phiMasses_.size()) { mass_ [2].push_back(phiMasses_[0]*sqrt(1.+2.*double(ix))); width_[2].push_back(gammaPhi_*mass_[2].back()); } // input for lower ones else { mass_ [2].push_back(phiMasses_[ix]); width_[2].push_back(phiWidths_[ix]); } } for(unsigned int ix=0;ix en; vector re0,im0; vector re1,im1; Energy mK = getParticleData(ParticleID::Kplus)->mass(); - Energy2 step = (sqr(eMax_)-sqr(2.*mK))/nMax_; + Energy maxE = eMax_>ZERO ? eMax_ : 10.*GeV; + Energy2 step = (sqr(maxE)-sqr(2.*mK))/nMax_; Energy2 Q2 = sqr(2.*mK); for(unsigned int ix=0;ix0 && icharge !=0)) return false; // check the total isospin if(flavour.I!=IsoSpin::IUnknown) { if(flavour.I!=IsoSpin::IOne && flavour.I!=IsoSpin::IZero ) return false; } // check I_3 if(flavour.I3!=IsoSpin::I3Unknown&&flavour.I==IsoSpin::IOne) { switch(flavour.I3) { case IsoSpin::I3Zero: if(imode==0) return false; break; case IsoSpin::I3One: if(imode!=0 || icharge ==-3) return false; break; case IsoSpin::I3MinusOne: if(imode!=0 || icharge ==3) return false; break; default: return false; } } if(imode==0 && (flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero )) return false; if(imode!=0 && (flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero and flavour.strange != Strangeness::ssbar )) return false; if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return false; if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return false; // make sure that the decays are kinematically allowed tPDPtr part[2]; if(imode==0) { part[0]=getParticleData(ParticleID::Kplus); part[1]=getParticleData(ParticleID::Kbar0); } else if(imode==1|| imode==2) { part[0]=getParticleData(ParticleID::K_S0); part[1]=getParticleData(ParticleID::K_L0); } else { part[0]=getParticleData(ParticleID::Kplus); part[1]=getParticleData(ParticleID::Kminus); } Energy min(part[0]->massMin()+part[1]->massMin()); if(min>upp) return false; - eMax_=upp; + eMax_=max(upp,eMax_); // set the resonances vector res; if(icharge==0) { res.push_back(getParticleData(113 )); res.push_back(getParticleData(100113)); res.push_back(getParticleData(30113 )); res.push_back(getParticleData( 223)); res.push_back(getParticleData( 333)); } else { res.push_back(getParticleData(213 )); res.push_back(getParticleData(100213)); res.push_back(getParticleData(30213 )); if(icharge==-3) { for(unsigned int ix=0;ix<3;++ix) { if(res[ix]&&res[ix]->CC()) res[ix]=res[ix]->CC(); } } } // create the channels for(unsigned int ix=0;ix=3) continue; PhaseSpaceChannel newChannel((PhaseSpaceChannel(phase),ires,res[ix],ires+1,iloc+1,ires+1,iloc+2)); mode->addChannel(newChannel); } // reset the masses in the intergrators for(unsigned int ix=0;ix<3;++ix) { if(ixresetIntermediate(res[ix],rhoMasses_[ix],rhoWidths_[ix]); } } if(res.size()>3) { mode->resetIntermediate(res[3],omegaMasses_[0],omegaWidths_[0]); mode->resetIntermediate(res[4],phiMasses_ [0], phiWidths_[0]); } // return if successful return true; } // the particles produced by the current tPDVector TwoKaonCzyzCurrent::particles(int icharge, unsigned int imode, int,int) { tPDVector output(2); if(imode==0) { output[0]=getParticleData(ParticleID::Kplus); output[1]=getParticleData(ParticleID::K0); } else if(imode==1||imode==2) { output[0]=getParticleData(ParticleID::K_S0); output[1]=getParticleData(ParticleID::K_L0); } else { output[0]=getParticleData(ParticleID::Kplus ); output[1]=getParticleData(ParticleID::Kminus); } if(icharge==-3) { for(unsigned int ix=0;ixCC()) output[ix]=output[ix]->CC(); } } return output; } // hadronic current vector TwoKaonCzyzCurrent::current(tcPDPtr resonance, FlavourInfo flavour, const int imode, const int ichan,Energy & scale, const tPDVector & outgoing, const vector & momenta, DecayIntegrator::MEOption) const { useMe(); // check the total isospin if(flavour.I!=IsoSpin::IUnknown) { if(flavour.I!=IsoSpin::IOne && flavour.I!=IsoSpin::IZero ) return vector(); } // check I_3 int icharge = outgoing[0]->iCharge()+outgoing[1]->iCharge(); if(flavour.I3!=IsoSpin::I3Unknown&&flavour.I==IsoSpin::IOne) { switch(flavour.I3) { case IsoSpin::I3Zero: if(imode==0) return vector(); break; case IsoSpin::I3One: if(imode!=0 || icharge ==-3) return vector(); break; case IsoSpin::I3MinusOne: if(imode!=0 || icharge ==3) return vector(); break; default: return vector(); } } // momentum difference and sum of the mesons Lorentz5Momentum pdiff(momenta[0]-momenta[1]); Lorentz5Momentum psum (momenta[0]+momenta[1]); psum.rescaleMass(); scale=psum.mass(); // mass2 of vector intermediate state Energy2 q2(psum.m2()); double dot(psum*pdiff/q2); psum *=dot; // calculate the current Complex FK = Fkaon(q2,imode,ichan, flavour.I,flavour.strange,resonance, momenta[0].mass(),momenta[1].mass()); // compute the current pdiff -= psum; return vector(1,FK*pdiff); } bool TwoKaonCzyzCurrent::accept(vector id) { // check there are only two particles if(id.size()!=2) return false; // pion modes if((id[0]==ParticleID::Kminus && id[1]==ParticleID::K0) || (id[0]==ParticleID::K0 && id[1]==ParticleID::Kminus) || (id[0]==ParticleID::Kplus && id[1]==ParticleID::Kbar0) || (id[0]==ParticleID::Kbar0 && id[1]==ParticleID::Kplus)) return true; else if((id[0]==ParticleID::Kminus && id[1]==ParticleID::Kplus) || (id[0]==ParticleID::Kplus && id[1]==ParticleID::Kminus)) return true; else if((id[0]==ParticleID::K_S0 && id[1]==ParticleID::K_L0) || (id[0]==ParticleID::K_L0 && id[1]==ParticleID::K_S0)) return true; else return false; } // the decay mode unsigned int TwoKaonCzyzCurrent::decayMode(vector idout) { unsigned int nk0(0),nkp(0); for(unsigned int ix=0;ix=0) { if(ichan<3) { on[1]=on[2]=false; imin = ichan; imax = ichan+1; } else if(ichan==3) { on[0]=on[2]=false; imin=0; imax=1; } else if(ichan==4) { on[0]=on[1]=false; imin=0; imax=1; } else assert(false); } if(resonance) { switch(resonance->id()%1000) { case 223: imin=0; on[0]=on[2]=false; break; case 333: imin=0; on[0]=on[1]=false; break; case 113: switch(resonance->id()/1000) { case 0: imin=0; break; case 100: imin = 1; break; case 30 : imin = 2; break; default : assert(false); } on[1]=on[2]=false; break; default: assert(false); } imax = imin+1; } // calculate the form factor Complex FK(0.); for(unsigned int ix=imin;ix> beta_ >> omegaWgt_ >> omegaMag_ >> omegaPhase_ >> iunit(omegaMass_,GeV) >> iunit(omegaWidth_,GeV) >> rhoWgt_ >> rhoMag_ >> rhoPhase_ >> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV) >> iunit(mass_,GeV) >> iunit(width_,GeV) >> coup_ >> dh_ >> iunit(hres_,GeV2) >> iunit(h0_,GeV2) >> nMax_ >> iunit(eMax_,GeV) >> fpiRe_ >> fpiIm_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigTwoPionCzyzCurrent("Herwig::TwoPionCzyzCurrent", "HwWeakCurrents.so"); void TwoPionCzyzCurrent::Init() { static ClassDocumentation documentation ("The TwoPionCzyzCurrent class uses the currents from " "PRD 81 094014 for the weak current with two pions", "The current for two pions from \\cite{Czyz:2010hj} was used.", "%\\cite{Czyz:2010hj}\n" "\\bibitem{Czyz:2010hj}\n" "H.~Czyz, A.~Grzelinska and J.~H.~Kuhn,\n" "%``Narrow resonances studies with the radiative return method,''\n" "Phys.\\ Rev.\\ D {\\bf 81} (2010) 094014\n" "doi:10.1103/PhysRevD.81.094014\n" "[arXiv:1002.0279 [hep-ph]].\n" "%%CITATION = doi:10.1103/PhysRevD.81.094014;%%\n" "%28 citations counted in INSPIRE as of 30 Jul 2018\n"); static ParVector interfaceRhoMasses ("RhoMasses", "The masses of the different rho resonances for the pi pi channel", &TwoPionCzyzCurrent::rhoMasses_, MeV, -1, 775.8*MeV, ZERO, 10000.*MeV, false, false, true); static ParVector interfaceRhoWidths ("RhoWidths", "The widths of the different rho resonances for the pi pi channel", &TwoPionCzyzCurrent::rhoWidths_, MeV, -1, 150.3*MeV, ZERO, 1000.*MeV, false, false, true); static ParVector interfaceRhoMagnitude ("RhoMagnitude", "Magnitude of the weight of the different resonances for the pi pi channel", &TwoPionCzyzCurrent::rhoMag_, -1, 0., 0, 0, false, false, Interface::nolimits); static ParVector interfaceRhoPhase ("RhoPhase", "Phase of the weight of the different resonances for the pi pi channel", &TwoPionCzyzCurrent::rhoPhase_, -1, 0., 0, 0, false, false, Interface::nolimits); static Parameter interfacenMax ("nMax", "The maximum number of resonances to include in the sum," " should be approx infinity", &TwoPionCzyzCurrent::nMax_, 1000, 10, 10000, false, false, Interface::limited); static Parameter interfacebeta ("beta", "The beta parameter for the couplings", &TwoPionCzyzCurrent::beta_, 2.148, 0.0, 100., false, false, Interface::limited); static Parameter interfaceOmegaMass ("OmegaMass", "The mass of the omega meson", &TwoPionCzyzCurrent::omegaMass_, MeV,782.4*MeV, 0.0*MeV, 1000.0*MeV, false, false, Interface::limited); static Parameter interfaceOmegaWidth ("OmegaWidth", "The mass of the omega meson", &TwoPionCzyzCurrent::omegaWidth_, MeV, 8.33*MeV, 0.0*MeV, 500.0*MeV, false, false, Interface::limited); static Parameter interfaceOmegaMagnitude ("OmegaMagnitude", "The magnitude of the omega couplings", &TwoPionCzyzCurrent::omegaMag_, 18.7e-4, 0.0, 10.0, false, false, Interface::limited); static Parameter interfaceOmegaPhase ("OmegaPhase", "The magnitude of the omega couplings", &TwoPionCzyzCurrent::omegaPhase_, 0.106, 0.0, 2.*Constants::pi, false, false, Interface::limited); } void TwoPionCzyzCurrent::doinit() { WeakCurrent::doinit(); // check consistency of parametrers if(rhoMasses_.size()!=rhoWidths_.size()) throw InitException() << "Inconsistent parameters in TwoPionCzyzCurrent" << "::doinit()" << Exception::abortnow; // weights for the rho channels if(rhoMag_.size()!=rhoPhase_.size()) throw InitException() << "The vectors containing the weights and phase for the" << " rho channel must be the same size in " << "TwoPionCzyzCurrent::doinit()" << Exception::runerror; Complex rhoSum(0.); for(unsigned int ix=0;ix0) rhoSum +=rhoWgt_.back(); } omegaWgt_ = omegaMag_*(cos(omegaPhase_)+Complex(0.,1.)*sin(omegaPhase_)); // set up the masses and widths of the rho resonances double gamB(tgamma(2.-beta_)); Complex cwgt(0.); Energy mpi(getParticleData(ParticleID::piplus)->mass()); for(unsigned int ix=0;ix0) { gamB *= ((1.-beta_+double(ix)))/double(ix); } Complex c_n = tgamma(beta_-0.5) /(0.5+double(ix)) / sqrt(Constants::pi) * sin(Constants::pi*(beta_-1.-double(ix)))/Constants::pi*gamB; if(ix%2!=0) c_n *= -1.; // set the masses and widths // calc for higher resonances if(ix>=rhoMasses_.size()) { mass_ .push_back(rhoMasses_[0]*sqrt(1.+2.*double(ix))); width_.push_back(rhoWidths_[0]/rhoMasses_[0]*mass_.back()); } // input for lower ones else { mass_ .push_back(rhoMasses_[ix]); width_.push_back(rhoWidths_[ix]); if(ix>0) cwgt += c_n; } // parameters for the gs propagators hres_.push_back(Resonance::Hhat(sqr(mass_.back()),mass_.back(),width_.back(),mpi,mpi)); dh_ .push_back(Resonance::dHhatds(mass_.back(),width_.back(),mpi,mpi)); h0_.push_back(Resonance::H(ZERO,mass_.back(),width_.back(),mpi,mpi,dh_.back(),hres_.back())); coup_.push_back(c_n); } // fix up the early weights for(unsigned int ix=1;ixmass()); vector en; vector re,im; - Energy step = (eMax_-2.*mpi)/nMax_; + Energy maxE = eMax_>ZERO ? eMax_ : 10.*GeV; + Energy step = (maxE-2.*mpi)/nMax_; Energy Q = 2.*mpi; for(unsigned int ix=0;ix0 && icharge !=0)) return false; // check the total isospin if(flavour.I!=IsoSpin::IUnknown) { if(flavour.I!=IsoSpin::IOne) return false; } // check I_3 if(flavour.I3!=IsoSpin::I3Unknown) { switch(flavour.I3) { case IsoSpin::I3Zero: if(imode==0) return false; break; case IsoSpin::I3One: if(imode==1 || icharge ==-3) return false; break; case IsoSpin::I3MinusOne: if(imode==1 || icharge ==3) return false; break; default: return false; } } if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero ) return false; if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return false; if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return false; // make sure that the decays are kinematically allowed tPDPtr part[2]; if(imode==0) { part[0]=getParticleData(ParticleID::piplus); part[1]=getParticleData(ParticleID::pi0); } else { part[0]=getParticleData(ParticleID::piplus); part[1]=getParticleData(ParticleID::piminus); } Energy min(part[0]->massMin()+part[1]->massMin()); if(min>upp) return false; - eMax_=upp; + eMax_=max(upp,eMax_); // set up the resonances tPDPtr res[3]; if(icharge==0) { res[0] =getParticleData(113); res[1] =getParticleData(100113); res[2] =getParticleData(30113); } else { res[0] =getParticleData(213); res[1] =getParticleData(100213); res[2] =getParticleData(30213); if(icharge==-3) { for(unsigned int ix=0;ix<3;++ix) { if(res[ix]&&res[ix]->CC()) res[ix]=res[ix]->CC(); } } } // create the channels for(unsigned int ix=0;ix<3;++ix) { if(!res[ix]) continue; if(resonance && resonance != res[ix]) continue; mode->addChannel((PhaseSpaceChannel(phase),ires,res[ix],ires+1,iloc+1,ires+1,iloc+2)); } // reset the masses in the intergrators for(unsigned int ix=0;ix<3;++ix) { if(ixresetIntermediate(res[ix],rhoMasses_[ix],rhoWidths_[ix]); } } return true; } // the particles produced by the current tPDVector TwoPionCzyzCurrent::particles(int icharge, unsigned int imode, int,int) { tPDVector output(2); if(imode==0) { output[0]=getParticleData(ParticleID::piplus); output[1]=getParticleData(ParticleID::pi0); if(icharge==-3) { for(unsigned int ix=0;ixCC()) output[ix]=output[ix]->CC(); } } } else { output[0]=getParticleData(ParticleID::piplus); output[1]=getParticleData(ParticleID::piminus); } return output; } // hadronic current vector TwoPionCzyzCurrent::current(tcPDPtr resonance, FlavourInfo flavour, const int imode, const int ichan,Energy & scale, const tPDVector & outgoing, const vector & momenta, DecayIntegrator::MEOption) const { useMe(); // check the isospin if(flavour.I!=IsoSpin::IUnknown && flavour.I!=IsoSpin::IOne) return vector(); int icharge = outgoing[0]->iCharge()+outgoing[1]->iCharge(); // check I_3 if(flavour.I3!=IsoSpin::I3Unknown) { switch(flavour.I3) { case IsoSpin::I3Zero: if(imode==0) return vector(); break; case IsoSpin::I3One: if(imode==1 || icharge ==-3) return vector(); break; case IsoSpin::I3MinusOne: if(imode==1 || icharge ==3) return vector(); break; default: return vector(); } } if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return vector(); if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return vector(); if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return vector(); // momentum difference and sum of the mesons Lorentz5Momentum pdiff(momenta[0]-momenta[1]); Lorentz5Momentum psum (momenta[0]+momenta[1]); psum.rescaleMass(); scale=psum.mass(); // mass2 of vector intermediate state Energy2 q2(psum.m2()); double dot(psum*pdiff/q2); psum *=dot; // compute the form factor Complex FPI=Fpi(q2,imode,ichan,resonance,momenta[0].mass(),momenta[1].mass()); // calculate the current pdiff -= psum; return vector(1,FPI*pdiff); } bool TwoPionCzyzCurrent::accept(vector id) { // check there are only two particles if(id.size()!=2) return false; // pion modes if((abs(id[0])==ParticleID::piplus && id[1] ==ParticleID::pi0 ) || ( id[0] ==ParticleID::pi0 && abs(id[1])==ParticleID::piplus)) return true; else if((id[0]==ParticleID::piminus && id[1]==ParticleID::piplus) || (id[0]==ParticleID::piplus && id[1]==ParticleID::piminus)) return true; else return false; } // the decay mode unsigned int TwoPionCzyzCurrent::decayMode(vector idout) { unsigned int npi(0); for(unsigned int ix=0;ix0) { imin = ichan; imax = ichan+1; } if(resonance) { switch(resonance->id()/1000) { case 0: imax = 1; break; case 100: imin = 1; imax = 2; break; case 30 : imin = 2; imax = 3; break; default: assert(false); } } for(unsigned int ix=imin;ix