diff --git a/Decay/ScalarMeson/CLEODptoKmPipPip.cc b/Decay/ScalarMeson/CLEODptoKmPipPip.cc --- a/Decay/ScalarMeson/CLEODptoKmPipPip.cc +++ b/Decay/ScalarMeson/CLEODptoKmPipPip.cc @@ -1,114 +1,112 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the CLEODptoKmPipPip class. // #include "CLEODptoKmPipPip.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" using namespace Herwig; CLEODptoKmPipPip::CLEODptoKmPipPip() : WeakDalitzDecay(5./GeV,1.5/GeV,false) {} IBPtr CLEODptoKmPipPip::clone() const { return new_ptr(*this); } IBPtr CLEODptoKmPipPip::fullclone() const { return new_ptr(*this); } void CLEODptoKmPipPip::doinit() { WeakDalitzDecay::doinit(); static const double degtorad = Constants::pi/180.; // create the resonances addResonance(DalitzResonance(getParticleData(-313 ), 896 *MeV, 50.3*MeV,0,1,2,-1. , 0. )); addResonance(DalitzResonance(getParticleData(-313 ), 896 *MeV, 50.3*MeV,0,2,1,-1. , 0. )); addResonance(DalitzResonance(getParticleData(-10311 ),1463 *MeV,163.8*MeV,0,1,2,3. , 49.7*degtorad)); addResonance(DalitzResonance(getParticleData(-10311 ),1463 *MeV,163.8*MeV,0,2,1,3. , 49.7*degtorad)); addResonance(DalitzResonance(getParticleData(-315 ),1432.4*MeV, 109*MeV,0,1,2,0.962, -29.9*degtorad)); addResonance(DalitzResonance(getParticleData(-315 ),1432.4*MeV, 109*MeV,0,2,1,0.962, -29.9*degtorad)); addResonance(DalitzResonance(getParticleData(-30313 ),1717 *MeV, 322*MeV,0,1,2,-6.5 , 29.0*degtorad)); addResonance(DalitzResonance(getParticleData(-30313 ),1717 *MeV, 322*MeV,0,2,1,-6.5 , 29.0*degtorad)); addResonance(DalitzResonance(getParticleData(-9000311), 809 *MeV, 470*MeV,0,1,2,5.01 ,-163.7*degtorad)); addResonance(DalitzResonance(getParticleData(-9000311), 809 *MeV, 470*MeV,0,2,1,5.01 ,-163.7*degtorad)); // D+ -> K- pi+ pi+ createMode(getParticleData(ParticleID::Dplus), {getParticleData(ParticleID::Kminus), getParticleData(ParticleID::piplus), getParticleData(ParticleID::piplus)}); } void CLEODptoKmPipPip::doinitrun() { WeakDalitzDecay::doinitrun(); } void CLEODptoKmPipPip::persistentOutput(PersistentOStream & os) const { } void CLEODptoKmPipPip::persistentInput(PersistentIStream & is, int) { } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigCLEODptoKmPipPip("Herwig::CLEODptoKmPipPip", "HwSMDecay.so"); void CLEODptoKmPipPip::Init() { static ClassDocumentation documentation ("There is no documentation for the CLEODptoKmPipPip class"); } int CLEODptoKmPipPip::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { int id0(parent->id()); // incoming particle must be D+ if(abs(id0)!=ParticleID::Dplus) return -1; cc = id0<0; // must be three decay products if(children.size()!=3) return -1; unsigned int npip(0),npim(0),nk(0); for(tPDPtr child : children) { long id= child->id(); if(id ==ParticleID::piplus) ++npip; else if(id ==ParticleID::piminus) ++npim; else if(abs(id)==ParticleID::Kplus) ++nk; } if((id0==ParticleID::Dplus &&(nk==1&&npip==2))|| (id0==ParticleID::Dminus&&(nk==1&&npim==2))) return 0; else return -1; } Complex CLEODptoKmPipPip::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 K- pi+ pi+ + // createMode(getParticleData(ParticleID::Dplus), + // {getParticleData(ParticleID::Kminus), + // getParticleData(ParticleID::piplus), + // getParticleData(ParticleID::piplus)}); +} + +void FOCUSDptoKmPipPip::doinitrun() { + WeakDalitzDecay::doinitrun(); +} + +void FOCUSDptoKmPipPip::persistentOutput(PersistentOStream & os) const { os << ounit(beta_,GeV) << theta_ << gamma_ << ounit(g1_,GeV) << ounit(g2_,GeV) << c1_ << c2_ << c3_ - << KHalf_ << KThreeHalf_ << maxWgt_ << weights_ - << ounit(rD0_,1./GeV) << ounit(rRes_,1./GeV) - << ounit(mRes_,GeV) << ounit(wRes_,GeV); + << KHalf_ << KThreeHalf_; } -void DtoKPiPiFOCUS::persistentInput(PersistentIStream & is, int) { +void FOCUSDptoKmPipPip::persistentInput(PersistentIStream & is, int) { is >> iunit(beta_,GeV) >> theta_ >> gamma_ >> iunit(g1_,GeV) >> iunit(g2_,GeV) >> c1_ >> c2_ >> c3_ - >> KHalf_ >> KThreeHalf_ >> maxWgt_ >> weights_ - >> iunit(rD0_,1./GeV) >> iunit(rRes_,1./GeV) - >> iunit(mRes_,GeV) >> iunit(wRes_,GeV); + >> KHalf_ >> KThreeHalf_; } - // The following static variable is needed for the type // description system in ThePEG. -DescribeClass -describeHerwigDtoKPiPiFOCUS("Herwig::DtoKPiPiFOCUS", "HwSMDecay.so"); +DescribeClass +describeHerwigFOCUSDptoKmPipPip("Herwig::FOCUSDptoKmPipPip", "HwSMDecay.so"); -void DtoKPiPiFOCUS::Init() { +void FOCUSDptoKmPipPip::Init() { - static ClassDocumentation documentation - ("There is no documentation for the DtoKPiPiFOCUS class"); - - static Reference interfaceKHalf + static ClassDocumentation documentation + ("The FOCUSDptoKmPipPip class implements the fits of FOCUS, " + "Phys.Lett. B653 (2007) 1-11 for the decay D+ -> K-pi+pi-.", + "The decay $D^+\\to K^-\\pi^+\\pi^+$ was simulated using the" + " model of FOCUS \\cite{Pennington:2007se}.", + "\\bibitem{Pennington:2007se}\n" + "J.~M.~Link {\\it et al.} [FOCUS Collaboration],\n" + "%``Dalitz plot analysis of the $D^{+} \\to K^{-} \\pi^{+} \\pi^{+}$ decay in the FOCUS experiment,''\n" + "Phys.\\ Lett.\\ B {\\bf 653} (2007) 1\n" + "doi:10.1016/j.physletb.2007.06.070\n" + "[arXiv:0705.2248 [hep-ex]].\n"); + + static Reference interfaceKHalf ("KHalf", "The K-matrix for the I=1/2 s-wave component", - &DtoKPiPiFOCUS::KHalf_, false, false, true, false, false); + &FOCUSDptoKmPipPip::KHalf_, false, false, true, false, false); - static Reference interfaceKThreeHalf + static Reference interfaceKThreeHalf ("KThreeHalf", "The K-matrix for the I=1/2 s-wave component", - &DtoKPiPiFOCUS::KThreeHalf_, false, false, true, false, false); + &FOCUSDptoKmPipPip::KThreeHalf_, false, false, true, false, false); - static Parameter interfaceg1 + static Parameter interfaceg1 ("g1", "The coupling of the first channel to K*0", - &DtoKPiPiFOCUS::g1_, 0.31072*GeV, -10.*GeV, 10.*GeV, + &FOCUSDptoKmPipPip::g1_, 0.31072*GeV, -10.*GeV, 10.*GeV, false, false, Interface::limited); - static Parameter interfaceg2 + static Parameter interfaceg2 ("g2", "The coupling of the first channel to K*0", - &DtoKPiPiFOCUS::g2_, -0.02323*GeV, -10.*GeV, 10.*GeV, + &FOCUSDptoKmPipPip::g2_, -0.02323*GeV, -10.*GeV, 10.*GeV, false, false, Interface::limited); - static Parameter interfacebeta + static Parameter interfacebeta ("beta", "The beta coefficient for the P-vector", - &DtoKPiPiFOCUS::beta_, 3.389*GeV, -10.*GeV, 10.*GeV, + &FOCUSDptoKmPipPip::beta_, 3.389*GeV, -10.*GeV, 10.*GeV, false, false, Interface::limited); - static Parameter interfacetheta + static Parameter interfacetheta ("theta", "The theta phase (in degrees) for the P-vector", - &DtoKPiPiFOCUS::theta_, 286, 0.0, 360., + &FOCUSDptoKmPipPip::theta_, 286, 0.0, 360., false, false, Interface::limited); - static ParVector interfacegamma + static ParVector interfacegamma ("gamma", "The gamma phases in degress", - &DtoKPiPiFOCUS::gamma_, 3, 0., 0.0, 360.0, + &FOCUSDptoKmPipPip::gamma_, 3, 0., 0.0, 360.0, false, false, Interface::limited); - static ParVector interfacec1 + static ParVector interfacec1 ("c1", "The c1 coefficients for the P-vector", - &DtoKPiPiFOCUS::c1_, -1, 0., -100., 100., + &FOCUSDptoKmPipPip::c1_, -1, 0., -100., 100., false, false, Interface::limited); - static ParVector interfacec2 + static ParVector interfacec2 ("c2", "The c2 coefficients for the P-vector", - &DtoKPiPiFOCUS::c2_, -1, 0., -100., 100., + &FOCUSDptoKmPipPip::c2_, -1, 0., -100., 100., false, false, Interface::limited); - static ParVector interfacec3 + static ParVector interfacec3 ("c3", "The c3 coefficients for the P-vector", - &DtoKPiPiFOCUS::c3_, -1, 0., -100., 100., - false, false, Interface::limited); - - static Parameter interfaceDRadius - ("DRadius", - "The radius parameter for the Blatt-Weisskopf form-factor for the D", - &DtoKPiPiFOCUS::rD0_, 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", - &DtoKPiPiFOCUS::rRes_, 1./GeV, 1.5/GeV, ZERO, 10./GeV, - false, false, Interface::limited); - - static ParVector interfaceMasses - ("Masses", - "The masses of the resonances", - &DtoKPiPiFOCUS::mRes_, GeV, -1, 1.0*GeV, 0.0*GeV, 10.0*GeV, - false, false, Interface::limited); - - static ParVector interfaceWidths - ("Widths", - "The widths of the resonances", - &DtoKPiPiFOCUS::wRes_, GeV, -1, 1.0*GeV, 0.0*GeV, 10.0*GeV, + &FOCUSDptoKmPipPip::c3_, -1, 0., -100., 100., false, false, Interface::limited); } -int DtoKPiPiFOCUS::modeNumber(bool & cc,tcPDPtr parent, - const tPDVector & children) const { +int FOCUSDptoKmPipPip::modeNumber(bool & cc,tcPDPtr parent, + const tPDVector & children) const { int id0(parent->id()); // incoming particle must be D+ if(abs(id0)!=ParticleID::Dplus) return -1; cc = id0<0; // must be three decay products if(children.size()!=3) return -1; - unsigned int npip(0),npim(0),nkm(0),nk0(0),npi0(0); - int id; + unsigned int npip(0),npim(0),nk(0); for(tPDPtr child : children) { - id=child->id(); + long id= child->id(); if(id ==ParticleID::piplus) ++npip; - else if(id ==ParticleID::pi0) ++npi0; else if(id ==ParticleID::piminus) ++npim; - else if(abs(id)==ParticleID::K0) ++nk0; - else if(id ==ParticleID::K_L0) ++nk0; - else if(id ==ParticleID::K_S0) ++nk0; - else if(abs(id)==ParticleID::Kplus) ++nkm; + else if(abs(id)==ParticleID::Kplus) ++nk; } - if((id0==ParticleID::Dplus &&(nkm==1&&npip==2))|| - (id0==ParticleID::Dminus&&(nkm==1&&npim==2))) return 0; + if((id0==ParticleID::Dplus &&(nk==1&&npip==2))|| + (id0==ParticleID::Dminus&&(nk==1&&npim==2))) return 0; else return -1; } -void DtoKPiPiFOCUS:: -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); + +Complex FOCUSDptoKmPipPip::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 & 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); - } - Complex amp(0.); - for(unsigned int ipi=1;ipi<3;++ipi) { - unsigned int iother = ipi==1 ? 2 : 1; - Lorentz5Momentum pres=momenta[0]+momenta[ipi]; - pres.rescaleMass(); - // momentum for the D Blatt-Weisskopf factor - Energy pD = Kinematics::pstarTwoBodyDecay(part.mass(),pres.mass(),momenta[iother].mass()); - // momentum for the resonance Blatt-Weisskopf factor - Energy pAB = Kinematics::pstarTwoBodyDecay(pres.mass(),momenta[0].mass(),momenta[ipi].mass()); - // loop over the resonances - double rD = rD0_*pD, rR = rRes_*pAB; - // most resonaces vectors so calculate that now - double Fd = 1./sqrt(1.+sqr(rD)); - double Fr = 1./sqrt(1.+sqr(rR)); - unsigned int ipow=3; - for(unsigned int ires=0;ires<4;++ires) { - if(ichan>=0 && ichan+1!=int(ires*2+ipi)) continue; - Energy pR = Kinematics::pstarTwoBodyDecay(mRes_[ires],momenta[0].mass(),momenta[ipi].mass()); - // last resonance is the K_2 - if(ires==3) { - Fd = 1./sqrt(9.+3.*sqr(rD)+Math::powi(rD,4)); - Fr = 1./sqrt(9.+3.*sqr(rR)+Math::powi(rR,4)); - ipow = 5; - } - double rB = rRes_*pR; - double Fb = ires<3 ? 1./sqrt(1.+sqr(rB)) : 1./sqrt(9.+3.*sqr(rB)+Math::powi(rB,4)); - Energy gamma = wRes_[ires]*mRes_[ires]/pres.mass()*sqr(Fr/Fb)*Math::powi(pAB/pR,ipow); - complex BW = Fd*Fr/(sqr(mRes_[ires])-pres.mass2()-Complex(0.,1.)*gamma*mRes_[ires]); - // // and the decay momenta - // - // Energy pR = Kinematics::pstarTwoBodyDecay(mres,mA,mB); - // double Fd(1.),Fr(1.),s(1.); - // switch(ispin) { - // case 0: - // // default values of parameters are correct - // break; - // case 1: - // Fr = sqrt((1.+sqr(_rres*pR ))/(1.+sqr(_rres*pAB ))); - // Fd = - // s = ((mAC-mBC)*(mAC+mBC)+(mD-mC)*(mD+mC)*(mB-mA)*(mB+mA)/sqr(mres))/GeV2; - // break; - // case 2: - // Fr = sqrt((9.+3.*sqr(_rres*pR )+Math::powi(_rres*pR ,4))/ - // (9.+3.*sqr(_rres*pAB )+Math::powi(_rres*pAB ,4))); - // Fd = - // s = sqr(((mBC-mAC)*(mBC+mAC)+(mD-mC)*(mD+mC)*(mA-mB)*(mA+mB)/sqr(mres))/GeV2) - // -(mAB*mAB-2.*mD*mD-2.*mC*mC+sqr((mD-mC)*(mD+mC))/sqr(mres))* - // (mAB*mAB-2.*mA*mA-2.*mB*mB+sqr((mA-mB)*(mA+mB))/sqr(mres))/3./GeV2/GeV2; - // break; - // default: - // throw Exception() << "D0toKPiPiCLEO::amplitude spin is too high ispin = " - // << ispin << Exception::runerror; - // } - // } - } - if(ichan>=0 && ichan!=int(7+ipi)) continue; - // finally the scalar piece - static const Energy2 sc=2.*GeV2; - ublas::vector pVector(2); - double fact = 1.-pres.mass2()/KHalf_->poles()[0]; - Complex f1=exp(Complex(0.,1.)*gamma_[0]/180.*Constants::pi)*fact; - double shat = (pres.mass2()-sc)/GeV2; - pVector[0]=0.; - for(unsigned int ix=0;ix p3Vector(1); - p3Vector[0]=0.; - Complex f3=exp(Complex(0.,1.)*gamma_[2]/180.*Constants::pi); - for(unsigned int ix=0;ixpoles()[0]); - pVector[1] += beta_*phase*g2_/KHalf_->poles()[0]; - amp += KHalf_->amplitudes(pres.mass2(), pVector,true)[0]; - amp += KThreeHalf_->amplitudes(pres.mass2(),p3Vector,true)[0]; - } - // now compute the matrix element - (*ME())(0,0,0,0)=amp; - return norm(amp); -} - -void DtoKPiPiFOCUS::dataBaseOutput(ofstream & output, bool header) const { - if(header) output << "update decayers set parameters=\""; - // parameters for the DecayIntegrator base class - DecayIntegrator::dataBaseOutput(output,false); - // // parameters - // output << "newdef " << name() << ":KmPipPipNonResonantMagnitude " - // << _a1NR << "\n"; - // output << "newdef " << name() << ":KmPipPipNonResonantPhase " - // << _phi1NR << "\n"; - // output << "newdef " << name() << ":KmPipPipK892Magnitude " - // << _a1K892 << "\n"; - // output << "newdef " << name() << ":KmPipPipK892Phase " - // << _phi1K892 << "\n"; - // output << "newdef " << name() << ":KmPipPipK1430Magnitude " - // << _a1K1430 << "\n"; - // output << "newdef " << name() << ":KmPipPipK1430Phase " - // << _phi1K1430 << "\n"; - // output << "newdef " << name() << ":KmPipPipK1680Magnitude " - // << _a1K1680 << "\n"; - // output << "newdef " << name() << ":KmPipPipK1680Phase " - // << _phi1K1680 << "\n"; - // output << "newdef " << name() << ":KmPipPi0NonResonantMagnitude " - // << _a2NR << "\n"; - // output << "newdef " << name() << ":KmPipPi0NonResonantPhase " - // << _phi2NR << "\n"; - // output << "newdef " << name() << ":KmPipPi0K8920Magnitude " - // << _a2K8920 << "\n"; - // output << "newdef " << name() << ":KmPipPi0K8920Phase " - // << _phi2K8920 << "\n"; - // output << "newdef " << name() << ":KmPipPi0K892mMagnitude " - // << _a2K892m << "\n"; - // output << "newdef " << name() << ":KmPipPi0K892mPhase " - // << _phi2K892m << "\n"; - // output << "newdef " << name() << ":KmPipPi0RhoMagnitude " - // << _a2rho << "\n"; - // output << "newdef " << name() << ":KmPipPi0RhoPhase " - // << _phi2rho << "\n"; - // output << "newdef " << name() << ":K0PipPimNonResonantMagnitude " - // << _a3NR << "\n"; - // output << "newdef " << name() << ":K0PipPimNonResonantPhase " - // << _phi3NR << "\n"; - // output << "newdef " << name() << ":K0PipPimK892Magnitude " - // << _a3K892 << "\n"; - // output << "newdef " << name() << ":K0PipPimK892Phase " - // << _phi3K892 << "\n"; - // output << "newdef " << name() << ":K0PipPimRhoMagnitude " - // << _a3rho << "\n"; - // output << "newdef " << name() << ":K0PipPimRhoPhase " - // << _phi3rho << "\n"; - // output << "newdef " << name() << ":LocalParameters " << _localparameters << "\n"; - // output << "newdef " << name() << ":K8920Mass " << _mK8920/GeV << "\n"; - // output << "newdef " << name() << ":K8920Width " << _wK8920/GeV << "\n"; - // output << "newdef " << name() << ":K892MinusMass " << _mK892m/GeV << "\n"; - // output << "newdef " << name() << ":K892MinusWidth " << _wK892m/GeV << "\n"; - // output << "newdef " << name() << ":K1680Mass " << _mK1680/GeV << "\n"; - // output << "newdef " << name() << ":K1680Width " << _wK1680/GeV << "\n"; - // output << "newdef " << name() << ":K1430Mass " << _mK1430/GeV << "\n"; - // output << "newdef " << name() << ":K1430Width " << _wK1430/GeV << "\n"; - // output << "newdef " << name() << ":Rho0Mass " << _mrho0 /GeV << "\n"; - // output << "newdef " << name() << ":Rho0Width " << _wrho0 /GeV << "\n"; - // output << "newdef " << name() << ":RhoPlusMass " << _mrhop /GeV << "\n"; - // output << "newdef " << name() << ":RhoPlusWidth " << _wrhop /GeV << "\n"; - // for(unsigned int ix=0;ix<_maxwgt.size();++ix) { - // output << "insert " << name() << ":MaximumWeights " - // << ix << " " << _maxwgt[ix] << "\n"; - // } - // for(unsigned int ix=0;ix<_weights.size();++ix) { - // output << "insert " << name() << ":Weights " - // << ix << " " << _weights[ix] << "\n"; - // } - if(header) { - output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; - } -} - -void DtoKPiPiFOCUS::doinit() { - DecayIntegrator::doinit(); - // set up the phase space - // intermediate resonances - tPDPtr resonances[5]={getParticleData(-313 ), - getParticleData(-100313), - getParticleData(-30313 ), - getParticleData(-315 ), - getParticleData(-10311 )}; - // D+ -> K-pi+pi+ - tPDPtr in = getParticleData(ParticleID::Dplus); - tPDVector out= {getParticleData(ParticleID::Kminus), - getParticleData(ParticleID::piplus), - getParticleData(ParticleID::piplus)}; - PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,maxWgt_)); - vector channels; - for(tPDPtr res : resonances) { - channels.push_back((PhaseSpaceChannel(mode),0,res,0,3,1,1,1,2)); - channels.push_back((PhaseSpaceChannel(mode),0,res,0,2,1,1,1,3)); - } - // add the mode - for(unsigned int ix=0;ixaddChannel(channels[ix]); - } - addMode(mode); -} - -void DtoKPiPiFOCUS::doinitrun() { - DecayIntegrator::doinitrun(); -} diff --git a/Decay/ScalarMeson/FOCUSDptoKmPipPip.h b/Decay/ScalarMeson/FOCUSDptoKmPipPip.h --- a/Decay/ScalarMeson/FOCUSDptoKmPipPip.h +++ b/Decay/ScalarMeson/FOCUSDptoKmPipPip.h @@ -1,246 +1,176 @@ // -*- C++ -*- -#ifndef Herwig_DtoKPiPiFOCUS_H -#define Herwig_DtoKPiPiFOCUS_H +#ifndef Herwig_FOCUSDptoKmPipPip_H +#define Herwig_FOCUSDptoKmPipPip_H // -// This is the declaration of the DtoKPiPiFOCUS class. +// This is the declaration of the FOCUSDptoKmPipPip class. // -#include "Herwig/Decay/DecayIntegrator.h" #include "Herwig/Decay/FormFactors/KMatrix.h" +#include "WeakDalitzDecay.h" namespace Herwig { using namespace ThePEG; /** - * Here is the documentation of the DtoKPiPiFOCUS class. + * The FOCUSDptoKmPipPip class implements the fits of FOCUS, + * Phys.Lett. B653 (2007) 1-11 for the decay \f$D^+\toK^-\pi^+\pi^-\f$ * - * @see \ref DtoKPiPiFOCUSInterfaces "The interfaces" - * defined for DtoKPiPiFOCUS. + * @see \ref FOCUSDptoKmPipPipInterfaces "The interfaces" + * defined for FOCUSDptoKmPipPip. */ -class DtoKPiPiFOCUS: public DecayIntegrator { +class FOCUSDptoKmPipPip: public WeakDalitzDecay { public: /** * The default constructor. */ - DtoKPiPiFOCUS(); + FOCUSDptoKmPipPip(); /** * Which of the possible decays is required * @param cc Is this mode the charge conjugate * @param parent The decaying particle * @param children The decay products */ virtual int modeNumber(bool & cc, tcPDPtr parent, const tPDVector & children) const; - /** - * Return the matrix element squared for a given mode and phase-space channel. - * @param ichan The channel we are calculating the matrix element for. - * @param part The decaying Particle. - * @param outgoing The particles produced in the decay - * @param momenta The momenta of the particles produced in the decay - * @param meopt Option for the calculation of the matrix element - * @return The matrix element squared for the phase-space configuration. - */ - double me2(const int ichan,const Particle & part, - const tPDVector & outgoing, - const vector & momenta, - MEOption meopt) const; - - /** - * Construct the SpinInfos for the particles produced in the decay - */ - virtual void constructSpinInfo(const Particle & part, - ParticleVector outgoing) const; - - /** - * Output the setup information for the particle database - * @param os The stream to output the information to - * @param header Whether or not to output the information for MySQL - */ - virtual void dataBaseOutput(ofstream & os,bool header) const; - public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); - + +protected: + + /** @name Clone Methods. */ + //@{ + /** + * Make a simple clone of this object. + * @return a pointer to the new object. + */ + virtual IBPtr clone() const; + + /** Make a clone of this object, possibly modifying the cloned object + * to make it sane. + * @return a pointer to the new object. + */ + virtual IBPtr fullclone() const; + //@} + protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving 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: - /** @name Clone Methods. */ - //@{ /** - * Make a simple clone of this object. - * @return a pointer to the new object. + * Calculate the amplitude */ - 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; - //@} + 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. */ - DtoKPiPiFOCUS & operator=(const DtoKPiPiFOCUS &); + FOCUSDptoKmPipPip & operator=(const FOCUSDptoKmPipPip &) = delete; private: /** - * Parameters for the Blatt-Weisskopf form-factors - */ - //@{ - /** - * Radial size for the \f$D^0\f$ - */ - InvEnergy rD0_; - - /** - * Radial size for the light resonances - */ - InvEnergy rRes_; - //@} - - /** * The K-matrices for the s-wave component */ //@{ /** * The \f$I=\frac12\f$ component */ KMatrixPtr KHalf_; /** * The \f$I=\frac32\f$ component */ KMatrixPtr KThreeHalf_; //@} /** * Parameters for the f$pf$-vector */ //@{ /** * Pole couplings */ Energy g1_,g2_; /** * \f$\beta\f$ */ Energy beta_; /** * \f$\theta\f$ */ double theta_; /** * \f$\gamma\f$ */ vector gamma_; /** * \f$c_{1i}\f$ */ vector c1_; /** * \f$c_{2i}\f$ */ vector c2_; /** * \f$c_{3i}\f$ */ vector c3_; //@} - - /** - * Parameters for the phase-space integration - */ - //@{ - /** - * Maximum weights for the various modes - */ - double maxWgt_; - - /** - * Weights for the different integration channels - */ - vector weights_; - //@} - - /** - * Masses and widths of the resonances - */ - //@{ - /** - * Masses - */ - vector mRes_; - - /** - * Widths - */ - vector wRes_; - //@} - - /** - * Spin density matrix - */ - mutable RhoDMatrix rho_; - }; } -#endif /* Herwig_DtoKPiPiFOCUS_H */ +#endif /* Herwig_FOCUSDptoKmPipPip_H */