ThreeMesonDefaultCurrent.h

// -*- C++ -*-
//
// ThreeMesonDefaultCurrent.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_ThreeMesonDefaultCurrent_H
#define HERWIG_ThreeMesonDefaultCurrent_H
//
// This is the declaration of the ThreeMesonDefaultCurrent class.
//
#include "ThreeMesonCurrentBase.h"
#include "Herwig++/Utilities/Interpolator.h"
#include "Herwig++/Utilities/Kinematics.h"
#include "ThePEG/StandardModel/StandardModelBase.h"

namespace Herwig {
using namespace ThePEG;

class ThreeMesonDefaultCurrent: public ThreeMesonCurrentBase {

  friend class Defaulta1MatrixElement;

public:

  ThreeMesonDefaultCurrent();

  void persistentOutput(PersistentOStream & os) const;

  void persistentInput(PersistentIStream & is, int version);

  static void Init();

public:

  
  virtual bool createMode(int icharge,unsigned int imode,DecayPhaseSpaceModePtr mode,
                          unsigned int iloc,unsigned int ires,
                          DecayPhaseSpaceChannelPtr phase,Energy upp);

  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;

protected:

  virtual bool acceptMode(int) const;

  virtual FormFactors calculateFormFactors(const int ichan, const int imode,
                                           Energy2 q2,
                                           Energy2 s1, Energy2 s2, Energy2 s3) const;

protected:

  virtual IBPtr clone() const {return new_ptr(*this);}

  virtual IBPtr fullclone() const {return new_ptr(*this);}

protected:

  virtual void doinit();

  virtual void doinitrun();

  virtual void doupdate();

private:

  ThreeMesonDefaultCurrent & operator=(const ThreeMesonDefaultCurrent &);

private:
  
  Complex BrhoF123(Energy2 q2,int ires) const {
    Complex output(0.),norm(0.);
    for(unsigned int ix=0,N=min(3,int(_rhoF123wgts.size()));ix<N;++ix) {
      norm+=_rhoF123wgts[ix];
    }
    if(ires<0) {
      for(unsigned int ix=0,N=min(3,int(_rhoF123wgts.size()));ix<N;++ix) {
        output+=_rhoF123wgts[ix]*rhoKBreitWigner(q2,0,ix);
      }
    }
    else {
      unsigned int temp(ires);
      if(temp<_rhoF123wgts.size()&&temp<3)
        output=_rhoF123wgts[temp]*rhoKBreitWigner(q2,0,temp);
      else
        output=0.;
    }
    return output/norm;
  }

  Complex BrhoF5(Energy2 q2,int ires) const {
    Complex output(0.),norm(0.);
    for(unsigned int ix=0,N=min(3,int(_rhoF5wgts.size()));ix<N;++ix) {
      norm+=_rhoF5wgts[ix];
    }
    if(ires<0) {
      for(unsigned int ix=0,N=min(3,int(_rhoF5wgts.size()));ix<N;++ix) {
        output+=_rhoF5wgts[ix]*rhoKBreitWigner(q2,1,ix);
      }
    }
    else {
      unsigned int temp(ires);
      if(temp<_rhoF5wgts.size()&&temp<3) {
        output=_rhoF5wgts[temp]*rhoKBreitWigner(q2,1,temp);
      }
    }
    return output/norm;
  }

  Complex BKstarF123(Energy2 q2,int ires) const {
    Complex output(0.),norm(0.);
    for(unsigned int ix=0,N=min(3,int(_kstarF123wgts.size()));ix<N;++ix) {
      norm+=_kstarF123wgts[ix];
    }
    if(ires<0) {
      for(unsigned int ix=0,N=min(3,int(_kstarF123wgts.size()));ix<N;++ix) {
        output+=_kstarF123wgts[ix]*rhoKBreitWigner(q2,2,ix);
      }
    }
    else {
      unsigned int temp(ires);
      if(temp<_kstarF123wgts.size()&&temp<3) {
        output=_kstarF123wgts[temp]*rhoKBreitWigner(q2,2,temp);
      }
    }
    return output/norm;
  }

  Complex BKstarF5(Energy2 q2,int ires) const {
    Complex output(0.),norm(0.);
    for(unsigned int ix=0,N=min(3,int(_kstarF5wgts.size()));ix<N;++ix) {
      norm+=_kstarF5wgts[ix];
    }
    if(ires<0) {
      for(unsigned int ix=0,N=min(3,int(_kstarF5wgts.size()));ix<N;++ix) {
        output+=_kstarF5wgts[ix]*rhoKBreitWigner(q2,3,ix);
      }
    }
    else {
      unsigned int temp(ires);
      if(temp<_kstarF5wgts.size()&&temp<3) {
        output=_kstarF5wgts[ires]*rhoKBreitWigner(q2,3,temp);
      }
    }
    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);}
    output /=(1.+_rhoKstarwgt);
    return output;
  }
  
  Complex a1BreitWigner(Energy2 q2) const  {
    Complex ii(0.,1.);
    Energy2 m2(_a1mass*_a1mass);
    Energy  q(sqrt(q2));
    return m2/(m2-q2-ii*q*a1Width(q2));
  }
  
  Complex K1BreitWigner(Energy2 q2) const {
    Energy2 m2 = sqr(_k1mass);
    Complex ii(0.,1.);
    complex<Energy2> fact(m2 - ii*_k1mass*_k1width);
    return fact/(fact-q2);
  }
  
  Energy a1Width(Energy2 q2) const {
    Energy output;
    if(!_a1opt) output = _a1mass*_a1width*g(q2)/g(sqr(_a1mass))/sqrt(q2);
    else        output = (*_a1runinter)(q2);
    return output;
  }
  
  double g(Energy2 q2) const {
    double output;
    if(q2 < 9.*sqr(_mpi)) {
      output=0.;
    }
    else if(q2 < sqr(_rhoF123masses[0]+_mpi)) {
      double diff = (q2-9.*sqr(_mpi))/GeV2;
      
      output = 4.1*sqr(diff)*diff*(1.-3.3*diff+5.8*sqr(diff));
    }
    else {
      double ratio = q2/GeV2;
      output = ratio*(1.623+10.38/ratio-9.32/sqr(ratio)+0.65/(ratio*sqr(ratio)));
    }
    return output;
  }

  void inita1Width(int iopt);

  Complex rhoKBreitWigner(Energy2 q2,unsigned int itype,unsigned int ires) const;

private:
  
  vector<double> _rhoF123wgts;

  vector<double> _kstarF123wgts;
  
  vector<double> _rhoF5wgts;

  vector<double> _kstarF5wgts;
  
  double _rhoKstarwgt;
  
  vector<Energy>  _a1runwidth;

  vector<Energy2> _a1runq2;


  Interpolator<Energy,Energy2>::Ptr _a1runinter;

  bool _initializea1;
  
  Energy _a1mass;

  Energy _a1width;

  Energy _k1mass;

  Energy _k1width;

  Energy _fpi;

  Energy _mpi;

  Energy _mK;

  bool _rhoparameters;

  vector<Energy> _rhoF123masses;

  vector<Energy> _rhoF5masses;

  vector<Energy> _rhoF123widths;

  vector<Energy> _rhoF5widths;
  
  bool _kstarparameters;

  vector<Energy> _kstarF123masses;

  vector<Energy> _kstarF5masses;

  vector<Energy> _kstarF123widths;

  vector<Energy> _kstarF5widths;
  
  bool _a1parameters;
  
  bool _k1parameters;

  bool _a1opt;

  Energy _maxmass;

  Energy _maxcalc;
  
};

}

#endif /* THEPEG_ThreeMesonDefaultCurrent_H */