Page MenuHomeHEPForge
Files F19252023

No OneTemporary
Actions

Size
18 KB
Referenced Files
None
Subscribers
None

View Options

diff --git a/Helicity/Correlations/DecayMatrixElement.cc b/Helicity/Correlations/DecayMatrixElement.cc
new file mode 100644
--- /dev/null
+++ b/Helicity/Correlations/DecayMatrixElement.cc
@@ -0,0 +1,144 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the DecayMatrixElement class.
+//
+// Author: Peter Richardson
+//
+
+#include "DecayMatrixElement.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+
+#ifdef ThePEG_TEMPLATES_IN_CC_FILE
+// #include "DecayMatrixElement.tcc"
+#endif
+
+namespace Herwig {
+namespace Helicity {
+
+using namespace ThePEG;
+
+DecayMatrixElement::~DecayMatrixElement() {}
+
+NoPIOClassDescription<DecayMatrixElement> DecayMatrixElement::initDecayMatrixElement;
+// Definition of the static class description member.
+
+void DecayMatrixElement::Init() {
+
+ static ClassDocumentation<DecayMatrixElement> documentation
+ ("The\\classname{DecayMatrixElement} class is designed to store the helicity "
+ "amplitude expression for the matrix element of a decay.");
+
+}
+
+// calculate the decay matrix for this decay
+RhoDMatrix DecayMatrixElement::calculateDMatrix(vector<RhoDMatrix> rhoout)
+{
+ // vectors for the helicities
+ vector<int> ihel1(_outspin.size()+1),ihel2(_outspin.size()+1);
+ // rhomatrix to be returned
+ RhoDMatrix output(_inspin);
+ // make sure that this is zeroed
+ for(int ix=0;ix<_inspin;++ix)
+ {
+ int ixa=ix-_inspin/2;if(_inspin%2==0&&ixa>=0){++ixa;}
+ for(int iy=0;iy<_inspin;++iy)
+ {
+ int iya=iy-_inspin/2;if(_inspin%2==0&&iya>=0){++iya;}
+ output(ixa,iya)=0.;
+ }
+ }
+ // loop over all helicity components of the matrix element
+ // outer loop
+ Complex temp;
+ for(unsigned int ix=0;ix<_matrixelement.size();++ix)
+ {
+ // map the vector index to the helicities
+ for(int ixa=_outspin.size();ixa>0;--ixa)
+ {
+ ihel1[ixa]=(ix%_constants[ixa])/_constants[ixa+1]-int(_outspin[ixa-1]/2);
+ if(_outspin[ixa-1]%2==0&&ihel1[ixa]>=0){++ihel1[ixa];}
+ }
+ ihel1[0]=(ix%_constants[0])/_constants[1]-int(_inspin/2);
+ if(_inspin%2==0&&ihel1[0]>=0){++ihel1[0];}
+ // inner loop
+ for(unsigned int iy=0;iy<_matrixelement.size();++iy)
+ {
+ // map the vector index to the helicities
+ for(unsigned int iya=_outspin.size();iya>0;--iya)
+ {
+ ihel2[iya]=(iy%_constants[iya])/_constants[iya+1]
+ -int(_outspin[iya-1]/2);
+ if(_outspin[iya-1]%2==0&&ihel2[iya]>=0){++ihel2[iya];}
+ }
+ ihel2[0]=(iy%_constants[0])/_constants[1]-int(_inspin/2);
+ if(_inspin%2==0&&ihel2[0]>=0){++ihel2[0];}
+ // matrix element piece
+ temp=_matrixelement[ix]*conj(_matrixelement[iy]);
+ // spin density matrices for the outgoing particles
+ for(unsigned int iz=0;iz<_outspin.size();++iz)
+ {temp*=rhoout[iz](ihel1[iz+1],ihel2[iz+1]);}
+ output(ihel1[0],ihel2[0])+=temp;
+ }
+ }
+ // ensure unit trace for the matrix
+ output.normalize();
+ // return the answer
+ return output;
+}
+
+// calculate the rho matrix for a given outgoing particle
+RhoDMatrix DecayMatrixElement::calculateRhoMatrix(int id,RhoDMatrix rhoin,
+ vector<RhoDMatrix>rhoout)
+{
+ // vectors for the helicities
+ vector<int> ihel1(_outspin.size()+1),ihel2(_outspin.size()+1);
+ // rhomatrix to be returned
+ RhoDMatrix output(_outspin[id]); output.zero();
+ // loop over all helicity components of the matrix element
+ // outer loop
+ Complex temp;
+ for(unsigned int ix=0;ix<_matrixelement.size();++ix)
+ {
+ // map the vector index to the helicities
+ for(unsigned int ixa=_outspin.size();ixa>0;--ixa)
+ {
+ ihel1[ixa]=(ix%_constants[ixa])/_constants[ixa+1]-int(_outspin[ixa-1]/2);
+ if(_outspin[ixa-1]%2==0&&ihel1[ixa]>=0){++ihel1[ixa];}
+ }
+ ihel1[0]=(ix%_constants[0])/_constants[1]-int(_inspin/2);
+ if(_inspin%2==0&&ihel1[0]>=0){++ihel1[0];}
+ // inner loop
+ for(unsigned int iy=0;iy<_matrixelement.size();++iy)
+ {
+ // map the vector index to the helicities
+ for(unsigned int iya=_outspin.size();iya>0;--iya)
+ {
+ ihel2[iya]=(iy%_constants[iya])/_constants[iya+1]
+ -int(_outspin[iya-1]/2);
+ if(_outspin[iya-1]%2==0&&ihel2[iya]>=0){++ihel2[iya];}
+ }
+ ihel2[0]=(iy%_constants[0])/_constants[1]-int(_inspin/2);
+ if(_inspin%2==0&&ihel2[0]>=0){++ihel2[0];}
+ // matrix element piece
+ temp=_matrixelement[ix]*conj(_matrixelement[iy]);
+ // spin denisty matrix for the incoming particle
+ temp*=rhoin(ihel1[0],ihel2[0]);
+ // spin density matrix for the outgoing particles
+ for(unsigned int iz=0;iz<_outspin.size()-1;++iz)
+ {
+ if(int(iz)<id){temp*=rhoout[iz](ihel1[iz+1],ihel2[iz+1]);}
+ else{temp*=rhoout[iz](ihel1[iz+2],ihel2[iz+2]);}
+ }
+ // add to the rho matrix
+ output(ihel1[id+1],ihel2[id+1])+=temp;
+ }
+ }
+ // return the answer
+ output.normalize();
+ return output;
+}
+
+}
+}
+
diff --git a/Helicity/Correlations/DecayMatrixElement.h b/Helicity/Correlations/DecayMatrixElement.h
new file mode 100644
--- /dev/null
+++ b/Helicity/Correlations/DecayMatrixElement.h
@@ -0,0 +1,173 @@
+// -*- C++ -*-
+#ifndef HERWIG_DecayMatrixElement_H
+#define HERWIG_DecayMatrixElement_H
+//
+// This is the declaration of the <!id>DecayMatrixElement<!!id> class.
+//
+// CLASSDOC SUBSECTION Description:
+//
+// Implementation of the complex matrix element for a decay
+// an arbitary number of external particles are supported.
+//
+// CLASSDOC SUBSECTION See also:
+//
+// <a href="RhoDMatrix.html">RhoDMatrix.h</a>,
+// <a href="ProductionMatrixElement.html">ProductionMatrixElement.h</a>,
+// <a href="DecayVertex.html">DecayVertex.h</a>.
+//
+// Author: Peter Richardson
+//
+
+#include <ThePEG/Config/ThePEG.h>
+#include <ThePEG/Utilities/ClassDescription.h>
+#include <ThePEG/Helicity/RhoDMatrix.h>
+// #include "DecayMatrixElement.fh"
+// #include "DecayMatrixElement.xh"
+
+namespace Herwig {
+namespace Helicity {
+
+using namespace ThePEG;
+using ThePEG::Helicity::RhoDMatrix;
+
+class DecayMatrixElement: public Base {
+
+public:
+
+ inline DecayMatrixElement();
+
+ inline DecayMatrixElement(int,int,int);
+ // constructor for two body decay
+
+ inline DecayMatrixElement(int,int,int,int);
+ // constructor for three body decay
+
+ inline DecayMatrixElement(int,int,int,int,int);
+ // constructor for four body decay
+
+ inline DecayMatrixElement(int,int,int,int,int,int);
+ // constructor for five body decay
+
+ inline DecayMatrixElement(int,int,int,int,int,int,int);
+ // constructor for six body decay
+
+ inline DecayMatrixElement(int,vector<int>);
+ inline DecayMatrixElement(vector<int>);
+ // constructors for arbitray body decay
+
+ inline DecayMatrixElement(const DecayMatrixElement &);
+ virtual ~DecayMatrixElement();
+ // Standard ctors and dtor.
+
+public:
+
+ inline int inspin();
+ // get the spin of the incoming particle
+
+ inline vector<int> outspin();
+ // get the spins of the outgoing particles
+
+public:
+ RhoDMatrix calculateDMatrix(vector<RhoDMatrix>);
+ // calculate the decay matrix for this decay
+
+ RhoDMatrix calculateRhoMatrix(int,RhoDMatrix,vector<RhoDMatrix>);
+ // calculate the rho matrix for a given outgoing particle
+
+ inline Complex contract(RhoDMatrix &);
+ // contract the matrix element with the rho matrix of the incoming particle
+
+public:
+
+ // access to the individual helicity components
+ inline Complex operator () (int,int,int) const;
+ inline Complex & operator () (int,int,int);
+ // 2 body decay
+ inline Complex operator () (int,int,int,int) const;
+ inline Complex & operator () (int,int,int,int);
+ // 3 body decay
+ inline Complex operator () (int,int,int,int,int) const;
+ inline Complex & operator () (int,int,int,int,int);
+ // 4 body decay
+ inline Complex operator () (int,int,int,int,int,int) const;
+ inline Complex & operator () (int,int,int,int,int,int);
+ // 5 body decay
+ inline Complex operator () (int,int,int,int,int,int,int) const;
+ inline Complex & operator () (int,int,int,int,int,int,int);
+ // 6 body decay
+ inline Complex operator () (vector<int>) const;
+ inline Complex & operator () (vector<int>);
+ // n body decay
+
+public:
+
+ inline void reset(const DecayMatrixElement &) const;
+ // reset the matrix element
+
+public:
+
+ static void Init();
+ // Standard Init function used to initialize the interfaces.
+
+private:
+
+ static NoPIOClassDescription<DecayMatrixElement> initDecayMatrixElement;
+ // Describe a concrete class without persistent data.
+
+ DecayMatrixElement & operator=(const DecayMatrixElement &);
+ // Private and non-existent assignment operator.
+
+private:
+
+ inline void setMESize();
+ // set the size of the vector containing the matrix element
+
+private:
+
+ mutable int _nout;
+ // number of outgoing particles
+ mutable int _inspin;
+ // spin of the incoming particle as 2s+1
+ mutable vector<int> _outspin;
+ // spins of the outgonig particles
+ mutable vector<Complex> _matrixelement;
+ // storage of the matrix element, a vector is better for memory usage
+ mutable vector<int> _constants;
+ // constants needed to map the index of the vector to a helicity structure
+};
+
+}
+}
+
+// CLASSDOC OFF
+
+namespace ThePEG {
+
+// The following template specialization informs ThePEG about the
+// base class of DecayMatrixElement.
+template <>
+struct BaseClassTrait<Herwig::Helicity::DecayMatrixElement,1> {
+ typedef Base NthBase;
+};
+
+// The following template specialization informs ThePEG about the
+// name of this class and the shared object where it is defined.
+template <>
+struct ClassTraits<Herwig::Helicity::DecayMatrixElement>
+ : public ClassTraitsBase<Herwig::Helicity::DecayMatrixElement> {
+ static string className() { return "/Herwig++/Helicity/DecayMatrixElement"; }
+ // Return the class name.
+ static string library() { return "hwCorrelations.so"; }
+ // Return the name of the shared library to be loaded to get
+ // access to this class and every other class it uses
+ // (except the base class).
+};
+
+}
+
+#include "DecayMatrixElement.icc"
+#ifndef ThePEG_TEMPLATES_IN_CC_FILE
+// #include "DecayMatrixElement.tcc"
+#endif
+
+#endif /* HERWIG_DecayMatrixElement_H */
diff --git a/Helicity/Correlations/DecayMatrixElement.icc b/Helicity/Correlations/DecayMatrixElement.icc
new file mode 100644
--- /dev/null
+++ b/Helicity/Correlations/DecayMatrixElement.icc
@@ -0,0 +1,267 @@
+// -*- C++ -*-
+//
+// This is the implementation of the inlined member functions of
+// the DecayMatrixElement class.
+//
+// Author: Peter Richardson
+//
+
+namespace Herwig {
+namespace Helicity {
+
+using namespace ThePEG;
+
+// default constructor
+inline DecayMatrixElement::DecayMatrixElement() {}
+
+// copy constructor
+inline DecayMatrixElement::DecayMatrixElement(const DecayMatrixElement & x)
+ : Base(x), _nout(x._nout) ,_inspin(x._inspin), _outspin(x._outspin),
+ _matrixelement(x._matrixelement), _constants(x._constants) {}
+
+// constructor for two body decay
+inline DecayMatrixElement::DecayMatrixElement(int in,int out1,int out2)
+{
+ _nout=2; _inspin=in;_outspin.push_back(out1);_outspin.push_back(out2);
+ setMESize();
+}
+
+// constructor for three body decay
+inline DecayMatrixElement::DecayMatrixElement(int in,int out1,int out2,int out3)
+{
+ _nout=3; _inspin=in;_outspin.push_back(out1);_outspin.push_back(out2);
+ _outspin.push_back(out3);
+ setMESize();
+}
+
+// constructor for four body decay
+inline DecayMatrixElement::DecayMatrixElement(int in,int out1,int out2,int out3,
+ int out4)
+{
+ _nout=4; _inspin=in;_outspin.push_back(out1);_outspin.push_back(out2);
+ _outspin.push_back(out3);_outspin.push_back(out4);
+ setMESize();
+}
+
+// constructor for five body decay
+inline DecayMatrixElement::DecayMatrixElement(int in,int out1,int out2,int out3,
+ int out4,int out5)
+{
+ _nout=5; _inspin=in;_outspin.push_back(out1);_outspin.push_back(out2);
+ _outspin.push_back(out3);_outspin.push_back(out4);_outspin.push_back(out5);
+ setMESize();
+}
+
+// constructor for six body decay
+inline DecayMatrixElement::DecayMatrixElement(int in,int out1,int out2,int out3,
+ int out4,int out5,int out6)
+{
+ _nout=6; _inspin=in;_outspin.push_back(out1);_outspin.push_back(out2);
+ _outspin.push_back(out3);_outspin.push_back(out4);_outspin.push_back(out5);
+ _outspin.push_back(out6);
+ setMESize();
+}
+
+// constructor for arbitary body decay
+inline DecayMatrixElement::DecayMatrixElement(int in,vector<int> out)
+{_nout=out.size(); _inspin=in;_outspin=out;setMESize();}
+
+// constructor for arbitray body decay
+inline DecayMatrixElement::DecayMatrixElement(vector<int> out)
+{
+ _nout=out.size();
+ _inspin=out[0];
+ _outspin=vector<int>(out.begin()+1,out.end());
+ setMESize();
+}
+
+// set the size of the vector containing the matrix element
+inline void DecayMatrixElement::setMESize()
+{
+ int isize=_inspin;
+ for(unsigned int ix=0;ix<_outspin.size();++ix)
+ {
+ isize*=_outspin[ix];
+ }
+ // zero the matrix element
+ _matrixelement.resize(isize,0.);
+ // set up the constants for the mapping of helicity to vectro index
+ _constants.resize(_outspin.size()+2);
+ int temp=1;
+ for(unsigned int ix=_outspin.size();ix>0;--ix)
+ {
+ temp*=_outspin[ix-1];_constants[ix]=temp;
+ }
+ temp*=_inspin;_constants[0]=temp;
+ _constants[_outspin.size()+1]=1;
+}
+
+// reset the matrix element
+inline void DecayMatrixElement::reset(const DecayMatrixElement & x) const
+{
+ _nout = x._nout;
+ _inspin = x._inspin;
+ _outspin =x._outspin;
+ _matrixelement=x._matrixelement;
+ _constants=x._constants;
+}
+
+// get the spins of the incoming particles
+inline int DecayMatrixElement::inspin(){return _inspin;}
+
+// get the spins of the outgoing particles
+inline vector<int> DecayMatrixElement::outspin(){return _outspin;}
+
+// access a component for a two body decay
+inline Complex DecayMatrixElement::operator() (int ia,int ib,int ic) const
+{
+ vector<int> itemp(3);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;
+ return (*this)(itemp);
+}
+
+// access a component of a three body decay
+inline Complex DecayMatrixElement::operator() (int ia,int ib,int ic,int id) const
+{
+ vector<int> itemp(4);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;itemp[3]=id;
+ return (*this)(itemp);
+}
+
+// access a component of a four body decay
+inline Complex DecayMatrixElement::operator() (int ia,int ib,int ic,int id,
+ int ie) const
+{
+ vector<int> itemp(5);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;itemp[3]=id;itemp[4]=ie;
+ return (*this)(itemp);
+}
+
+// access a component of a five body decay
+inline Complex DecayMatrixElement::operator() (int ia,int ib,int ic,int id,
+ int ie,int ig) const
+{
+ vector<int> itemp(6);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;
+ itemp[3]=id;itemp[4]=ie;itemp[5]=ig;
+ return (*this)(itemp);
+}
+
+// access a component of a six body decay
+inline Complex DecayMatrixElement::operator() (int ia,int ib,int ic,int id,
+ int ie,int ig, int ih) const
+{
+ vector<int> itemp(7);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;
+ itemp[3]=id;itemp[4]=ie;itemp[5]=ig;itemp[6]=ih;
+ return (*this)(itemp);
+}
+
+// access a component
+inline Complex DecayMatrixElement::operator () (vector<int> in) const
+{
+ if(in.size()!=_outspin.size()+1)
+ {cerr << "Requested the components of a " << in.size()-1
+ << " body decay but this decay is " << _outspin.size() << endl; return 0.;}
+ unsigned int iloc=0,itemp;
+ // contribution for the incoming particle
+ itemp = _inspin/2+ in[0]-int( _inspin%2==0&& in[0]>0);
+ iloc+=itemp*_constants[1];
+ // contributions for the outgoing particles
+ for(unsigned int ix=1;ix<in.size();++ix)
+ {
+ itemp = (_outspin[ix-1]+1)/2+in[ix]-int(_outspin[ix-1]%2==0&&in[ix]>0);
+ iloc+=itemp*_constants[ix+1];
+ }
+ if(iloc>=_matrixelement.size())
+ {cerr << "Invalid component of decay matrix element requested " << iloc << endl;
+ return 0.;}
+ return _matrixelement[iloc];
+}
+
+// set a component for a two body decay
+inline Complex & DecayMatrixElement::operator() (int ia,int ib,int ic)
+{
+ vector<int> itemp(3);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;
+ return (*this)(itemp);
+}
+
+// set a component of a three body decay
+inline Complex & DecayMatrixElement::operator() (int ia,int ib,int ic,int id)
+{
+ vector<int> itemp(4);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;itemp[3]=id;
+ return (*this)(itemp);
+}
+
+// set a component of a four body decay
+inline Complex & DecayMatrixElement::operator() (int ia,int ib,int ic,int id,
+ int ie)
+{
+ vector<int> itemp(5);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;itemp[3]=id;itemp[4]=ie;
+ return (*this)(itemp);
+}
+
+// set a component of a five body decay
+inline Complex & DecayMatrixElement::operator() (int ia,int ib,int ic,int id,
+ int ie,int ig)
+{
+ vector<int> itemp(6);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;
+ itemp[3]=id;itemp[4]=ie;itemp[5]=ig;
+ return (*this)(itemp);
+}
+
+// set a component of a six body decay
+inline Complex & DecayMatrixElement::operator() (int ia,int ib,int ic,int id,
+ int ie,int ig, int ih)
+{
+ vector<int> itemp(7);itemp[0]=ia;itemp[1]=ib;itemp[2]=ic;
+ itemp[3]=id;itemp[4]=ie;itemp[5]=ig;itemp[6]=ih;
+ return (*this)(itemp);
+}
+
+// set a component
+inline Complex & DecayMatrixElement::operator () (vector<int> in)
+{
+ static Complex dummy=0.;
+ if(in.size()!=_outspin.size()+1)
+ {cerr << "Requested the components of a " << in.size()-1
+ << " body decay but this decay is " << _outspin.size() << endl; return dummy;}
+ unsigned int iloc=0,itemp;
+ // contribution for the incoming particle
+ itemp = _inspin/2+ in[0]-int( _inspin%2==0&& in[0]>0);
+ iloc+=itemp*_constants[1];
+ // contributions for the outgoing particles
+ for(unsigned int ix=1;ix<in.size();++ix)
+ {
+ itemp = _outspin[ix-1]/2+in[ix]-int(_outspin[ix-1]%2==0&&in[ix]>0);
+ iloc+=itemp*_constants[ix+1];
+ }
+ if(iloc>=_matrixelement.size())
+ {cerr << "Invalid component of decay matrix element requested " << iloc << endl;
+ return dummy;}
+ return _matrixelement[iloc];
+}
+
+// contract the matrix element with the rho matrix of the incoming particle
+inline Complex DecayMatrixElement::contract(RhoDMatrix &in)
+{
+ Complex me=0.; int in1,in2;
+ for(int ix=0;ix<_constants[1];++ix)
+ {
+ in1=0;
+ for(int inhel1=-_inspin/2;inhel1<_inspin/2+1;++inhel1)
+ {
+ if(inhel1==0&&_inspin%2==0){++inhel1;}
+ in2=0;
+ for(int inhel2=-_inspin/2;inhel2<_inspin/2+1;++inhel2)
+ {
+ if(inhel2==0&&_inspin%2==0){++inhel2;}
+ // compute the term
+ me+=_matrixelement[in1*_constants[1]+ix]*
+ conj(_matrixelement[in2*_constants[1]+ix])*in(inhel1,inhel2);
+ ++in2;
+ }
+ ++in1;
+ }
+ }
+ return me;
+}
+
+}
+}
+

File Metadata

Mime Type
text/x-diff
Expires
Tue, Sep 30, 6:13 AM (5 h, 48 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6566476
Default Alt Text
(18 KB)

Event Timeline

HEPForge.org · .