TwoMesonRhoKStarCurrent.h

// -*- C++ -*-
//
// TwoMesonRhoKStarCurrent.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_TwoMesonRhoKStarCurrent_H
#define HERWIG_TwoMesonRhoKStarCurrent_H
// This is the declaration of the TwoMesonRhoKStarCurrent class.

#include "WeakDecayCurrent.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "Herwig++/Utilities/Kinematics.h"
#include "ThePEG/StandardModel/StandardModelBase.h"

namespace Herwig {
using namespace ThePEG;

class TwoMesonRhoKStarCurrent: public WeakDecayCurrent {

public:

  TwoMesonRhoKStarCurrent();

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

  virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);

  virtual vector<LorentzPolarizationVectorE> 
  current(const int imode, const int ichan,Energy & scale,  
          const ParticleVector & decay, 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:

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

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

protected:


  virtual void doinit();

private:

  static ClassDescription<TwoMesonRhoKStarCurrent> initTwoMesonRhoKStarCurrent;

  TwoMesonRhoKStarCurrent & operator=(const TwoMesonRhoKStarCurrent &);

private:

  Complex BreitWigner(Energy2 q2, unsigned int imodel,unsigned int itype,
                      unsigned int ires) const {
    // workout the index of the resonace
    unsigned int iy(3*itype+ires);
    // off shell mass
    Energy q(sqrt(q2));
    // and the running width
    Energy pq(pcm(iy,q));
    double ratio(pow<3,1>(pq/_mom[iy]));
    Energy gamrun(_width[iy]*_mass[iy]*ratio/q);
    // work out the denominator
    complex<Energy2> denom, numer;
    Complex ii(0.,1.);
    if(imodel==0) {
      denom=q2-_mass2[iy]+ii*_mass[iy]*gamrun;
      numer=-_mass2[iy];
    }
    else if(imodel==1) {
      denom = q2 - _mass2[iy] + ii*_mass[iy]*gamrun
        - ( sqr(pq) * (GSModelhFunction(iy,q)-_hm2[iy])
            - sqr(_mom[iy]) * (q2-_mass2[iy]) * _dhdq2[iy]);
      
      numer = -(_mass2[iy]+_dparam[iy]*_mass[iy]*_width[iy]);
    }
    else assert(false);
    // return the answer
    return numer/denom;
  }
  
  double GSModelDParameter(const unsigned int ires)const {
    Energy mpi(0.5*(_massa[ires]+_massb[ires]));
    using Constants::pi;
    return 3.*mpi*mpi/pi/sqr(_mom[ires])*log(0.5*(_mass[ires]+2.*_mom[ires])/mpi)
      +0.5*_mass[ires]/pi/_mom[ires]-mpi*mpi*_mass[ires]/pi/pow<3,1>(_mom[ires]);
  }
  
  InvEnergy2 GSModeldhdq2Parameter(const unsigned int ires)const  {
    Energy mpi = 0.5 * (_massa[ires] + _massb[ires]);
    return _width[ires] / Constants::pi / pow<3,1>(_mom[ires]) *
      (sqr(mpi) / _mass[ires] / _mom[ires] * 
       log(0.5 * (_mass[ires] + 2.*_mom[ires])/mpi) + 0.5);
  }
  
  double GSModelhFunction(const unsigned int ires, const Energy q) const  {
    Energy pq(pcm(ires,q));
    return _width[ires]*_mass2[ires]/pow<3,1>(_mom[ires])
      *2.*pq/Constants::pi/q*log((q+2.*pq)/(_massa[ires]+_massb[ires]));
  }
  
  Energy pcm(const unsigned int ires, const Energy q) const {
    Energy2 q2(q*q);
    if(q <= _massa[ires]+_massb[ires]) return ZERO;
    return 0.5/q*sqrt((q2-sqr(_massa[ires]+_massb[ires]))*
                      (q2-sqr(_massa[ires]-_massb[ires])));
  }

private:

  vector<Complex> _piwgt;

  vector<double> _pimag;

  vector<double> _piphase;

  vector<Complex> _kwgt;

  vector<double> _kmag;

  vector<double> _kphase;

  int _pimodel;

  int _kmodel;

  bool _rhoparameters;

  vector<Energy> _rhomasses;

  vector<Energy> _rhowidths;
  
  bool _kstarparameters;

  vector<Energy> _kstarmasses;

  vector<Energy> _kstarwidths;

  vector<Energy> _mass;

  vector<Energy> _width;

  vector<Energy2> _mass2;

  vector<Energy2> _massw;

  vector<Energy> _massa,_massb;

  vector<Energy>  _mom;

  vector<InvEnergy2> _dhdq2;

  vector<double> _hm2;

  vector<double> _dparam;
};

}


#include "ThePEG/Utilities/ClassTraits.h"

namespace ThePEG {

template <>
struct BaseClassTrait<Herwig::TwoMesonRhoKStarCurrent,1> {
  typedef Herwig::WeakDecayCurrent NthBase;
};

template <>
struct ClassTraits<Herwig::TwoMesonRhoKStarCurrent>
  : public ClassTraitsBase<Herwig::TwoMesonRhoKStarCurrent> {
  static string className() { return "Herwig::TwoMesonRhoKStarCurrent"; }
  static string library() { return "HwWeakCurrents.so"; }

};

}

#endif /* HERWIG_TwoMesonRhoKStarCurrent_H */