a1SimpleDecayer.h

// -*- C++ -*-
//
// a1SimpleDecayer.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_a1SimpleDecayer_H
#define HERWIG_a1SimpleDecayer_H
//
// This is the declaration of the a1SimpleDecayer class.
//
 
#include "Herwig/Decay/DecayIntegrator.h"
#include "Herwig/Decay/PhaseSpaceMode.h"
#include "ThePEG/Helicity/LorentzPolarizationVector.h"
#include "Herwig/Decay/ResonanceHelpers.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
class a1SimpleDecayer: public DecayIntegrator {
 
public:
 
  a1SimpleDecayer();
 
  virtual int modeNumber(bool & cc, tcPDPtr parent, 
                         const tPDVector & children) const;
 
  double me2(const int ichan,const Particle & part,
             const tPDVector & outgoing,
             const vector<Lorentz5Momentum> & momenta,
             MEOption meopt) const;
 
  virtual void constructSpinInfo(const Particle & part,
                                 ParticleVector outgoing) const;
 
  virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
 
  virtual 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;
 
  virtual void dataBaseOutput(ofstream & os,bool header) 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();
 
  virtual void doinitrun();
 
private:
  
  Complex rhoFormFactor(Energy2 q2,int ires) const {
    Complex output(0.),norm(0.);
    for(unsigned int ix=0;ix<3;++ix) norm += _rhowgts[ix];
    if(ires<0) {
      for(unsigned int ix=0;ix<3;++ix)
        output+=_rhowgts[ix]*Resonance::BreitWignerPWave(q2,_rhomass[ix],_rhowidth[ix],_mpi,_mpi);
    }
    else {
      assert(ires<3);
      output=_rhowgts[ires]*Resonance::BreitWignerPWave(q2,_rhomass[ires],_rhowidth[ires],_mpi,_mpi);
    }
    return output/norm;
  }
 
private:
 
  a1SimpleDecayer & operator=(const a1SimpleDecayer &) = delete;
 
private:
 
  vector<Energy> _rhomass;
  
  vector<Energy> _rhowidth;
 
  vector<double> _rhowgts;
 
  bool _localparameters;
 
  InvEnergy _coupling;
 
  double _onemax;
 
  double _twomax;
 
  double _threemax;
  
  vector<double> _onewgts;
  
  vector<double> _twowgts;
  
  vector<double> _threewgts;
 
  Energy _mpi;
 
  mutable RhoDMatrix _rho;
 
  mutable vector<Helicity::LorentzPolarizationVector> _vectors;
 
};
 
}
 
#endif /* HERWIG_a1SimpleDecayer_H */