Page MenuHomeHEPForge
Files F19251274

No OneTemporary
Actions

Size
53 KB
Referenced Files
None
Subscribers
None

View Options

diff --git a/Decay/ScalarMeson/ScalarMesonFactorizedDecayer.cc b/Decay/ScalarMeson/ScalarMesonFactorizedDecayer.cc
--- a/Decay/ScalarMeson/ScalarMesonFactorizedDecayer.cc
+++ b/Decay/ScalarMeson/ScalarMesonFactorizedDecayer.cc
@@ -1,720 +1,686 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ScalarMesonFactorizedDecayer class.
//
#include "ScalarMesonFactorizedDecayer.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/epsilon.h"
#include "ThePEG/Helicity/WaveFunction/TensorWaveFunction.h"
-namespace Herwig {
-using namespace ThePEG;
+using namespace Herwig;
using namespace ThePEG::Helicity;
-using Helicity::ScalarWaveFunction;
-using Helicity::VectorWaveFunction;
-using Helicity::TensorWaveFunction;
-using Helicity::incoming;
-using Helicity::outgoing;
+using Herwig::Helicity::ScalarWaveFunction;
+using Herwig::Helicity::VectorWaveFunction;
+using Herwig::Helicity::TensorWaveFunction;
+using Herwig::Helicity::incoming;
+using Herwig::Helicity::outgoing;
inline void ScalarMesonFactorizedDecayer::doinit() throw(InitException) {
DecayIntegrator::doinit();
unsigned int ix,iy,iz,icurr,iform;
// get the CKM matrix (unsquared for interference)
Complex ckmmat[3][3];
vector< vector<Complex > > CKM(_theCKM->getUnsquaredMatrix(SM().families()));
for(ix=0;ix<3;++ix){for(iy=0;iy<3;++iy){ckmmat[ix][iy]=CKM[ix][iy];}}
int id0,id1,Wcharge,iq,ia,jspin,spect,inq,outq;
// make sure the currents and form factors got initialised
for(ix=0;ix<_current.size();++ix){_current[ix]->init();}
for(ix=0;ix<_form.size();++ix){_form[ix]->init();}
// find all the possible modes
vector<unsigned int> tformmap[2],tcurrmap[2];
vector<int> inquark,outquark,currq,curra;
vector<PDVector> particles;
PDVector extpart,ptemp;
Energy min,minb;
// loop over the modes in the form factors and currents
- for(iform=0;iform<_form.size();++iform)
- {
- for(ix=0;ix<_form[iform]->numberOfFactors();++ix)
- {
- // particles from the form-factor
+ for(iform=0;iform<_form.size();++iform) {
+ for(ix=0;ix<_form[iform]->numberOfFactors();++ix) {
+ // particles from the form-factor
+ extpart.resize(2);
+ _form[iform]->particleID(ix,id0,id1);
+ _form[iform]->formFactorInfo(ix,jspin,spect,inq,outq);
+ // particles from the form factor
+ extpart[0]=getParticleData(id0);
+ extpart[1]=getParticleData(id1);
+ // charge of the decay products
+ Wcharge =extpart[0]->iCharge()-extpart[1]->iCharge();
+ // max mass for the particles in the current
+ min = extpart[0]->massMax()-extpart[1]->massMin();
+ for(icurr=0;icurr<_current.size();++icurr) {
+ for(iy=0;iy<_current[icurr]->numberOfModes();++iy) {
extpart.resize(2);
- _form[iform]->particleID(ix,id0,id1);
- _form[iform]->formFactorInfo(ix,jspin,spect,inq,outq);
- // particles from the form factor
- extpart[0]=getParticleData(id0);
- extpart[1]=getParticleData(id1);
- // charge of the decay products
- Wcharge =extpart[0]->iCharge()-extpart[1]->iCharge();
- // max mass for the particles in the current
- min = extpart[0]->massMax()-extpart[1]->massMin();
- for(icurr=0;icurr<_current.size();++icurr)
- {
- for(iy=0;iy<_current[icurr]->numberOfModes();++iy)
- {
- extpart.resize(2);
- // get the particles from the current
- _current[icurr]->decayModeInfo(iy,iq,ia);
- ptemp=_current[icurr]->particles(Wcharge,iy,iq,ia);
- minb=0.*MeV;
- for(iz=0;iz<ptemp.size();++iz)
- {
- extpart.push_back(ptemp[iz]);
- minb+=ptemp[iz]->massMin();
- }
- // add this mode to the list
- if(extpart.size()>2&&minb<min&&
- (Wcharge!=0||(Wcharge==0&&((inq>0&&inq%2!=iq%2)||
- (inq<0&&abs(inq)%2!=abs(ia)%2)))))
- {
- tformmap[0].push_back(iform);tformmap[1].push_back(ix);
- tcurrmap[0].push_back(icurr);tcurrmap[1].push_back(iy);
- particles.push_back(extpart);
- inquark.push_back(inq);outquark.push_back(outq);
- currq.push_back(iq);curra.push_back(ia);
- }
- // if the meson in the current is neutral try the CC mode
- if(Wcharge==0&&iq!=-ia&&((inq>0&&inq%2!=iq%2)||
- (inq<0&&abs(inq)%2!=abs(ia)%2)))
- {
- extpart.resize(2);
- // get the particles from the current
- ptemp=_current[icurr]->particles(Wcharge,iy,-ia,-iq);
- minb=0.*MeV;
- for(iz=0;iz<ptemp.size();++iz)
- {
- extpart.push_back(ptemp[iz]);
- minb+=ptemp[iz]->massMin();
- }
- if(extpart.size()>2&&minb<min)
- {
- tformmap[0].push_back(iform);tformmap[1].push_back(ix);
- tcurrmap[0].push_back(icurr);tcurrmap[1].push_back(iy);
- particles.push_back(extpart);
- inquark.push_back(inq);outquark.push_back(outq);
- currq.push_back(-ia);curra.push_back(-iq);
- }
- }
- }
+ // get the particles from the current
+ _current[icurr]->decayModeInfo(iy,iq,ia);
+ ptemp=_current[icurr]->particles(Wcharge,iy,iq,ia);
+ minb=0.*MeV;
+ for(iz=0;iz<ptemp.size();++iz) {
+ extpart.push_back(ptemp[iz]);
+ minb+=ptemp[iz]->massMin();
+ }
+ // add this mode to the list
+ if(extpart.size()>2&&minb<min&&
+ (Wcharge!=0||(Wcharge==0&&((inq>0&&inq%2!=iq%2)||
+ (inq<0&&abs(inq)%2!=abs(ia)%2))))) {
+ tformmap[0].push_back(iform);tformmap[1].push_back(ix);
+ tcurrmap[0].push_back(icurr);tcurrmap[1].push_back(iy);
+ particles.push_back(extpart);
+ inquark.push_back(inq);outquark.push_back(outq);
+ currq.push_back(iq);curra.push_back(ia);
+ }
+ // if the meson in the current is neutral try the CC mode
+ if(Wcharge==0&&iq!=-ia&&((inq>0&&inq%2!=iq%2)||
+ (inq<0&&abs(inq)%2!=abs(ia)%2))) {
+ extpart.resize(2);
+ // get the particles from the current
+ ptemp=_current[icurr]->particles(Wcharge,iy,-ia,-iq);
+ minb=0.*MeV;
+ for(iz=0;iz<ptemp.size();++iz) {
+ extpart.push_back(ptemp[iz]);
+ minb+=ptemp[iz]->massMin();
}
+ if(extpart.size()>2&&minb<min) {
+ tformmap[0].push_back(iform);tformmap[1].push_back(ix);
+ tcurrmap[0].push_back(icurr);tcurrmap[1].push_back(iy);
+ particles.push_back(extpart);
+ inquark.push_back(inq);outquark.push_back(outq);
+ currq.push_back(-ia);curra.push_back(-iq);
+ }
+ }
}
+ }
}
+ }
// loop over the modes and find the dupliciates
vector<bool> modecc; vector<unsigned int> modeloc,tformpart,ttform[2],ttcurr[2];
vector<Complex> tCKM; Complex ckm;
DecayPhaseSpaceModePtr mode;
DecayPhaseSpaceChannelPtr channel;
bool done;
int id,idbar;
double maxweight;
vector<double> channelwgts;
unsigned int isize;double ort(sqrt(0.5));
vector<double>::iterator start,end;
- for(ix=0;ix<particles.size();++ix)
- {
- while (true)
- {
- if(particles[ix].size()==0){break;}
- findModes(ix,particles,modeloc,modecc);
- // if more than three particles only allow one diagram
- if(particles[ix].size()>3&&modeloc.size()!=0){break;}
- // create the mode and set the particles as for the first instance
- mode=new_ptr(DecayPhaseSpaceMode(particles[ix],this));
- channel = new_ptr(DecayPhaseSpaceChannel(mode));
- channel->addIntermediate(particles[ix][0],0,0.0,1,-1);
- min = particles[ix][0]->massMax()-particles[ix][1]->massMin();
- Wcharge = particles[ix][0]->iCharge()-particles[ix][1]->iCharge();
- done=_current[tcurrmap[0][ix]]->createMode(Wcharge,tcurrmap[1][ix],
- mode,2,1,channel,min);
- if(!done){throw InitException() << "Failed to construct mode in "
- << "ScalarMesonFactorizedDecayer::doinit()."
- << Exception::abortnow;}
- // set the parameters for the additional modes
- ttform[0].resize(0);ttform[1].resize(0);
- ttcurr[0].resize(0);ttcurr[1].resize(0);
- ttform[0].push_back(tformmap[0][ix]);ttform[1].push_back(tformmap[1][ix]);
- ttcurr[0].push_back(tcurrmap[0][ix]);ttcurr[1].push_back(tcurrmap[1][ix]);
- tformpart.resize(0);tformpart.push_back(0);
- id=particles[ix][1]->id();
- if(particles[ix][1]->CC()){idbar=particles[ix][1]->CC()->id();}
- else{idbar=id;}
- for(iy=0;iy<modeloc.size();++iy)
- {
- ttform[0].push_back(tformmap[0][modeloc[iy]]);
- ttform[1].push_back(tformmap[1][modeloc[iy]]);
- ttcurr[0].push_back(tcurrmap[0][modeloc[iy]]);
- ttcurr[1].push_back(tcurrmap[1][modeloc[iy]]);
- iz=1;
- do
- {
- if(( modecc[iy]&&particles[modeloc[iy]][iz]->id()==idbar)||
- (!modecc[iy]&&particles[modeloc[iy]][iz]->id()==id))
- {tformpart.push_back(iz-1);}
- ++iz;
- }
- while(tformpart.size()!=iy+2&&iz<3);
- }
- // calculate ckm factors
- tCKM.resize(0);
- for(iy=0;iy<ttcurr[0].size();++iy)
- {
- // get the quarks involved in the process
- if(iy==0)
- {
- iq=currq[ix];ia=curra[ix];
- inq=inquark[ix];outq=outquark[ix];
- }
- else
- {
- if(!modecc[iy-1])
- {
- iq=currq[modeloc[iy-1]];ia=curra[modeloc[iy-1]];
- inq=inquark[modeloc[iy-1]];outq=outquark[modeloc[iy-1]];
- }
- else
- {
- ia=-currq[modeloc[iy-1]];iq=-curra[modeloc[iy-1]];
- inq=-inquark[modeloc[iy-1]];outq=-outquark[modeloc[iy-1]];
- }
- }
- _form[ttform[0][iy]]->particleID(ttform[1][iy],id0,id1);
- Wcharge = getParticleData(id0)->iCharge()-getParticleData(id1)->iCharge();
- ckm=1.;
- if(Wcharge!=0)
- {
- ckm=1.;
- if(abs(iq)%2==0){ckm *= conj(ckmmat[abs(iq)/2-1][(abs(ia)-1)/2]);}
- else{ckm *= conj(ckmmat[abs(ia)/2-1][(abs(iq)-1)/2]);}
- if(abs(inq)%2==0){ckm *= ckmmat[abs(inq)/2-1][(abs(outq)-1)/2];}
- else{ckm *= ckmmat[abs(outq)/2-1][(abs(inq)-1)/2];}
- if(abs(inq)==5){ckm*=_a1b;}
- else{ckm*=_a1c;}
- }
- else
- {
- ckm=1.;
- if(inq>0)
- {
- if(abs(inq)%2==0){ckm *= ckmmat[abs(inq)/2-1][(abs(iq)-1)/2];}
- else{ckm *= ckmmat[abs(iq)/2-1][(abs(inq)-1)/2];}
- if(abs(outq)%2==0)
- {ckm *= conj(ckmmat[abs(outq)/2-1][(abs(ia)-1)/2]);}
- else{ckm *= conj(ckmmat[abs(ia)/2-1][(abs(outq)-1)/2]);}
- }
- else
- {
- if(abs(inq)%2==0){ckm *= ckmmat[abs(inq)/2-1][(abs(ia)-1)/2];}
- else{ckm *= ckmmat[abs(ia)/2-1][(abs(inq)-1)/2];}
- if(abs(outq)%2==0)
- {ckm *= conj(ckmmat[abs(outq)/2-1][(abs(iq)-1)/2]);}
- else{ckm *= conj(ckmmat[abs(iq)/2-1][(abs(outq)-1)/2]);}
- }
- if(abs(inq)==5){ckm*=_a2b;}
- else{ckm*=_a2c;}
- }
- if((abs(inq)%2==0&&inq<0)||(abs(inq)%2!=0&&inq>0)){ckm=conj(ckm);}
- tCKM.push_back(ckm);
- }
- // special if the particles are idential add additional modes and
- // identical particle factors
- if(particles[ix][1]->id()==particles[ix][2]->id()&&particles[ix].size()==3)
- {
- isize=ttcurr[0].size();
- for(iy=0;iy<isize;++iy)
- {
- ttcurr[0].push_back(ttcurr[0][iy]);ttcurr[1].push_back(ttcurr[1][iy]);
- ttform[0].push_back(ttform[0][iy]);ttform[1].push_back(ttform[1][iy]);
- if(tformpart[ix]==0){tformpart.push_back(1);}
- else{tformpart.push_back(0);}
- tCKM[iy]*=ort;tCKM.push_back(tCKM[iy]);
- }
- }
- // add the parameters for the mode to the list
- _currentmapA.push_back(ttcurr[0]);_currentmapB.push_back(ttcurr[1]);
- _formmapA.push_back(ttform[0]);_formmapB.push_back(ttform[1]);
- _formpart.push_back(tformpart);
- _CKMfact.push_back(tCKM);
- // add the mode to the list
- if(_wgtmax.size()>numberModes()){maxweight=_wgtmax[numberModes()];}
- else{maxweight=0.;}
- // the weights for the channel
- if(_wgtloc.size()>numberModes()&&
- _wgtloc[numberModes()]+mode->numberChannels()<=_weights.size())
- {
- start=_weights.begin()+_wgtloc[numberModes()];
- end = start+mode->numberChannels();
- channelwgts=vector<double>(start,end);
- }
- else
- {channelwgts.resize(mode->numberChannels(),1./(mode->numberChannels()));}
- // don't need channels for two body decays
- if(particles[ix].size()==3)
- {
- channelwgts.resize(0);
- mode=new_ptr(DecayPhaseSpaceMode(particles[ix],this));
- }
- addMode(mode,maxweight,channelwgts);
- // resize the duplicate modes to remove them
- for(iy=0;iy<modeloc.size();++iy){particles[modeloc[iy]]=PDVector(0);}
- break;
+ for(ix=0;ix<particles.size();++ix) {
+ while (true) {
+ if(particles[ix].size()==0){break;}
+ findModes(ix,particles,modeloc,modecc);
+ // if more than three particles only allow one diagram
+ if(particles[ix].size()>3&&modeloc.size()!=0){break;}
+ // create the mode and set the particles as for the first instance
+ mode=new_ptr(DecayPhaseSpaceMode(particles[ix],this));
+ channel = new_ptr(DecayPhaseSpaceChannel(mode));
+ channel->addIntermediate(particles[ix][0],0,0.0,1,-1);
+ min = particles[ix][0]->massMax()-particles[ix][1]->massMin();
+ Wcharge = particles[ix][0]->iCharge()-particles[ix][1]->iCharge();
+ done=_current[tcurrmap[0][ix]]->createMode(Wcharge,tcurrmap[1][ix],
+ mode,2,1,channel,min);
+ if(!done) throw InitException() << "Failed to construct mode in "
+ << "ScalarMesonFactorizedDecayer::doinit()."
+ << Exception::abortnow;
+ // set the parameters for the additional modes
+ ttform[0].resize(0);ttform[1].resize(0);
+ ttcurr[0].resize(0);ttcurr[1].resize(0);
+ ttform[0].push_back(tformmap[0][ix]);ttform[1].push_back(tformmap[1][ix]);
+ ttcurr[0].push_back(tcurrmap[0][ix]);ttcurr[1].push_back(tcurrmap[1][ix]);
+ tformpart.resize(0);tformpart.push_back(0);
+ id=particles[ix][1]->id();
+ if(particles[ix][1]->CC()){idbar=particles[ix][1]->CC()->id();}
+ else{idbar=id;}
+ for(iy=0;iy<modeloc.size();++iy) {
+ ttform[0].push_back(tformmap[0][modeloc[iy]]);
+ ttform[1].push_back(tformmap[1][modeloc[iy]]);
+ ttcurr[0].push_back(tcurrmap[0][modeloc[iy]]);
+ ttcurr[1].push_back(tcurrmap[1][modeloc[iy]]);
+ iz=1;
+ do {
+ if(( modecc[iy]&&particles[modeloc[iy]][iz]->id()==idbar)||
+ (!modecc[iy]&&particles[modeloc[iy]][iz]->id()==id))
+ tformpart.push_back(iz-1);
+ ++iz;
}
+ while(tformpart.size()!=iy+2&&iz<3);
+ }
+ // calculate ckm factors
+ tCKM.resize(0);
+ for(iy=0;iy<ttcurr[0].size();++iy) {
+ // get the quarks involved in the process
+ if(iy==0) {
+ iq=currq[ix];ia=curra[ix];
+ inq=inquark[ix];outq=outquark[ix];
+ }
+ else {
+ if(!modecc[iy-1]) {
+ iq=currq[modeloc[iy-1]];ia=curra[modeloc[iy-1]];
+ inq=inquark[modeloc[iy-1]];outq=outquark[modeloc[iy-1]];
+ }
+ else {
+ ia=-currq[modeloc[iy-1]];iq=-curra[modeloc[iy-1]];
+ inq=-inquark[modeloc[iy-1]];outq=-outquark[modeloc[iy-1]];
+ }
+ }
+ _form[ttform[0][iy]]->particleID(ttform[1][iy],id0,id1);
+ Wcharge = getParticleData(id0)->iCharge()-getParticleData(id1)->iCharge();
+ ckm=1.;
+ if(Wcharge!=0) {
+ ckm=1.;
+ if(abs(iq)%2==0) ckm *= conj(ckmmat[abs(iq)/2-1][(abs(ia)-1)/2]);
+ else ckm *= conj(ckmmat[abs(ia)/2-1][(abs(iq)-1)/2]);
+ if(abs(inq)%2==0) ckm *= ckmmat[abs(inq)/2-1][(abs(outq)-1)/2];
+ else ckm *= ckmmat[abs(outq)/2-1][(abs(inq)-1)/2];
+ if(abs(inq)==5) ckm*=_a1b;
+ else ckm*=_a1c;
+ }
+ else {
+ ckm=1.;
+ if(inq>0) {
+ if(abs(inq)%2==0) ckm *= ckmmat[abs(inq)/2-1][(abs(iq)-1)/2];
+ else ckm *= ckmmat[abs(iq)/2-1][(abs(inq)-1)/2];
+ if(abs(outq)%2==0) ckm *= conj(ckmmat[abs(outq)/2-1][(abs(ia)-1)/2]);
+ else ckm *= conj(ckmmat[abs(ia)/2-1][(abs(outq)-1)/2]);
+ }
+ else {
+ if(abs(inq)%2==0) ckm *= ckmmat[abs(inq)/2-1][(abs(ia)-1)/2];
+ else ckm *= ckmmat[abs(ia)/2-1][(abs(inq)-1)/2];
+ if(abs(outq)%2==0) ckm *= conj(ckmmat[abs(outq)/2-1][(abs(iq)-1)/2]);
+ else ckm *= conj(ckmmat[abs(iq)/2-1][(abs(outq)-1)/2]);
+ }
+ if(abs(inq)==5) ckm*=_a2b;
+ else ckm*=_a2c;
+ }
+ if((abs(inq)%2==0&&inq<0)||(abs(inq)%2!=0&&inq>0)){ckm=conj(ckm);}
+ tCKM.push_back(ckm);
+ }
+ // special if the particles are idential add additional modes and
+ // identical particle factors
+ if(particles[ix][1]->id()==particles[ix][2]->id()&&particles[ix].size()==3) {
+ isize=ttcurr[0].size();
+ for(iy=0;iy<isize;++iy) {
+ ttcurr[0].push_back(ttcurr[0][iy]);ttcurr[1].push_back(ttcurr[1][iy]);
+ ttform[0].push_back(ttform[0][iy]);ttform[1].push_back(ttform[1][iy]);
+ if(tformpart[ix]==0){tformpart.push_back(1);}
+ else{tformpart.push_back(0);}
+ tCKM[iy]*=ort;tCKM.push_back(tCKM[iy]);
+ }
+ }
+ // add the parameters for the mode to the list
+ _currentmapA.push_back(ttcurr[0]);_currentmapB.push_back(ttcurr[1]);
+ _formmapA.push_back(ttform[0]);_formmapB.push_back(ttform[1]);
+ _formpart.push_back(tformpart);
+ _CKMfact.push_back(tCKM);
+ // add the mode to the list
+ if(_wgtmax.size()>numberModes()){maxweight=_wgtmax[numberModes()];}
+ else{maxweight=0.;}
+ // the weights for the channel
+ if(_wgtloc.size()>numberModes()&&
+ _wgtloc[numberModes()]+mode->numberChannels()<=_weights.size()) {
+ start=_weights.begin()+_wgtloc[numberModes()];
+ end = start+mode->numberChannels();
+ channelwgts=vector<double>(start,end);
+ }
+ else {
+ channelwgts.resize(mode->numberChannels(),1./(mode->numberChannels()));
+ }
+ // don't need channels for two body decays
+ if(particles[ix].size()==3) {
+ channelwgts.resize(0);
+ mode=new_ptr(DecayPhaseSpaceMode(particles[ix],this));
+ }
+ addMode(mode,maxweight,channelwgts);
+ // resize the duplicate modes to remove them
+ for(iy=0;iy<modeloc.size();++iy) particles[modeloc[iy]]=PDVector(0);
+ break;
}
+ }
}
bool ScalarMesonFactorizedDecayer::accept(const DecayMode & dm) const {
// N.B. this is a necessary but not sufficient test
bool allowed(false),dummy;
// find the ids of the particles
ParticleMSet::const_iterator pit = dm.products().begin();
ParticleMSet::const_iterator pend = dm.products().end();
vector<int> ids,idcurr;
int id(dm.parent()->id());
for( ; pit!=pend;++pit){ids.push_back((**pit).id());}
// loop over the possible particles in the formfactor
unsigned int ipart(0),iform,icurr,ix;
- do
- {
- idcurr.resize(0);
- for(ix=0;ix<ids.size();++ix){if(ix!=ipart){idcurr.push_back(ids[ix]);}}
- iform=0;
- do
- {
- // check if possible from the form factor
- if(_form[iform]->formFactorNumber(id,ids[ipart],dummy)>=0)
- {
- // check if possible from the current
- icurr=0;
- do{allowed=_current[icurr]->accept(idcurr);++icurr;}
- while(!allowed&&icurr<_current.size());
- }
- ++iform;
+ do {
+ idcurr.resize(0);
+ for(ix=0;ix<ids.size();++ix){if(ix!=ipart){idcurr.push_back(ids[ix]);}}
+ iform=0;
+ do {
+ // check if possible from the form factor
+ if(_form[iform]->formFactorNumber(id,ids[ipart],dummy)>=0) {
+ // check if possible from the current
+ icurr=0;
+ do {
+ allowed=_current[icurr]->accept(idcurr);
+ ++icurr;
}
- while(!allowed&&iform<_form.size());
- ++ipart;
+ while(!allowed&&icurr<_current.size());
+ }
+ ++iform;
}
+ while(!allowed&&iform<_form.size());
+ ++ipart;
+ }
while(!allowed&&ipart<ids.size());
return allowed;
}
int ScalarMesonFactorizedDecayer::modeNumber(bool & cc,const DecayMode & dm) const {
int imode(-1);
// id's of the particles and CC
// of the parent
int id0(dm.parent()->id()),id0bar(id0);
if(dm.parent()->CC()){id0bar=dm.parent()->CC()->id();}
// of the products
vector<int> ids,idbars;
ParticleMSet::const_iterator pit = dm.products().begin();
ParticleMSet::const_iterator pend = dm.products().end();
for( ;pit!=pend;++pit) {
ids.push_back((**pit).id());
- if((**pit).CC()){idbars.push_back((**pit).CC()->id());}
- else{idbars.push_back(ids.back());}
+ if((**pit).CC()) idbars.push_back((**pit).CC()->id());
+ else idbars.push_back(ids.back());
}
// loop over the modes
vector<bool> done(ids.size(),false);
unsigned int nfound,ix,iy,iz;
int idtemp;
bool found;
cc=false;
ix=0;
do {
// particle mode
if(id0==mode(ix)->externalParticles(0)->id()&&
ids.size()+1==mode(ix)->numberofParticles()) {
nfound=0;
for(iy=0;iy<ids.size();++iy){done[iy]=false;}
for(iy=0;iy<ids.size();++iy) {
idtemp=mode(ix)->externalParticles(iy+1)->id();
iz=0;found=false;
do{if(idtemp==ids[iz]&&!done[iz]){done[iz]=true;found=true;}++iz;}
while(iz<ids.size()&&!found);
if(found){++nfound;}
}
if(nfound==ids.size()){cc=false;imode=ix;}
}
// CC mode
if(id0bar==mode(ix)->externalParticles(0)->id()&&
ids.size()+1==mode(ix)->numberofParticles()) {
nfound=0;
for(iy=0;iy<idbars.size();++iy) done[iy]=false;
for(iy=0;iy<idbars.size();++iy) {
idtemp=mode(ix)->externalParticles(iy+1)->id();
iz=0;found=false;
do {
if(idtemp==idbars[iz]&&!done[iz]) {
done[iz]=true;
found=true;
}
++iz;
}
while(iz<idbars.size()&&!found);
if(found) ++nfound;
}
if(nfound==idbars.size()) {
cc=true;
imode=ix;
}
}
++ix;
}
while(imode<0&&ix<numberModes());
if(imode<0){throw DecayIntegratorError() << "Unable to find the mode in "
<< "ScalarMesonFactorizedDecayer::decay()"
<< Exception::abortnow;}
return imode;
}
void ScalarMesonFactorizedDecayer::persistentOutput(PersistentOStream & os) const {
os << _current << _form << _theCKM << ounit(_GF,1/GeV2)
<< _a1b << _a2b << _a1c << _a2c
<< _currentmapA << _currentmapB
<< _formmapA << _formmapB << _formpart << _wgtloc
<< _wgtmax << _weights << _CKMfact ;
}
void ScalarMesonFactorizedDecayer::persistentInput(PersistentIStream & is, int) {
is >> _current >> _form >> _theCKM >> iunit(_GF ,1/GeV2)
>> _a1b >> _a2b >> _a1c >> _a2c
>> _currentmapA >> _currentmapB
>> _formmapA >> _formmapB >> _formpart >> _wgtloc
>> _wgtmax >> _weights >> _CKMfact;
}
ClassDescription<ScalarMesonFactorizedDecayer>
ScalarMesonFactorizedDecayer::initScalarMesonFactorizedDecayer;
// Definition of the static class description member.
void ScalarMesonFactorizedDecayer::Init() {
static ClassDocumentation<ScalarMesonFactorizedDecayer> documentation
("The ScalarMesonFactorizedDecayer class is designed for the weak decay of"
" scalar mesons using the factorization approximation.");
static RefVector<ScalarMesonFactorizedDecayer,WeakDecayCurrent> interfaceCurrents
("Currents",
"A vector of references to the currents",
&ScalarMesonFactorizedDecayer::_current, -1, false, false, true, false, false);
static RefVector<ScalarMesonFactorizedDecayer,ScalarFormFactor> interfaceFormFactors
("FormFactors",
"A vector of references to the form-factors",
&ScalarMesonFactorizedDecayer::_form, -1, false, false, true, false, false);
static Parameter<ScalarMesonFactorizedDecayer,InvEnergy2> interfaceGFermi
("GFermi",
"The Fermi coupling constant",
&ScalarMesonFactorizedDecayer::_GF,
1./GeV2, 1.16639E-5/GeV2, 0./GeV2, 1.0e-4/GeV2,
false, false, false);
static Parameter<ScalarMesonFactorizedDecayer,double> interfacea1Bottom
("a1Bottom",
"The factorization paramter a_1 for decays of bottom baryons",
&ScalarMesonFactorizedDecayer::_a1b, 1.1, -10.0, 10.0,
false, false, true);
static Parameter<ScalarMesonFactorizedDecayer,double> interfacea2Bottom
("a2Bottom",
"The factorization paramter a_2 for decays of bottom baryons",
&ScalarMesonFactorizedDecayer::_a2b, -0.24, -10.0, 10.0,
false, false, true);
static Parameter<ScalarMesonFactorizedDecayer,double> interfacea1Charm
("a1Charm",
"The factorization paramter a_1 for decays of charm baryons",
&ScalarMesonFactorizedDecayer::_a1c, 1.3, -10.0, 10.0,
false, false, true);
static Parameter<ScalarMesonFactorizedDecayer,double> interfacea2Charm
("a2Charm",
"The factorization paramter a_2 for decays of charm baryons",
&ScalarMesonFactorizedDecayer::_a2c, -0.55, -10.0, 10.0,
false, false, true);
static Reference<ScalarMesonFactorizedDecayer,StandardCKM> interfaceCKM
("CKM",
"Reference to the Standard Model object",
&ScalarMesonFactorizedDecayer::_theCKM, false, false, true, false);
static ParVector<ScalarMesonFactorizedDecayer,int> interfaceWeightLocation
("WeightLocation",
"The locations of the weights for a given channel in the vector",
&ScalarMesonFactorizedDecayer::_wgtloc,
0, 0, 0, 0, 10000, false, false, true);
static ParVector<ScalarMesonFactorizedDecayer,double> interfaceWeightMax
("MaximumWeight",
"The maximum weight for a given channel.",
&ScalarMesonFactorizedDecayer::_wgtmax,
0, 0, 0, 0., 100., false, false, true);
static ParVector<ScalarMesonFactorizedDecayer,double> interfaceWeights
("Weights",
"The weights for the integration.",
&ScalarMesonFactorizedDecayer::_weights,
0, 0, 0, 0., 1., false, false, true);
}
double ScalarMesonFactorizedDecayer::me2(bool vertex, const int ichan,
const Particle & part,
- const ParticleVector & decay) const
-{
+ const ParticleVector & decay) const {
// find the mode
unsigned int mode(imode()),ix,iy,chel,fhel;
int id0(part.id()),id1;
Complex ii(0.,1.);
-
// workaround for gcc 3.2.3 bug
// spin info for the decaying particle
//ALB ScalarWaveFunction(const_ptr_cast<tPPtr>(&part),incoming,true,vertex);
tPPtr mytempPart = const_ptr_cast<tPPtr>(&part);
ScalarWaveFunction(mytempPart,incoming,true,vertex);
-
vector<unsigned int> ihel(decay.size());
// get the wavefunctions of the decay products
vector<vector<LorentzPolarizationVector> > vecwave(decay.size());
vector<vector<LorentzTensor<double> > > tenwave(decay.size());
- for(ix=0;ix<decay.size();++ix)
- {
- iy=decay[ix]->dataPtr()->iSpin();
-
- // workaround for gcc 3.2.3 bug
- //ALB if(iy==1){ScalarWaveFunction(decay[ix],outgoing,true,vertex);}
- if(iy==1){PPtr mytemp=decay[ix]; ScalarWaveFunction(mytemp,outgoing,true,vertex);}
- else if(iy==3)
- //ALB {VectorWaveFunction(vecwave[ix],decay[ix],outgoing,true,false,vertex);}
- {
- vector<LorentzPolarizationVector> mytempLPV ;
- VectorWaveFunction(mytempLPV,decay[ix],outgoing,true,false,vertex);
- vecwave[ix]=mytempLPV ;
- }
- else if(iy==5)
- //ALB {TensorWaveFunction(tenwave[ix],decay[ix],outgoing,true,false,vertex);}
- {
- vector<LorentzTensor<double> > mytempLT;
- TensorWaveFunction(mytempLT,decay[ix],outgoing,true,false,vertex);
- tenwave[ix]=mytempLT;
- }
+ for(ix=0;ix<decay.size();++ix) {
+ iy=decay[ix]->dataPtr()->iSpin();
+ // workaround for gcc 3.2.3 bug
+ //ALB if(iy==1){ScalarWaveFunction(decay[ix],outgoing,true,vertex);}
+ if(iy==1){PPtr mytemp=decay[ix]; ScalarWaveFunction(mytemp,outgoing,true,vertex);}
+ else if(iy==3) {
+ //ALB {VectorWaveFunction(vecwave[ix],decay[ix],outgoing,true,false,vertex);}
+ vector<LorentzPolarizationVector> mytempLPV ;
+ VectorWaveFunction(mytempLPV,decay[ix],outgoing,true,false,vertex);
+ vecwave[ix]=mytempLPV ;
}
+ else if(iy==5) {
+ //ALB {TensorWaveFunction(tenwave[ix],decay[ix],outgoing,true,false,vertex);}
+ vector<LorentzTensor<double> > mytempLT;
+ TensorWaveFunction(mytempLT,decay[ix],outgoing,true,false,vertex);
+ tenwave[ix]=mytempLT;
+ }
+ }
// create the matrix element
vector<PDT::Spin> spin;
for(ix=0;ix<decay.size();++ix){spin.push_back(decay[ix]->dataPtr()->iSpin());}
DecayMatrixElement newME(PDT::Spin0,spin);
// loop over the different diagrams
vector<LorentzPolarizationVectorE> form;
Complex fp,f0,A0,A1,A2,A3,V,k;
complex<InvEnergy2> h,bp,bm;
// complex<Energy2> dot;
Lorentz5Momentum q,sum;
Energy2 q2;
Energy MP(part.mass()),MV,msum,mdiff,scale;
LorentzPolarizationVectorE dotv;
double pre;
ParticleVector cpart;
- for(unsigned int iy=0;iy<_CKMfact[mode].size();++iy)
- {
- MV=decay[_formpart[mode][iy]]->mass();
- id1=decay[_formpart[mode][iy]]->id();
- int id0t,id1t;
- _form[_formmapA[mode][iy]]->particleID(_formmapB[mode][iy],id0t,id1t);
- bool cc(id0t!=id0);
- // calculate the form-factor part
- form.resize(0);
- q = part.momentum()-decay[_formpart[mode][iy]]->momentum(); q.rescaleMass();
- sum = part.momentum()+decay[_formpart[mode][iy]]->momentum();sum.rescaleMass();
- q2=q.mass2();
- if(decay[_formpart[mode][iy]]->dataPtr()->iSpin()==1)
- {
- _form[_formmapA[mode][iy]]->ScalarScalarFormFactor(q2,_formmapB[mode][iy],
- id0,id1,MP,MV,f0,fp);
- pre=(MP*MP-MV*MV)/q2;
- form.push_back(fp*sum+pre*(f0-fp)*q);
- }
- else if(decay[_formpart[mode][iy]]->dataPtr()->iSpin()==3)
- {
- msum=MP+MV;mdiff=MP-MV;
- _form[_formmapA[mode][iy]]->ScalarVectorFormFactor(q2,_formmapB[mode][iy],id0,
- id1,MP,MV,A0,A1,A2,V);
- if(cc){V=-V;}
- A3 = 0.5/MV*(msum*A1-mdiff*A2);
- // compute the hadron currents
- for(ix=0;ix<3;++ix)
- {
- // dot product
- complex<Energy>
- dot = vecwave[_formpart[mode][iy]][ix]*part.momentum();
- // current
- form.push_back(-ii*msum*A1*vecwave[_formpart[mode][iy]][ix]
- +ii*A2/msum*dot*sum
- +2.*ii*MV/q2*(A3-A0)*dot*q
- +2.*V/msum*
- Helicity::
- epsilon(vecwave[_formpart[mode][iy]][ix],
- part.momentum(),
- decay[_formpart[mode][iy]]->momentum()));
- }
- }
- else if(decay[_formpart[mode][iy]]->dataPtr()->iSpin()==5)
- {
- _form[_formmapA[mode][iy]]->ScalarTensorFormFactor(q2,_formmapB[mode][iy],
- id0,id1,MP,MV,h,k,bp,bm);
- if(cc){h=-h;}
- // compute the hadron currents
- for(ix=0;ix<5;++ix)
- {
- dotv = tenwave[_formpart[mode][iy]][ix]*part.momentum();
- complex<Energy2> dot = dotv*part.momentum();
- form.push_back(ii*h*Helicity::epsilon(dotv,sum,q)-k*dotv
- -bp*dot*sum-bm*dot*q);
- }
- }
- // find the particles for the current
- cpart.resize(0);
- for(ix=0;ix<decay.size();++ix)
- {if(ix!=_formpart[mode][iy]){cpart.push_back(decay[ix]);}}
- ix=decay.size();
- vector<unsigned int> constants(decay.size()+1),ihel(decay.size()+1);
- int itemp(1);
- do
- {
- --ix;
- if(ix!=_formpart[mode][iy])
- {itemp*=decay[ix]->data().iSpin();constants[ix]=itemp;}
- }
- while(ix!=0);
- constants[decay.size()]=1;
- if(_formpart[mode][iy]!=decay.size())
- {constants[_formpart[mode][iy]]=constants[_formpart[mode][iy]+1];}
- // calculate the current
- vector<LorentzPolarizationVectorE>
- curr=_current[_currentmapA[mode][iy]]->current(vertex,_currentmapB[mode][iy],ichan,
+ for(unsigned int iy=0;iy<_CKMfact[mode].size();++iy) {
+ MV=decay[_formpart[mode][iy]]->mass();
+ id1=decay[_formpart[mode][iy]]->id();
+ int id0t,id1t;
+ _form[_formmapA[mode][iy]]->particleID(_formmapB[mode][iy],id0t,id1t);
+ bool cc(id0t!=id0);
+ // calculate the form-factor part
+ form.resize(0);
+ q = part.momentum()-decay[_formpart[mode][iy]]->momentum(); q.rescaleMass();
+ sum = part.momentum()+decay[_formpart[mode][iy]]->momentum();sum.rescaleMass();
+ q2=q.mass2();
+ if(decay[_formpart[mode][iy]]->dataPtr()->iSpin()==1) {
+ _form[_formmapA[mode][iy]]->ScalarScalarFormFactor(q2,_formmapB[mode][iy],
+ id0,id1,MP,MV,f0,fp);
+ pre=(MP*MP-MV*MV)/q2;
+ form.push_back(fp*sum+pre*(f0-fp)*q);
+ }
+ else if(decay[_formpart[mode][iy]]->dataPtr()->iSpin()==3) {
+ msum=MP+MV;mdiff=MP-MV;
+ _form[_formmapA[mode][iy]]->ScalarVectorFormFactor(q2,_formmapB[mode][iy],id0,
+ id1,MP,MV,A0,A1,A2,V);
+ if(cc){V=-V;}
+ A3 = 0.5/MV*(msum*A1-mdiff*A2);
+ // compute the hadron currents
+ for(ix=0;ix<3;++ix) {
+ // dot product
+ complex<Energy> dot = vecwave[_formpart[mode][iy]][ix]*part.momentum();
+ // current
+ form.push_back(-ii*msum*A1*vecwave[_formpart[mode][iy]][ix]
+ +ii*A2/msum*dot*sum
+ +2.*ii*MV/q2*(A3-A0)*dot*q
+ +2.*V/msum*
+ Helicity::
+ epsilon(vecwave[_formpart[mode][iy]][ix],
+ part.momentum(),
+ decay[_formpart[mode][iy]]->momentum()));
+ }
+ }
+ else if(decay[_formpart[mode][iy]]->dataPtr()->iSpin()==5) {
+ _form[_formmapA[mode][iy]]->ScalarTensorFormFactor(q2,_formmapB[mode][iy],
+ id0,id1,MP,MV,h,k,bp,bm);
+ if(cc){h=-h;}
+ // compute the hadron currents
+ for(ix=0;ix<5;++ix) {
+ dotv = tenwave[_formpart[mode][iy]][ix]*part.momentum();
+ complex<Energy2> dot = dotv*part.momentum();
+ form.push_back(ii*h*Helicity::epsilon(dotv,sum,q)-k*dotv
+ -bp*dot*sum-bm*dot*q);
+ }
+ }
+ // find the particles for the current
+ cpart.resize(0);
+ for(ix=0;ix<decay.size();++ix)
+ {if(ix!=_formpart[mode][iy]){cpart.push_back(decay[ix]);}}
+ ix=decay.size();
+ vector<unsigned int> constants(decay.size()+1),ihel(decay.size()+1);
+ int itemp(1);
+ do {
+ --ix;
+ if(ix!=_formpart[mode][iy]) {
+ itemp*=decay[ix]->data().iSpin();
+ constants[ix]=itemp;
+ }
+ }
+ while(ix!=0);
+ constants[decay.size()]=1;
+ if(_formpart[mode][iy]!=decay.size())
+ constants[_formpart[mode][iy]]=constants[_formpart[mode][iy]+1];
+ // calculate the current
+ vector<LorentzPolarizationVectorE>
+ curr=_current[_currentmapA[mode][iy]]->current(vertex,_currentmapB[mode][iy],ichan,
scale,cpart);
- pre = (pow(part.mass()/scale,int(cpart.size()-2)));
- // loop over the helicities to calculate the matrix element
- ihel[0]=0;
- for(chel=0;chel<curr.size();++chel)
- {
- for(ix=decay.size();ix>0;--ix)
- {if(ix!=_formpart[mode][iy]+1)
- {ihel[ix]=(chel%constants[ix-1])/constants[ix];}}
- for(fhel=0;fhel<form.size();++fhel)
- {
- ihel[_formpart[mode][iy]+1]=fhel;
- newME(ihel)+=pre*_CKMfact[mode][iy]*form[fhel].dot(curr[chel])*_GF;
- }
- }
+ pre = (pow(part.mass()/scale,int(cpart.size()-2)));
+ // loop over the helicities to calculate the matrix element
+ ihel[0]=0;
+ for(chel=0;chel<curr.size();++chel) {
+ for(ix=decay.size();ix>0;--ix) {
+ if(ix!=_formpart[mode][iy]+1)
+ ihel[ix]=(chel%constants[ix-1])/constants[ix];
+ }
+ for(fhel=0;fhel<form.size();++fhel) {
+ ihel[_formpart[mode][iy]+1]=fhel;
+ newME(ihel)+=pre*_CKMfact[mode][iy]*form[fhel].dot(curr[chel])*_GF;
+ }
}
+ }
// store the matrix element
ME(newME);
// perform the contraction
RhoDMatrix rhoin(PDT::Spin0);rhoin.average();
return 0.5*(newME.contract(rhoin)).real();
}
-
+
void ScalarMesonFactorizedDecayer::findModes(unsigned int imode,
vector<PDVector> & particles,
vector<unsigned int> & loc,
- vector<bool> & cc)
-{
+ vector<bool> & cc) {
unsigned int ix,iy,nfound,iz;
// resize the vectors
loc.resize(0);cc.resize(0);
// get the id's for the mode
vector<int> id,idbar;
int idtemp; bool found;
- for(ix=0;ix<particles[imode].size();++ix)
- {
- id.push_back(particles[imode][ix]->id());
- if(particles[imode][ix]->CC()){idbar.push_back(particles[imode][ix]->CC()->id());}
- else{idbar.push_back(id[ix]);}
- }
+ for(ix=0;ix<particles[imode].size();++ix) {
+ id.push_back(particles[imode][ix]->id());
+ if(particles[imode][ix]->CC()) idbar.push_back(particles[imode][ix]->CC()->id());
+ else idbar.push_back(id[ix]);
+ }
vector<bool> done(id.size(),false);
// loop over the modes
- for(ix=0;ix<particles.size();++ix)
- {
- if(ix!=imode)
- {
- // the particle mode
- if(particles[ix][0]->id()==id[0]&&particles[ix].size()==id.size())
- {
- nfound=1;
- for(iy=0;iy<id.size();++iy){done[iy]=false;}
- for(iy=1;iy<id.size();++iy)
- {
- idtemp=particles[ix][iy]->id();
- iz=1;found=false;
- do{if(idtemp==id[iz]&&!done[iz]){done[iz]=true;found=true;}++iz;}
- while(iz<id.size()&&!found);
- if(found){++nfound;}
- }
- if(nfound==id.size()){cc.push_back(false);loc.push_back(ix);}
+ for(ix=0;ix<particles.size();++ix) {
+ if(ix!=imode) {
+ // the particle mode
+ if(particles[ix][0]->id()==id[0]&&particles[ix].size()==id.size()) {
+ nfound=1;
+ for(iy=0;iy<id.size();++iy){done[iy]=false;}
+ for(iy=1;iy<id.size();++iy) {
+ idtemp=particles[ix][iy]->id();
+ iz=1;
+ found=false;
+ do {
+ if(idtemp==id[iz]&&!done[iz]) {
+ done[iz]=true;
+ found=true;
}
- // the charge conjugate mode
- if(particles[ix][0]->id()==idbar[0]&&particles[ix].size()==idbar.size())
- {
- nfound=1;
- for(iy=0;iy<idbar.size();++iy){done[iy]=false;}
- for(iy=1;iy<idbar.size();++iy)
- {
- idtemp=particles[ix][iy]->id();
- iz=1;found=false;
- do{if(idtemp==idbar[iz]&&!done[iz]){done[iz]=true;found=true;}++iz;}
- while(iz<idbar.size()&&!found);
- if(found){++nfound;}
- }
- if(nfound==idbar.size()){cc.push_back(false);loc.push_back(ix);}
+ ++iz;
+ }
+ while(iz<id.size()&&!found);
+ if(found) ++nfound;
+ }
+ if(nfound==id.size()) {
+ cc.push_back(false);
+ loc.push_back(ix);
+ }
+ }
+ // the charge conjugate mode
+ if(particles[ix][0]->id()==idbar[0]&&particles[ix].size()==idbar.size()) {
+ nfound=1;
+ for(iy=0;iy<idbar.size();++iy) done[iy]=false;
+ for(iy=1;iy<idbar.size();++iy) {
+ idtemp=particles[ix][iy]->id();
+ iz=1;
+ found=false;
+ do {
+ if(idtemp==idbar[iz]&&!done[iz]) {
+ done[iz]=true;
+ found=true;
}
+ ++iz;
+ }
+ while(iz<idbar.size()&&!found);
+ if(found){++nfound;}
}
+ if(nfound==idbar.size()){cc.push_back(false);loc.push_back(ix);}
+ }
}
+ }
}
void ScalarMesonFactorizedDecayer::dataBaseOutput(ofstream & output,
- bool header) const
-{
+ bool header) const {
unsigned int ix;
- if(header){output << "update decayers set parameters=\"";}
+ if(header) output << "update decayers set parameters=\"";
DecayIntegrator::dataBaseOutput(output,false);
output << "set " << fullName() << ":GFermi " << _GF*GeV2 << "\n";
output << "set " << fullName() << ":a1Bottom " << _a1b << "\n";
output << "set " << fullName() << ":a2Bottom " << _a2b << "\n";
output << "set " << fullName() << ":a1Charm " << _a1c << "\n";
output << "set " << fullName() << ":a2Charm " << _a2c << "\n";
output << "set " << fullName() << ":CKM " << _theCKM->fullName() << " \n";
- for(ix=0;ix<_current.size();++ix)
- {
- _current[ix]->dataBaseOutput(output,false,true);
- output << "insert " << fullName() << ":Currents " << ix << " "
- << _current[ix]->fullName() << " \n";
- }
- for(ix=0;ix<_form.size();++ix)
- {
- _form[ix]->dataBaseOutput(output,false,true);
- output << "insert " << fullName() << ":FormFactors " << ix << " "
- << _form[ix]->fullName() << " \n";
- }
- for(ix=0;ix<_wgtloc.size();++ix)
- {output << "insert " << fullName() << ":WeightLocation " << ix << " "
- << _wgtloc[ix] << "\n";}
- for(ix=0;ix<_wgtmax.size();++ix)
- {output << "insert " << fullName() << ":MaximumWeight " << ix << " "
- << _wgtmax[ix] << "\n";}
- for(ix=0;ix<_weights.size();++ix)
- {output << "insert " << fullName() << ":Weights " << ix << " "
- << _weights[ix] << "\n";}
- if(header){output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;}
+ for(ix=0;ix<_current.size();++ix) {
+ _current[ix]->dataBaseOutput(output,false,true);
+ output << "insert " << fullName() << ":Currents " << ix << " "
+ << _current[ix]->fullName() << " \n";
+ }
+ for(ix=0;ix<_form.size();++ix) {
+ _form[ix]->dataBaseOutput(output,false,true);
+ output << "insert " << fullName() << ":FormFactors " << ix << " "
+ << _form[ix]->fullName() << " \n";
+ }
+ for(ix=0;ix<_wgtloc.size();++ix) {
+ output << "insert " << fullName() << ":WeightLocation " << ix << " "
+ << _wgtloc[ix] << "\n";
+ }
+ for(ix=0;ix<_wgtmax.size();++ix) {
+ output << "insert " << fullName() << ":MaximumWeight " << ix << " "
+ << _wgtmax[ix] << "\n";
+ }
+ for(ix=0;ix<_weights.size();++ix) {
+ output << "insert " << fullName() << ":Weights " << ix << " "
+ << _weights[ix] << "\n";
+ }
+ if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}
-}
diff --git a/Decay/ScalarMeson/SemiLeptonicScalarDecayer.cc b/Decay/ScalarMeson/SemiLeptonicScalarDecayer.cc
--- a/Decay/ScalarMeson/SemiLeptonicScalarDecayer.cc
+++ b/Decay/ScalarMeson/SemiLeptonicScalarDecayer.cc
@@ -1,300 +1,300 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SemiLeptonicScalarDecayer class.
//
#include "SemiLeptonicScalarDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Helicity/LorentzPolarizationVector.h"
#include "ThePEG/Helicity/epsilon.h"
#include "ThePEG/Helicity/LorentzTensor.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/TensorWaveFunction.h"
using namespace Herwig;
using namespace ThePEG;
using namespace ThePEG::Helicity;
using namespace ThePEG::Helicity;
void SemiLeptonicScalarDecayer::doinit() throw(InitException) {
DecayIntegrator::doinit();
// make sure the current got initialised
_current->init();
// and the form factors
_form->init();
// the channels
PDVector extpart,ptemp;
_modemap.resize(0);
double maxweight;
vector<double> channelwgts(1,1.);
int id0(0),id1(0),Wcharge(0);
Energy min;
DecayPhaseSpaceModePtr mode;
DecayPhaseSpaceChannelPtr channel;
int iq,ia; unsigned int ix,iy,iz;
bool done;
for(ix=0;ix<_form->numberOfFactors();++ix) {
// get the external particles for this mode
extpart.resize(2);
_form->particleID(ix,id0,id1);
extpart[0]=getParticleData(id0);
extpart[1]=getParticleData(id1);
_modemap.push_back(numberModes());
if(!extpart[0]||!extpart[1]) continue;
Wcharge =(extpart[0]->iCharge()-extpart[1]->iCharge());
min = extpart[0]->mass()+extpart[0]->widthUpCut()
-extpart[1]->mass()+extpart[1]->widthLoCut();
for(iy=0;iy<_current->numberOfModes();++iy) {
extpart.resize(2);
_current->decayModeInfo(iy,iq,ia);
ptemp=_current->particles(Wcharge,iy,iq,ia);
for(iz=0;iz<ptemp.size();++iz) extpart.push_back(ptemp[iz]);
// create the mode
mode=new_ptr(DecayPhaseSpaceMode(extpart,this));
// create the first piece of the channel
channel = new_ptr(DecayPhaseSpaceChannel(mode));
channel->addIntermediate(extpart[0],0,0.0,1,-1);
done=_current->createMode(Wcharge,iy,mode,2,1,channel,min);
if(done) {
// the maximum weight
if(_maxwgt.size()>numberModes()) maxweight=_maxwgt[numberModes()];
else maxweight=2.;
channelwgts.resize(mode->numberChannels(),1./(mode->numberChannels()));
addMode(mode,maxweight,channelwgts);
}
}
}
- id0=511;
- id1=-411;
- Energy2 q2=5.22436*GeV2;
- bool cc(false);
- Energy m0(getParticleData(id0)->mass()),m1(getParticleData(id1)->mass());
- unsigned int iloc(_form->formFactorNumber(id0,id1,cc));
- Complex f0,fp;
- _form->ScalarScalarFormFactor(q2,iloc,id0,id1,m0,m1,f0,fp);
+// id0=511;
+// id1=-411;
+// Energy2 q2=5.22436*GeV2;
+// bool cc(false);
+// Energy m0(getParticleData(id0)->mass()),m1(getParticleData(id1)->mass());
+// unsigned int iloc(_form->formFactorNumber(id0,id1,cc));
+// Complex f0,fp;
+// _form->ScalarScalarFormFactor(q2,iloc,id0,id1,m0,m1,f0,fp);
//_form->ScalarVectorFormFactor(q2,iloc,id0,id1,MP,MV,A0,A1,A2,V);
}
bool SemiLeptonicScalarDecayer::accept(const DecayMode & dm) const {
// find the non-lepton
int imes(0),idtemp,idin(dm.parent()->id());
vector<int> idother; bool dummy;
ParticleMSet::const_iterator pit = dm.products().begin();
ParticleMSet::const_iterator pend = dm.products().end();
for( ; pit!=pend;++pit) {
idtemp=(**pit).id();
if(abs(idtemp)>16) imes=idtemp;
else idother.push_back(idtemp);
}
// check that the form factor exists
if(_form->formFactorNumber(idin,imes,dummy)<0) return false;
// and the current
return _current->accept(idother);
}
int SemiLeptonicScalarDecayer::modeNumber(bool & cc,const DecayMode & dm) const {
// find the ids of the particles for the decay current
ParticleMSet::const_iterator pit = dm.products().begin();
ParticleMSet::const_iterator pend = dm.products().end();
int idtemp,imes(0),idin(dm.parent()->id());
vector<int> idother;
cc=false;
for( ; pit!=pend;++pit) {
idtemp=(**pit).id();
if(abs(idtemp)>16) imes=idtemp;
else idother.push_back(idtemp);
}
return _modemap[_form->formFactorNumber(idin,imes,cc)]
+_current->decayMode(idother);
}
void SemiLeptonicScalarDecayer::persistentOutput(PersistentOStream & os) const {
os << _current << _form << _maxwgt << _modemap << ounit(_GF,1/GeV2);
}
void SemiLeptonicScalarDecayer::persistentInput(PersistentIStream & is, int) {
is >> _current >> _form >> _maxwgt >> _modemap >> iunit(_GF,1/GeV2);
}
ClassDescription<SemiLeptonicScalarDecayer> SemiLeptonicScalarDecayer::initSemiLeptonicScalarDecayer;
// Definition of the static class description member.
void SemiLeptonicScalarDecayer::Init() {
static ClassDocumentation<SemiLeptonicScalarDecayer> documentation
("The SemiLeptonicScalarDecayer class is designed for the"
"semi-leptonic decay of a (pseudo)-scalar meson.");
static Parameter<SemiLeptonicScalarDecayer,InvEnergy2> interfaceGFermi
("GFermi",
"The Fermi coupling constant",
&SemiLeptonicScalarDecayer::_GF, 1./GeV2, 1.16639E-5/GeV2,
0./GeV2, 1.e-4/GeV2,
false, false, false);
static Reference<SemiLeptonicScalarDecayer,LeptonNeutrinoCurrent> interfaceCurrent
("Current",
"The current for the leptons produced in the decay.",
&SemiLeptonicScalarDecayer::_current, true, true, true, false, false);
static Reference<SemiLeptonicScalarDecayer,ScalarFormFactor> interfaceFormFactor
("FormFactor",
"The form factor",
&SemiLeptonicScalarDecayer::_form, true, true, true, false, false);
static ParVector<SemiLeptonicScalarDecayer,double> interfaceMaximumWeight
("MaximumWeight",
"The maximum weights for the decays",
&SemiLeptonicScalarDecayer::_maxwgt,
0, 0, 0, 0, 100., false, false, true);
}
// combine the currents and form-factors to give the matrix element
double SemiLeptonicScalarDecayer::me2(bool vertex, const int ichan,
const Particle & inpart,
const ParticleVector & decay) const {
// workaround for gcc 3.2.3 bug
// spin info for the decaying particle
//ALB ScalarWaveFunction(const_ptr_cast<tPPtr>(&inpart),incoming,true,vertex);
tPPtr mytempInpart = const_ptr_cast<tPPtr>(&inpart);
ScalarWaveFunction(mytempInpart,incoming,true,vertex);
// get the information on the form-factor
int jspin(0),id0(inpart.id()),id1(decay[0]->id());
bool cc=false;
unsigned int iloc(_form->formFactorNumber(id0,id1,cc));
int spect,iq,ia;
_form->formFactorInfo(iloc,jspin,spect,iq,ia);
// work out the value of q and calculate the form factors
Lorentz5Momentum q(inpart.momentum()-decay[0]->momentum());
q.rescaleMass();
Energy2 q2(q.mass2());
Lorentz5Momentum sum(inpart.momentum()+decay[0]->momentum());
// calculate the hadronic current for the decay
Complex ii(0.,1.);
vector<LorentzPolarizationVectorE> hadron;
if(jspin==0) {
// workaround for gcc 3.2.3 bug
//ALB ScalarWaveFunction(decay[0],outgoing,true,vertex);
PPtr mytemp = decay[0];
ScalarWaveFunction(mytemp,outgoing,true,vertex);
Complex fp,f0;
_form->ScalarScalarFormFactor(q2,iloc,id0,id1,inpart.mass(),decay[0]->mass(),
f0,fp);
Complex pre((sqr(inpart.mass())-sqr(decay[0]->mass()))/q2*(f0-fp));
hadron.push_back(fp*sum+(pre*q));
}
else if(jspin==1) {
vector<LorentzPolarizationVector> vwave;
VectorWaveFunction(vwave,decay[0],outgoing,true,false,vertex);
Complex A0,A1,A2,A3,V;
complex<Energy> dot;
Energy MP(inpart.mass()),MV(decay[0]->mass()),msum(MP+MV),mdiff(MP-MV);
_form->ScalarVectorFormFactor(q2,iloc,id0,id1,MP,MV,A0,A1,A2,V);
A3 = 0.5/MV*(msum*A1-mdiff*A2);
if(cc) V*=-1.;
// compute the hadron currents
for(unsigned int ix=0;ix<3;++ix) {
// dot product
dot = vwave[ix]*inpart.momentum();
// current
hadron.push_back(-ii*msum*A1*vwave[ix]
+ii*A2/msum*dot*sum
+2.*ii*MV/q2*(A3-A0)*dot*q
- +2.*V/msum*Helicity::epsilon(vwave[ix],inpart.momentum(),
+ -2.*V/msum*Helicity::epsilon(vwave[ix],inpart.momentum(),
decay[0]->momentum()));
}
}
else if(jspin==2) {
vector<LorentzTensor<double> > twave;
TensorWaveFunction(twave,decay[0],outgoing,true,false,vertex);
complex<InvEnergy2> h,bp,bm;
complex<double> k;
complex<Energy2> dot;
_form->ScalarTensorFormFactor(q2,iloc,id0,id1,inpart.mass(),decay[0]->mass(),
h,k,bp,bm);
if(!cc) h*=-1.;
LorentzPolarizationVectorE dotv;
// compute the hadron currents
for(unsigned int ix=0;ix<5;++ix) {
dotv = twave[ix]*inpart.momentum();
dot = dotv*inpart.momentum();
hadron.push_back(ii*h*Helicity::epsilon(dotv,sum,q)
-k*dotv-bp*dot*sum-bm*dot*q);
}
}
int mode=(abs(decay[1]->id())-11)/2;
// construct the lepton current
Energy scale;
ParticleVector leptons;
leptons.push_back(decay[decay.size()-2]);
leptons.push_back(decay[decay.size()-1]);
vector<LorentzPolarizationVectorE>
lepton(_current->current(vertex,mode,ichan,scale,leptons));
// work out the mapping for the lepton vector
vector<unsigned int> constants(decay.size()+1),ihel(decay.size()+1);
vector<PDT::Spin> ispin(decay.size());
unsigned int itemp(1),imes(0);
for(int ix=int(decay.size()-1);ix>=0;--ix) {
ispin[ix]=decay[ix]->data().iSpin();
if(abs(decay[ix]->id())<=16) {
itemp*=ispin[ix];
constants[ix]=itemp;
}
else imes=ix;
}
constants[decay.size()]=1;
constants[imes]=constants[imes+1];
DecayMatrixElement newME(PDT::Spin0,ispin);
for(unsigned int mhel=0;mhel<hadron.size();++mhel) {
for(unsigned int lhel=0;lhel<lepton.size();++lhel) {
// map the index for the leptons to a helicity state
for(unsigned int ix=decay.size();ix>0;--ix) {
if(ix-1!=imes) ihel[ix]=(lhel%constants[ix-1])/constants[ix];
}
// helicities of mesons
ihel[0]=0;
ihel[imes+1]=mhel;
newME(ihel)= lepton[lhel].dot(hadron[mhel])*_GF;
}
}
RhoDMatrix temp(PDT::Spin0); temp.average();
// store the matrix element
ME(newME);
double ckm(1.);
if(iq<=6) {
if(iq%2==0) ckm = SM().CKM(abs(iq)/2-1,(abs(ia)-1)/2);
else ckm = SM().CKM(abs(ia)/2-1,(abs(iq)-1)/2);
}
// return the answer
return 0.5*(newME.contract(temp)).real()*ckm;
}
// output the setup information for the particle database
void SemiLeptonicScalarDecayer::dataBaseOutput(ofstream & output,
bool header) const {
if(header) output << "update decayers set parameters=\"";
DecayIntegrator::dataBaseOutput(output,false);
output << "set " << fullName() << ":GFermi " << _GF*GeV2 << "\n";
for(unsigned int ix=0;ix<_maxwgt.size();++ix) {
output << "insert " << fullName() << ":MaximumWeight " << ix << " "
<< _maxwgt[ix] << "\n";
}
_current->dataBaseOutput(output,false,true);
output << "set " << fullName() << ":Current " << _current->fullName() << " \n";
_form->dataBaseOutput(output,false,true);
output << "set " << fullName() << ":FormFactor " << _form->fullName() << " \n";
if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}
diff --git a/Decay/ScalarMeson/SemiLeptonicScalarDecayer.icc b/Decay/ScalarMeson/SemiLeptonicScalarDecayer.icc
--- a/Decay/ScalarMeson/SemiLeptonicScalarDecayer.icc
+++ b/Decay/ScalarMeson/SemiLeptonicScalarDecayer.icc
@@ -1,37 +1,35 @@
// -*- C++ -*-
//
// This is the implementation of the inlined member functions of
// the SemiLeptonicScalarDecayer class.
//
namespace Herwig {
using namespace ThePEG;
-inline SemiLeptonicScalarDecayer::SemiLeptonicScalarDecayer() {
- // default value of the fermi constant taken from PDG 2002
- _GF = 1.16639E-5/GeV2;
+inline SemiLeptonicScalarDecayer::SemiLeptonicScalarDecayer() : _GF(1.16639E-5/GeV2) {
// intermediates
generateIntermediates(true);
}
inline IBPtr SemiLeptonicScalarDecayer::clone() const {
return new_ptr(*this);
}
inline IBPtr SemiLeptonicScalarDecayer::fullclone() const {
return new_ptr(*this);
}
inline void SemiLeptonicScalarDecayer::doinitrun() {
_current->initrun();
_form->initrun();
DecayIntegrator::doinitrun();
if(initialize()) {
_maxwgt.resize(0);
for(unsigned int ix=0;ix<numberModes();++ix) {
_maxwgt.push_back(mode(ix)->maxWeight());
}
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Sep 30, 5:49 AM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6544063
Default Alt Text
(53 KB)

Event Timeline

HEPForge.org · .