herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
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
29namespace Herwig {
30
31using namespace ThePEG;
32
41
42public:
43
48
53
54public:
55
63 bool recursive, const DecayMode & dm) const;
64
65public:
66
74
80 void persistentInput(PersistentIStream & is, int version);
82
89 static void Init();
90
91protected:
92
99 EvtParticle *EvtGenParticle(const Particle & part) const;
100
105 ParticleVector decayProducts(EvtParticle* evtpart, bool boost=true) const;
106
111 EvtVector4R EvtGenMomentum(const Lorentz5Momentum & mom) const {
112 return EvtVector4R(mom.t()/GeV,mom.x()/GeV,mom.y()/GeV,mom.z()/GeV);
113 }
114
120 EvtId EvtGenID(int id,bool exception=true) const;
121
127 EvtDiracSpinor EvtGenSpinor(const LorentzSpinor<SqrtEnergy> & sp) const {
128 InvSqrtEnergy norm(sqrt(0.5)/sqrt(GeV));
129 EvtDiracSpinor output;
130 output.set(EvtGenComplex(-norm*( sp.s1()+sp.s3())),
131 EvtGenComplex(-norm*( sp.s2()+sp.s4())),
132 EvtGenComplex(-norm*(-sp.s1()+sp.s3())),
133 EvtGenComplex(-norm*(-sp.s2()+sp.s4())));
134 return output;
135 }
136
141 EvtVector4C EvtGenPolarization(const LorentzPolarizationVector & eps) const {
142 return EvtVector4C(EvtGenComplex(eps.t()),EvtGenComplex(eps.x()),
143 EvtGenComplex(eps.y()),EvtGenComplex(eps.z()));
144 }
145
150 EvtRaritaSchwinger EvtGenRSSpinor(const LorentzRSSpinor<SqrtEnergy> & sp) const {
151 InvSqrtEnergy norm(sqrt(0.5)/sqrt(GeV));
152 complex<double> out[4][4];
153 for(unsigned int ix=0;ix<4;++ix) {
154 out[ix][0] = -Complex(norm*( sp(ix,0)+sp(ix,2)));
155 out[ix][1] = -Complex(norm*( sp(ix,1)+sp(ix,3)));
156 out[ix][2] = -Complex(norm*(-sp(ix,0)+sp(ix,2)));
157 out[ix][3] = -Complex(norm*(-sp(ix,1)+sp(ix,3)));
158 }
159 EvtRaritaSchwinger output;
160 unsigned int ix,iy;
161 // remember we have vec,spin and evtgen spin,vec
162 for(ix=0;ix<4;++ix) {
163 for(iy=0;iy<4;++iy) output.set(ix,iy,EvtGenComplex(out[iy][ix]));
164 }
165 return output;
166 }
167
172 EvtTensor4C EvtGenTensor(const LorentzTensor<double> & ten) const {
173 EvtTensor4C output;
174 unsigned int ix,iy;
175 for(ix=0;ix<4;++ix){
176 for(iy=0;iy<4;++iy) output.set(ix,iy,EvtGenComplex(ten(ix,iy)));
177 }
178 return output;
179 }
180
185 EvtSpinDensity EvtGenSpinDensity(const RhoDMatrix & rho) const {
186 EvtSpinDensity rhoout;
187 unsigned int ix,iy,ispin(rho.iSpin());
188 rhoout.setDim(ispin);
189 for(ix=0;ix<ispin;++ix) {
190 for(iy=0;iy<ispin;++iy)
191 rhoout.set(ix,iy,EvtGenComplex(rho(ix,iy)));
192 }
193 return rhoout;
194 }
195
199 EvtComplex EvtGenComplex(Complex c) const {
200 return EvtComplex(c.real(),c.imag());
201 }
203
214 PPtr ThePEGParticle(EvtParticle *part, tcPDPtr pd,bool spin=true) const {
215 PPtr output(new_ptr(Particle(pd)));
216 output->set5Momentum(ThePEGMomentum(part->getP4(),part->mass()));
217 if(spin) ThePEGSpin(output,part);
218 // make the daughters
219 ParticleVector daug(decayProducts(part,spin));
220 for(unsigned int ix=0;ix<daug.size();++ix) output->addChild(daug[ix]);
221 return output;
222 }
223
229 void ThePEGSpin(PPtr pegpart,EvtParticle *evtpart) const;
230
236 int ThePEGID(EvtId id,bool exception=true) const;
237
243 Lorentz5Momentum ThePEGMomentum(const EvtVector4R & mom,double mass) const {
244 return Lorentz5Momentum(mom.get(1)*GeV,mom.get(2)*GeV,
245 mom.get(3)*GeV,mom.get(0)*GeV,mass*GeV);
246 }
250 Complex ThePEGComplex(EvtComplex c) const {
251 return Complex(real(c),imag(c));
252 }
253
259 RhoDMatrix ThePEGSpinDensity(const EvtSpinDensity & rho, int id) const;
260
266 LorentzSpinor<SqrtEnergy> ThePEGSpinor(const EvtDiracSpinor & sp) const {
267 SqrtEnergy norm(sqrt(0.5)*sqrt(GeV));
268 vector<complex<SqrtEnergy> > evtSpin(4);
269 for(unsigned int ix=0;ix<4;++ix) evtSpin[ix] = -norm*ThePEGComplex(sp.get_spinor(ix));
270 return LorentzSpinor<SqrtEnergy>(evtSpin[0]-evtSpin[2],evtSpin[1]-evtSpin[3],
271 evtSpin[0]+evtSpin[2],evtSpin[1]+evtSpin[3]);
272 }
273
278 LorentzPolarizationVector ThePEGPolarization(const EvtVector4C & eps) const {
279 return LorentzPolarizationVector(conj(ThePEGComplex(eps.get(1))),
280 conj(ThePEGComplex(eps.get(2))),
281 conj(ThePEGComplex(eps.get(3))),
282 conj(ThePEGComplex(eps.get(0))));
283 }
284
289 LorentzRSSpinor<SqrtEnergy> ThePEGRSSpinor(const EvtRaritaSchwinger & sp) const {
290 complex<SqrtEnergy> evtSpin[4][4];
291 SqrtEnergy norm(sqrt(0.5)*sqrt(GeV));
292 // normalisation and swap vec,spin order
293 for(unsigned int ix=0;ix<4;++ix) {
294 for(unsigned int iy=0;iy<4;++iy) evtSpin[ix][iy]=-norm*ThePEGComplex(sp.get(iy,ix));
295 }
297 for(unsigned int ix=0;ix<4;++ix) {
298 output(ix,0) = evtSpin[ix][0] - evtSpin[ix][2];
299 output(ix,1) = evtSpin[ix][1] - evtSpin[ix][3];
300 output(ix,2) = evtSpin[ix][0] + evtSpin[ix][2];
301 output(ix,3) = evtSpin[ix][1] + evtSpin[ix][3];
302 }
303 // output.changeRep(Helicity::defaultDRep);
304 return output;
305 }
306
311 LorentzTensor<double> ThePEGTensor(const EvtTensor4C & ten) const {
313 unsigned int ix,iy;
314 for(ix=0;ix<4;++ix) {
315 for(iy=0;iy<4;++iy)output(ix,iy)=conj(ThePEGComplex(ten.get(ix,iy)));
316 }
317 return output;
318 }
320
324 void checkConversion() const;
325
330 void outputEvtGenDecays(long id) const;
331
336 int EvtGenChannel(const DecayMode &dm) const;
337
342 tSpinPtr getSpinInfo(const Particle &part) const {
343 // return spin info if exists
344 if(part.spinInfo()) {
345 return dynamic_ptr_cast<tSpinPtr>(const_ptr_cast<tPPtr>(&part)->spinInfo());
346 }
347 // otherwise make it
348 tPPtr ptemp(const_ptr_cast<tPPtr>(&part));
349 PDT::Spin spin(part.dataPtr()->iSpin());
350 SpinPtr pspin;
351 if(spin==PDT::Spin0) pspin=new_ptr(ScalarSpinInfo(part.momentum(),true));
352 else if(spin==PDT::Spin1Half) pspin=new_ptr(FermionSpinInfo(part.momentum(),true));
353 else if(spin==PDT::Spin1) pspin=new_ptr(VectorSpinInfo(part.momentum(),true));
354 else if(spin==PDT::Spin3Half) pspin=new_ptr(RSFermionSpinInfo(part.momentum(),true));
355 else if(spin==PDT::Spin2) pspin=new_ptr(TensorSpinInfo(part.momentum(),true));
356 else throw Exception() << "Can't create spinInfo for decaying particle in "
357 << "EvtGen::checkSpinInfo for spin " << spin << "particle "
359 ptemp->spinInfo(pspin);
360 return pspin;
361 }
362
370 void constructVertex(const Particle & parent,ParticleVector products,
371 EvtDecayAmp* damp) const;
372
373protected:
374
381 virtual IBPtr clone() const;
382
387 virtual IBPtr fullclone() const;
389
390protected:
391
398 virtual void doinitrun();
399
404 virtual void dofinish();
406
407private:
408
414
415private:
416
425
429 string pdtName_;
430
434 vector<string> userDecays_;
436
441
446
450 vector<long> convID_;
451
455 string p8Data_;
456
457private:
458
462 EvtRandomEngine * evtrnd_;
463
467 EvtGen * evtgen_;
468
472 mutable ofstream logFile_;
473
474};
475
476}
477
478#endif /* Herwig_EvtGenInterface_H */
The EvtGenInterface class is the main class for the use of the EvtGen decay package with Herwig.
EvtRaritaSchwinger EvtGenRSSpinor(const LorentzRSSpinor< SqrtEnergy > &sp) const
Convert our Rarita-Schwinger spinor to the EvtGen one.
RhoDMatrix ThePEGSpinDensity(const EvtSpinDensity &rho, int id) const
Convert a spin density to a ThePEG one from an EvtGen one.
Complex ThePEGComplex(EvtComplex c) const
Convert from EvtGen complex to ours.
ParticleVector decayProducts(EvtParticle *evtpart, bool boost=true) const
Return the decay products of an EvtGen particle in as ThePEG particles.
void ThePEGSpin(PPtr pegpart, EvtParticle *evtpart) const
Set the SpinInfo for a ThePEG particle using an EvtGen particle.
EvtId EvtGenID(int id, bool exception=true) const
Convert a PDG code from ThePEG into an EvtGen particle id.
EvtGen * evtgen_
Main EvtGen object.
virtual void doinitrun()
Initialize this object.
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
virtual IBPtr clone() const
Make a simple clone of this object.
int EvtGenChannel(const DecayMode &dm) const
Find the location in the EvtGen list of decay channels for a given decay mode.
void checkConversion() const
Check the conversion of particles between Herwig and EvtGen.
virtual void dofinish()
Finalize this object.
vector< string > userDecays_
Names of addition user specified decays.
string pdtName_
The name of the file containing the particle data.
EvtTensor4C EvtGenTensor(const LorentzTensor< double > &ten) const
Convert our tensor to the EvtGen one.
LorentzRSSpinor< SqrtEnergy > ThePEGRSSpinor(const EvtRaritaSchwinger &sp) const
Convert an EvtGen Rarita-Schwinger spinor to ours.
LorentzPolarizationVector ThePEGPolarization(const EvtVector4C &eps) const
Convert an EvtGen complex 4-vector to a LorentzPolarizationVector.
EvtVector4C EvtGenPolarization(const LorentzPolarizationVector &eps) const
Convert a LorentzPolarizationVector to a complex EvtGen 4-vector.
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.
EvtGenInterface()
The default constructor.
EvtRandomEngine * evtrnd_
Pointer to the random number generator for EvtGen.
LorentzSpinor< SqrtEnergy > ThePEGSpinor(const EvtDiracSpinor &sp) const
Convert an EvtDiracSpinor a LorentzSpinor.
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
EvtGenInterface & operator=(const EvtGenInterface &)=delete
The assignment operator is private and must never be called.
string p8Data_
Location of the PYTHIA8 data directory.
int ThePEGID(EvtId id, bool exception=true) const
Convert an EvtGen EvtId to a PDG code in our conventions.
EvtParticle * EvtGenParticle(const Particle &part) const
Convert a particle to an EvtGen particle.
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
ofstream logFile_
File to output the log info to.
bool reDirect_
Whether or not to redirect cout and cerr when EvtGen is running.
void constructVertex(const Particle &parent, ParticleVector products, EvtDecayAmp *damp) const
Construct the DecayVertex for Herwig using the information from EvtGen.
EvtComplex EvtGenComplex(Complex c) const
Convert from our complex to the EvtGen one.
string decayName_
Names of the various EvtGen parameter files.
tSpinPtr getSpinInfo(const Particle &part) const
Check the particle has SpinInfo and if not create it.
vector< long > convID_
Particles for which to output the EvtGen decays.
static void Init()
The standard Init function used to initialize the interfaces.
Lorentz5Momentum ThePEGMomentum(const EvtVector4R &mom, double mass) const
Convert from EvtGen momentum to Lorentz5Momentum.
ParticleVector decay(const Particle &parent, bool recursive, const DecayMode &dm) const
Use EvtGen to perform a decay.
void outputEvtGenDecays(long id) const
Output the EvtGen decay modes for a given particle.
EvtGenInterface(const EvtGenInterface &)
The copy constructor (explicit as cannot copy streams)
bool checkConv_
Check the conversion of the particles.
EvtVector4R EvtGenMomentum(const Lorentz5Momentum &mom) const
Convert a Lorentz5Momentum to a real EvtGen 4-vector.
EvtDiracSpinor EvtGenSpinor(const LorentzSpinor< SqrtEnergy > &sp) const
Convert a LorentzSpinor to an EvtGen one.
PPtr ThePEGParticle(EvtParticle *part, tcPDPtr pd, bool spin=true) const
Functions to convert between EvtGen and Herwig classes.
complex< Value > s2() const
complex< Value > s3() const
complex< Value > s1() const
complex< Value > s4() const
tcSpinPtr spinInfo() const
const Lorentz5Momentum & momentum() const
tcEventPDPtr dataPtr() const
PDT::Spin iSpin() const
-*- C++ -*-
LorentzVector< complex< double > > LorentzPolarizationVector
ThePEG::Ptr< Particle >::transient_pointer tPPtr
ThePEG::Ptr< SpinInfo >::pointer SpinPtr
double sqrt(int x)
std::complex< double > Complex
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
ThePEG::Ptr< Particle >::pointer PPtr
vector< PPtr > ParticleVector
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
ThePEG::Ptr< SpinInfo >::transient_pointer tSpinPtr