diff --git a/Decay/WeakCurrents/WeakBaryonCurrent.cc b/Decay/WeakCurrents/WeakBaryonCurrent.cc
--- a/Decay/WeakCurrents/WeakBaryonCurrent.cc
+++ b/Decay/WeakCurrents/WeakBaryonCurrent.cc
@@ -1,214 +1,213 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the WeakBaryonCurrent class.
//
#include "WeakBaryonCurrent.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "ThePEG/Helicity/HelicityFunctions.h"
using namespace Herwig;
WeakBaryonCurrent::WeakBaryonCurrent() {}
IBPtr WeakBaryonCurrent::clone() const {
return new_ptr(*this);
}
IBPtr WeakBaryonCurrent::fullclone() const {
return new_ptr(*this);
}
void WeakBaryonCurrent::persistentOutput(PersistentOStream & os) const {
os << formFactor_;
}
void WeakBaryonCurrent::persistentInput(PersistentIStream & is, int) {
is >> formFactor_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<WeakBaryonCurrent,WeakCurrent>
describeHerwigWeakBaryonCurrent("Herwig::WeakBaryonCurrent", "Herwig.so");
void WeakBaryonCurrent::Init() {
static ClassDocumentation<WeakBaryonCurrent> documentation
("The WeakBaryonCurrent class is a wrapper for the BaryonFormFactor"
" so it can be used as a WeakCurrent");
static Reference<WeakBaryonCurrent,BaryonFormFactor> interfaceFormFactor
("FormFactor",
"The baryon form factor",
&WeakBaryonCurrent::formFactor_, false, false, true, false, false);
}
void WeakBaryonCurrent::doinit() {
// initialize the form factor
formFactor_->init();
// set the modes
for(unsigned int iloc=0;iloc<formFactor_->numberOfFactors();++iloc) {
int ispin(0), ospin(0), spect1(0), spect2(0), inquark(0), outquark(0);
formFactor_->formFactorInfo(iloc,ispin,ospin,spect1,spect2,inquark,outquark);
addDecayMode(outquark,-inquark);
}
setInitialModes(formFactor_->numberOfFactors());
WeakCurrent::doinit();
}
// complete the construction of the decay mode for integration
bool WeakBaryonCurrent::createMode(int icharge, tcPDPtr ,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
unsigned int imode,PhaseSpaceModePtr mode,
unsigned int iloc,int ires,
PhaseSpaceChannel phase, Energy upp ) {
// todo isospin in the form factors
// no isospin here
if(Itotal!=IsoSpin::IUnknown || i3 !=IsoSpin::I3Unknown) return false;
unsigned int iq(0),ia(0);
tPDVector out = particles(icharge,imode,iq,ia);
// make sure the the decays are kinematically allowed
Energy min =out[0]->massMin()+out[1]->massMin();
if(min>=upp) return false;
// set the resonances and check charge
tPDPtr res;
if(icharge==3) res=getParticleData(ParticleID::Wplus );
else if(icharge==-3) res=getParticleData(ParticleID::Wminus);
else res=getParticleData(ParticleID::gamma );
// create the channel
mode->addChannel((PhaseSpaceChannel(phase),ires,res,ires+1,iloc+1,ires+1,iloc+2));
// return if successful
return true;
}
// the particles produced by the current
tPDVector WeakBaryonCurrent::particles(int icharge, unsigned int imode, int , int ) {
tPDVector extpart(2);
int id0(0),id1(0);
formFactor_->particleID(imode,id0,id1);
extpart[0] = getParticleData(id0);
if(extpart[0]->CC()) extpart[0]=extpart[0]->CC();
extpart[1] = getParticleData(id1);
int charge = extpart[0]->iCharge()+extpart[1]->iCharge();
if(charge==icharge)
return extpart;
else if(charge==-icharge) {
for(unsigned int ix=0;ix<2;++ix)
if(extpart[ix]->CC()) extpart[ix]=extpart[ix]->CC();
return extpart;
}
else
return tPDVector();
}
void WeakBaryonCurrent::constructSpinInfo(ParticleVector decay) const {
- assert(false);
- // if(decay[0]->id()>0) {
- // SpinorWaveFunction ::constructSpinInfo(wave_ ,decay[1],outgoing,true);
- // SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
- // }
- // else {
- // SpinorWaveFunction ::constructSpinInfo( wave_,decay[0],outgoing,true);
- // SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[1],outgoing,true);
- // }
+ if(decay[0]->id()>0) {
+ SpinorWaveFunction ::constructSpinInfo(wave_ ,decay[1],outgoing,true);
+ SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
+ }
+ else {
+ SpinorWaveFunction ::constructSpinInfo( wave_,decay[0],outgoing,true);
+ SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[1],outgoing,true);
+ }
}
// hadronic current
vector<LorentzPolarizationVectorE>
WeakBaryonCurrent::current(tcPDPtr ,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
const int, const int, Energy & scale,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
DecayIntegrator::MEOption) const {
// no isospin here
if(Itotal!=IsoSpin::IUnknown || i3 !=IsoSpin::I3Unknown) return vector<LorentzPolarizationVectorE>();
useMe();
Lorentz5Momentum q = momenta[0]+momenta[1];
q.rescaleMass();
scale=q.mass();
int in = abs(outgoing[0]->id());
int out = abs(outgoing[1]->id());
Energy m1 = outgoing[0]->mass();
Energy m2 = outgoing[1]->mass();
bool cc = false;
unsigned int imode = formFactor_->formFactorNumber(in,out,cc);
// todo generalize to spin != 1/2
assert(outgoing[0]->iSpin()==PDT::Spin1Half &&
outgoing[1]->iSpin()==PDT::Spin1Half );
- vector<LorentzSpinor <SqrtEnergy> > f2(2);
- vector<LorentzSpinorBar<SqrtEnergy> > a2(2);
+ wave_.resize(2);
+ wavebar_.resize(2);
for(unsigned int ix=0;ix<2;++ix) {
- a2[ix] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ix,Helicity::outgoing);
- f2[ix] = HelicityFunctions::dimensionedSpinor (-momenta[1],ix,Helicity::outgoing);
+ wavebar_[ix] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ix,Helicity::outgoing);
+ wave_[ix] = HelicityFunctions::dimensionedSpinor (-momenta[1],ix,Helicity::outgoing);
}
// get the form factors
Complex f1v(0.),f2v(0.),f3v(0.),f1a(0.),f2a(0.),f3a(0.);
formFactor_->SpinHalfSpinHalfFormFactor(sqr(scale),imode,in,out,m1,m2,
f1v,f2v,f3v,f1a,f2a,f3a,
BaryonFormFactor::TimeLike);
Complex left = f1v - f1a + f2v -double((m1-m2)/(m1+m2))*f2a;
Complex right = f1v + f1a + f2v +double((m1-m2)/(m1+m2))*f2a;
vector<LorentzPolarizationVectorE> baryon;
Lorentz5Momentum diff = momenta[0]-momenta[1];
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
LorentzPolarizationVectorE
- vtemp = f2[ohel2].generalCurrent(a2[ohel1],left,right);
- complex<Energy> vspin=f2[ohel2].scalar (a2[ohel1]);
- complex<Energy> aspin=f2[ohel2].pseudoScalar(a2[ohel1]);
+ vtemp = wave_[ohel2].generalCurrent(wavebar_[ohel1],left,right);
+ complex<Energy> vspin=wave_[ohel2].scalar (wavebar_[ohel1]);
+ complex<Energy> aspin=wave_[ohel2].pseudoScalar(wavebar_[ohel1]);
vtemp-= (f2v*vspin+f2a*aspin)/(m1+m2)*diff;
vtemp+= (f3v*vspin+f3a*aspin)/(m1+m2)*q;
baryon.push_back(vtemp);
}
}
// return the answer
return baryon;
}
bool WeakBaryonCurrent::accept(vector<int> id) {
assert(false);
// bool allowed(false);
// if(id.size()!=2) return false;
// if(abs(id[0])%2==0) {
// if((id[0]> 10&&id[0]< 18&&id[1]==-id[0]+1)||
// (id[0]<-10&&id[0]>-18&&id[1]==-id[0]-1)) allowed=true;
// }
// else {
// if((id[1]> 10&&id[1]< 18&&id[0]==-id[1]+1)||
// (id[1]<-10&&id[1]>-18&&id[0]==-id[1]-1)) allowed=true;
// }
// return allowed;
}
// the decay mode
unsigned int WeakBaryonCurrent::decayMode(vector<int> idout) {
assert(idout.size()==2);
int itemp[2] = {idout[0],idout[1]};
for(unsigned int ix=0;ix<2;++ix)
if(itemp[ix]<0) itemp[ix]=-itemp[ix];
bool cc = false;
return formFactor_->formFactorNumber(itemp[0],itemp[1],cc);
}
// output the information for the database
void WeakBaryonCurrent::dataBaseOutput(ofstream & output,bool header,
bool create) const {
if(header) output << "update decayers set parameters=\"";
if(create) output << "create Herwig::WeakBaryonCurrent " << name() << "\n";
output << "newdef " << name() << ":FormFactor " << formFactor_ << "\n";
WeakCurrent::dataBaseOutput(output,false,false);
if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}
diff --git a/Decay/WeakCurrents/WeakBaryonCurrent.h b/Decay/WeakCurrents/WeakBaryonCurrent.h
--- a/Decay/WeakCurrents/WeakBaryonCurrent.h
+++ b/Decay/WeakCurrents/WeakBaryonCurrent.h
@@ -1,194 +1,207 @@
// -*- C++ -*-
#ifndef Herwig_WeakBaryonCurrent_H
#define Herwig_WeakBaryonCurrent_H
//
// This is the declaration of the WeakBaryonCurrent class.
//
#include "WeakCurrent.h"
#include "Herwig/Decay/FormFactors/BaryonFormFactor.h"
+#include "ThePEG/Helicity/LorentzSpinorBar.h"
namespace Herwig {
using namespace ThePEG;
/**
* The WeakBaryonCurrent class is a wrapper around the BaryonFormFactor so that
* they can be used as a current.
*
* @see \ref WeakBaryonCurrentInterfaces "The interfaces"
* defined for WeakBaryonCurrent.
*/
class WeakBaryonCurrent: public WeakCurrent {
public:
/**
* The default constructor.
*/
WeakBaryonCurrent();
public:
/** @name Methods for the construction of the phase space integrator. */
//@{
/**
* Complete the construction of the decay mode for integration.classes inheriting
* from this one.
* This method is purely virtual and must be implemented in the classes inheriting
* from WeakCurrent.
* @param icharge The total charge of the outgoing particles in the current.
* @param resonance If specified only include terms with this particle
* @param Itotal If specified the total isospin of the current
* @param I3 If specified the thrid component of isospin
* @param imode The mode in the current being asked for.
* @param mode The phase space mode for the integration
* @param iloc The location of the of the first particle from the current in
* the list of outgoing particles.
* @param ires The location of the first intermediate for the current.
* @param phase The prototype phase space channel for the integration.
* @param upp The maximum possible mass the particles in the current are
* allowed to have.
* @return Whether the current was sucessfully constructed.
*/
virtual bool createMode(int icharge, tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
unsigned int imode,PhaseSpaceModePtr mode,
unsigned int iloc,int ires,
PhaseSpaceChannel phase, Energy upp );
/**
* The particles produced by the current. This just returns the leptons.
* @param icharge The total charge of the particles in the current.
* @param imode The mode for which the particles are being requested
* @param iq The PDG code for the quark
* @param ia The PDG code for the antiquark
* @return The external particles for the current.
*/
virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
//@}
/**
* Hadronic current. This method is purely virtual and must be implemented in
* all classes inheriting from this one.
* @param resonance If specified only include terms with this particle
* @param Itotal If specified the total isospin of the current
* @param I3 If specified the thrid component of isospin
* @param imode The mode
* @param ichan The phase-space channel the current is needed for.
* @param scale The invariant mass of the particles in the current.
* @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 current.
*/
virtual vector<LorentzPolarizationVectorE>
current(tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
const int imode, const int ichan,Energy & scale,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
DecayIntegrator::MEOption meopt) const;
/**
* Construct the SpinInfo for the decay products
*/
void constructSpinInfo(ParticleVector decay) const;
/**
* Accept the decay. Checks that this is one of the allowed modes.
* @param id The id's of the particles in the current.
* @return Can this current have the external particles specified.
*/
virtual bool accept(vector<int> id);
/**
* Returns the decay mode number for a given set of particles in the current.
* @param id The id's of the particles in the current.
* @return The number of the mode
*/
virtual unsigned int decayMode(vector<int> id);
/**
* Output the setup information for the particle database
* @param os The stream to output the information to
* @param header Whether or not to output the information for MySQL
* @param create Whether or not to add a statement creating the object
*/
virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @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();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- WeakBaryonCurrent & operator=(const WeakBaryonCurrent &);
+ WeakBaryonCurrent & operator=(const WeakBaryonCurrent &) = delete;
private:
/**
* The baryon form factor
*/
BaryonFormFactorPtr formFactor_;
+
+private:
+
+ /**
+ * Spinors for the decay products
+ */
+ mutable vector<Helicity::LorentzSpinor <SqrtEnergy> > wave_;
+
+ /**
+ * barred spinors for the decay products
+ */
+ mutable vector<Helicity::LorentzSpinorBar<SqrtEnergy> > wavebar_;
};
}
#endif /* Herwig_WeakBaryonCurrent_H */