KPiCurrent.h

// -*- C++ -*-
//
// KPiCurrent.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_KPiCurrent_H
#define HERWIG_KPiCurrent_H
//
// This is the declaration of the KPiCurrent class.
//
 
#include "WeakCurrent.h"
#include "Herwig/Utilities/Kinematics.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
class KPiCurrent: public WeakCurrent {
 
public:
 
  KPiCurrent();
 
 
  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 sWaveBreitWigner(Energy2 q2,unsigned int ires) const {
    Energy q=sqrt(q2),gam(ZERO);
    Energy2 m2=sqr(_scamass[ires]);
    if(q>_mK+_mpi) {
      Energy pX=Kinematics::pstarTwoBodyDecay(_scamass[ires],_mK,_mpi);
      Energy p =Kinematics::pstarTwoBodyDecay( q            ,_mK,_mpi);
      gam = _scawidth[ires]*m2/q2*p/pX;
    }
    return m2/(m2-q2-Complex(0.,1.)*q*gam);
  }
  
  Complex pWaveBreitWigner(Energy2 q2,unsigned int ires) const {
    Energy q=sqrt(q2),gam(ZERO);
    Energy2 m2=sqr(_vecmass[ires]);
    if(q>_mK+_mpi) {
      Energy pX=Kinematics::pstarTwoBodyDecay(_vecmass[ires],_mK,_mpi);
      Energy p =Kinematics::pstarTwoBodyDecay( q            ,_mK,_mpi);
      double ratio=p/pX;
      gam = _vecwidth[ires]*m2/q2*ratio*sqr(ratio);
    }
    return m2/(m2-q2-Complex(0.,1.)*q*gam);
  }
 
protected:
 
  virtual IBPtr clone() const {return new_ptr(*this);}
  
  virtual IBPtr fullclone() const {return new_ptr(*this);}
 
protected:
 
 
  virtual void doinit();
 
private:
 
  KPiCurrent & operator=(const KPiCurrent &) = delete;
 
private:
 
  bool _localparameters;
 
  bool _transverse;
 
  double _cV;
 
  double _cS;
 
  vector<double> _vecmag;
 
  vector<double> _vecphase;
 
  vector<Complex> _vecwgt;
 
  vector<Energy> _vecmass;
 
  vector<Energy> _vecwidth;
 
  vector<double> _scamag;
 
  vector<double> _scaphase;
 
  vector<Complex> _scawgt;
 
  vector<Energy> _scamass;
 
  vector<Energy> _scawidth;
 
  Energy _mpi;
 
  Energy _mK;
 
  vector<int> _resmap;
};
 
}
 
#endif /* HERWIG_KPiCurrent_H */