Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-08-06 09:24:04

0001 // -*- C++ -*-
0002 //
0003 // FivePionCurrent.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef HERWIG_FivePionCurrent_H
0010 #define HERWIG_FivePionCurrent_H
0011 //
0012 // This is the declaration of the FivePionCurrent class.
0013 //
0014 
0015 #include "WeakCurrent.h"
0016 #include "ThePEG/Helicity/epsilon.h"
0017 
0018 namespace Herwig {
0019 
0020 using namespace ThePEG;
0021 
0022 /**
0023  * Here is the documentation of the FivePionCurrent class.
0024  *
0025  * @see \ref FivePionCurrentInterfaces "The interfaces"
0026  * defined for FivePionCurrent.
0027  */
0028 class FivePionCurrent: public WeakCurrent {
0029 
0030 public:
0031 
0032   /**
0033    * The default constructor.
0034    */
0035   FivePionCurrent();
0036 
0037   /** @name Methods for the construction of the phase space integrator. */
0038   //@{
0039   /**
0040    * Complete the construction of the decay mode for integration.classes inheriting
0041    * from this one.
0042    * This method is purely virtual and must be implemented in the classes inheriting
0043    * from WeakCurrent.
0044    * @param icharge   The total charge of the outgoing particles in the current.
0045    * @param resonance If specified only include terms with this particle
0046    * @param flavour Information on the required flavours of the quarks
0047    * @param imode     The mode in the current being asked for.
0048    * @param mode      The phase space mode for the integration
0049    * @param iloc      The location of the of the first particle from the current in
0050    *                  the list of outgoing particles.
0051    * @param ires      The location of the first intermediate for the current.
0052    * @param phase     The prototype phase space channel for the integration.
0053    * @param upp       The maximum possible mass the particles in the current are
0054    *                  allowed to have.
0055    * @return Whether the current was sucessfully constructed.
0056    */
0057   virtual bool createMode(int icharge, tcPDPtr resonance,
0058               FlavourInfo flavour,
0059               unsigned int imode,PhaseSpaceModePtr mode,
0060               unsigned int iloc,int ires,
0061               PhaseSpaceChannel phase, Energy upp );
0062 
0063   /**
0064    * The particles produced by the current.
0065    * @param icharge The total charge of the particles in the current.
0066    * @param imode The mode for which the particles are being requested
0067    * @param iq The PDG code for the quark
0068    * @param ia The PDG code for the antiquark
0069    * @return The external particles for the current.
0070    */
0071   virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
0072   //@}
0073 
0074   /**
0075    * Hadronic current. This method is purely virtual and must be implemented in
0076    * all classes inheriting from this one.
0077    * @param resonance If specified only include terms with this particle
0078    * @param flavour Information on the required flavours of the quarks
0079    * @param imode The mode
0080    * @param ichan The phase-space channel the current is needed for.
0081    * @param scale The invariant mass of the particles in the current.
0082    * @param outgoing The particles produced in the decay
0083    * @param momenta  The momenta of the particles produced in the decay
0084    * @param meopt Option for the calculation of the matrix element
0085    * @return The current. 
0086    */
0087   virtual vector<LorentzPolarizationVectorE> 
0088   current(tcPDPtr resonance,
0089       FlavourInfo flavour,
0090       const int imode, const int ichan,Energy & scale,
0091       const tPDVector & outgoing,
0092       const vector<Lorentz5Momentum> & momenta,
0093       DecayIntegrator::MEOption meopt) const;
0094 
0095   /**
0096    * Accept the decay. 
0097    * @param id The id's of the particles in the current.
0098    * @return Can this current have the external particles specified.
0099    */
0100   virtual bool accept(vector<int> id);
0101 
0102   /**
0103    * Return the decay mode number for a given set of particles in the current. 
0104    * @param id The id's of the particles in the current.
0105    * @return The number of the mode
0106    */
0107   virtual unsigned int decayMode(vector<int> id);
0108 
0109   /**
0110    * Output the setup information for the particle database
0111    * @param os The stream to output the information to
0112    * @param header Whether or not to output the information for MySQL
0113    * @param create Whether or not to add a statement creating the object
0114    */
0115   virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
0116 
0117 public:
0118 
0119   /** @name Functions used by the persistent I/O system. */
0120   //@{
0121   /**
0122    * Function used to write out object persistently.
0123    * @param os the persistent output stream written to.
0124    */
0125   void persistentOutput(PersistentOStream & os) const;
0126 
0127   /**
0128    * Function used to read in object persistently.
0129    * @param is the persistent input stream read from.
0130    * @param version the version number of the object when written.
0131    */
0132   void persistentInput(PersistentIStream & is, int version);
0133   //@}
0134 
0135   /**
0136    * The standard Init function used to initialize the interfaces.
0137    * Called exactly once for each class by the class description system
0138    * before the main function starts or
0139    * when this class is dynamically loaded.
0140    */
0141   static void Init();
0142 
0143 protected:
0144 
0145   /**
0146    *  Methods to calculate the Breit-Wigner distributions for the various
0147    *  mesons.
0148    */
0149   //@{
0150   /**
0151    * Breit-Wigner for the \f$\rho\f$.
0152    * @param scale The virtual mass
0153    */
0154   Complex rhoBreitWigner(Energy2 scale) const {
0155     Energy2 m2=sqr(_rhomass);
0156     return m2/(m2-scale-Complex(0.,1.)*_rhomass*_rhowidth);
0157   }
0158 
0159   /**
0160    * Breit-Wigner for the \f$a_1\f$.
0161    * @param scale The virtual mass
0162    */
0163   Complex a1BreitWigner(Energy2 scale) const {
0164     Energy2 m2=sqr(_a1mass);
0165     return m2/(m2-scale-Complex(0.,1.)*_a1mass*_a1width);
0166   }
0167   
0168   /**
0169    * Breit-Wigner for the \f$\omega\f$.
0170    * @param scale The virtual mass
0171    */
0172   Complex omegaBreitWigner(Energy2 scale) const {
0173     Energy2 m2=sqr(_omegamass);
0174     return m2/(m2-scale-Complex(0.,1.)*_omegamass*_omegawidth);
0175   }
0176 
0177   /**
0178    * Breit-Wigner for the \f$\sigma\f$.
0179    * @param scale The virtual mass
0180    */
0181   Complex sigmaBreitWigner(Energy2 scale) const {
0182     Energy2 m2=sqr(_sigmamass);
0183     return m2/(m2-scale-Complex(0.,1.)*_sigmamass*_sigmawidth);
0184   }
0185   //@}
0186 
0187   /**
0188    *  Currents for the different channels
0189    */
0190   //@{
0191   /**
0192    * The \f$\rho\omega\f$ current
0193    * @param iopt Option for the inclusion of \f$\rho\f$ Breit-Wigner terms in the 
0194    * \f$\omega\f$ decay piece
0195    * @param Q The total momentum for the current
0196    * @param q1 The first momentum
0197    * @param q2 The first momentum
0198    * @param q3 The first momentum
0199    * @param q4 The first momentum
0200    * @param q5 The first momentum
0201    */
0202   LorentzVector<complex<InvEnergy2> >
0203   rhoOmegaCurrent(unsigned int iopt,
0204           const Lorentz5Momentum & Q,
0205           const Lorentz5Momentum & q1,
0206           const Lorentz5Momentum & q2,
0207           const Lorentz5Momentum & q3,
0208           const Lorentz5Momentum & q4,
0209           const Lorentz5Momentum & q5) const {
0210     // prefactor
0211     complex<InvEnergy7> pre(_preomega*a1BreitWigner(Q.m2())*
0212                 omegaBreitWigner((q1+q2+q3).m2())*
0213                 rhoBreitWigner((q4+q5).m2()));
0214     // omega piece
0215     Complex omega(-1.);
0216     if(_rhoomega) {
0217       if(iopt==1)      omega=rhoBreitWigner((q2+q3).m2());
0218       else if(iopt==2) omega=rhoBreitWigner((q1+q3).m2());
0219       else if(iopt==3) omega=rhoBreitWigner((q1+q2).m2());
0220       else            
0221     omega=rhoBreitWigner((q2+q3).m2())+rhoBreitWigner((q1+q3).m2())+
0222       rhoBreitWigner((q1+q2).m2());
0223     }
0224     LorentzVector<complex<Energy3> > omegacurrent(Helicity::epsilon(q1,q2,q3));
0225     LorentzVector<complex<InvEnergy2> > output =
0226       pre * omega * Helicity::epsilon(q4-q5,omegacurrent,Q);
0227     return output;
0228   }
0229   
0230   /**
0231    *  The \f$a_1\sigma\f$ current
0232    * @param iopt Option for the inclusion of \f$\rho\f$ Breit-Wigner terms in the 
0233    * \f$a_1\f$ decay piece
0234    * @param Q The total momentum for the current
0235    * @param q1 The first momentum
0236    * @param q2 The first momentum
0237    * @param q3 The first momentum
0238    * @param q4 The first momentum
0239    * @param q5 The first momentum
0240    */
0241   LorentzVector<complex<InvEnergy2> > 
0242   a1SigmaCurrent(unsigned int iopt,
0243          const Lorentz5Momentum & Q,
0244          const Lorentz5Momentum & q1,
0245          const Lorentz5Momentum & q2,
0246          const Lorentz5Momentum & q3,
0247          const Lorentz5Momentum & q4,
0248          const Lorentz5Momentum & q5) const {
0249     Lorentz5Momentum pa1(q1+q2+q3);pa1.rescaleMass();
0250     Energy2 ma12(pa1.m2());
0251     complex<InvEnergy3> pre(_presigma*a1BreitWigner(Q.m2())*a1BreitWigner(ma12)*
0252                 sigmaBreitWigner((q4+q5).m2()));
0253     Energy2 pdot[2]={q2*(q1-q3),q1*(q2-q3)};
0254     LorentzPolarizationVectorE rho[2] = 
0255       {(pdot[0]/ma12*pa1-q1+q3)*rhoBreitWigner((q1+q3).m2()),
0256        (pdot[1]/ma12*pa1-q2+q3)*rhoBreitWigner((q2+q3).m2())};
0257     LorentzPolarizationVectorE total;
0258     if(iopt==1)      total = rho[0];
0259     else if(iopt==2) total = rho[1];
0260     else             total = rho[0]+rho[1];
0261     Complex qdot = total * Q / Q.m2();
0262     LorentzPolarizationVectorE cq(Q);
0263     cq = cq * qdot;
0264     cq -= total;
0265     return pre * cq;
0266   }
0267   //@}
0268 
0269 protected:
0270 
0271   /** @name Clone Methods. */
0272   //@{
0273   /**
0274    * Make a simple clone of this object.
0275    * @return a pointer to the new object.
0276    */
0277   virtual IBPtr clone() const {return new_ptr(*this);}
0278 
0279   /** Make a clone of this object, possibly modifying the cloned object
0280    * to make it sane.
0281    * @return a pointer to the new object.
0282    */
0283   virtual IBPtr fullclone() const {return new_ptr(*this);}
0284   //@}
0285 
0286 protected:
0287 
0288   /** @name Standard Interfaced functions. */
0289   //@{
0290   /**
0291    * Initialize this object after the setup phase before saving an
0292    * EventGenerator to disk.
0293    * @throws InitException if object could not be initialized properly.
0294    */
0295   virtual void doinit();
0296   //@}
0297 
0298 private:
0299 
0300   /**
0301    * The assignment operator is private and must never be called.
0302    * In fact, it should not even be implemented.
0303    */
0304   FivePionCurrent & operator=(const FivePionCurrent &) = delete;
0305 
0306 private:
0307 
0308   /**
0309    * The masses and widths of the intermediate particles 
0310    */
0311   //@{
0312   /**
0313    * The mass of the \f$\rho\f$ for the current.
0314    */
0315   Energy _rhomass;
0316 
0317   /**
0318    * The mass of the \f$a_1\f$ for the current.
0319    */
0320   Energy _a1mass;
0321 
0322   /**
0323    * The mass of the \f$\omega\f$ for the current.
0324    */
0325   Energy _omegamass;
0326 
0327   /**
0328    * The mass of the \f$\sigma\f$ for the current.
0329    */
0330   Energy _sigmamass;
0331 
0332   /**
0333    * The width for the \f$\rho\f$.
0334    */
0335   Energy _rhowidth;
0336 
0337   /**
0338    *  The \f$a_1\f$ width
0339    */
0340   Energy _a1width;
0341 
0342   /**
0343    *  The \f$\omega\f$ width.
0344    */
0345   Energy _omegawidth;
0346 
0347   /**
0348    *  The \f$\sigma\f$ width.
0349    */
0350   Energy _sigmawidth;
0351   //@}
0352 
0353   /**
0354    * use local values of the particle masses
0355    */
0356   bool _localparameters;
0357 
0358   /**
0359    *  Option for the treatment of \f$\rho\f$ Breit-Wigners in \f$\omega\f$ decay
0360    */
0361   bool _rhoomega;
0362 
0363   /**
0364    *  Normalisation parameters for the different currents
0365    */
0366   //@{
0367   /**
0368    *  The \f$c\f$ parameter
0369    */
0370   Energy2 _c;
0371 
0372   /**
0373    *  The \f$c_0\f$ parameter
0374    */
0375   double _c0;
0376 
0377   /**
0378    *  The \f$f_{\omega\rho\pi}\f$ parameter
0379    */
0380   InvEnergy _fomegarhopi;
0381 
0382   /**
0383    * The \f$g_{\rho\pi\pi}\f$ parameter
0384    */
0385   double _grhopipi;
0386 
0387   /**
0388    * The \f$G_{a\rho\pi}\f$ parameter
0389    */
0390   Energy _garhopi;
0391 
0392   /**
0393    *  The \f$f_{aaf}\f$ parameter
0394    */
0395   Energy _faaf;
0396 
0397   /**
0398    *  The \f$f_{f\pi\pi}\f$ parameter
0399    */
0400   Energy _ffpipi;
0401   //@}
0402 
0403   /**
0404    *  Values cached to avoid unnessacary calculations
0405    */
0406   //@{
0407   /**
0408    *  Prefactor for the \f$\rho\omega\f$ current
0409    */
0410   InvEnergy7 _preomega;
0411 
0412   /**
0413    *  Prefactor for the \f$a_1\sigma\f$ current
0414    */
0415   InvEnergy3 _presigma;
0416   //@}
0417 };
0418 
0419 }
0420 
0421 #endif /* HERWIG_FivePionCurrent_H */