FFDipole.h

// -*- C++ -*-
//
// FFDipole.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_FFDipole_H
#define HERWIG_FFDipole_H
//
// This is the declaration of the FFDipole class.
//
 
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig/Decay/DecayIntegrator.fh"
#include "Herwig/Utilities/Kinematics.h"
#include "Herwig/Utilities/Maths.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Vectors/Lorentz5Vector.h"
#include "ThePEG/Interface/Interfaced.h"
#include "FFDipole.fh"
 
namespace Herwig {
using namespace ThePEG;
using ThePEG::Constants::pi;
 
class FFDipole: public Interfaced {
 
public:
 
  FFDipole() :
    _emin(1.e-6*MeV), _eminrest(100*MeV), _eminlab(100*MeV), _emax(),
    _multiplicity(), _m(3), _charge(), _qdrf(2),
    _qnewdrf(2), _qprf(2), _qnewprf(2), _qlab(2), _qnewlab(2), _dipolewgt(),
    _yfswgt(), _jacobianwgt(), _mewgt(), _maxwgt(7.0), _mode(1), _maxtry(500),
    _energyopt(1), _betaopt(4), _dipoleopt(), _nweight(0), _wgtsum(0.), _wgtsq(0.),
    _weightOutput(false) {}
 
public:
 
  virtual ParticleVector generatePhotons(const Particle & p,
                                         ParticleVector children,
                                         tDecayIntegratorPtr decayer);
 
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:
 
  double photon(double beta1,double ombeta1, double beta2, double ombeta2);
 
  double exactDipoleWeight(double beta1,double ombeta1,
                           double beta2,double ombeta2,unsigned int iphot) {
    double opbc,ombc;
    // if cos is greater than zero use result accurate as cos->1
    if(_cosphot[iphot]>0) {
      opbc=1.+beta2*_cosphot[iphot];
      ombc=ombeta1+beta1*sqr(_sinphot[iphot])/(1.+_cosphot[iphot]);
    }
    // if cos is less    than zero use result accurate as cos->-1
    else {
      opbc=ombeta2+beta2*sqr(_sinphot[iphot])/(1.-_cosphot[iphot]);
      ombc=1.-beta1*_cosphot[iphot];
    }
    return 0.5/(opbc*ombc)*(1.+beta1*beta2
                            -0.5*ombeta1*(1.+beta1)*opbc/ombc            
                            -0.5*ombeta2*(1.+beta2)*ombc/opbc);
  }
  
  double jacobianWeight() {
    Energy pcm1=Kinematics::pstarTwoBodyDecay(_m[0],_m[1],_m[2]);
    Energy m12 =sqrt((_qnewdrf[0]+_qnewdrf[1]).m2())            ;
    Energy pcm2=Kinematics::pstarTwoBodyDecay(m12,_m[1],_m[2])  ;
    double betaprobeta = pcm2*_m[0]/pcm1/m12   ;
    double spros       = sqr(m12/_m[0])        ;
    double deltafn     = m12/(m12+_bigLdrf.e());
    return betaprobeta*spros*deltafn           ;
  }
 
  double meWeight(const ParticleVector & children);
 
  double makePhotons(const Boost & boost, 
                     const ParticleVector & children);
 
  bool boostMomenta(const Boost & boost);
 
  unsigned int removePhotons();
 
  double collinearWeight(const ParticleVector & children);
 
  double virtualWeight(const ParticleVector & children);
 
protected:
 
  virtual void dofinish();
 
private:
 
  FFDipole & operator=(const FFDipole &) = delete;
 
private:
 
  void printDebugInfo(const Particle & p,
                      const ParticleVector & children,
                      double wgt) const;
 
private:
 
  Energy _emin;
 
  Energy _eminrest;
 
  Energy _eminlab;
 
  Energy _emax;
 
  unsigned int _multiplicity;
 
  vector<Energy> _m;
 
  double _charge;
 
  vector<Lorentz5Momentum> _qdrf;
 
  vector<Lorentz5Momentum> _qnewdrf;
 
  vector<Lorentz5Momentum> _ldrf;
 
  Lorentz5Momentum _bigLdrf;
 
  vector<Lorentz5Momentum> _qprf;
 
  vector<Lorentz5Momentum> _qnewprf;
 
  vector<Lorentz5Momentum> _lprf;
 
  Lorentz5Momentum _bigLprf;
 
  vector<Lorentz5Momentum> _qlab;
  
  vector<Lorentz5Momentum> _qnewlab;
 
  vector<Lorentz5Momentum> _llab;
 
  Lorentz5Momentum _bigLlab;
 
 
  double _dipolewgt;
 
  double _yfswgt;
 
  double _jacobianwgt;
 
  double _mewgt;
 
  double _maxwgt;
 
  vector<double> _cosphot;
 
  vector<double> _sinphot;
 
  vector<double> _photonwgt;
 
  vector<bool> _photcut;
 
  unsigned int _mode;
 
  unsigned int _maxtry;
 
  unsigned int _energyopt;
 
  unsigned int _betaopt;
 
  unsigned int _dipoleopt;
 
  tDecayIntegratorPtr _decayer;
 
  tPPtr _parent;
 
  unsigned int _nweight;
 
  double _wgtsum;
 
  double _wgtsq;
 
  bool _weightOutput;
};
 
}
 
#endif /* HERWIG_FFDipole_H */