TwoPionPhotonCurrent.h

// -*- C++ -*-
//
// TwoPionPhotonCurrent.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_TwoPionPhotonCurrent_H
#define HERWIG_TwoPionPhotonCurrent_H
//
// This is the declaration of the TwoPionPhotonCurrent class.
//
#include "WeakCurrent.h"
 
namespace Herwig {
using namespace ThePEG;
 
class TwoPionPhotonCurrent: public WeakCurrent {
 
public:
 
  TwoPionPhotonCurrent();
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
public:
 
  virtual bool createMode(int icharge, tcPDPtr resonance,
                          FlavourInfo flavour,
                          unsigned int imode,PhaseSpaceModePtr mode,
                          unsigned int iloc,int ires,
                          PhaseSpaceChannel phase, Energy upp );
 
  virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
 
  virtual vector<LorentzPolarizationVectorE> 
  current(tcPDPtr resonance,
          FlavourInfo flavour,
          const int imode, const int ichan,Energy & scale,
          const tPDVector & outgoing,
          const vector<Lorentz5Momentum> & momenta,
          DecayIntegrator::MEOption meopt) const;
 
  virtual void constructSpinInfo(ParticleVector decay) const;
 
  virtual bool accept(vector<int> id);
 
  virtual unsigned int decayMode(vector<int> id);
 
  virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
 
protected:
 
  virtual IBPtr clone() const {return new_ptr(*this);}
 
  virtual IBPtr fullclone() const {return new_ptr(*this);}
 
private:
 
  TwoPionPhotonCurrent & operator=(const TwoPionPhotonCurrent &) = delete;
 
private:
  
  complex<InvEnergy> FFunction(Energy2 q2) const {
    complex<InvEnergy2> output(ZERO);
    for(unsigned int ix=0; ix<_resweights.size();++ix) {
      output -= _resweights[ix]*BreitWigner(q2,ix);
    }
    return output*_grho*_grhoomegapi*sqrt(2.);
  }
 
  complex<InvEnergy2> BreitWigner(Energy2 q2,unsigned int ires) const {
    static const Complex ii(0.,1.);
    complex<Energy2> denom;
    if(ires<_rhomasses.size()) {
      denom = q2-_rhomasses[ires]*_rhomasses[ires]+ii*_rhomasses[ires]*_rhowidths[ires];
    }
    else if(ires==10) {
      denom = q2-_omegamass*_omegamass+ii*_omegamass*_omegawidth;
    }
    else assert(false);
    return 1./denom;
  }
  
private:
  
  Energy2 _grho;
  
  InvEnergy _grhoomegapi;
  
  vector<double> _resweights;
 
  vector<Energy> _rhomasses; 
 
  vector<Energy> _rhowidths;
 
  Energy _omegamass;
 
  Energy _omegawidth;
 
  Energy _intmass;
 
  Energy _intwidth; 
 
};
 
}
 
 
#endif /* HERWIG_TwoPionPhotonCurrent_H */