#include "BEAM/Main/Beam_Spectra_Handler.H"
#include "BEAM/Main/Beam_Base.H"
#include "BEAM/Main/Monochromatic.H"
#include "BEAM/Main/Spectrum_Reader.H"
#include "BEAM/Main/Laser_Backscattering.H"
#include "BEAM/Main/EPA.H"
#include "ATOOLS/Org/Run_Parameter.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Org/Scoped_Settings.H"
#include <stdio.h>
#include "ATOOLS/Org/Exception.H"
#include "ATOOLS/Org/My_MPI.H"
#include "ATOOLS/Phys/KF_Table.H"


using namespace ATOOLS;
using namespace BEAM;
using namespace std;


Beam_Spectra_Handler::Beam_Spectra_Handler():
  p_BeamBase(NULL)
{
  RegisterDefaults();

  p_BeamBase = new Beam_Base*[2];
  for (short int i=0;i<2;i++) p_BeamBase[i] = NULL;

  if (!(SpecifySpectra() && InitKinematics())) {
    msg_Error()<<"Error in Beam_Spectra_Handler::Beam_Spectra_Handler :"<<endl
	       <<"    Could not init spectra or kinematics. Abort program."<<endl;
    ATOOLS::Abort();
  }

  m_mode = 0;
  m_polarisation = 0;
  for (short int i=0;i<2;i++) {
    if (p_BeamBase[i]->On()) m_mode += i+1;
    if (p_BeamBase[i]->PolarisationOn()) m_polarisation += i+1;
  }
  ATOOLS::rpa->gen.SetBeam1(p_BeamBase[0]->Beam());
  ATOOLS::rpa->gen.SetBeam2(p_BeamBase[1]->Beam());
  ATOOLS::rpa->gen.SetPBeam(0,p_BeamBase[0]->InMomentum());
  ATOOLS::rpa->gen.SetPBeam(1,p_BeamBase[1]->InMomentum());
}

void Beam_Spectra_Handler::RegisterDefaults()
{
  Settings& s = Settings::GetMainSettings();

  // NOTE: for backwards-compatibility we allow the use of <setting_name>_i
  // with i=1,2 as alternatives for BEAMS, BEAM_SPECTRA, and BEAM_ENERGIES. We
  // do not advertise this in the manual, it's only to make conversion of run
  // cards less error-prone.

  const auto defbeam = 0;
  const auto beam1 = s["BEAM_1"].SetDefault(defbeam).Get<int>();
  const auto beam2 = s["BEAM_2"].SetDefault(defbeam).Get<int>();
  s["BEAMS"].SetDefault({beam1, beam2});

  std::string defspectrum {"Monochromatic"};
  const auto spectrum1
    = s["BEAM_SPECTRUM_1"].SetDefault(defspectrum).Get<std::string>();
  const auto spectrum2
    = s["BEAM_SPECTRUM_2"].SetDefault(defspectrum).Get<std::string>();
  s["BEAM_SPECTRA"].SetDefault({spectrum1, spectrum2});

  const auto defenergy = 0.0;
  const auto energy1 = s["BEAM_ENERGY_1"].SetDefault(defenergy).Get<double>();
  const auto energy2 = s["BEAM_ENERGY_2"].SetDefault(defenergy).Get<double>();
  s["BEAM_ENERGIES"].SetDefault({energy1, energy2});

  s["BEAM_POLARIZATIONS"].SetDefault(0.0);
  s["BEAM_SMIN"].SetDefault(1e-10);
  s["BEAM_SMAX"].SetDefault(1.0);
}

void Beam_Spectra_Handler::RegisterLaserDefaults(const int& lasermode)
{
  Settings& s = Settings::GetMainSettings();
  RegisterLaserEnergyAndPolarizationDefaults();
  s["LASER_MODE"].SetDefault(lasermode);
  s["LASER_ANGLES"].SetDefault(false);
  s["LASER_NONLINEARITY"].SetDefault(false);
}

void Beam_Spectra_Handler::RegisterSpectrumReaderDefaults()
{
  Settings& s = Settings::GetMainSettings();
  RegisterLaserEnergyAndPolarizationDefaults();
  s["SPECTRUM_FILES"].SetDefault("");
}

void Beam_Spectra_Handler::RegisterLaserEnergyAndPolarizationDefaults()
{
  Settings& s = Settings::GetMainSettings();
  s["E_LASER"].SetDefault(0.0);
  s["P_LASER"].SetDefault(0.0);
}

Beam_Spectra_Handler::~Beam_Spectra_Handler() { 
  for (short int i=0;i<2;i++) {
    if (p_BeamBase[i]) { delete p_BeamBase[i]; p_BeamBase[i] = NULL; }
  }
  if (p_BeamBase) { delete [] p_BeamBase; p_BeamBase = NULL; }
}

bool Beam_Spectra_Handler::Init()
{
  bool init(p_BeamBase[0]->Init());
  if (!p_BeamBase[1]->Init()) init=false;
  return init;
}

bool Beam_Spectra_Handler::SpecifySpectra()
{
  bool okay(true);
  char help[20];
  Beam_Type::code      beam_spec;
  std::vector<std::string> beam_spectra{
    Settings::GetMainSettings()["BEAM_SPECTRA"].GetVector<std::string>() };
  if (beam_spectra.size() == 0 || beam_spectra.size() > 2)
    THROW(fatal_error, "Specify either one or two values for `BEAM_SPECTRA'.");
  for (short int num=0;num<2;num++) {
    sprintf(help,"%i",num+1);
    std::string number   = string(help); 
    std::string bs{ (num == 0) ? beam_spectra.front() : beam_spectra.back() };
    if      (bs=="Monochromatic")        beam_spec=Beam_Type::Monochromatic;
    else if (bs=="Gaussian")             beam_spec=Beam_Type::Gaussian;
    else if (bs=="Laser_Backscattering") beam_spec=Beam_Type::Laser_Back;
    else if (bs=="Simple_Compton")       beam_spec=Beam_Type::Simple_Compton;
    else if (bs=="Spectrum_Reader")      beam_spec=Beam_Type::Spec_Read;
    else if (bs=="EPA")                  beam_spec=Beam_Type::EPA;
    else                                 beam_spec=Beam_Type::Unknown;
    switch (beam_spec) {
    case Beam_Type::Monochromatic :
      okay = okay&&InitializeMonochromatic(num);
      break;
    case Beam_Type::Gaussian :
      msg_Error()<<"Error in Beam_Initialization::SpecifySpectra :"<<endl
		 <<"   Gaussian beam spectra still have to be implemented."<<endl 
		 <<"   Will read in parameters, check the procedure and abort later."<<endl;
      okay = 0;
      break;
    case Beam_Type::Simple_Compton :
      RegisterLaserDefaults(-1);
      okay = okay&&InitializeLaserBackscattering(num);
      break;
    case Beam_Type::Laser_Back :
      RegisterLaserDefaults(0);
      okay = okay&&InitializeLaserBackscattering(num);
      break;
    case Beam_Type::EPA :
      okay = okay&&InitializeEPA(num);
      break;
    case Beam_Type::Spec_Read :
      RegisterSpectrumReaderDefaults();
      okay = okay&&InitializeSpectrumReader(num);
      break;
    default :
      msg_Error()<<"Warning in Beam_Initialization::SpecifySpectra :"<<endl
		 <<"   No beam spectrum specified for beam "<<num+1<<endl
		 <<"   Will initialize monochromatic beam."<<endl;
      okay = okay&&InitializeMonochromatic(num);
      break;
    }
  }
  return okay;
}

bool Beam_Spectra_Handler::InitializeLaserBackscattering(int num)
{
  Settings& s = Settings::GetMainSettings();
  std::vector<int> beam{ s["BEAMS"].GetVector<int>() };
  if (beam.size() != 1 && beam.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAMS'.");
  int flav{ (num == 0) ? beam.front() : beam.back() };
  InitializeFlav((kf_code)abs(flav));
  Flavour beam_particle     = Flavour((kf_code)abs(flav));
  if (flav<0) beam_particle = beam_particle.Bar();

  std::vector<double> beam_energies{ s["BEAM_ENERGIES"].GetVector<double>() };
  if (beam_energies.size() != 1 && beam_energies.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAM_ENERGIES'.");
  double beam_energy{ (num == 0) ? beam_energies.front() : beam_energies.back() };

  std::vector<double> beam_polarizations{ s["BEAM_POLARIZATIONS"].GetVector<double>() };
  if (beam_polarizations.size() != 1 && beam_polarizations.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAM_POLARIZATIONS'.");
  double beam_polarization{ (num == 0) ? beam_polarizations.front() : beam_polarizations.back() };

  if ((beam_particle!=Flavour(kf_e)) && (beam_particle!=Flavour(kf_e).Bar())) {
    msg_Error()<<"Error in Beam_Initialization::SpecifySpectra :"<<endl
	       <<"   Tried to initialize Laser_Backscattering for "
	       <<beam_particle<<"."<<endl
	       <<"   This option is not available. "
	       <<"Result will be to terminate program."<<endl;
    return false;
  }

  std::vector<double> laser_energies{ s["E_LASER"].GetVector<double>() };
  if (laser_energies.size() != 1 && laser_energies.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `E_LASER'.");
  double Laser_energy{ (num == 0) ? laser_energies.front() : laser_energies.back() };

  std::vector<double> laser_polarizations{ s["P_LASER"].GetVector<double>() };
  if (laser_polarizations.size() != 1 && laser_polarizations.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `P_LASER'.");
  double Laser_polarization{ (num == 0) ? laser_polarizations.front() : laser_polarizations.back() };

  int mode = s["LASER_MODE"].Get<int>();
  bool angles = s["LASER_ANGLES"].Get<bool>();
  bool nonlin = s["LASER_NONLINEARITY"].Get<bool>();

  p_BeamBase[num]= new Laser_Backscattering(
      beam_particle,beam_energy,beam_polarization,
      Laser_energy,Laser_polarization,
      mode,(int)angles,(int)nonlin,1-2*num);

  return true;
}

bool Beam_Spectra_Handler::InitializeSpectrumReader(int num)
{
  Settings& s = Settings::GetMainSettings();
  std::vector<int> beam{ s["BEAMS"].GetVector<int>() };
  if (beam.size() != 1 && beam.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAMS'.");

  int flav{ (num == 0) ? beam.front() : beam.back() };
  InitializeFlav((kf_code)abs(flav));
  Flavour beam_particle     = Flavour((kf_code)abs(flav));
  if (flav<0) beam_particle = beam_particle.Bar();

  std::vector<double> beam_energies{ s["BEAM_ENERGIES"].GetVector<double>() };
  if (beam_energies.size() != 1 && beam_energies.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAM_ENERGIES'.");
  double beam_energy{ (num == 0) ? beam_energies.front() : beam_energies.back() };

  std::vector<double> beam_polarizations{ s["BEAM_POLARIZATIONS"].GetVector<double>() };
  if (beam_polarizations.size() != 1 && beam_polarizations.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAM_POLARIZATIONS'.");
  double beam_polarization{ (num == 0) ? beam_polarizations.front() : beam_polarizations.back() };

  std::vector<double> laser_energies{ s["E_LASER"].GetVector<double>() };
  if (laser_energies.size() != 1 && laser_energies.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `E_LASER'.");
  double laser_energy{ (num == 0) ? laser_energies.front() : laser_energies.back() };

  std::vector<double> laser_polarizations{ s["P_LASER"].GetVector<double>() };
  if (laser_polarizations.size() != 1 && laser_polarizations.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `P_LASER'.");
  double laser_polarization{ (num == 0) ? laser_polarizations.front() : laser_polarizations.back() };

  std::vector<std::string> spectrum_files{ s["SPECTRUM_FILES"].GetVector<std::string>() };
  if (spectrum_files.size() != 1 && spectrum_files.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `SPECTRUM_FILES'.");
  std::string spectrumfile{ (num == 0) ? spectrum_files.front() : spectrum_files.back() };

  p_BeamBase[num] = new Spectrum_Reader(beam_particle,beam_energy,beam_polarization,
					laser_energy, laser_polarization,
					spectrumfile,1-2*num);

  return true;
}

bool Beam_Spectra_Handler::InitializeMonochromatic(int num)
{
  Settings& s = Settings::GetMainSettings();
  std::vector<int> beam{ s["BEAMS"].GetVector<int>() };
  if (beam.size() != 1 && beam.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAMS'.");

  int flav{ (num == 0) ? beam.front() : beam.back() };
  InitializeFlav((kf_code)abs(flav));
  Flavour beam_particle     = Flavour((kf_code)abs(flav));
  if (flav<0) beam_particle = beam_particle.Bar();

  std::vector<double> beam_energies{ s["BEAM_ENERGIES"].GetVector<double>() };
  if (beam_energies.size() != 1 && beam_energies.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAM_ENERGIES'.");
  double beam_energy{ (num == 0) ? beam_energies.front() : beam_energies.back() };

  std::vector<double> beam_polarizations{ s["BEAM_POLARIZATIONS"].GetVector<double>() };
  if (beam_polarizations.size() != 1 && beam_polarizations.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAM_POLARIZATIONS'.");
  double beam_polarization{ (num == 0) ? beam_polarizations.front() : beam_polarizations.back() };

  p_BeamBase[num] = new Monochromatic(beam_particle,beam_energy,beam_polarization,1-2*num);

  return true;
}


bool Beam_Spectra_Handler::InitializeEPA(int num)
{
  Settings& s = Settings::GetMainSettings();
  std::vector<int> beam{ s["BEAMS"].GetVector<int>() };
  if (beam.size() != 1 && beam.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAMS'.");

  int flav{ (num == 0) ? beam.front() : beam.back() };
  InitializeFlav((kf_code)abs(flav));
  Flavour beam_particle     = Flavour((kf_code)abs(flav));
  if (flav<0) beam_particle = beam_particle.Bar();

  if (!(abs(flav)==kf_p_plus || abs(flav)==kf_e) &&
      !beam_particle.IsIon()) {
    msg_Error()<<"Error in Beam_Initialization::SpecifySpectra :"<<endl
               <<"   Tried to initialize EPA for "<<beam_particle<<"."<<endl
	       <<"   This option is not available. "
	       <<"Result will be to terminate program."<<endl;
    return false;
  }

  std::vector<double> beam_energies{ s["BEAM_ENERGIES"].GetVector<double>() };
  if (beam_energies.size() != 1 && beam_energies.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAM_ENERGIES'.");
  double beam_energy{ (num == 0) ? beam_energies.front() : beam_energies.back() };

  if (beam_particle.IsIon()) { 
    beam_energy *= beam_particle.GetAtomicNumber();
    // for ions the energy is specified as nucleon energy and not as energy of the
    // whole beam particle
  }
  msg_Tracking() << "InitializeEPA: Beam energy " << beam_energy << std::endl;

  std::vector<double> beam_polarizations{ s["BEAM_POLARIZATIONS"].GetVector<double>() };
  if (beam_polarizations.size() != 1 && beam_polarizations.size() != 2)
    THROW(fatal_error, "Specify either one or two values for `BEAM_POLARIZATIONS'.");
  double beam_polarization{ (num == 0) ? beam_polarizations.front() : beam_polarizations.back() };

  double  beam_mass         = beam_particle.Mass(true);
  double  beam_charge       = beam_particle.Charge();
  p_BeamBase[num]           = new EPA(beam_particle,beam_mass,beam_charge,
				      beam_energy,beam_polarization,
				      1-2*num);

  return true;
}

bool Beam_Spectra_Handler::InitKinematics()
{
  Settings& settings = Settings::GetMainSettings();
  // cms system from beam momenta - this is for potential asymmetric collisions.
  Vec4D  P      = p_BeamBase[0]->InMomentum()+p_BeamBase[1]->InMomentum();
  double s      = P.Abs2();
  double E      = sqrt(s);
  rpa->gen.SetEcms(E);
  settings.AddGlobalTag("E_CMS", ToString(rpa->gen.Ecms()));

  m_splimits[0] = s*settings["BEAM_SMIN"].Get<double>();
  m_splimits[1] = s*ATOOLS::Min(settings["BEAM_SMAX"].Get<double>(),Upper1()*Upper2());
  m_splimits[2] = s;
  m_ylimits[0]  = -10.;
  m_ylimits[1]  = 10.;
  m_exponent[0] = .5;
  m_exponent[1] = .98 * ( p_BeamBase[0]->Exponent() + p_BeamBase[1]->Exponent());
  m_mass12      = sqr(p_BeamBase[0]->Bunch().Mass());
  m_mass22      = sqr(p_BeamBase[1]->Bunch().Mass());
  double x      = 1./2.+(m_mass12-m_mass22)/(2.*E*E);
  double E1     = x*E;
  double E2     = E-E1;
  m_fiXVECs[0]  = Vec4D(E1,0.,0., sqrt(sqr(E1)-m_mass12));
  m_fiXVECs[1]  = Vec4D(E2,0.,0.,-sqrt(sqr(E1)-m_mass12));
  m_asymmetric  = 0;
  if ((dabs((m_fiXVECs[0]+(-1.)*p_BeamBase[0]->InMomentum()).Abs2())>0.0000001) ||
      (dabs((m_fiXVECs[1]+(-1.)*p_BeamBase[1]->InMomentum()).Abs2())>0.0000001) ) m_asymmetric = 1;


  m_type = p_BeamBase[0]->Type() + std::string("*") + p_BeamBase[1]->Type();
  return true;
}


void Beam_Spectra_Handler::Output() {
  msg_Info()<<"Beam_Spectra_Handler : "<<endl
	   <<"   type = "<<m_type<<endl
	   <<"   for    "<<p_BeamBase[0]->Beam()<<"  ("<<p_BeamBase[0]->InMomentum()<<")"<<endl
	   <<"   and    "<<p_BeamBase[1]->Beam()<<"  ("<<p_BeamBase[1]->InMomentum()<<")"<<endl;
}


bool Beam_Spectra_Handler::CheckConsistency(ATOOLS::Flavour * _beams,
					    ATOOLS::Flavour * _bunches) {
  bool fit = 1;
  for (int i=0;i<2;i++) {
    if ((_beams[i]!=GetBeam(i)->Beam()) || (_bunches[i]!=GetBeam(i)->Bunch())) {
      fit = 0;
      break;
    }
    /*
      if (p_BeamBase[i]->Type() == string("Laser_Backscattering")) {
      if (! ( ((_beams[i]==Flavour(kf_e)) || (_beams[i]==Flavour(kf_e).Bar())) &&
      (_bunches[i]==Flavour(kf_photon))         ) ) {
      fit = 0;
      break;
      }
      }
      if (p_BeamBase[i]->Type() == string("Beam_Strahlung")) {
      if (! ( ((_beams[i] == Flavour(kf_e)) || (_beams[i] == Flavour(kf_e).Bar())) &&
      (_beams[i] == _bunches[i])         ) ) {
      fit = 0;
      break;
      }
      }
      if (p_BeamBase[i]->Type() == string("Monochromatic") ||
      p_BeamBase[i]->Type() == string("Gaussian") ) {
      if (_bunches[i]!=_beams[i]) {
      fit = 0;
      break;
      }
      }
    */
  }
  return fit;
}

bool Beam_Spectra_Handler::CheckConsistency(ATOOLS::Flavour * _bunches) {
  bool fit = 1;
  for (int i=0;i<2;i++) {
    if (_bunches[i]!=GetBeam(i)->Bunch()) {
      fit = 0;
      break;
    }
    /*
      if (p_BeamBase[i]->Type() == string("Laser_Backscattering")) {
      if (_bunches[i]!=Flavour(kf_photon)) {
      fit = 0;
      break;
      }
      }
      if (p_BeamBase[i]->Type() == string("Beam_Strahlung")) {
      if ((_bunches[i]!=Flavour(kf_e) && _bunches[i]!=Flavour(kf_e).Bar()) ||
      (_bunches[i]!=GetBeam(i)->Bunch()) ){
      fit = 0;
      break;
      }
      }
      if (p_BeamBase[i]->Type() == string("Monochromatic") ||
      p_BeamBase[i]->Type() == string("Gaussian") ) {
      if (_bunches[i]!=GetBeam(i)->Bunch()) {
      fit = 0;
      break;
      }
      }
    */
  }
  return fit;
}


bool Beam_Spectra_Handler::MakeBeams(Vec4D * p) 
{
  double sprime=m_spkey[3], y=m_ykey[2];
  if (m_mode==0) {
    m_x1 = m_x2 = 1.;
    p[0] = m_fiXVECs[0];
    p[1] = m_fiXVECs[1];
    return true;
  }
  else {
    if ( (sprime<m_splimits[0]) || (sprime>m_splimits[1]) || m_splimits[0]==m_splimits[1] ) {
      return false;
    }

    double E      = sqrt(m_splimits[2]);
    double Eprime = sqrt(sprime);
    double x      = 1./2.+(m_mass12-m_mass22)/(2.*sprime);
    double E1     = x*Eprime;
    double E2     = Eprime-E1;
    
    p[0]          = Vec4D(E1,0.,0.,sqrt(sqr(E1)-m_mass12));
    p[1]          = Vec4D(E2,(-1.)*Vec3D(p[0]));
    E1            = exp(y);  
    E2            = exp(-y);  


    m_CMSBoost    = Poincare(Vec4D(E1+E2,0.,0.,E1-E2));
    
    Vec4D p1      = p[0];
    Vec4D p2      = p[1];
    m_CMSBoost.BoostBack(p1);
    m_CMSBoost.BoostBack(p2);
    m_x1          = 2.*p1[0]/E;
    m_x2          = 2.*p2[0]/E;

    if (m_mode==1) m_x2 = 1.;
    if (m_mode==2) m_x1 = 1.;

    p_BeamBase[0]->SetX(m_x1);
    p_BeamBase[1]->SetX(m_x2);

    return true;
  }
}

void Beam_Spectra_Handler::InitializeFlav(kf_code flav)
{
  if (s_kftable.find(flav)==s_kftable.end()) {
    bool initialize_diquarks(false);
    if (flav==kf_p_plus) {
      s_kftable[flav]=
        new Particle_Info(kf_p_plus,0.938272,0,3,1,1,1,"P+","P^{+}");
      initialize_diquarks=true;
    }
    else if (flav==kf_n) {
      s_kftable[flav]=
        new Particle_Info(kf_n,0.939566,7.424e-28,0,0,1,1,"n","n");
      initialize_diquarks=true;
    }
    else if (flav==kf_e) {
      s_kftable[flav]=
        new Particle_Info(kf_e,0.000511,.0,-3,0,1,0,1,1,0,"e-","e+","e^{-}","e^{+}");
    }
    else if (flav==kf_photon) {
      s_kftable[flav]=
        new Particle_Info(22,.0,.0,0,0,2,-1,1,1,0,"P","P","P","P");
    }
    else if (flav==1000822080) {
      s_kftable[flav]=
        new Particle_Info(1000822080, 193.75, 246, 0, 0, "Pb208", "Pb208");
    }
    else if (flav==1000822070) {
      s_kftable[flav]=
        new Particle_Info(1000822070, 192.82, 246, -1, 2, "Pb207", "Pb207");
    }
    else if (flav==1000822060) {
      s_kftable[flav]=
        new Particle_Info(1000822060, 192.82, 246, 0, 2, "Pb206", "Pb206");
    }
    else if (flav==1000791970) {
      s_kftable[flav]=
        new Particle_Info(1000791970, 183.5, 237, 3, 2, "Au197", "Au197");
    }
    else if (flav==1000200400) {
      s_kftable[flav]=
        new Particle_Info(1000200400, 37.26, 60, 0, 2, "Ca40", "Ca40");
    }
    else {
      THROW(fatal_error,"You specified a beam particle "+ToString(flav)+
            "which is not contained in your chosen model. Will abort.");
    }
      

    if (initialize_diquarks) {
      s_kftable[1103]=new Particle_Info(1103,0.77133,0,-2,-3,2,0,0,1,1,"dd_1","dd_1b","dd_1b","dd_1b");
      s_kftable[2101]=new Particle_Info(2101,0.57933,0,1,-3,0,0,0,1,1,"ud_0","ud_0b","ud_0b","ud_0b");
      s_kftable[2103]=new Particle_Info(2103,0.77133,0,1,-3,2,0,0,1,1,"ud_1","ud_1b","ud_1b","ud_1b");
      s_kftable[2203]=new Particle_Info(2203,0.77133,0,4,-3,2,0,0,1,1,"uu_1","uu_1b","uu_1b","uu_1b");
      s_kftable[3101]=new Particle_Info(3101,0.80473,0,-2,-3,0,0,0,1,1,"sd_0","sd_0b","sd_0b","sd_0b");
      s_kftable[3103]=new Particle_Info(3103,0.92953,0,-2,-3,2,0,0,1,1,"sd_1","sd_1b","sd_1b","sd_1b");
      s_kftable[3201]=new Particle_Info(3201,0.80473,0,1,-3,0,0,0,1,1,"su_0","su_0b","su_0b","su_0b");
      s_kftable[3203]=new Particle_Info(3203,0.92953,0,1,-3,2,0,0,1,1,"su_1","su_1b","su_1b","su_1b");
      s_kftable[3303]=new Particle_Info(3303,1.09361,0,-2,-3,2,0,0,1,1,"ss_1","ss_1b","ss_1b","ss_1b");
      s_kftable[4101]=new Particle_Info(4101,1.96908,0,1,-3,0,0,0,1,1,"cd_0","cd_0b","cd_0b","cd_0b");
      s_kftable[4103]=new Particle_Info(4103,2.00808,0,1,-3,2,0,0,1,1,"cd_1","cd_1b","cd_1b","cd_1b");
      s_kftable[4201]=new Particle_Info(4201,1.96908,0,4,-3,0,0,0,1,1,"cu_0","cu_0b","cu_0b","cu_0b");
      s_kftable[4203]=new Particle_Info(4203,2.00808,0,4,-3,2,0,0,1,1,"cu_1","cu_1b","cu_1b","cu_1b");
      s_kftable[4301]=new Particle_Info(4301,2.15432,0,1,-3,0,0,0,1,1,"cs_0","cs_0b","cs_0b","cs_0b");
      s_kftable[4303]=new Particle_Info(4303,2.17967,0,1,-3,2,0,0,1,1,"cs_1","cs_1b","cs_1b","cs_1b");
      s_kftable[4403]=new Particle_Info(4403,3.27531,0,4,-3,2,0,0,1,1,"cc_1","cc_1b","cc_1b","cc_1b");
      s_kftable[5101]=new Particle_Info(5101,5.38897,0,-2,-3,0,0,0,1,1,"bd_0","bd_0b","bd_0b","bd_0b");
      s_kftable[5103]=new Particle_Info(5103,5.40145,0,-2,-3,2,0,0,1,1,"bd_1","bd_1b","bd_1b","bd_1b");
      s_kftable[5201]=new Particle_Info(5201,5.38897,0,1,-3,0,0,0,1,1,"bu_0","bu_0b","bu_0b","bu_0b");
      s_kftable[5203]=new Particle_Info(5203,5.40145,0,1,-3,2,0,0,1,1,"bu_1","bu_1b","bu_1b","bu_1b");
      s_kftable[5301]=new Particle_Info(5301,5.56725,0,-2,-3,0,0,0,1,1,"bs_0","bs_0b","bs_0b","bs_0b");
      s_kftable[5303]=new Particle_Info(5303,5.57536,0,-2,-3,2,0,0,1,1,"bs_1","bs_1b","bs_1b","bs_1b");
      s_kftable[5401]=new Particle_Info(5401,6.67143,0,1,-3,0,0,0,1,1,"bc_0","bc_0b","bc_0b","bc_0b");
      s_kftable[5403]=new Particle_Info(5403,6.67397,0,1,-3,2,0,0,1,1,"bc_1","bc_1b","bc_1b","bc_1b");
      s_kftable[5503]=new Particle_Info(5503,10.07354,0,-2,-3,2,0,0,1,1,"bb_1","bb_1b","bb_1b","bb_1b");
    }
  }
}

/* ----------------------------------------------------------------

   Weight calculation 

   ---------------------------------------------------------------- */


bool Beam_Spectra_Handler::CalculateWeight(double scale) 
{
  switch (m_mode) {
  case 3 :
    if ( (p_BeamBase[0]->CalculateWeight(m_x1,scale)) && 
         (p_BeamBase[1]->CalculateWeight(m_x2,scale)) ) return true;
    break;
  case 2 :
    if (p_BeamBase[1]->CalculateWeight(m_x2,scale))     return true;
    break;
  case 1 :
    if (p_BeamBase[0]->CalculateWeight(m_x1,scale))     return true;
    break;
  }
  return false;
}


double Beam_Spectra_Handler::Weight(Flavour * flin)
{
  if (flin==NULL) return (p_BeamBase[0]->Weight() * p_BeamBase[1]->Weight());
  return (p_BeamBase[0]->Weight(flin[0]) * p_BeamBase[1]->Weight(flin[1]));
}


double Beam_Spectra_Handler::Weight(int * pol_types, double *dofs)
{
  double weight = 1.;
  for (int i=0;i<2;++i) {
    if (p_BeamBase[i]->PolarisationOn()) {
      if (pol_types[i]!=99) {
	double hel=(double)pol_types[i];
	double pol=p_BeamBase[i]->Polarisation();
	double dof=dofs[i];
	if (hel*pol>0.) 
	  weight*=(1.+dabs(pol)*(dof-1.))/dof;
	else
	  weight*=(1.-dabs(pol))/dof;

	//assuming 2 degrees of freedom
	//	weight*=dabs(hel+pol)/2.;
      }
      else {
	msg_Out()<<"ERROR: unpolarised cross section for polarised beam!! "<<endl;
      } 
    }
  }
  return weight; 
}

/* ----------------------------------------------------------------

   Boosts

   ---------------------------------------------------------------- */


void  Beam_Spectra_Handler::BoostInCMS(Vec4D* p,int n) {
  for (int i=0; i<n; ++i) m_CMSBoost.Boost(p[i]);
}

void  Beam_Spectra_Handler::BoostInLab(Vec4D* p,int n) {
  for (int i=0; i<n; ++i) m_CMSBoost.BoostBack(p[i]);
}

void   Beam_Spectra_Handler::SetSprimeMin(double _spl)      
{ 
  m_splimits[0]  = Max(m_splimits[0],_spl);
  if (m_splimits[0]>m_splimits[1])  m_splimits[0]=m_splimits[1];
}

void   Beam_Spectra_Handler::SetSprimeMax(double _spl)      { m_splimits[1]  = Min(m_splimits[1],_spl); }

void Beam_Spectra_Handler::AssignKeys(ATOOLS::Integration_Info *const info)
{
  m_spkey.Assign("s' beam",4,0,info);
  m_ykey.Assign("y beam",3,0,info);
  m_xkey.Assign("x beam",5,0,info);
}

void Beam_Spectra_Handler::SetLimits() 
{
  for (short int i=0;i<2;++i) {
    m_spkey[i]=m_splimits[i];
    m_ykey[i]=m_ylimits[i];
  }
  m_spkey[2]=ATOOLS::sqr(ATOOLS::rpa->gen.Ecms());
  m_xkey[0]=-std::numeric_limits<double>::max();
  m_xkey[2]=-std::numeric_limits<double>::max();
  m_xkey[1]=log(Upper1());
  m_xkey[3]=log(Upper2());
}