ThreeBodyAllOnCalculator.h

// -*- C++ -*-
//
// ThreeBodyAllOnCalculator.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_ThreeBodyAllOnCalculator_H
#define HERWIG_ThreeBodyAllOnCalculator_H
// This is the declaration of the ThreeBodyAllOnCalculator class.
 
#include "WidthCalculatorBase.h"
#include "Herwig/Utilities/GSLIntegrator.h"
 
namespace Herwig {
using namespace ThePEG;
 
template <class T>
class ThreeBodyAllOnCalculator;
 
 
template <class T>
class ThreeBodyAllOnCalculator: public WidthCalculatorBase {
 
 struct Outer {
 
   Outer(typename Ptr<Herwig::ThreeBodyAllOnCalculator<T> >::const_pointer in,
         double relerr)
    : _integrand(in), _integrator(1e-35,relerr,1000)
  {}
  
  Energy4 operator ()(double x) const {
    Energy2 low, upp;
    _integrand->outerVariables(x,low,upp);
    return _integrator.value(*_integrand,low,upp);
  }
  typedef double ArgType;
  typedef Energy4 ValType;
 
  typename Ptr<Herwig::ThreeBodyAllOnCalculator<T> >::const_pointer _integrand;
  
  GSLIntegrator _integrator;
};
 
public:
 
  friend struct ThreeBodyAllOnOuter;
 
public:
 
  ThreeBodyAllOnCalculator(vector<double> inweights,
                           vector<int> intype,
                           vector<Energy> inmass,
                           vector<Energy> inwidth,
                           vector<double> inpow,
                           T inme, int mode,
                           Energy m1,Energy m2,Energy m3,
                           double relerr=1e-3)
    : _channelweights(inweights),_channeltype(intype),_channelmass(inmass),
      _channelwidth(inwidth),_channelpower(inpow),_theME(inme),_mode(mode),
      _thechannel(0),_mapping(inweights.size(),0),_souter(ZERO),
      _integrator(1e-35,relerr,1000),_relerr(relerr) {
    _m.resize(4);
    _m[1]=m1;_m[2]=m2;_m[3]=m3;
    _m2.resize(4);
    for(int ix=1;ix<4;++ix) {
      _m2[ix]=sqr(_m[ix]);
    }
  }
 
  Energy partialWidth(Energy2 q2) const;
 
  void resetMass(int imass,Energy mass) {
    assert(imass<4);
    _m[imass]=mass;
    _m2[imass]=mass*mass;
  }
 
  Energy getMass(const int imass) const {
    assert(imass>=0&&imass<4);
    return _m[imass];
  }
 
  Energy otherMass(const int imass) const {
    assert(imass>0&&imass<4);
    if(imass==1)      return _m[2]+_m[3];
    else if(imass==2) return _m[1]+_m[3];
    else              return _m[1]+_m[2];
  }
 
  Energy2 operator ()(Energy2 argument) const;
  typedef Energy2 ArgType;
  typedef Energy2 ValType;
 
 
protected:
 
  void outerVariables(double x, Energy2 & low, Energy2 & upp) const;
 
private:
 
  ThreeBodyAllOnCalculator & operator=(const ThreeBodyAllOnCalculator &) = delete;
 
private:
  
  vector<double> _channelweights;
 
  vector<int> _channeltype;
 
  vector<Energy> _channelmass;
 
  vector<Energy> _channelwidth;
 
  vector<double> _channelpower;
 
  T _theME;
 
  int _mode;
 
  mutable int _thechannel;
 
  mutable vector<int> _mapping;
 
  mutable Energy2 _souter;
 
  mutable vector<Energy>  _m;
 
  mutable vector<Energy2> _m2;
 
  GSLIntegrator _integrator;
 
  double _relerr;
};
}
 
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "ThreeBodyAllOnCalculator.tcc"
#endif
 
#endif /* HERWIG_ThreeBodyAllOnCalculator_H */