FivePionCurrent.h

// -*- C++ -*-
//
// FivePionCurrent.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_FivePionCurrent_H
#define HERWIG_FivePionCurrent_H
//
// This is the declaration of the FivePionCurrent class.
//
 
#include "WeakCurrent.h"
#include "ThePEG/Helicity/epsilon.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
class FivePionCurrent: public WeakCurrent {
 
public:
 
  FivePionCurrent();
 
  virtual bool createMode(int icharge, tcPDPtr resonance,
                          FlavourInfo flavour,
                          unsigned int imode,PhaseSpaceModePtr mode,
                          unsigned int iloc,int ires,
                          PhaseSpaceChannel phase, Energy upp );
 
  virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
 
  virtual vector<LorentzPolarizationVectorE> 
  current(tcPDPtr resonance,
          FlavourInfo flavour,
          const int imode, const int ichan,Energy & scale,
          const tPDVector & outgoing,
          const vector<Lorentz5Momentum> & momenta,
          DecayIntegrator::MEOption meopt) const;
 
  virtual bool accept(vector<int> id);
 
  virtual unsigned int decayMode(vector<int> id);
 
  virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  Complex rhoBreitWigner(Energy2 scale) const {
    Energy2 m2=sqr(_rhomass);
    return m2/(m2-scale-Complex(0.,1.)*_rhomass*_rhowidth);
  }
 
  Complex a1BreitWigner(Energy2 scale) const {
    Energy2 m2=sqr(_a1mass);
    return m2/(m2-scale-Complex(0.,1.)*_a1mass*_a1width);
  }
  
  Complex omegaBreitWigner(Energy2 scale) const {
    Energy2 m2=sqr(_omegamass);
    return m2/(m2-scale-Complex(0.,1.)*_omegamass*_omegawidth);
  }
 
  Complex sigmaBreitWigner(Energy2 scale) const {
    Energy2 m2=sqr(_sigmamass);
    return m2/(m2-scale-Complex(0.,1.)*_sigmamass*_sigmawidth);
  }
 
  LorentzVector<complex<InvEnergy2> >
  rhoOmegaCurrent(unsigned int iopt,
                  const Lorentz5Momentum & Q,
                  const Lorentz5Momentum & q1,
                  const Lorentz5Momentum & q2,
                  const Lorentz5Momentum & q3,
                  const Lorentz5Momentum & q4,
                  const Lorentz5Momentum & q5) const {
    // prefactor
    complex<InvEnergy7> pre(_preomega*a1BreitWigner(Q.m2())*
                            omegaBreitWigner((q1+q2+q3).m2())*
                            rhoBreitWigner((q4+q5).m2()));
    // omega piece
    Complex omega(-1.);
    if(_rhoomega) {
      if(iopt==1)      omega=rhoBreitWigner((q2+q3).m2());
      else if(iopt==2) omega=rhoBreitWigner((q1+q3).m2());
      else if(iopt==3) omega=rhoBreitWigner((q1+q2).m2());
      else            
        omega=rhoBreitWigner((q2+q3).m2())+rhoBreitWigner((q1+q3).m2())+
          rhoBreitWigner((q1+q2).m2());
    }
    LorentzVector<complex<Energy3> > omegacurrent(Helicity::epsilon(q1,q2,q3));
    LorentzVector<complex<InvEnergy2> > output =
      pre * omega * Helicity::epsilon(q4-q5,omegacurrent,Q);
    return output;
  }
  
  LorentzVector<complex<InvEnergy2> > 
  a1SigmaCurrent(unsigned int iopt,
                 const Lorentz5Momentum & Q,
                 const Lorentz5Momentum & q1,
                 const Lorentz5Momentum & q2,
                 const Lorentz5Momentum & q3,
                 const Lorentz5Momentum & q4,
                 const Lorentz5Momentum & q5) const {
    Lorentz5Momentum pa1(q1+q2+q3);pa1.rescaleMass();
    Energy2 ma12(pa1.m2());
    complex<InvEnergy3> pre(_presigma*a1BreitWigner(Q.m2())*a1BreitWigner(ma12)*
                            sigmaBreitWigner((q4+q5).m2()));
    Energy2 pdot[2]={q2*(q1-q3),q1*(q2-q3)};
    LorentzPolarizationVectorE rho[2] = 
      {(pdot[0]/ma12*pa1-q1+q3)*rhoBreitWigner((q1+q3).m2()),
       (pdot[1]/ma12*pa1-q2+q3)*rhoBreitWigner((q2+q3).m2())};
    LorentzPolarizationVectorE total;
    if(iopt==1)      total = rho[0];
    else if(iopt==2) total = rho[1];
    else             total = rho[0]+rho[1];
    Complex qdot = total * Q / Q.m2();
    LorentzPolarizationVectorE cq(Q);
    cq = cq * qdot;
    cq -= total;
    return pre * cq;
  }
 
protected:
 
  virtual IBPtr clone() const {return new_ptr(*this);}
 
  virtual IBPtr fullclone() const {return new_ptr(*this);}
 
protected:
 
  virtual void doinit();
 
private:
 
  FivePionCurrent & operator=(const FivePionCurrent &) = delete;
 
private:
 
  Energy _rhomass;
 
  Energy _a1mass;
 
  Energy _omegamass;
 
  Energy _sigmamass;
 
  Energy _rhowidth;
 
  Energy _a1width;
 
  Energy _omegawidth;
 
  Energy _sigmawidth;
 
  bool _localparameters;
 
  bool _rhoomega;
 
  Energy2 _c;
 
  double _c0;
 
  InvEnergy _fomegarhopi;
 
  double _grhopipi;
 
  Energy _garhopi;
 
  Energy _faaf;
 
  Energy _ffpipi;
 
  InvEnergy7 _preomega;
 
  InvEnergy3 _presigma;
};
 
}
 
#endif /* HERWIG_FivePionCurrent_H */