Beam.h

// Header file for information on incoming beams, including PDF's.

#ifndef Beam_H
#define Beam_H

#include <iostream>
#include <vector>
#include <cmath>
#include <ctime>
#include "Pythia7/Config/Pythia7.h"
#include "ThePEG/PDF/PDFBase.h"

namespace Pythia7 {
namespace Shower {

using namespace std;

//**************************************************************************

class PDF {

public:
  PDF(long idBeamin = 2212)
    : idBeam(idBeamin), xSav(-1.0), Q2Sav(-1.0), xu(0.0), xd(0.0), xubar(0.0),
      xdbar(0.0), xs(0.0), xc(0.0), xb(0.0), xg(0.0), xlepton(0.0),
      xgamma(0.0) {}
  virtual ~PDF() {}
  double xfx(long, double, double);  
  virtual bool isParton(long) const; 
  virtual bool isValence(long) const; 
  
protected:
  long idBeam;
  double xSav, Q2Sav;
  double xu, xd, xubar, xdbar, xs, xc, xb, xg, xlepton, xgamma;
  virtual void xfUpdate(double x, double Q2) = 0; 

};

//*********


class GRV94L : public PDF {

public:
  GRV94L(long idBeamin = 2212) : PDF(idBeamin) {;}

private:
  void xfUpdate(double x, double Q2);
  double grvv (double x, double n, double ak, double bk, double a, 
    double b, double c, double d);
  double grvw (double x, double s, double al, double be, double ak, 
    double bk, double a, double b, double c, double d, double e, double es);
  double grvs (double x, double s, double sth, double al, double be, 
    double ak, double ag, double b, double d, double e, double es);
};

//*********

 
class Lepton : public PDF {

public:
  Lepton(long idBeamin = 11) : PDF(idBeamin) {;}

private:
  void xfUpdate(double x, double Q2);
};

//**************************************************************************


class BeamParticle : public Particle {

public:
  BeamParticle() : Particle() {pdfbm = 0;}  
  BeamParticle(long idin, Vec4 pin = Vec4(0.,0.,0.,0.), double min = 0., 
    PDF* pdfin = 0) :  Particle(idin, -9, -1, -1, 0, 0, pin, min, 0.) 
    {pdf(pdfin);} 
  BeamParticle(const BeamParticle& pt) : Particle(pt) {pdfbm = pt.pdfbm;}  
  BeamParticle& operator=(const BeamParticle& pt) {if (this != &pt) {
    pdfbm = pt.pdfbm; return *this;} } 
  BeamParticle(const Particle& pt) : Particle(pt) {pdfbm = 0;}  
  BeamParticle& operator=(const Particle& pt) {Particle::operator=(pt);   
    pdfbm = 0; return *this; }  


  double xfx(long id, double x, double Q2) {return pdfbm->xfx(id, x, Q2);}
  void pdf(PDF* pdfin) {pdfbm = pdfin;}
  PDF* pdf() const {return pdfbm;}
  bool isParton(long id) const {return pdfbm->isParton(id);} 
  bool isValence(long id) const {return pdfbm->isValence(id);} 
  bool isHadron() const {if (abs(id()) > 100) return true; return false;}
  bool isLepton() const {if (abs(id()) == 11 || abs(id()) == 13
    || abs(id()) == 15) return true; return false;}


private: 
  PDF* pdfbm;
};

class ThePEGPDF: public PDF {

public:

  ThePEGPDF(tcPDFPtr pdf, tcPDPtr parent);

  virtual void xfUpdate(double x, double Q2);

  void set(long id, double xf);

private:

  tcPDFPtr thePDF;
  tcPDPtr theParent;
  map<long, tcPDPtr> thePartons;
  long theLeptonID;
  

};



//**************************************************************************
}
}

#endif // Beam_H

---