TwoKaonOnePionCurrent.h

// -*- C++ -*-
//
// TwoKaonOnePionCurrent.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_TwoKaonOnePionCurrent_H
#define HERWIG_TwoKaonOnePionCurrent_H
//
// This is the declaration of the TwoKaonOnePionCurrent class.
//
 
#include "WeakCurrent.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
class TwoKaonOnePionCurrent: public WeakCurrent {
 
public:
 
  TwoKaonOnePionCurrent();
 
  virtual bool createMode(int icharge, tcPDPtr resonance,
                          FlavourInfo flavour,
                          unsigned int imode,PhaseSpaceModePtr mode,
                          unsigned int iloc,int ires,
                          PhaseSpaceChannel phase, Energy upp );
 
 
  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 tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
  
  virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
  
  double threeBodyMatrixElement(const int imode,  const Energy2 q2,
                                const Energy2 s3, const Energy2 s2, 
                                const Energy2 s1, const Energy  m1, 
                                const Energy  m2, const Energy  m3) const;
  
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
protected:
 
  virtual void doinit();
 
  virtual void doinitrun();
 
  virtual void doupdate();
 
private:
 
  TwoKaonOnePionCurrent & operator=(const TwoKaonOnePionCurrent &) = delete;
 
private:
 
  Complex Trho1(Energy2 q2,int ires) const {
    if(ires>=int(_rho1wgts.size())) return 0.;
    double norm = std::accumulate(_rho1wgts.begin(),_rho1wgts.end(),0.);
    unsigned int imin=0,imax=_rho1wgts.size();
    if(ires>0) {
      imin=ires;
      imax=imin+1;
    }
    Complex output(0.);
    for(unsigned int ix=imin;ix<imax;++ix)
      output+=_rho1wgts[ix]*
        Resonance::BreitWignerPWave(q2,_rho1mass[ix],_rho1width[ix],_mpi,_mpi);
    return output/norm;
  }
 
  Complex Trho2(Energy2 q2,int ires) const {
    if(ires>=int(_rho2wgts.size())) return 0.;
    double norm = std::accumulate(_rho2wgts.begin(),_rho2wgts.end(),0.);
    unsigned int imin=0,imax=_rho2wgts.size();
    if(ires>0) {
      imin=ires;
      imax=imin+1;
    }
    Complex output(0.);
    for(unsigned int ix=imin;ix<imax;++ix)
      output+=_rho2wgts[ix]*
        Resonance::BreitWignerPWave(q2,_rho2mass[ix],_rho2width[ix],_mpi,_mpi);
    return output/norm;
  }
 
  Complex TKstar1(Energy2 q2,int ires) const {
    if(ires>=int(_kstar1wgts.size())) return 0.;
    double norm = std::accumulate(_kstar1wgts.begin(),_kstar1wgts.end(),0.);
    unsigned int imin=0,imax=_kstar1wgts.size();
    if(ires>0) {
      imin=ires;
      imax=imin+1;
    }
    Complex output(0.);
    for(unsigned int ix=imin;ix<imax;++ix)
      output+=_kstar1wgts[ix]*
        Resonance::BreitWignerPWave(q2,_kstar1mass[ix],_kstar1width[ix],_mK,_mpi);
    return output/norm;
  }
 
  Complex a1BreitWigner(Energy2 q2) const {
    Complex ii(0.,1.);
    Energy2 m2(_a1mass*_a1mass);
    Energy  q(sqrt(q2));
    Energy gamma = !_a1opt ?
      _a1mass*_a1width*Resonance::ga1(q2)/Resonance::ga1(_a1mass*_a1mass)/sqrt(q2) : (*_a1runinter)(q2);
    return m2/(m2-q2-ii*q*gamma);
  }
 
  void inita1Width(int iopt);
 
  Complex Tomega(Energy2 q2, int ires) const;
 
  Complex OmegaPhiBreitWigner(Energy2 q2, unsigned int ires) const {
    Energy2 m2,mg;
    if(ires==0) {
      m2=sqr(_omegamass);
      mg=_omegamass*_omegawidth;
    }
    else {
      m2=sqr(_phimass);
      mg=_phimass*_phiwidth;
    }
    return (-m2+Complex(0.,1.)*mg)/(q2-m2+Complex(0.,1.)*mg);
  }
 
  Complex TOmegaKStar(Energy2 s1,Energy2 s2,int ires) const;
 
private:
 
  vector<double> _rho1wgts;
 
  vector<Energy> _rho1mass;
 
  vector<Energy> _rho1width;
 
  vector<double> _rho2wgts;
 
  vector<Energy> _rho2mass;
 
  vector<Energy> _rho2width;
 
  vector<double> _kstar1wgts;
 
  vector<Energy> _kstar1mass;
 
  vector<Energy> _kstar1width;
 
  Energy _a1mass;
  
  Energy _a1width;
 
  vector<Energy>  _a1runwidth;
 
  vector<Energy2> _a1runq2;
 
  Interpolator<Energy,Energy2>::Ptr _a1runinter;
 
  double _epsomega;
 
  Energy _omegamass;
 
  Energy _omegawidth;
 
  Energy _phimass;
 
  Energy _phiwidth;
 
  double _omegaKstarwgt;
 
  Energy _fpi;
 
  Energy _mpi;
 
  Energy _mK;
 
  bool _initializea1;
 
  bool _a1opt;
 
  Energy _maxmass;
 
  Energy _maxcalc;
};
 
}
 
#endif /* HERWIG_TwoKaonOnePionCurrent_H */