#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());
}