herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
EvtGenInterface.h
1 // -*- C++ -*-
2 #ifndef Herwig_EvtGenInterface_H
3 #define Herwig_EvtGenInterface_H
4 //
5 // This is the declaration of the EvtGenInterface class.
6 //
7 
8 #include "EvtGenInterface.fh"
9 #include "ThePEG/Interface/Interfaced.h"
10 #include "ThePEG/Vectors/Lorentz5Vector.h"
11 #include "ThePEG/Helicity/ScalarSpinInfo.h"
12 #include "ThePEG/Helicity/FermionSpinInfo.h"
13 #include "ThePEG/Helicity/VectorSpinInfo.h"
14 #include "ThePEG/Helicity/RSFermionSpinInfo.h"
15 #include "ThePEG/Helicity/TensorSpinInfo.h"
16 #include "ThePEG/EventRecord/Particle.h"
17 
18 #include "EvtGenRandom.h"
19 #include "EvtGen/EvtGen.hh"
20 #include "EvtGenBase/EvtParticle.hh"
21 #include "EvtGenBase/EvtSpinDensity.hh"
22 #include "EvtGenBase/EvtVector4R.hh"
23 #include "EvtGenBase/EvtVector4C.hh"
24 #include "EvtGenBase/EvtTensor4C.hh"
25 #include "EvtGenBase/EvtDiracSpinor.hh"
26 #include "EvtGenBase/EvtRaritaSchwinger.hh"
27 #include "EvtGenBase/EvtDecayAmp.hh"
28 
29 namespace Herwig {
30 
31 using namespace ThePEG;
32 
40 class EvtGenInterface: public Interfaced {
41 
42 public:
43 
50 
55 
59  virtual ~EvtGenInterface();
61 
62 public:
63 
70  ParticleVector decay(const Particle &parent,
71  bool recursive, const DecayMode & dm) const;
72 
73 public:
74 
81  void persistentOutput(PersistentOStream & os) const;
82 
88  void persistentInput(PersistentIStream & is, int version);
90 
97  static void Init();
98 
99 protected:
100 
107  EvtParticle *EvtGenParticle(const Particle & part) const;
108 
113  ParticleVector decayProducts(EvtParticle* evtpart, bool boost=true) const;
114 
119  EvtVector4R EvtGenMomentum(const Lorentz5Momentum & mom) const {
120  return EvtVector4R(mom.t()/GeV,mom.x()/GeV,mom.y()/GeV,mom.z()/GeV);
121  }
122 
128  EvtId EvtGenID(int id,bool exception=true) const;
129 
135  EvtDiracSpinor EvtGenSpinor(const LorentzSpinor<SqrtEnergy> & sp) const {
136  InvSqrtEnergy norm(sqrt(0.5)/sqrt(GeV));
137  EvtDiracSpinor output;
138  output.set(EvtGenComplex(-norm*( sp.s1()+sp.s3())),
139  EvtGenComplex(-norm*( sp.s2()+sp.s4())),
140  EvtGenComplex(-norm*(-sp.s1()+sp.s3())),
141  EvtGenComplex(-norm*(-sp.s2()+sp.s4())));
142  return output;
143  }
144 
149  EvtVector4C EvtGenPolarization(const LorentzPolarizationVector & eps) const {
150  return EvtVector4C(EvtGenComplex(eps.t()),EvtGenComplex(eps.x()),
151  EvtGenComplex(eps.y()),EvtGenComplex(eps.z()));
152  }
153 
158  EvtRaritaSchwinger EvtGenRSSpinor(const LorentzRSSpinor<SqrtEnergy> & sp) const {
159  InvSqrtEnergy norm(sqrt(0.5)/sqrt(GeV));
160  complex<double> out[4][4];
161  for(unsigned int ix=0;ix<4;++ix) {
162  out[ix][0] = -Complex(norm*( sp(ix,0)+sp(ix,2)));
163  out[ix][1] = -Complex(norm*( sp(ix,1)+sp(ix,3)));
164  out[ix][2] = -Complex(norm*(-sp(ix,0)+sp(ix,2)));
165  out[ix][3] = -Complex(norm*(-sp(ix,1)+sp(ix,3)));
166  }
167  EvtRaritaSchwinger output;
168  unsigned int ix,iy;
169  // remember we have vec,spin and evtgen spin,vec
170  for(ix=0;ix<4;++ix) {
171  for(iy=0;iy<4;++iy) output.set(ix,iy,EvtGenComplex(out[iy][ix]));
172  }
173  return output;
174  }
175 
180  EvtTensor4C EvtGenTensor(const LorentzTensor<double> & ten) const {
181  EvtTensor4C output;
182  unsigned int ix,iy;
183  for(ix=0;ix<4;++ix){
184  for(iy=0;iy<4;++iy) output.set(ix,iy,EvtGenComplex(ten(ix,iy)));
185  }
186  return output;
187  }
188 
193  EvtSpinDensity EvtGenSpinDensity(const RhoDMatrix & rho) const {
194  EvtSpinDensity rhoout;
195  unsigned int ix,iy,ispin(rho.iSpin());
196  rhoout.setDim(ispin);
197  for(ix=0;ix<ispin;++ix) {
198  for(iy=0;iy<ispin;++iy)
199  rhoout.set(ix,iy,EvtGenComplex(rho(ix,iy)));
200  }
201  return rhoout;
202  }
203 
207  EvtComplex EvtGenComplex(Complex c) const {
208  return EvtComplex(c.real(),c.imag());
209  }
211 
222  PPtr ThePEGParticle(EvtParticle *part, tcPDPtr pd,bool spin=true) const {
223  PPtr output(new_ptr(Particle(pd)));
224  output->set5Momentum(ThePEGMomentum(part->getP4(),part->mass()));
225  if(spin) ThePEGSpin(output,part);
226  // make the daughters
227  ParticleVector daug(decayProducts(part,spin));
228  for(unsigned int ix=0;ix<daug.size();++ix) output->addChild(daug[ix]);
229  return output;
230  }
231 
237  void ThePEGSpin(PPtr pegpart,EvtParticle *evtpart) const;
238 
244  int ThePEGID(EvtId id,bool exception=true) const;
245 
251  Lorentz5Momentum ThePEGMomentum(const EvtVector4R & mom,double mass) const {
252  return Lorentz5Momentum(mom.get(1)*GeV,mom.get(2)*GeV,
253  mom.get(3)*GeV,mom.get(0)*GeV,mass*GeV);
254  }
258  Complex ThePEGComplex(EvtComplex c) const {
259  return Complex(real(c),imag(c));
260  }
261 
267  RhoDMatrix ThePEGSpinDensity(const EvtSpinDensity & rho, int id) const;
268 
274  LorentzSpinor<SqrtEnergy> ThePEGSpinor(const EvtDiracSpinor & sp) const {
275  SqrtEnergy norm(sqrt(0.5)*sqrt(GeV));
276  vector<complex<SqrtEnergy> > evtSpin(4);
277  for(unsigned int ix=0;ix<4;++ix) evtSpin[ix] = -norm*ThePEGComplex(sp.get_spinor(ix));
278  return LorentzSpinor<SqrtEnergy>(evtSpin[0]-evtSpin[2],evtSpin[1]-evtSpin[3],
279  evtSpin[0]+evtSpin[2],evtSpin[1]+evtSpin[3]);
280  }
281 
286  LorentzPolarizationVector ThePEGPolarization(const EvtVector4C & eps) const {
287  return LorentzPolarizationVector(conj(ThePEGComplex(eps.get(1))),
288  conj(ThePEGComplex(eps.get(2))),
289  conj(ThePEGComplex(eps.get(3))),
290  conj(ThePEGComplex(eps.get(0))));
291  }
292 
297  LorentzRSSpinor<SqrtEnergy> ThePEGRSSpinor(const EvtRaritaSchwinger & sp) const {
298  complex<SqrtEnergy> evtSpin[4][4];
299  SqrtEnergy norm(sqrt(0.5)*sqrt(GeV));
300  // normalisation and swap vec,spin order
301  for(unsigned int ix=0;ix<4;++ix) {
302  for(unsigned int iy=0;iy<4;++iy) evtSpin[ix][iy]=-norm*ThePEGComplex(sp.get(iy,ix));
303  }
305  for(unsigned int ix=0;ix<4;++ix) {
306  output(ix,0) = evtSpin[ix][0] - evtSpin[ix][2];
307  output(ix,1) = evtSpin[ix][1] - evtSpin[ix][3];
308  output(ix,2) = evtSpin[ix][0] + evtSpin[ix][2];
309  output(ix,3) = evtSpin[ix][1] + evtSpin[ix][3];
310  }
311  // output.changeRep(Helicity::defaultDRep);
312  return output;
313  }
314 
319  LorentzTensor<double> ThePEGTensor(const EvtTensor4C & ten) const {
320  LorentzTensor<double> output;
321  unsigned int ix,iy;
322  for(ix=0;ix<4;++ix) {
323  for(iy=0;iy<4;++iy)output(ix,iy)=conj(ThePEGComplex(ten.get(ix,iy)));
324  }
325  return output;
326  }
328 
332  void checkConversion() const;
333 
338  void outputEvtGenDecays(long id) const;
339 
344  int EvtGenChannel(const DecayMode &dm) const;
345 
350  tSpinPtr getSpinInfo(const Particle &part) const {
351  // return spin info if exists
352  if(part.spinInfo()) {
353  return dynamic_ptr_cast<tSpinPtr>(const_ptr_cast<tPPtr>(&part)->spinInfo());
354  }
355  // otherwise make it
356  tPPtr ptemp(const_ptr_cast<tPPtr>(&part));
357  PDT::Spin spin(part.dataPtr()->iSpin());
358  SpinPtr pspin;
359  if(spin==PDT::Spin0) pspin=new_ptr(ScalarSpinInfo(part.momentum(),true));
360  else if(spin==PDT::Spin1Half) pspin=new_ptr(FermionSpinInfo(part.momentum(),true));
361  else if(spin==PDT::Spin1) pspin=new_ptr(VectorSpinInfo(part.momentum(),true));
362  else if(spin==PDT::Spin3Half) pspin=new_ptr(RSFermionSpinInfo(part.momentum(),true));
363  else if(spin==PDT::Spin2) pspin=new_ptr(TensorSpinInfo(part.momentum(),true));
364  else throw Exception() << "Can't create spinInfo for decaying particle in "
365  << "EvtGen::checkSpinInfo for spin " << spin << "particle "
367  ptemp->spinInfo(pspin);
368  return pspin;
369  }
370 
378  void constructVertex(const Particle & parent,ParticleVector products,
379  EvtDecayAmp* damp) const;
380 
381 protected:
382 
389  virtual IBPtr clone() const;
390 
395  virtual IBPtr fullclone() const;
397 
398 protected:
399 
406  virtual void doinitrun();
407 
412  virtual void dofinish();
414 
415 private:
416 
421  EvtGenInterface & operator=(const EvtGenInterface &) = delete;
422 
423 private:
424 
432  string decayName_;
433 
437  string pdtName_;
438 
442  vector<string> userDecays_;
444 
448  bool reDirect_;
449 
454 
458  vector<long> convID_;
459 
463  string p8Data_;
464 
465 private:
466 
470  EvtRandomEngine * evtrnd_;
471 
475  EvtGen * evtgen_;
476 
480  mutable ofstream logFile_;
481 
482 };
483 
484 }
485 
486 #endif /* Herwig_EvtGenInterface_H */
complex< Value > s2() const
ThePEG::Ptr< Particle >::pointer PPtr
string pdtName_
The name of the file containing the particle data.
bool reDirect_
Whether or not to redirect cout and cerr when EvtGen is running.
double sqrt(int x)
vector< long > convID_
Particles for which to output the EvtGen decays.
std::complex< double > Complex
ThePEG::Ptr< SpinInfo >::transient_pointer tSpinPtr
vector< string > userDecays_
Names of addition user specified decays.
Definition: expm-1.h:45
LorentzSpinor< SqrtEnergy > ThePEGSpinor(const EvtDiracSpinor &sp) const
Convert an EvtDiracSpinor a LorentzSpinor.
EvtComplex EvtGenComplex(Complex c) const
Convert from our complex to the EvtGen one.
tSpinPtr getSpinInfo(const Particle &part) const
Check the particle has SpinInfo and if not create it.
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
LorentzRSSpinor< SqrtEnergy > ThePEGRSSpinor(const EvtRaritaSchwinger &sp) const
Convert an EvtGen Rarita-Schwinger spinor to ours.
EvtRandomEngine * evtrnd_
Pointer to the random number generator for EvtGen.
ofstream logFile_
File to output the log info to.
EvtTensor4C EvtGenTensor(const LorentzTensor< double > &ten) const
Convert our tensor to the EvtGen one.
complex< Value > s4() const
complex< Value > s3() const
complex< Value > s1() const
const Lorentz5Momentum & momentum() const
LorentzPolarizationVector ThePEGPolarization(const EvtVector4C &eps) const
Convert an EvtGen complex 4-vector to a LorentzPolarizationVector.
EvtRaritaSchwinger EvtGenRSSpinor(const LorentzRSSpinor< SqrtEnergy > &sp) const
Convert our Rarita-Schwinger spinor to the EvtGen one.
PPtr ThePEGParticle(EvtParticle *part, tcPDPtr pd, bool spin=true) const
Functions to convert between EvtGen and Herwig classes.
EvtSpinDensity EvtGenSpinDensity(const RhoDMatrix &rho) const
Convert a spin density matrix to an EvtGen spin density matrix.
LorentzTensor< double > ThePEGTensor(const EvtTensor4C &ten) const
Convert an EvtGen tensor to ThePEG.
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
Lorentz5Momentum ThePEGMomentum(const EvtVector4R &mom, double mass) const
Convert from EvtGen momentum to Lorentz5Momentum.
bool checkConv_
Check the conversion of the particles.
string decayName_
Names of the various EvtGen parameter files.
ThePEG::Ptr< Particle >::transient_pointer tPPtr
EvtVector4R EvtGenMomentum(const Lorentz5Momentum &mom) const
Convert a Lorentz5Momentum to a real EvtGen 4-vector.
Complex ThePEGComplex(EvtComplex c) const
Convert from EvtGen complex to ours.
ThePEG::Ptr< SpinInfo >::pointer SpinPtr
EvtVector4C EvtGenPolarization(const LorentzPolarizationVector &eps) const
Convert a LorentzPolarizationVector to a complex EvtGen 4-vector.
tcEventPDPtr dataPtr() const
tcSpinPtr spinInfo() const
EvtDiracSpinor EvtGenSpinor(const LorentzSpinor< SqrtEnergy > &sp) const
Convert a LorentzSpinor to an EvtGen one.
-*- C++ -*-
EvtGen * evtgen_
Main EvtGen object.
vector< PPtr > ParticleVector
string p8Data_
Location of the PYTHIA8 data directory.
PDT::Spin iSpin() const
The EvtGenInterface class is the main class for the use of the EvtGen decay package with Herwig...