diff --git a/Decay/ScalarMeson/CLEOD0toK0PipPim.cc b/Decay/ScalarMeson/CLEOD0toK0PipPim.cc --- a/Decay/ScalarMeson/CLEOD0toK0PipPim.cc +++ b/Decay/ScalarMeson/CLEOD0toK0PipPim.cc @@ -1,131 +1,447 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the CLEOD0toK0PipPim class. // #include "CLEOD0toK0PipPim.h" #include "ThePEG/Interface/ClassDocumentation.h" +#include "ThePEG/Interface/Parameter.h" +#include "ThePEG/Interface/Switch.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" using namespace Herwig; -CLEOD0toK0PipPim::CLEOD0toK0PipPim() : WeakDalitzDecay(5./GeV,1.5/GeV,true) -{} +CLEOD0toK0PipPim::CLEOD0toK0PipPim() : WeakDalitzDecay(5./GeV,1.5/GeV,true) { + // amplitudes and phases for D0 -> K0pi+pi- + aKstarp_ = 0.11 ; phiKstarp_ = 321; + arho_ = 1.00 ; phirho_ = 0; + aomega_ = 0.037 ; phiomega_ = 114; + aKstarm_ = 1.56 ; phiKstarm_ = 150; + af980_ = 0.34*GeV2; phif980_ = 188; + af2_ = 0.7/GeV2 ; phif2_ = 308; + af1370_ = 1.8*GeV2 ; phif1370_ = 85; + aK14300_ = 2.0*GeV2 ; phiK14300_ = 3; + aK14302_ = 1.0/GeV2 ; phiK14302_ = 335; + aK1680_ = 5.6 ; phiK1680_ = 174; + aNR_ = 1.1 ; phiNR_ = 340; + // masses and widths + f0opt_ = true; + gpi_ = 0.09; + gK_ = 0.02; + mK892_ = 891.66*MeV; wK892_ = 50.8 *MeV; + mrho_ = 769.3 *MeV; wrho_ = 150.2 *MeV; + momega_ = 782.57*MeV; womega_ = 8.44*MeV; + mf980_ = 977.00*MeV; wf980_ = 50. *MeV; + mf2_ = 1275.4 *MeV; wf2_ = 185.1 *MeV; + mf1370_ = 1310 *MeV; wf1370_ = 272.0 *MeV; + mK14300_ = 1412 *MeV; wK14300_ = 294 *MeV; + mK14302_ = 1425.6 *MeV; wK14302_ = 98.5 *MeV; + mK1680_ = 1717 *MeV; wK1680_ = 322 *MeV; + // intermediates + generateIntermediates(true); +} IBPtr CLEOD0toK0PipPim::clone() const { return new_ptr(*this); } IBPtr CLEOD0toK0PipPim::fullclone() const { return new_ptr(*this); } -void CLEOD0toK0PipPim::persistentOutput(PersistentOStream & ) const { +void CLEOD0toK0PipPim::persistentOutput(PersistentOStream & os) const { + os << ounit(mrho_,GeV) << ounit(wrho_,GeV) << ounit(momega_,GeV) << ounit(womega_,GeV) + << ounit(mf980_,GeV) << ounit(wf980_,GeV) << gpi_ << gK_ << f0opt_ + << ounit(mf2_,GeV) << ounit(wf2_,GeV) << ounit(mf1370_,GeV) << ounit(wf1370_,GeV) + << ounit(mK892_,GeV) << ounit(wK892_,GeV) << ounit(mK14300_,GeV) << ounit(wK14300_,GeV) + << ounit(mK14302_,GeV) << ounit(wK14302_,GeV) << ounit(mK1680_,GeV) << ounit(wK1680_,GeV) + << aKstarp_ << phiKstarp_ << arho_ << phirho_ << aomega_ << phiomega_ + << aKstarm_ << phiKstarm_ << ounit(af980_,GeV2) << phif980_ << ounit(af2_,1./GeV2) << phif2_ + << ounit(af1370_,GeV2) << phif1370_ << ounit(aK14300_,GeV2) << phiK14300_ + << ounit(aK14302_,1./GeV2) << phiK14302_ + << aK1680_ << phiK1680_ << aNR_ << phiNR_ << cNR_; } -void CLEOD0toK0PipPim::persistentInput(PersistentIStream & , int) { +void CLEOD0toK0PipPim::persistentInput(PersistentIStream & is, int) { + is >> iunit(mrho_,GeV) >> iunit(wrho_,GeV) >> iunit(momega_,GeV) >> iunit(womega_,GeV) + >> iunit(mf980_,GeV) >> iunit(wf980_,GeV) >> gpi_ >> gK_ >> f0opt_ + >> iunit(mf2_,GeV) >> iunit(wf2_,GeV) >> iunit(mf1370_,GeV) >> iunit(wf1370_,GeV) + >> iunit(mK892_,GeV) >> iunit(wK892_,GeV) >> iunit(mK14300_,GeV) >> iunit(wK14300_,GeV) + >> iunit(mK14302_,GeV) >> iunit(wK14302_,GeV) >> iunit(mK1680_,GeV) >> iunit(wK1680_,GeV) + >> aKstarp_ >> phiKstarp_ >> arho_ >> phirho_ >> aomega_ >> phiomega_ + >> aKstarm_ >> phiKstarm_ >> iunit(af980_,GeV2) >> phif980_ >> iunit(af2_,1./GeV2) >> phif2_ + >> iunit(af1370_,GeV2) >> phif1370_ >> iunit(aK14300_,GeV2) >> phiK14300_ + >> iunit(aK14302_,1./GeV2) >> phiK14302_ + >> aK1680_ >> phiK1680_ >> aNR_ >> phiNR_ >> cNR_; } - // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigCLEOD0toK0PipPim("Herwig::CLEOD0toK0PipPim", "CLEOD0toK0PipPim.so"); void CLEOD0toK0PipPim::Init() { static ClassDocumentation documentation ("The CLEOD0toK0PipPim class implements the model of CLEO for" " D0 -> Kbar0 pi+pi-", "The CLEO fit of \\cite{Muramatsu:2002jp} was used for the decay $D^0\\to\\bar{K}^0\\pi^+\\pi^-$", "\\bibitem{Muramatsu:2002jp} H.~Muramatsu {\\it et al.} " "[CLEO Collaboration],Phys.\\ Rev.\\ Lett.\\ {\\bf 89} (2002) 251802" "[Erratum-ibid.\\ {\\bf 90} (2003) 059901] [arXiv:hep-ex/0207067].\n"); + static Parameter interfacerhoMass + ("RhoMass", + "The mass of the rho0 meson", + &CLEOD0toK0PipPim::mrho_, MeV, 769.3 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacerhoWidth + ("RhoWidth", + "The width of the rho0 meson in D", + &CLEOD0toK0PipPim::wrho_, MeV, 150.2 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceOmegaMass + ("OmegaMass", + "The mass of the omega meson", + &CLEOD0toK0PipPim::momega_, MeV, 782.57*MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceOmegaWidth + ("OmegaWidth", + "The width of the omega meson", + &CLEOD0toK0PipPim::womega_, MeV, 8.44*MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacef980Mass + ("f980Mass", + "The mass of the f_0(980) meson", + &CLEOD0toK0PipPim::mf980_, MeV, 977.00*MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacef980Width + ("f980Width", + "The width of the f_0(980) meson", + &CLEOD0toK0PipPim::wf980_, MeV, 50. *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacegPi + ("gPi", + "The g_pi coupling for the f_0(980) width", + &CLEOD0toK0PipPim::gpi_, 0.09, 0.0, 1., + false, false, Interface::limited); + + static Parameter interfacegK + ("gK", + "The g_K coupling for the f_0(980) width", + &CLEOD0toK0PipPim::gK_, 0.02, 0.0, 1., + false, false, Interface::limited); + + static Switch interfacef0Option + ("f0Option", + "Option for the treatment of the f_0(980) width", + &CLEOD0toK0PipPim::f0opt_, true, false, false); + static SwitchOption interfacef0OptionCoupled + (interfacef0Option, + "Coupled", + "Use the coupling pion and kaon channels", + true); + static SwitchOption interfacef0OptionSWave + (interfacef0Option, + "SWave", + "Use a simple s-wave running width", + false); + + static Parameter interfacef_2Mass + ("f_2Mass", + "The mass of the f_2 meson", + &CLEOD0toK0PipPim::mf2_, MeV, 1275.4 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacef_2Width + ("f_2Width", + "The width of the f_2 meson", + &CLEOD0toK0PipPim::wf2_, MeV, 185.1 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacef1370Mass + ("f1370Mass", + "The mass of the f_0(1370) meson", + &CLEOD0toK0PipPim::mf1370_, MeV, 1310 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacef1370Width + ("f1370Width", + "The width of the f_0(1370) meson", + &CLEOD0toK0PipPim::wf1370_, MeV, 272.0 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstarMass + ("KstarMass", + "The mass of the K*+(892) meson", + &CLEOD0toK0PipPim::mK892_, MeV, 891.66*MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstarWidth + ("KstarWidth", + "The width of the K*+(892) meson", + &CLEOD0toK0PipPim::wK892_, MeV, 50.8 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceK_01430Mass + ("K_01430Mass", + "The mass of the K_0(1430) meson", + &CLEOD0toK0PipPim::mK14300_, MeV, 1412 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceK_01430Width + ("K_01430Width", + "The width of the K_0(1430) meson", + &CLEOD0toK0PipPim::wK14300_, MeV, 294 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceK_21430Mass + ("K_21430Mass", + "The mass of the K_2(1430) meson", + &CLEOD0toK0PipPim::mK14302_, MeV, 1425.6 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + static Parameter interfaceK_21430Width + ("K_21430Width", + "The width of the K_2(1430) meson", + &CLEOD0toK0PipPim::wK14302_, MeV, 98.5 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstar1680Mass + ("Kstar1680Mass", + "The mass of the K*(1680) meson", + &CLEOD0toK0PipPim::mK1680_, MeV, 1717 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstar1680Width + ("Kstar1680Width", + "The width of the K*(1680) meson", + &CLEOD0toK0PipPim::wK1680_, MeV, 322 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKStarPlusAmplitude + ("KStarPlusAmplitude", + "Amplitude for the K*(892)+ component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::aKstarp_, 0.11, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceKStarPlusPhase + ("KStarPlusPhase", + "Phase for the K*(892)+ component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phiKstarp_, 321., 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfaceRhoAmplitude + ("RhoAmplitude", + "Amplitude for the rho0 component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::arho_, 1.0, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceRhoPhase + ("RhoPhase", + "Phase for the rho0 component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phirho_, 0.0, 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfaceOmegaAmplitude + ("OmegaAmplitude", + "Amplitude for the omega component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::aomega_, 0.037, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceOmegaPhase + ("OmegaPhase", + "Phase for the omega component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phiomega_, 114., 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfaceKStarMinusAmplitude + ("KStarMinusAmplitude", + "Amplitude for the K*(892)- component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::aKstarm_, 1.56, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceKStarMinusPhase + ("KStarMinusPhase", + "Phase for the K*(892)- component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phiKstarm_, 150., 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfacef980Amplitude + ("f980Amplitude", + "Amplitude for the f_0(980) component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::af980_, GeV2, 0.34*GeV2, ZERO, 10.0*GeV2, + false, false, Interface::limited); + + static Parameter interfacef980Phase + ("f980Phase", + "Phase for the f_0(980) component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phif980_, 188., 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfacef2Amplitude + ("f2Amplitude", + "Amplitude for the f_2 component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::af2_, 1./GeV2, 0.7/GeV2, ZERO, 10.0/GeV2, + false, false, Interface::limited); + + static Parameter interfacef2Phase + ("f2Phase", + "Phase for the f_0(2) component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phif2_, 308., 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfacef1370Amplitude + ("f1370Amplitude", + "Amplitude for the f_0(1370) component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::af1370_, GeV2, 1.8*GeV2, ZERO, 10.0*GeV2, + false, false, Interface::limited); + + static Parameter interfacef1370Phase + ("f1370Phase", + "Phase for the f_0(1370) component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phif1370_, 85., 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfaceKK_0MinusAmplitude + ("KK_0MinusAmplitude", + "Amplitude for the K_0(1430)- component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::aK14300_, GeV2, 2.0*GeV2, ZERO, 10.0*GeV2, + false, false, Interface::limited); + + static Parameter interfaceKK_0MinusPhase + ("KK_0MinusPhase", + "Phase for the K_0(1430)- component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phiK14300_, 3, 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfaceKK_2MinusAmplitude + ("KK_2MinusAmplitude", + "Amplitude for the K_2(1430)- component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::aK14302_, 1./GeV2, 1.0/GeV2, ZERO, 10.0/GeV2, + false, false, Interface::limited); + + static Parameter interfaceKK_2MinusPhase + ("KK_2MinusPhase", + "Phase for the K_2(1430)- component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phiK14302_, 335, 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfaceK1680MinusAmplitude + ("K1680MinusAmplitude", + "Amplitude for the K*(892)- component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::aK1680_, 5.6, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceK1680MinusPhase + ("K1680MinusPhase", + "Phase for the K*(892)- component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phiK1680_, 174, 0.0, 360.0, + false, false, Interface::limited); + + static Parameter interfaceNonResonantAmplitude + ("NonResonantAmplitude", + "Amplitude for the non-resonant component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::aNR_, 1.1, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceNonResonantPhase + ("NonResonantPhase", + "Phase for the non-resonant component for D0 -> Kbar0 pi+ pi-", + &CLEOD0toK0PipPim::phiNR_, 340, 0.0, 360.0, + false, false, Interface::limited); } void CLEOD0toK0PipPim::doinit() { WeakDalitzDecay::doinit(); static const double degtorad = Constants::pi/180.; - // create the resonances - addResonance(DalitzResonance(getParticleData( 323),0.89166*GeV,0.0508*GeV,0,1,2,-0.11 , 321*degtorad)); - addResonance(DalitzResonance(getParticleData( 113),0.7693 *GeV,0.1502*GeV,1,2,0,-1. , 0 )); - addResonance(DalitzResonance(getParticleData( 223),0.78257*GeV, 8.44*MeV,1,2,0,-0.037 , 114*degtorad)); - addResonance(DalitzResonance(getParticleData( -323),0.89166*GeV,0.0508*GeV,0,2,1,-1.56 , 150*degtorad)); - addResonance(DalitzResonance(getParticleData(9010221),0.977 *GeV,0.050 *GeV,1,2,0, 0.34 , 188*degtorad)); - addResonance(DalitzResonance(getParticleData( 225),1.2754 *GeV,0.1851*GeV,1,2,0, 0.7 , 308*degtorad)); - addResonance(DalitzResonance(getParticleData( 10211),1.310 *GeV,0.2720*GeV,1,2,0, 1.8 , 85*degtorad)); - addResonance(DalitzResonance(getParticleData( -10321),1.412 *GeV,0.294 *GeV,0,2,1, 2.0 , 3*degtorad)); - addResonance(DalitzResonance(getParticleData( -325),1.4256 *GeV,0.0985*GeV,0,2,1, 1.0 , 335*degtorad)); - addResonance(DalitzResonance(getParticleData( -30323),1.717 *GeV,0.322 *GeV,0,2,1,-5.6 , 174*degtorad)); + // non-resonant amplitude + cNR_ = aNR_*Complex(cos(phiNR_*degtorad),sin(phiNR_*degtorad)); + // resonances + addResonance(DalitzResonance(getParticleData( 323),mK892_ , wK892_ ,0,1,2,-aKstarp_ , phiKstarp_*degtorad)); + addResonance(DalitzResonance(getParticleData( 113),mrho_ , wrho_ ,1,2,0,-arho_ , phirho_ *degtorad)); + addResonance(DalitzResonance(getParticleData( 223),momega_ , womega_ ,1,2,0,-aomega_ , phiomega_ *degtorad)); + addResonance(DalitzResonance(getParticleData( -323),mK892_ , wK892_ ,0,2,1,-aKstarm_ , phiKstarm_*degtorad)); + addResonance(DalitzResonance(getParticleData(9010221),mf980_ , wf980_ ,1,2,0, af980_/GeV2 , phif980_ *degtorad)); + addResonance(DalitzResonance(getParticleData( 225),mf2_ , wf2_ ,1,2,0, af2_*GeV2 , phif2_ *degtorad)); + addResonance(DalitzResonance(getParticleData( 10221),mf1370_ , wf1370_ ,1,2,0, af1370_/GeV2 , phif1370_ *degtorad)); + addResonance(DalitzResonance(getParticleData( -10321),mK14300_, wK14300_,0,2,1, aK14300_/GeV2, phiK14300_*degtorad)); + addResonance(DalitzResonance(getParticleData( -325),mK14302_, wK14302_,0,2,1, aK14302_*GeV2, phiK14302_*degtorad)); + addResonance(DalitzResonance(getParticleData( -30323),mK1680_ , wK1680_ ,0,2,1,-aK1680_ , phiK1680_ *degtorad)); // D0 -> K- pi+ pi0 createMode(getParticleData(ParticleID::D0), {getParticleData(ParticleID::Kbar0), getParticleData(ParticleID::piplus), getParticleData(ParticleID::piminus)}); } void CLEOD0toK0PipPim::doinitrun() { WeakDalitzDecay::doinitrun(); } int CLEOD0toK0PipPim::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { int id0(parent->id()); // incoming particle must be D0 if(abs(id0)!=ParticleID::D0) return -1; cc = id0==ParticleID::Dbar0; // must be three decay products if(children.size()!=3) return -1; tPDVector::const_iterator pit = children.begin(); unsigned int npip(0),npim(0),nkm(0),nk0(0),npi0(0); for( ;pit!=children.end();++pit) { id0=(**pit).id(); if(id0 ==ParticleID::piplus) ++npip; else if(id0 ==ParticleID::pi0) ++npi0; else if(id0 ==ParticleID::piminus) ++npim; else if(abs(id0)==ParticleID::K0) ++nk0; else if(id0 ==ParticleID::K_L0) ++nk0; else if(id0 ==ParticleID::K_S0) ++nk0; else if(abs(id0)==ParticleID::Kplus) ++nkm; } if(npim==1&&npip==1&&nk0==1) return 0; else return -1; } Complex CLEOD0toK0PipPim::amplitude(int ichan) const { Complex output(0.); static const Complex ii(0.,1.); unsigned int imin=0, imax(resonances().size()); if(ichan>=0) { imin=ichan; imax=imin+1; } for(unsigned int ix=imin;ixid()!=9010221) { output += resAmp(ix); } // special treatment (Flatte) for f0(980) else { - // output += resAmp(ix); - static const double gpi=0.09, gK=0.02; - const Energy & mAB = mInv(resonances()[ix].daughter1,resonances()[ix].daughter2); - Energy Gamma_pi = gpi*sqrt(0.25*sqr(mAB)-sqr(mOut(resonances()[ix].daughter1))); - Energy2 arg = 0.25*sqr(mAB)-sqr(mOut(resonances()[ix].spectator)); - complex Gamma_K = arg>=ZERO ? gK*sqrt(arg) : gK*ii*sqrt(-arg); - output += resonances()[ix].amp*GeV2/ - (sqr(resonances()[ix].mass)-sqr(mAB)-ii*resonances()[ix].mass*(Gamma_pi+Gamma_K)); + if(!f0opt_) { + output += resAmp(ix); + } + else { + const Energy & mAB = mInv(resonances()[ix].daughter1,resonances()[ix].daughter2); + Energy Gamma_pi = gpi_*sqrt(0.25*sqr(mAB)-sqr(mOut(resonances()[ix].daughter1))); + Energy2 arg = 0.25*sqr(mAB)-sqr(mOut(resonances()[ix].spectator)); + complex Gamma_K = arg>=ZERO ? gK_*sqrt(arg) : gK_*ii*sqrt(-arg); + output += resonances()[ix].amp*GeV2/ + (sqr(resonances()[ix].mass)-sqr(mAB)-ii*resonances()[ix].mass*(Gamma_pi+Gamma_K)); + } } } - if(ichan<0) output += 1.1*Complex(cos(5.93411945678072),sin(5.93411945678072)); + if(ichan<0) output += cNR_; return output; } diff --git a/Decay/ScalarMeson/CLEOD0toK0PipPim.h b/Decay/ScalarMeson/CLEOD0toK0PipPim.h --- a/Decay/ScalarMeson/CLEOD0toK0PipPim.h +++ b/Decay/ScalarMeson/CLEOD0toK0PipPim.h @@ -1,120 +1,351 @@ // -*- C++ -*- #ifndef Herwig_CLEOD0toK0PipPim_H #define Herwig_CLEOD0toK0PipPim_H // // This is the declaration of the CLEOD0toK0PipPim class. // #include "WeakDalitzDecay.h" namespace Herwig { using namespace ThePEG; /** * The CLEOD0toK0PipPim class implements the Dalitz plot fit * of the CLEO collaboration for \f$D^0\to\bar{K}^0\pi^+\pi^-\f$, * Phys. Rev. Lett. 89 (2002) 251802 * * @see \ref CLEOD0toK0PipPimInterfaces "The interfaces" * defined for CLEOD0toK0PipPim. */ class CLEOD0toK0PipPim: public WeakDalitzDecay { public: /** * The default constructor. */ CLEOD0toK0PipPim(); /** * Which of the possible decays is required * @param cc Is this mode the charge conjugate * @param parent The decaying particle * @param children The decay products */ virtual int modeNumber(bool & cc, tcPDPtr parent, const tPDVector & children) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); //@} protected: /** * Calculate the amplitude */ virtual Complex amplitude(int ichan) const; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - CLEOD0toK0PipPim & operator=(const CLEOD0toK0PipPim &); + CLEOD0toK0PipPim & operator=(const CLEOD0toK0PipPim &) = delete; +private: + + /** + * Mass, Widths and related parameters + */ + //@{ + /** + * Mass of the \f$\rho(770)\f$ + */ + Energy mrho_; + + /** + * Width of the \f$\rho(770)\f$ + */ + Energy wrho_; + + /** + * Mass of the \f$\omega\f$ + */ + Energy momega_; + + /** + * Width of the \f$\omega\f$ + */ + Energy womega_; + + /** + * Mass of the \f$f_0(980)\f$ + */ + Energy mf980_; + + /** + * Width of the \f$f_0(980)\f$ + */ + Energy wf980_; + + /** + * \f$g_\pi\f$ coupling for the \f$f_0(980)\f$ width + */ + double gpi_; + + /** + *\f$g_K\f$ coupling for the \f$f_0(980)\f$ width + */ + double gK_; + + /** + * Option for handling the width of the \f$f_0(980)\f$ + */ + bool f0opt_; + + /** + * Mass of the \f$f_2(1270)\f$ + */ + Energy mf2_; + + /** + * Width of the \f$f_2(1270)\f$ + */ + Energy wf2_; + + /** + * Mass of the \f$f_0(1370)\f$ + */ + Energy mf1370_; + + /** + * Width of the \f$f_0(1370)\f$ + */ + Energy wf1370_; + + /** + * Mass of the \f$K^{*+}(892)\f$ + */ + Energy mK892_; + + /** + * Width of the \f$K^{*+}(892)\f$ + */ + Energy wK892_; + /** + * Mass of the \f$K_0^*(1430)\f$ + */ + Energy mK14300_; + + /** + * Width of the \f$K_0^*(1430)\f$ + */ + Energy wK14300_; + + /** + * Mass of the \f$K_2^*(1430)\f$ + */ + Energy mK14302_; + + /** + * Width of the \f$K_2^*(1430)\f$ + */ + Energy wK14302_; + + /** + * Mass of the \f$K^*(1680)\f$ + */ + Energy mK1680_; + + /** + * Width of the \f$K^*(1680)\f$ + */ + Energy wK1680_; + //@} + + /** + * Magnitudes and phases of the amplitudes for \f$D^0\to \bar{K}^0\pi^+\pi^-\f$ + */ + //@{ + /** + * Amplitude for the \f$K^{*+}\f$ + */ + double aKstarp_; + + /** + * Phase for the \f$K^{*+}\f$ + */ + double phiKstarp_; + + /** + * Amplitude for the \f$\rho^0(770)\f$ + */ + double arho_; + + /** + * Phase for the \f$\rho^0(770)\f$ + */ + double phirho_; + + /** + * Amplitude for the \f$\omega\f$ + */ + double aomega_; + + /** + * Phase for the \f$\omega\f$ + */ + double phiomega_; + + /** + * Amplitude for the \f$K^{*-}\f$ + */ + double aKstarm_; + + /** + * Phase for the \f$K^{*-}\f$ + */ + double phiKstarm_; + + /** + * Amplitude for the \f$f_0(980)\f$ + */ + Energy2 af980_; + + /** + * Phase for the \f$f_0(980)\f$ + */ + double phif980_; + + /** + * Amplitude for the \f$f_2(1270)\f$ + */ + InvEnergy2 af2_; + + /** + * Phase for the \f$f_2(1270)\f$ + */ + double phif2_; + + /** + * Amplitude for the \f$f_0(1370)\f$ + */ + Energy2 af1370_; + + /** + * Phase for the \f$f_0(1370)\f$ + */ + double phif1370_; + + /** + * Amplitude for the \f$K^*_0(1430)^-\f$ + */ + Energy2 aK14300_; + + /** + * Phase for the \f$K^*_0(1430)^-\f$ + */ + double phiK14300_; + + /** + * Amplitude for the \f$K^*_2(1430)^-\f$ + */ + InvEnergy2 aK14302_; + + /** + * Phase for the \f$K^*_2(1430)^-\f$ + */ + double phiK14302_; + + /** + * Amplitude for the \f$K^*(1680)^-\f$ + */ + double aK1680_; + + /** + * Phase for the \f$K^*(1680)^-\f$ + */ + double phiK1680_; + + /** + * Amplitude of the non-resonant component + */ + double aNR_; + + /** + * Phase of the non=resonant component + */ + double phiNR_; + + /** + * Complex amplitude of the non-resonant component + */ + Complex cNR_; + //@} + }; } #endif /* Herwig_CLEOD0toK0PipPim_H */ diff --git a/Decay/ScalarMeson/CLEOD0toKmPipPi0.cc b/Decay/ScalarMeson/CLEOD0toKmPipPi0.cc --- a/Decay/ScalarMeson/CLEOD0toKmPipPi0.cc +++ b/Decay/ScalarMeson/CLEOD0toKmPipPi0.cc @@ -1,104 +1,308 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the CLEOD0toKmPipPi0 class. // #include "CLEOD0toKmPipPi0.h" #include "ThePEG/Interface/ClassDocumentation.h" +#include "ThePEG/Interface/Parameter.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" using namespace Herwig; -CLEOD0toKmPipPi0::CLEOD0toKmPipPi0() : WeakDalitzDecay(5./GeV,1.5/GeV,true) -{} +CLEOD0toKmPipPi0::CLEOD0toKmPipPi0() : WeakDalitzDecay(5./GeV,1.5/GeV,true) { + // masses and widths + mrho_ = 770 *MeV; wrho_ = 150.7 *MeV; + mK892_ = 891.5 *MeV; wK892_ = 50 *MeV; + mK8920_ = 896.1 *MeV; wK8920_ = 50.5 *MeV; + mK14300_ = 1412 *MeV; wK14300_ = 294 *MeV; + mrho1700_ = 1700 *MeV; wrho1700_ = 240 *MeV; + mK1680_ = 1717 *MeV; wK1680_ = 322 *MeV; + // amplitudes and phases for D0 -> K-pi+pi0 + aNR_ = 1.75 ; phiNR_ = 31.2; + arho_ = 1.00 ; phirho_ = 0. ; + aKstarm_ = 0.44 ; phiKstarm_ = 163 ; + aKstar0_ = 0.39 ; phiKstar0_ = -0.2; + aK1430m_ = 0.77*GeV2; phiK1430m_ = 55.5; + aK14300_ = 0.85*GeV2; phiK14300_ = 166 ; + arho1700_ = 2.50 ; phirho1700_ = 171 ; + aK1680_ = 2.50 ; phiK1680_ = 103 ; + // intermediates + generateIntermediates(true); +} IBPtr CLEOD0toKmPipPi0::clone() const { return new_ptr(*this); } IBPtr CLEOD0toKmPipPi0::fullclone() const { return new_ptr(*this); } // The following static variable is needed for the type // description system in ThePEG. -DescribeNoPIOClass +DescribeClass describeHerwigCLEOD0toKmPipPi0("Herwig::CLEOD0toKmPipPi0", "HwSMDecay.so"); +void CLEOD0toKmPipPi0::persistentOutput(PersistentOStream & os) const { + os << ounit(mK14300_,GeV) << ounit(wK14300_,GeV) << ounit(mK1680_,GeV) << ounit(wK1680_,GeV) + << ounit(mrho1700_,GeV) << ounit(wrho1700_,GeV) << ounit(mK8920_,GeV) << ounit(wK8920_,GeV) + << ounit(mK892_,GeV) << ounit(wK892_,GeV) << ounit(mrho_,GeV) << ounit(wrho_,GeV) + << aNR_ << phiNR_ << arho_ << phirho_ << aKstarm_ << phiKstarm_ << aKstar0_ << phiKstar0_ + << ounit(aK1430m_,GeV2) << phiK1430m_ << ounit(aK14300_,GeV2) << phiK14300_ + << arho1700_ << phirho1700_ << aK1680_ << phiK1680_ << cNR_; +} + +void CLEOD0toKmPipPi0::persistentInput(PersistentIStream & is, int) { + is >> iunit(mK14300_,GeV) >> iunit(wK14300_,GeV) >> iunit(mK1680_,GeV) >> iunit(wK1680_,GeV) + >> iunit(mrho1700_,GeV) >> iunit(wrho1700_,GeV) >> iunit(mK8920_,GeV) >> iunit(wK8920_,GeV) + >> iunit(mK892_,GeV) >> iunit(wK892_,GeV) >> iunit(mrho_,GeV) >> iunit(wrho_,GeV) + >> aNR_ >> phiNR_ >> arho_ >> phirho_ >> aKstarm_ >> phiKstarm_ >> aKstar0_ >> phiKstar0_ + >> iunit(aK1430m_,GeV2) >> phiK1430m_ >> iunit(aK14300_,GeV2) >> phiK14300_ + >> arho1700_ >> phirho1700_ >> aK1680_ >> phiK1680_ >> cNR_; +} + void CLEOD0toKmPipPi0::Init() { static ClassDocumentation documentation ("The CLEOD0toKmPipPi0 class implements the model of CLEO for " "D0 -> K- pi+ pi0, Phys. Rev. D63 (2001) 092001.", "The CLEO fit of \\cite{Kopp:2000gv} was" " used for the decay $D^0\\to K^-\\pi^+\\pi^0$.", "\\bibitem{Kopp:2000gv} S.~Kopp {\\it et al.} [CLEO Collaboration], " "Phys.\\ Rev.\\ D {\\bf 63} (2001) 092001 [arXiv:hep-ex/0011065]."); + static Parameter interfaceK_01430Mass + ("K_01430Mass", + "The mass of the K_0(1430) meson", + &CLEOD0toKmPipPi0::mK14300_, MeV, 1412 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceK_01430Width + ("K_01430Width", + "The width of the K_0(1430) meson", + &CLEOD0toKmPipPi0::wK14300_, MeV, 294 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstar1680Mass + ("Kstar1680Mass", + "The mass of the K*(1680) meson", + &CLEOD0toKmPipPi0::mK1680_, MeV, 1717 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstar1680Width + ("Kstar1680Width", + "The width of the K*(1680) meson", + &CLEOD0toKmPipPi0::wK1680_, MeV, 322 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacerho1700Mass + ("rho1700Mass", + "The mass of the rho(1700) meson", + &CLEOD0toKmPipPi0::mrho1700_, MeV, 1700 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacerho1700Width + ("rho1700Width", + "The width of the rho(1700) meson", + &CLEOD0toKmPipPi0::wrho1700_, MeV, 240 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstar0892Mass + ("Kstar0892Mass", + "The mass of the K*0(892) meson", + &CLEOD0toKmPipPi0::mK8920_, MeV, 896.1 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstar0892Width + ("Kstar0892Width", + "The width of the K*0(892) meson", + &CLEOD0toKmPipPi0::wK8920_, MeV, 50.5 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstarPlus892Mass + ("KstarPlus892Mass", + "The mass of the K*+(892) meson", + &CLEOD0toKmPipPi0::mK892_, MeV, 891.5 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceKstarPlus892Width + ("KstarPlus892Width", + "The width of the K*+(892) meson in D0 -> K-pi+pi0", + &CLEOD0toKmPipPi0::wK892_, MeV, 50 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacerhoMass + ("RhoMass", + "The mass of the rho+ meson", + &CLEOD0toKmPipPi0::mrho_, MeV, 770 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfacerhoWidth + ("RhoWidth", + "The width of the rho+ meson", + &CLEOD0toKmPipPi0::wrho_, MeV, 150.7 *MeV, ZERO, 10000.0*MeV, + false, false, Interface::limited); + + static Parameter interfaceChargedNonResonantAmplitude + ("ChargedNonResonantAmplitude", + "Amplitude for the non-resonant component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::aNR_, 1.75, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceChargedNonResonantPhase + ("ChargedNonResonantPhase", + "Phase for the non-resonant component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::phiNR_, 31.2, -180.0, 180.0, + false, false, Interface::limited); + + static Parameter interfaceChargedRhoAmplitude + ("ChargedRhoAmplitude", + "Amplitude for the rho+ component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::arho_, 1.0, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceChargedRhoPhase + ("ChargedRhoPhase", + "Phase for the rho+ component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::phirho_, 0.0, -180.0, 180.0, + false, false, Interface::limited); + + static Parameter interfaceChargedKStarMinusAmplitude + ("ChargedKStarMinusAmplitude", + "Amplitude for the K*(892)- component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::aKstarm_, 0.44, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceChargedKStarMinusPhase + ("ChargedKStarMinusPhase", + "Phase for the K*(892)- component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::phiKstarm_, 163, -180.0, 180.0, + false, false, Interface::limited); + + static Parameter interfaceChargedKStar0Amplitude + ("ChargedKStar0Amplitude", + "Amplitude for the K*(892)0 component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::aKstar0_, 0.39, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceChargedKStar0Phase + ("ChargedKStar0Phase", + "Phase for the K*(892)0 component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::phiKstar0_, -0.2, -180.0, 180.0, + false, false, Interface::limited); + + static Parameter interfaceChargedK_0MinusAmplitude + ("ChargedK_0MinusAmplitude", + "Amplitude for the K_0(1430)- component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::aK1430m_, GeV2, 0.77*GeV2, ZERO, 10.0*GeV2, + false, false, Interface::limited); + + static Parameter interfaceChargedK_0MinusPhase + ("ChargedK_0MinusPhase", + "Phase for the K_0(1430)- component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::phiK1430m_, 55.5, -180.0, 180.0, + false, false, Interface::limited); + + static Parameter interfaceChargedK_00Amplitude + ("ChargedK_00Amplitude", + "Amplitude for the K_0(1430)0 component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::aK14300_, GeV2, 0.85*GeV2, ZERO, 10.0*GeV2, + false, false, Interface::limited); + + static Parameter interfaceChargedK_00Phase + ("ChargedK_00Phase", + "Phase for the K_0(1430)0 component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::phiK14300_, 166, -180.0, 180.0, + false, false, Interface::limited); + + static Parameter interfaceChargedRho1700Amplitude + ("ChargedRho1700Amplitude", + "Amplitude for the rho1700+ component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::arho1700_, 2.5, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceChargedRho1700Phase + ("ChargedRho1700Phase", + "Phase for the rho1700+ component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::phirho1700_, 171., -180.0, 180.0, + false, false, Interface::limited); + + static Parameter interfaceChargedK1680MinusAmplitude + ("ChargedK1680MinusAmplitude", + "Amplitude for the K*(1680)- component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::aK1680_, 2.5, 0.0, 10.0, + false, false, Interface::limited); + + static Parameter interfaceChargedK1680MinusPhase + ("ChargedK1680MinusPhase", + "Phase for the K*(1680)- component for D0 -> K- pi+ pi0", + &CLEOD0toKmPipPi0::phiK1680_, 103, -180.0, 180.0, + false, false, Interface::limited); } void CLEOD0toKmPipPi0::doinit() { WeakDalitzDecay::doinit(); - static const double degtorad = Constants::pi/180.; + static const double degtorad = Constants::pi/180.; + // non-resonant amplitude + cNR_ = aNR_*Complex(cos(phiNR_*degtorad),sin(phiNR_*degtorad)); // create the resonances - addResonance(DalitzResonance(getParticleData( 213) , 0.770*GeV,0.1507*GeV,1,2,0,-1. , 0. )); - addResonance(DalitzResonance(getParticleData(-323) ,0.8915*GeV,0.050 *GeV,0,2,1,-0.44, 163 *degtorad)); - addResonance(DalitzResonance(getParticleData(-313) ,0.8961*GeV,0.0505*GeV,0,1,2,-0.39,- 0.2*degtorad)); - addResonance(DalitzResonance(getParticleData(-10321),1.412 *GeV,0.294*GeV ,0,2,1, 0.77, 55.5*degtorad)); - addResonance(DalitzResonance(getParticleData(-10311),1.412 *GeV,0.294*GeV ,0,1,2, 0.85, 166. *degtorad)); - addResonance(DalitzResonance(getParticleData( 30213),1.717 *GeV,0.322*GeV ,1,2,0,-2.50, 171 *degtorad)); - addResonance(DalitzResonance(getParticleData(-30323),1.717 *GeV,0.322*GeV ,0,2,1,-2.50, 103 *degtorad)); + addResonance(DalitzResonance(getParticleData( 213) , mrho_ , wrho_ ,1,2,0,-arho_ , phirho_ *degtorad)); + addResonance(DalitzResonance(getParticleData(-323) , mK892_ , wK892_ ,0,2,1,-aKstarm_ , phiKstarm_ *degtorad)); + addResonance(DalitzResonance(getParticleData(-313) ,mK8920_ , wK8920_ ,0,1,2,-aKstar0_ , phiKstar0_ *degtorad)); + addResonance(DalitzResonance(getParticleData(-10321),mK14300_ , wK14300_ ,0,2,1, aK1430m_/GeV2, phiK1430m_ *degtorad)); + addResonance(DalitzResonance(getParticleData(-10311),mK14300_ , wK14300_ ,0,1,2, aK14300_/GeV2, phiK14300_ *degtorad)); + addResonance(DalitzResonance(getParticleData( 30213),mrho1700_, wrho1700_,1,2,0,-arho1700_ , phirho1700_*degtorad)); + addResonance(DalitzResonance(getParticleData(-30323),mK1680_ , wK1680_ ,0,2,1,-aK1680_ , phiK1680_ *degtorad)); // D+ -> K- pi+ pi+ createMode(getParticleData(ParticleID::D0), {getParticleData(ParticleID::Kminus), getParticleData(ParticleID::piplus), getParticleData(ParticleID::pi0)}); } void CLEOD0toKmPipPi0::doinitrun() { WeakDalitzDecay::doinitrun(); } int CLEOD0toKmPipPi0::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { int id0(parent->id()); // incoming particle must be D0 or D+ if(abs(id0)!=ParticleID::D0) return -1; cc = id0<0; // must be three decay products if(children.size()!=3) return -1; tPDVector::const_iterator pit = children.begin(); unsigned int npip(0),npim(0),nkm(0),npi0(0); int id; for( ;pit!=children.end();++pit) { id=(**pit).id(); if(id ==ParticleID::piplus) ++npip; else if(id ==ParticleID::pi0) ++npi0; else if(id ==ParticleID::piminus) ++npim; else if(abs(id)==ParticleID::Kplus) ++nkm; } if(nkm==1&&(npip+npim)==1&&npi0==1) return 0; else return -1; } Complex CLEOD0toKmPipPi0::amplitude(int ichan) const { Complex output(0.); unsigned int imin=0, imax(resonances().size()); if(ichan>=0) { imin=ichan; imax=imin+1; } for(unsigned int ix=imin;ix> resonances_ >> maxWgt_ >> weights_; + is >> resonances_ >> maxWgt_ >> weights_ + >> iunit(rParent_,1./GeV) >> iunit(rResonance_,1./GeV); } // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractClass describeHerwigWeakDalitzDecay("Herwig::WeakDalitzDecay", "HwSMDecay.so"); void WeakDalitzDecay::Init() { static ClassDocumentation documentation ("The WeakDalitzDecay class provides a base class for " "weak three-body decays of bottom and charm mesons"); + static Parameter interfaceDRadius + ("DRadius", + "The radius parameter for the Blatt-Weisskopf form-factor for the D", + &WeakDalitzDecay::rParent_, 1./GeV, 5./GeV, ZERO, 10./GeV, + false, false, Interface::limited); + + static Parameter interfaceResonanceRadius + ("ResonanceRadius", + "The radius parameter for the Blatt-Weisskopf form-factor for the" + "intermediate resonances", + &WeakDalitzDecay::rResonance_, 1./GeV, 1.5/GeV, ZERO, 10./GeV, + false, false, Interface::limited); + + static Parameter interfaceMaximumWeight + ("MaximumWeight", + "The maximum weight for the phase-space sampling", + &WeakDalitzDecay::maxWgt_, 1.0, 0.0, 1e10, + false, false, Interface::limited); + + static ParVector interfaceWeights + ("Weights", + "The weights for the different channels for the phase-space integration", + &WeakDalitzDecay::weights_, -1, 1.0, 0.0, 1.0, + false, false, Interface::limited); } void WeakDalitzDecay:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // set up the spin information for the decay products ScalarWaveFunction::constructSpinInfo(const_ptr_cast(&part), incoming,true); for(unsigned int ix=0;ix<3;++ix) ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true); } void WeakDalitzDecay::doinit() { DecayIntegrator::doinit(); } void WeakDalitzDecay::doinitrun() { DecayIntegrator::doinitrun(); weights_.resize(mode(0)->channels().size()); maxWgt_ = mode(0)->maxWeight(); for(unsigned int iz=0;izchannels().size();++iz) { weights_[iz]=mode(0)->channels()[iz].weight(); } } double WeakDalitzDecay::me2(const int ichan, const Particle & part, const tPDVector & , const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin0,PDT::Spin0))); useMe(); if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(_rho,const_ptr_cast(&part),incoming); } // set the kinematics mD_ = part.mass(); for(unsigned int ix=0;ix(resonances_.size(),1./double(resonances_.size())); } unsigned int ix=0; for(DalitzResonance res : resonances_) { mode->addChannel((PhaseSpaceChannel(mode),0,res.resonance,0,res.spectator+1,1,res.daughter1+1,1,res.daughter2+1)); resetIntermediate(res.resonance,res.mass,res.width); ++ix; } addMode(mode); } Complex WeakDalitzDecay::resAmp(unsigned int i, bool gauss) const { Complex output = resonances_[i].amp; // locations of the outgoing particles const unsigned int &d1 = resonances_[i].daughter1; const unsigned int &d2 = resonances_[i].daughter2; const unsigned int &sp = resonances_[i].spectator; // mass and width of the resonance const Energy & mR = resonances_[i].mass ; const Energy & wR = resonances_[i].width; // momenta for the resonance decay // off-shell Energy pAB=sqrt(0.25*sqr(sqr(m2_[d1][d2]) -sqr(mOut_[d1])-sqr(mOut_[d2])) - sqr(mOut_[d1]*mOut_[d2]))/m2_[d1][d2]; // on-shell Energy pR=sqrt(0.25*sqr( mR*mR -sqr(mOut_[d1])-sqr(mOut_[d2])) - sqr(mOut_[d1]*mOut_[d2]))/mR; // for the D decay Energy pD = sqrt(max(ZERO,(0.25*sqr(sqr(mD_)-sqr(mR)-sqr(mOut_[sp])) - sqr(mR*mOut_[sp]))/sqr(mD_))); Energy pDAB= sqrt( 0.25*sqr(sqr(mD_)-sqr(m2_[d1][d2])-sqr(mOut_[sp])) - sqr(m2_[d1][d2]*mOut_[sp]))/mD_; // Blatt-Weisskopf factors double fR=1, fD=1; unsigned int power(1); double r1A(rResonance_*pR),r1B(rResonance_*pAB ); double r2A(rParent_ *pD),r2B(rParent_ *pDAB); // mass for thre denominator Energy mDen = useResonanceMass_ ? mR : m2_[d1][d2]; // Blatt-Weisskopf factors and spin piece switch (resonances_[i].resonance->iSpin()) { case PDT::Spin0: if(gauss) { fR = exp(-(r1B-r1A)/12.); fD = exp(-(r2B-r2A)/12.); } break; case PDT::Spin1: fR=sqrt( (1. + sqr(r1A)) / (1. + sqr(r1B)) ); fD=sqrt( (1. + sqr(r2A)) / (1. + sqr(r2B)) ); power=3; output *= fR*fD*(sqr(m2_[d2][sp])-sqr(m2_[d1][sp]) + ( sqr(mD_)-sqr(mOut_[sp]))*(sqr(mOut_[d1])-sqr(mOut_[d2]))/sqr(mDen) )/GeV2; break; case PDT::Spin2: fR = sqrt( (9. + sqr(r1A)*(3.+sqr(r1A))) / (9. + sqr(r1B)*(3.+sqr(r1B)))); fD = sqrt( (9. + sqr(r2A)*(3.+sqr(r2A))) / (9. + sqr(r2B)*(3.+sqr(r2B)))); power=5; output *= fR*fD/GeV2/GeV2*( sqr( sqr(m2_[d2][sp]) - sqr(m2_[d1][sp]) +(sqr(mD_)-sqr(mOut_[sp]))*(sqr(mOut_[d1])-sqr(mOut_[d2]))/(sqr(mDen))) - (sqr(m2_[d1][d2])-2* sqr(mD_)-2*sqr(mOut_[sp]) + sqr((sqr( mD_) - sqr(mOut_[sp]))/mDen))* (sqr(m2_[d1][d2])-2*sqr(mOut_[d1])-2*sqr(mOut_[d2]) + sqr((sqr(mOut_[d1]) - sqr(mOut_[d2]))/mDen))/3.); break; default : assert(false); } // multiply by Breit-Wigner piece and return Energy gam = wR*pow(pAB/pR,power)*(mR/m2_[d1][d2])*fR*fR; return output*GeV2/(sqr(mR)-sqr(m2_[d1][d2])-mR*gam*Complex(0.,1.)); }