OneKaonTwoPionDefaultCurrent.h

// -*- C++ -*-
//
// OneKaonTwoPionDefaultCurrent.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_OneKaonTwoPionDefaultCurrent_H
#define HERWIG_OneKaonTwoPionDefaultCurrent_H
//
// This is the declaration of the OneKaonTwoPionDefaultCurrent class.
//
#include "WeakCurrent.h"
#include "Herwig/Utilities/Interpolator.h"
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "Herwig/Decay/ResonanceHelpers.h"
#include <numeric>
 
namespace Herwig {
using namespace ThePEG;
 
class OneKaonTwoPionDefaultCurrent: public WeakCurrent {
 
public:
 
  OneKaonTwoPionDefaultCurrent();
 
  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);
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
public:
 
  virtual bool createMode(int icharge, tcPDPtr resonance,
                          FlavourInfo flavour,
                          unsigned int imode,PhaseSpaceModePtr mode,
                          unsigned int iloc,int ires,
                          PhaseSpaceChannel phase, Energy upp );
 
  virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
  
protected:
 
  virtual IBPtr clone() const {return new_ptr(*this);}
 
  virtual IBPtr fullclone() const {return new_ptr(*this);}
 
protected:
 
  virtual void doinit();
 
private:
 
  OneKaonTwoPionDefaultCurrent & operator=(const OneKaonTwoPionDefaultCurrent &) = delete;
 
private:
  
  Complex BrhoF123(Energy2 q2,int ires) const {
    if(ires>=int(_rhoF123wgts.size())) return 0.;
    Complex output(0.);
    Complex norm = std::accumulate(_rhoF123wgts.begin(),
                                   _rhoF123wgts.end(),Complex(0.));
    unsigned int imin=0,imax=_rhoF123wgts.size();
    if(ires>0) {
      imin=ires;
      imax=imin+1;
    }
    for(unsigned int ix=imin;ix<imax;++ix)
      output+=_rhoF123wgts[ix]*Resonance::BreitWignerPWave(q2,_rhoF123masses[ix],
                                                           _rhoF123widths[ix],_mpi,_mpi);
    return output/norm;
  }
 
  Complex BKstarF123(Energy2 q2,int ires) const {
    if(ires>=int(_kstarF123wgts.size())) return 0.;
    Complex output(0.);
    Complex norm = std::accumulate(_kstarF123wgts.begin(),
                                   _kstarF123wgts.end(),Complex(0.));
    unsigned int imin=0,imax=_kstarF123wgts.size();
    if(ires>0) {
      imin=ires;
      imax=imin+1;
    }
    assert(imax<=_kstarF123wgts.size());
    for(unsigned int ix=imin;ix<imax;++ix)
      output+=_kstarF123wgts[ix]*Resonance::BreitWignerPWave(q2,_kstarF123masses[ix],
                                                             _kstarF123widths[ix],_mpi,_mK);
    return output/norm;
  }
  
  Complex FKrho(Energy2 si,Energy2 sj,int ires) const {
    Complex output;
    if(ires<0)
      output = _rhoKstarwgt*BKstarF123(si,-1)+BrhoF123(sj,-1);
    else if(ires%2==0)
      output= _rhoKstarwgt*BKstarF123(si,ires/2);
    else if(ires%2==1)
      output=BrhoF123(sj,ires/2);
    return output/(1.+_rhoKstarwgt);
  }
 
private:
  
  vector<double> _rhoF123wgts;
 
  vector<double> _kstarF123wgts;
 
  vector<double> _kstarF5wgts;
  
  double _rhoKstarwgt;
 
  Energy _k1mass;
 
  Energy _k1width;
 
  Energy _fpi;
 
  Energy _mpi;
 
  Energy _mK;
 
  vector<Energy> _rhoF123masses;
 
  vector<Energy> _rhoF123widths;
 
  vector<Energy> _kstarF123masses;
 
  vector<Energy> _kstarF5masses;
 
  vector<Energy> _kstarF123widths;
 
  vector<Energy> _kstarF5widths;
};
 
}
 
#endif /* HERWIG_OneKaonTwoPionDefaultCurrent_H */