Decay/EvtGen/EvtGenInterface.h

0001 // -*- C++ -*-
0002 #ifndef Herwig_EvtGenInterface_H
0003 #define Herwig_EvtGenInterface_H
0004 //
0005 // This is the declaration of the EvtGenInterface class.
0006 //
0007
0008 #include "EvtGenInterface.fh"
0009 #include "ThePEG/Interface/Interfaced.h"
0010 #include "ThePEG/Vectors/Lorentz5Vector.h"
0011 #include "ThePEG/Helicity/ScalarSpinInfo.h"
0012 #include "ThePEG/Helicity/FermionSpinInfo.h"
0013 #include "ThePEG/Helicity/VectorSpinInfo.h"
0014 #include "ThePEG/Helicity/RSFermionSpinInfo.h"
0015 #include "ThePEG/Helicity/TensorSpinInfo.h"
0016 #include "ThePEG/EventRecord/Particle.h"
0017
0018 #include "EvtGenRandom.h"
0019 #include "EvtGen/EvtGen.hh"
0020 #include "EvtGenBase/EvtParticle.hh"
0021 #include "EvtGenBase/EvtSpinDensity.hh"
0022 #include "EvtGenBase/EvtVector4R.hh"
0023 #include "EvtGenBase/EvtVector4C.hh"
0024 #include "EvtGenBase/EvtTensor4C.hh"
0025 #include "EvtGenBase/EvtDiracSpinor.hh"
0026 #include "EvtGenBase/EvtRaritaSchwinger.hh"
0027 #include "EvtGenBase/EvtDecayAmp.hh"
0028
0029 namespace Herwig {
0030
0031 using namespace ThePEG;
0032
0033 /**
0034  * The EvtGenInterface class is the main class for the use of the EvtGen decay
0035  * package with Herwig.
0036  *
0037  * @see \ref EvtGenInterfaceInterfaces "The interfaces"
0038  * defined for EvtGenInterface.
0039  */
0040 class EvtGenInterface: public Interfaced {
0041
0042 public:
0043
0044   /**
0045    * The default constructor.
0046    */
0047   EvtGenInterface();
0048
0049   /**
0050    * The copy constructor (explicit as cannot copy streams)
0051    */
0052   EvtGenInterface(const EvtGenInterface &);
0053
0054 public:
0055
0056   /**
0057    * Use EvtGen to perform a decay
0058    * @param parent The decaying particle
0059    * @param dm The decaymode
0060    * @return The decay products
0061    */
0062   ParticleVector decay(const Particle &parent,
0063                bool recursive, const DecayMode & dm) const;
0064
0065 public:
0066
0067   /** @name Functions used by the persistent I/O system. */
0068   //@{
0069   /**
0070    * Function used to write out object persistently.
0071    * @param os the persistent output stream written to.
0072    */
0073   void persistentOutput(PersistentOStream & os) const;
0074
0075   /**
0076    * Function used to read in object persistently.
0077    * @param is the persistent input stream read from.
0078    * @param version the version number of the object when written.
0079    */
0080   void persistentInput(PersistentIStream & is, int version);
0081   //@}
0082
0083   /**
0084    * The standard Init function used to initialize the interfaces.
0085    * Called exactly once for each class by the class description system
0086    * before the main function starts or
0087    * when this class is dynamically loaded.
0088    */
0089   static void Init();
0090
0091 protected:
0092
0093   /** @name Functions to convert between EvtGen and Herwig classes */
0094   //@{
0095   /**
0096    * Convert a particle to an EvtGen particle.
0097    * @param part The particle to be converted.
0098    */
0099   EvtParticle *EvtGenParticle(const Particle & part) const;
0100
0101   /**
0102    *  Return the decay products of an EvtGen particle in as ThePEG particles
0103    * @param evtpart The EvtGen particle
0104    */
0105   ParticleVector decayProducts(EvtParticle* evtpart, bool boost=true) const;
0106
0107   /**
0108    * Convert a Lorentz5Momentum to a real EvtGen 4-vector
0109    * @param mom The momentum to be converted
0110    */
0111   EvtVector4R EvtGenMomentum(const Lorentz5Momentum & mom) const {
0112     return EvtVector4R(mom.t()/GeV,mom.x()/GeV,mom.y()/GeV,mom.z()/GeV);
0113   }
0114
0115   /**
0116    * Convert a PDG code from ThePEG into an EvtGen particle id
0117    * @param id The PDG code
0118    * @param exception Whether or not to throw an Exception if fails
0119    */
0120   EvtId EvtGenID(int id,bool exception=true) const;
0121
0122   /**
0123    * Convert a LorentzSpinor to an EvtGen one. The spinor is converted to the
0124    * EvtGen Dirac representation/
0125    * @param sp The LorentzSpinor
0126    */
0127   EvtDiracSpinor EvtGenSpinor(const LorentzSpinor<SqrtEnergy> & sp) const {
0128     InvSqrtEnergy norm(sqrt(0.5)/sqrt(GeV));
0129     EvtDiracSpinor output;
0130     output.set(EvtGenComplex(-norm*( sp.s1()+sp.s3())),
0131            EvtGenComplex(-norm*( sp.s2()+sp.s4())),
0132            EvtGenComplex(-norm*(-sp.s1()+sp.s3())),
0133            EvtGenComplex(-norm*(-sp.s2()+sp.s4())));
0134     return output;
0135   }
0136
0137   /**
0138    * Convert a LorentzPolarizationVector to a complex EvtGen 4-vector
0139    * @param eps The polarization vector to be converted
0140    */
0141   EvtVector4C EvtGenPolarization(const LorentzPolarizationVector & eps) const {
0142     return EvtVector4C(EvtGenComplex(eps.t()),EvtGenComplex(eps.x()),
0143                EvtGenComplex(eps.y()),EvtGenComplex(eps.z()));
0144   }
0145
0146   /**
0147    * Convert our Rarita-Schwinger spinor to the EvtGen one
0148    * @param sp Our  RS Spinor
0149    */
0150   EvtRaritaSchwinger EvtGenRSSpinor(const LorentzRSSpinor<SqrtEnergy> & sp) const {
0151     InvSqrtEnergy norm(sqrt(0.5)/sqrt(GeV));
0152     complex<double> out[4][4];
0153     for(unsigned int ix=0;ix<4;++ix) {
0154       out[ix][0] = -Complex(norm*( sp(ix,0)+sp(ix,2)));
0155       out[ix][1] = -Complex(norm*( sp(ix,1)+sp(ix,3)));
0156       out[ix][2] = -Complex(norm*(-sp(ix,0)+sp(ix,2)));
0157       out[ix][3] = -Complex(norm*(-sp(ix,1)+sp(ix,3)));
0158     }
0159     EvtRaritaSchwinger output;
0160     unsigned int ix,iy;
0161     // remember we have vec,spin and evtgen spin,vec
0162     for(ix=0;ix<4;++ix) {
0163       for(iy=0;iy<4;++iy) output.set(ix,iy,EvtGenComplex(out[iy][ix]));
0164     }
0165     return output;
0166   }
0167
0168   /**
0169    * Convert our tensor to the EvtGen one.
0170    * @param ten Our tensor
0171    */
0172   EvtTensor4C EvtGenTensor(const LorentzTensor<double> & ten) const {
0173     EvtTensor4C output;
0174     unsigned int ix,iy;
0175     for(ix=0;ix<4;++ix){
0176       for(iy=0;iy<4;++iy) output.set(ix,iy,EvtGenComplex(ten(ix,iy)));
0177     }
0178     return output;
0179   }
0180
0181   /**
0182    * Convert a spin density matrix to an EvtGen spin density matrix.
0183    * @param rho The spin density matrix to be converted.
0184    */
0185   EvtSpinDensity EvtGenSpinDensity(const RhoDMatrix & rho) const {
0186     EvtSpinDensity rhoout;
0187     unsigned int ix,iy,ispin(rho.iSpin());
0188     rhoout.setDim(ispin);
0189     for(ix=0;ix<ispin;++ix) {
0190       for(iy=0;iy<ispin;++iy)
0191     rhoout.set(ix,iy,EvtGenComplex(rho(ix,iy)));
0192     }
0193     return rhoout;
0194   }
0195
0196   /**
0197    * Convert from our complex to the EvtGen one
0198    */
0199   EvtComplex EvtGenComplex(Complex c) const {
0200     return EvtComplex(c.real(),c.imag());
0201   }
0202   //@}
0203
0204   /**
0205    *  Functions to convert between EvtGen and Herwig classes
0206    */
0207   //@{
0208   /**
0209    * Convert a particle from an EvtGen one to ThePEG one.
0210    * @param part The EvtGen particle.
0211    * @param pd Pointer to the particle data object of ThePEG for the particle.
0212    * @param spin Convert the spin information as well
0213    */
0214   PPtr ThePEGParticle(EvtParticle *part, tcPDPtr pd,bool spin=true) const {
0215     PPtr output(new_ptr(Particle(pd)));
0216     output->set5Momentum(ThePEGMomentum(part->getP4(),part->mass()));
0217     if(spin) ThePEGSpin(output,part);
0218     // make the daughters
0219     ParticleVector daug(decayProducts(part,spin));
0220     for(unsigned int ix=0;ix<daug.size();++ix) output->addChild(daug[ix]);
0221     return output;
0222   }
0223
0224   /**
0225    * Set the SpinInfo for a ThePEG particle using an EvtGen particle
0226    * @param pegpart ThePEG particle.
0227    * @param evtpart The EvtGen particle.
0228    */
0229   void ThePEGSpin(PPtr pegpart,EvtParticle *evtpart) const;
0230
0231   /**
0232    * Convert an EvtGen EvtId to a PDG code in our conventions
0233    * @param id The EvtGen ID.
0234    * @param exception Whether or not to throw an Exception if fails
0235    */
0236   int ThePEGID(EvtId id,bool exception=true) const;
0237
0238   /**
0239    * Convert from EvtGen momentum to Lorentz5Momentum
0240    * @param mom The EvtGen 4-momentum
0241    * @param mass The mass
0242    */
0243   Lorentz5Momentum ThePEGMomentum(const EvtVector4R & mom,double mass) const  {
0244     return Lorentz5Momentum(mom.get(1)*GeV,mom.get(2)*GeV,
0245                 mom.get(3)*GeV,mom.get(0)*GeV,mass*GeV);
0246   }
0247   /**
0248    * Convert from EvtGen complex to ours
0249    */
0250   Complex ThePEGComplex(EvtComplex c) const {
0251     return Complex(real(c),imag(c));
0252   }
0253
0254   /**
0255    * Convert a spin density to a ThePEG one from an EvtGen one
0256    * @param rho The spin density matrix to be converted
0257    * @param id The PDG code of the particle to get special cases right.
0258    */
0259   RhoDMatrix ThePEGSpinDensity(const EvtSpinDensity & rho, int id) const;
0260
0261   /**
0262    * Convert an EvtDiracSpinor a LorentzSpinor. This spinor is converted to
0263    * the default Dirac matrix representation used by ThePEG.
0264    * @param sp The EvtDiracSpinor
0265    */
0266   LorentzSpinor<SqrtEnergy> ThePEGSpinor(const EvtDiracSpinor & sp) const {
0267     SqrtEnergy norm(sqrt(0.5)*sqrt(GeV));
0268     vector<complex<SqrtEnergy> > evtSpin(4);
0269     for(unsigned int ix=0;ix<4;++ix) evtSpin[ix] = -norm*ThePEGComplex(sp.get_spinor(ix));
0270     return LorentzSpinor<SqrtEnergy>(evtSpin[0]-evtSpin[2],evtSpin[1]-evtSpin[3],
0271                      evtSpin[0]+evtSpin[2],evtSpin[1]+evtSpin[3]);
0272   }
0273
0274   /**
0275    * Convert an EvtGen complex 4-vector to a LorentzPolarizationVector
0276    * @param eps The complex 4-vector to be converted.
0277    */
0278   LorentzPolarizationVector ThePEGPolarization(const EvtVector4C & eps) const {
0279     return LorentzPolarizationVector(conj(ThePEGComplex(eps.get(1))),
0280                      conj(ThePEGComplex(eps.get(2))),
0281                      conj(ThePEGComplex(eps.get(3))),
0282                      conj(ThePEGComplex(eps.get(0))));
0283   }
0284
0285   /**
0286    * Convert an EvtGen Rarita-Schwinger spinor to ours
0287    * @param sp The EvtGen RS spinor.
0288    */
0289   LorentzRSSpinor<SqrtEnergy> ThePEGRSSpinor(const EvtRaritaSchwinger & sp) const {
0290     complex<SqrtEnergy> evtSpin[4][4];
0291     SqrtEnergy norm(sqrt(0.5)*sqrt(GeV));
0292     // normalisation and swap vec,spin order
0293     for(unsigned int ix=0;ix<4;++ix) {
0294       for(unsigned int iy=0;iy<4;++iy) evtSpin[ix][iy]=-norm*ThePEGComplex(sp.get(iy,ix));
0295     }
0296     LorentzRSSpinor<SqrtEnergy> output;
0297     for(unsigned int ix=0;ix<4;++ix) {
0298       output(ix,0) = evtSpin[ix][0] - evtSpin[ix][2];
0299       output(ix,1) = evtSpin[ix][1] - evtSpin[ix][3];
0300       output(ix,2) = evtSpin[ix][0] + evtSpin[ix][2];
0301       output(ix,3) = evtSpin[ix][1] + evtSpin[ix][3];
0302     }
0303     // output.changeRep(Helicity::defaultDRep);
0304     return output;
0305   }
0306
0307   /**
0308    * Convert an EvtGen tensor to ThePEG
0309    * @param ten The EvtGen tensor
0310    */
0311   LorentzTensor<double> ThePEGTensor(const EvtTensor4C & ten) const {
0312     LorentzTensor<double> output;
0313     unsigned int ix,iy;
0314     for(ix=0;ix<4;++ix) {
0315       for(iy=0;iy<4;++iy)output(ix,iy)=conj(ThePEGComplex(ten.get(ix,iy)));
0316     }
0317     return output;
0318   }
0319   //@}
0320
0321   /**
0322    *  Check the conversion of particles between Herwig and EvtGen
0323    */
0324   void checkConversion() const;
0325
0326   /**
0327    * Output the EvtGen decay modes for a given particle
0328    * @param id The PDG code of the particle to output
0329    */
0330   void outputEvtGenDecays(long id) const;
0331
0332   /**
0333    *  Find the location in the EvtGen list of decay channels for
0334    *  a given decay mode.
0335    */
0336   int EvtGenChannel(const DecayMode &dm) const;
0337
0338   /**
0339    * Check the particle has SpinInfo and if not create it
0340    * @param part The particle
0341    */
0342   tSpinPtr getSpinInfo(const Particle &part) const {
0343     // return spin info if exists
0344     if(part.spinInfo()) {
0345       return dynamic_ptr_cast<tSpinPtr>(const_ptr_cast<tPPtr>(&part)->spinInfo());
0346     }
0347     // otherwise make it
0348     tPPtr ptemp(const_ptr_cast<tPPtr>(&part));
0349     PDT::Spin spin(part.dataPtr()->iSpin());
0350     SpinPtr pspin;
0351     if(spin==PDT::Spin0)          pspin=new_ptr(ScalarSpinInfo(part.momentum(),true));
0352     else if(spin==PDT::Spin1Half) pspin=new_ptr(FermionSpinInfo(part.momentum(),true));
0353     else if(spin==PDT::Spin1)     pspin=new_ptr(VectorSpinInfo(part.momentum(),true));
0354     else if(spin==PDT::Spin3Half) pspin=new_ptr(RSFermionSpinInfo(part.momentum(),true));
0355     else if(spin==PDT::Spin2)     pspin=new_ptr(TensorSpinInfo(part.momentum(),true));
0356     else throw Exception() << "Can't create spinInfo for decaying particle in "
0357                << "EvtGen::checkSpinInfo for spin " << spin << "particle "
0358                << Exception::eventerror;
0359     ptemp->spinInfo(pspin);
0360     return pspin;
0361   }
0362
0363   /**
0364    *  Construct the DecayVertex for Herwig using the information from
0365    *  EvtGen
0366    * @param parent The decaying particle
0367    * @param products The outgoing particles
0368    * @param damp Pointer to the EvtGen decayer
0369    */
0370   void constructVertex(const Particle & parent,ParticleVector products,
0371                EvtDecayAmp* damp) const;
0372
0373 protected:
0374
0375   /** @name Clone Methods. */
0376   //@{
0377   /**
0378    * Make a simple clone of this object.
0379    * @return a pointer to the new object.
0380    */
0381   virtual IBPtr clone() const;
0382
0383   /** Make a clone of this object, possibly modifying the cloned object
0384    * to make it sane.
0385    * @return a pointer to the new object.
0386    */
0387   virtual IBPtr fullclone() const;
0388   //@}
0389
0390 protected:
0391
0392   /** @name Standard Interfaced functions. */
0393   //@{
0394   /**
0395    * Initialize this object. Called in the run phase just before
0396    * a run begins.
0397    */
0398   virtual void doinitrun();
0399
0400   /**
0401    * Finalize this object. Called in the run phase just after a
0402    * run has ended. Used eg. to write out statistics.
0403    */
0404   virtual void dofinish();
0405   //@}
0406
0407 private:
0408
0409   /**
0410    * The assignment operator is private and must never be called.
0411    * In fact, it should not even be implemented.
0412    */
0413   EvtGenInterface & operator=(const EvtGenInterface &) = delete;
0414
0415 private:
0416
0417   /**
0418    *    Names of the various EvtGen parameter files
0419    */
0420   //@{
0421   /**
0422    *  The name of the file containing the decays
0423    */
0424   string decayName_;
0425
0426   /**
0427    *  The name of the file containing the particle data
0428    */
0429   string pdtName_;
0430
0431   /**
0432    *  Names of addition user specified decays
0433    */
0434   vector<string> userDecays_;
0435   //@}
0436
0437   /**
0438    *  Whether or not to redirect cout and cerr when EvtGen is running
0439    */
0440   bool reDirect_;
0441
0442   /**
0443    *  Check the conversion of the particles
0444    */
0445   bool checkConv_;
0446
0447   /**
0448    *  Particles for which to output the EvtGen decays
0449    */
0450   vector<long> convID_;
0451
0452   /**
0453    *  Location of the PYTHIA8 data directory
0454    */
0455   string p8Data_;
0456
0457 private:
0458
0459   /**
0460    * Pointer to the random number generator for EvtGen
0461    */
0462   EvtRandomEngine * evtrnd_;
0463
0464   /**
0465    * Main EvtGen object
0466    */
0467   EvtGen * evtgen_;
0468
0469   /**
0470    *  File to output the log info to
0471    */
0472   mutable ofstream logFile_;
0473
0474 };
0475
0476 }
0477
0478 #endif /* Herwig_EvtGenInterface_H */