file: ../SURF/event.cc
#include <LEDA/array.h>
#include <math.h>
#include <assert.h>
#include "event.h"
#include "unistd.h"
#include "f2c.h"
#include "Cfort.h"
//passing pointers to function members is bad.
event* global;
// A more proper way to do things would be to inherit parameter into a
// new class with the correct function, then make all the functions friends
// of class event
event::event():process_vars(num_procs,num_params),process_var_count(num_procs),
g_sum(num_procs),gg_sum(num_procs),in_lepton_momentum(2),
in_lepton_energy(2),lepton_mass(2),majorana_mass(2),h_function(num_params),
z_mass(5),w_mass(5),A_value(num_params),
process_cross_section(num_procs),process_cross_section_std_dev(num_procs),
generation_efficiency(num_procs),g_function_max(0,num_procs-1,false,true),
accepted(num_procs),events(num_procs),attempted_events(num_procs),
sqrt_s(2),letters(3),kinematic_min(num_params),kinematic_max(num_params),
branching_ratio(num_procs),cross_section(num_procs),cut_min(num_params),
cut_max(num_params),variable(num_params)
{
int i;
for(i=0;i<num_procs;i++){
g_sum[i]=gg_sum[i]=accepted[i]=events[i]=0;
attempted_events[i]=0;
g_function_max(i,true)=g_function_max(i,false)=0.;
process_cross_section[i]=generation_efficiency[i]=0;
}
//Note that order is important here, because the values of the kinematic
//min's and max's are dependent.
process_vars(nr1,0)=process_vars(nr2,0)=plt;
process_vars(nr1,1)=process_vars(nr2,1)=t;
process_vars(nr1,2)=process_vars(nr2,2)=yl;
process_vars(production,0)=full_t;
process_vars(decay,0)=time;
process_var_count[nr1]=process_var_count[nr2]=3;
process_var_count[production]=process_var_count[decay]=1;
event *temp=global;//I'm putting this in to not add to the sins of nonreentry
//in case someone ever finds out how to get around the non-reentry of pythia
global=this;//I do this to avoid having lots of context words in
//the friend functions of pointers
numover=0;
cut_count=0;
final_state=active;
hadronisation=true;
process_number=nr1;
current_process=parameter("current_process",nr1); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
h_function[yl]=omega_to_yl;
h_function[t]=cnst;
h_function[full_t]=cnst;
h_function[plt]=cnst;
h_function[time]=exp_neg_v_Av;
//things which will change the output of maximizer
cut_size=.001;
safety_factor=1.1;
//free parameters and immediate dependencies;
majorana_mass[value]=parameter("majorana_mass",150); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
majorana_mass[squared]=parameter("majorana_mass_squared",make_major_square); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
w_mass[value]=parameter("w_mass",ludat2_1.pmas[23]); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
w_mass[squared]=parameter("w_mass_squared",make_w_square); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
z_mass[value]=parameter("z_mass",ludat2_1.pmas[22]); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
z_mass[squared]=parameter("z_mass_squared",make_z_square); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
z_mass[maj2_minus_val2]=parameter
("majorana_squared_minus_"+z_mass[squared].name(),make_m2mz2);
z_mass[maj2_plus_2val2]=parameter
("majorana_squared_plus_2"+z_mass[squared].name(),make_m2p2z2);
z_mass[gamma]=parameter("gamma"+z_mass[value].name(),make_gamma_z); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
w_mass[maj2_minus_val2]=parameter
("majorana_squared_minus_"+w_mass[squared].name(),make_m2mw2);
w_mass[maj2_plus_2val2]=parameter
("majorana_squared_plus_2"+w_mass[squared].name(),make_m2p2w2);
w_mass[gamma]=parameter("gamma"+w_mass[value].name(),make_gamma_w); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
full_gamma=parameter("full_gamma",make_full_gamma); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
in_lepton_energy[value]=parameter("in_lepton_energy",100); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
in_lepton_energy[squared]=parameter("in_lepton_energy_squared",make_iles); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
lepton_mass[value]=parameter("in_lepton_mass",ulmass_(&eleven)); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
lepton_mass[squared]=parameter("in_lepton_mass_squared",make_ilmas); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
in_lepton_momentum[squared]=parameter("in_lepton_momentum_squared",make_ilms); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
in_lepton_momentum[value]=parameter("in_lepton_momentum",make_ilm); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
mixing=parameter("mixing constant",.1); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
Vzhi=parameter("Vzhi",.1); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
theta_w=parameter("theta_w",.50); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
//The different variables to play with
variable[plt]=parameter("plt",25.); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
variable[yl]=parameter("yl",0.); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
variable[t]=parameter("t",-1000); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
variable[full_t]=parameter("full_t",-1000); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
variable[time]=parameter("time",0.); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
direction=parameter("direction",make_direction); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
//dependent parameters
//note that the least dependent ones must be given first
sqrt_s[value]=parameter("sqrt_s",make_sqrt_s); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
sqrt_s[squared]=parameter("s",make_s); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
u=parameter("u",makeu); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
u3=parameter("u3",makeu3); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kNt=parameter("kNt",makekNt); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
E_nu_out=parameter("E_nu_out",make_E_nu_out); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
Pnt3=parameter("Pnt3",make_Pnt3); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kNl=parameter("kNl",makekNl); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
E_nu=parameter("E_nu",make_E_nu); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
mt=parameter("mt",makemt); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
En=parameter("En",makeEn); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
En3=parameter("En3",make_En3); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
Enu3=parameter("Enu3",make_Enu3); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
yN=parameter("yN",makeyN); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
lhs_limit=parameter("lhs_limit",make_lhs_limit); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
rhs_limit=parameter("rhs_limit",make_rhs_limit); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kinematic_min[yl]=parameter("kinematic_min_yl",make_kinematic_min_yl); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kinematic_max[yl]=parameter("kinematic_max_yl",make_kinematic_max_yl); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kNt_min=parameter("kNt_min",make_kNt_min); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kinematic_min[t]=parameter("kinematic_min_t",make_kinematic_min_t); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kinematic_max[t]=parameter("kinematic_max_t",make_kinematic_max_t); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kinematic_min[full_t]=parameter("kinematic_min_full_t", //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
make_kinematic_min_full_t);
kinematic_max[full_t]=parameter("kinematic_max_full_t",0.); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kinematic_min[plt]=parameter("kinematic_min_plt",1e-6); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kinematic_max[plt]=parameter("kinematic_max_plt",make_kinematic_max_plt); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kinematic_min[time]=parameter("kinematic_min_time",0.); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
kinematic_max[time]=parameter("kinematic_max_time",100000); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
//cuts;
cut_min[yl]= -10;
cut_min[t] = -40000;
cut_min[full_t] = -40000;
cut_min[plt] = 0;
cut_min[time]= 0;
cut_max[yl]= 10;
cut_max[t] = 0;
cut_max[full_t] = 0;
cut_max[plt]= 100;
cut_max[time]=kinematic_max[time];
//calculating branching ratios
branching_ratio[nr1]=parameter("br_1",make_br12); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
branching_ratio[nr2]=parameter("br_2",make_br12); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
branching_ratio[nr3]=parameter("br_3",make_br3); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
//A = a constant which you tune for efficiency on the simulation
A_value[yl]=parameter("A_yl",makeAyl); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
A_value[t]=parameter("A_t",makeAt); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
A_value[plt]=parameter("A_plt",makeAplt); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
A_value[time]=parameter("A_time",make_A_time); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
A_value[full_t]=parameter("A_full_t",1); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
out_lepton_energy=parameter("out_lepton_energy",make_out_lepton_energy); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
out_lepton_pl=parameter("out_lepton_pl",make_out_lepton_pl); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
QQ=parameter("QQ", makeQQ); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
QQ3=parameter("QQ3", makeQQ3); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
mix_theta=parameter("mix_theta",make_mix_theta); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
letters[A]=parameter("A",makeA); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
letters[B]=parameter("B",makeB); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
letters[C]=parameter("C",makeC); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
first_term=parameter("first_term",make_first_term); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
current_process.setvalue(nr1);
cross_section[nr1]=parameter("nr1_cross",make_nr1_cross); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
current_process.setvalue(nr2);
cross_section[nr2]=parameter("nr2_cross",make_nr2_cross); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
current_process.setvalue(production);
cross_section[production]=parameter("production_cross", //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
make_majorana_production_cross);
current_process.setvalue(nr3);
cross_section[decay]=parameter("decay_cross",make_decay_time_cross); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
g_1=parameter("g_1",true); //OVERLOAD CALL: parameter: parameter.cc(parameter), parameter.cc(parameter), parameter.cc(parameter)
initial=false;
ludat3_1.mdme[173] = 0;//shut off the top-decay channel for W
ludat3_1.mdme[152] = 0;//shut off the top-decay channel for Z0
/* -- define the majorana neutrino in JETSET 7.4 code. */
kc = 79;
s_copy(ludat4_1.chaf + (kc - 1 << 3), "nu_major", 8L, 8L);//name
ludat2_1.kchg[kc - 1] = 0;//charge
ludat2_1.kchg[kc + 499] = 0;//colour information
ludat2_1.kchg[kc + 999] = 0;//particle = antiparticle
ludat2_1.pmas[kc - 1] = majorana_mass[value];
ludat2_1.pmas[kc + 499] = (float)0.;//width
ludat2_1.pmas[kc + 999] = (float)0.;
ludat2_1.pmas[kc + 1499] = (float)0.;//average lifetime
ludat3_2.mdcy[kc - 1] = 0;//Majorana neutrino as stable particle
ludat3_2.mdcy[kc + 499] = 0;//entry point into the decay channel
ludat3_2.mdcy[kc + 999] = 0;//total number of decay channels
//a nonreintrant part of the code
reinitialize();//pythia
//end nonreintrant part
global=temp;
}
void event::set_singular_var()
{//this function sets variables to be along the singularity.
//value of plt comes from solving lhs_limit for plt
//m_N^2-m_W^2
//-----------
//2(P_Nt + lhs_limit*m_t)
if (g_1) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
variable[plt].setvalue(w_mass[maj2_minus_val2]*.5/(kNt+(1-cut_size)*mt));
else
variable[plt].setvalue(w_mass[maj2_minus_val2]*.5/(kNt+(1+cut_size)*mt));
//yl along the singularity used when maximizing differential crosssection
//to see the reasoning behind this math look in:
//Johan Rathsman's diploma thesis
if (current_process==nr1) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
if (g_1) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
variable[yl].setvalue(yN);
else
variable[yl].setvalue(yN+log(1+cut_size+sqrt(2*cut_size+pow(cut_size,2))));
else
if (g_1) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
variable[yl].setvalue(-yN);
else
variable[yl].setvalue(-yN+log(1+cut_size+sqrt(2*cut_size+pow(cut_size,2))));
double temp=cross_section[current_process];
if (!g_1){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
double tempyl=variable[yl];
if (current_process==nr1) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
variable[yl].setvalue(yN-log(1+cut_size+sqrt(2*cut_size+pow(cut_size,2))));
else
variable[yl].setvalue(yN-log(1+cut_size+sqrt(2*cut_size+pow(cut_size,2))));
if (temp>cross_section[current_process]) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
variable[yl].setvalue(tempyl);
}
}
void event::maximizer()
{
g_1.setvalue(true);
variable[t].setvalue(kinematic_min[full_t]+.1);
set_singular_var();
double temp=cross_section[current_process];
variable[t].setvalue(kinematic_max[full_t]-.1);
g_function_max(current_process,g_1)=safety_factor
*max(temp,cross_section[current_process]);
g_1.setvalue(false);
variable[t].setvalue(kinematic_min[full_t]+.1);
set_singular_var();
temp=cross_section[current_process];
variable[t].setvalue(kinematic_max[full_t]-.1);
g_function_max(current_process,g_1)=safety_factor
*max(temp,cross_section[current_process]);
}
void event::print()
{//This prints out the state of an event and associated variables.
cout << " printing out variables " << endl;
int i;
cout << " cuts: {plt,t,yl}" << endl;
cout << cut_min[0] << " plt " << cut_max[0] << endl;
cout << cut_min[1] << " t " << cut_max[1] << endl;
cout << cut_min[2] << " yl " << cut_max[2] << endl;
cout << endl;
for (i=0;i<num_params;i++){
variable[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
kinematic_min[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
kinematic_max[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
cout << " h_function[" << i << "]:" << h_function[i]
<< endl;
}
for(i=0;i<2;i++){
majorana_mass[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
in_lepton_momentum[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
in_lepton_energy[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
lepton_mass[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
z_mass[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
w_mass[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
sqrt_s[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
}
QQ.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
out_lepton_pl.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
out_lepton_energy.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
mixing.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
Vzhi.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
theta_w.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
cout << "safety factor = " << safety_factor << endl;
u.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
lhs_limit.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
rhs_limit.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
mix_theta.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
kNl.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
mt.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
kNt.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
En.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
yN.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
current_process.print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
for(i=0;i<num_procs;i++){
cout << "function max 1 "<< g_function_max(i,true) << " function max 2 "
<< g_function_max(i,false) << endl;
cout << process_cross_section[i] << " "
<< process_cross_section_std_dev[i] << " "
<< generation_efficiency[i] << endl;
if (i<3) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
branching_ratio[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
}
for(i=0;i<3;i++){
letters[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
A_value[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
}
cout << " final_state:"<< final_state << endl;
cout << " hadronisation:"<< hadronisation <<endl;
cout << " process_number:"<< process_number << endl;
}
void event::reset_pythia()
{//postcondition: the event record in pythia will be reset
//pythia stuff
//This doesn't look so pretty. There are some problems inherent in
//interfacing with fortran code. Fortran uses 2-d arrays differently.
//The Fortran codes have alot of number based flags.
int i,j;
for (i=0;i<lujets_1.n;i++)//clean the slate
for (j=0;j<5;j++){
lujets_1.k[i+4000*j]=0;
lujets_1.p[i+4000*j]=lujets_1.v[i+4000*j]=0.;
}
lujets_1.n=0;
//give compressed story
lujets_1.k[0]=lujets_1.k[1]=21;
lujets_1.k[4000] = 11; //electron
lujets_1.k[4001] = -11; //positron
lujets_1.k[8000]=0;
lujets_1.k[8001]=0;
lujets_1.p[8000] = (in_lepton_momentum[value]);
lujets_1.p[8001] = -in_lepton_momentum[value];
lujets_1.p[12000]=lujets_1.p[12001]= in_lepton_energy[value];
lujets_1.p[16000] = lepton_mass[value];
lujets_1.p[16001] = lepton_mass[value];
lujets_1.n = 2;
}
void event::reinitialize()
{//To be called before the first occur()
//postcondition: all parameters are internally self consistent
//initialize pythia
if (process_number==nr1 || process_number==mixed123){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
current_process.setvalue(nr1);
maximizer();
}
if (process_number==nr2 || process_number==mixed123){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
current_process.setvalue(nr2);
maximizer();
}
if (process_number==nr3 || process_number==mixed123){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
current_process.setvalue(production);
variable[full_t].setvalue(kinematic_min[full_t]);
double temp=cross_section[production];
variable[full_t].setvalue(kinematic_max[full_t]);
g_function_max(production,true)=g_function_max(production,false)
=safety_factor*max(temp,cross_section[production]);
current_process.setvalue(decay);
variable[time].setvalue(0.);
g_function_max(decay,true)=g_function_max(decay,false)=
safety_factor*cross_section[decay];
}
reset_pythia();
real roots = sqrt_s[value];
ludat1_1.mstu[11] = 0;// no header
pypars_1.mstp[121] = 0;
pyinit_("NONE", "e-", "e+", &roots, 4L, 2L, 2L);
ludat3_2.brat[ludat3_2.mdcy[533] + 9] = 0.;//set top decay to 0
ludatr_.mrlu[0]=getpid();
ludatr_.mrlu[1]=0;
//setting maximum spread of z mass
ludat2_1.pmas[1022]=min(10,(majorana_mass[value])-
(z_mass[value])-1);
initial=true;
}
void event::getpoint()
{//Gets a point in phase space, according to the distribution
int i,j;
bool redo=true;
double temp;
double omega;
events[current_process]++;
do
{
bool cut_active=false;
attempted_events[current_process]++;
for (i=0;i<process_var_count[current_process];i++)
switch (h_function[process_vars(current_process,i)]){
case cnst:
variable[process_vars(current_process,i)]
.setvalue(kinematic_min[process_vars(current_process,i)]
+rlu_(&zeroi)
*(kinematic_max[process_vars(current_process,i)]
-kinematic_min[process_vars(current_process,i)]));
break;
case v_inv:
variable[process_vars(current_process,i)]
.setvalue(kinematic_min[process_vars(current_process,i)]*
pow(kinematic_max[process_vars(current_process,i)]
/kinematic_min[process_vars(current_process,i)]
,rlu_(&zeroi)));
break;
case v_inv_v_inv:
variable[process_vars(current_process,i)]
.setvalue(kinematic_max[process_vars(current_process,i)]
*kinematic_min[process_vars(current_process,i)]
/(kinematic_max[process_vars(current_process,i)]
-rlu_(&zeroi)*
(kinematic_max[process_vars(current_process,i)]-
kinematic_min[process_vars(current_process,i)])));
break;
case exp_neg_v_Av:
variable[process_vars(current_process,i)]
.setvalue(log(exp(-A_value[process_vars(current_process,i)]
*kinematic_min[process_vars(current_process,i)])
-rlu_(&zeroi)
*(exp(A_value[process_vars(current_process,i)]
*kinematic_min[process_vars(current_process,i)])
-exp(-A_value[process_vars(current_process,i)]
*kinematic_max[process_vars(current_process,i)])))
/-A_value[process_vars(current_process,i)]);
break;
case v_exp_neg_v_v_Av:
variable[process_vars(current_process,i)]
.setvalue(sqrt(-1./A_value[process_vars(current_process,i)]*
log(exp(-A_value[process_vars(current_process,i)]
*pow(kinematic_min[process_vars(current_process,i)],2))
-rlu_(&zeroi)
*(exp(A_value[process_vars(current_process,i)]
*pow(kinematic_min[process_vars(current_process,i)],2))
-exp(-A_value[process_vars(current_process,i)]
*pow(kinematic_max[process_vars(current_process,i)],2))))
/A_value[process_vars(current_process,i)]));
break;
case omega_to_yl:
g_1.setvalue(lhs_limit<1);
//from major line:1303 I also recalculated from Rathsman's thesis
if (g_1) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
omega=-2*rhs_limit
/((rhs_limit-1)*sin(PI*(rlu_(&zeroi)-0.5))
-1-rhs_limit);
else
omega=-2.*lhs_limit*rhs_limit
/((rhs_limit-lhs_limit)
*sin(PI*(rlu_(&zeroi)-0.5))-lhs_limit-rhs_limit);
assert(!omega<1);
if (current_process==nr1) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
if(rlu_(&zeroi)>.5) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
variable[process_vars(current_process,i)]
.setvalue(yN+log(omega+sqrt(omega*omega-1)));
else
variable[process_vars(current_process,i)]
.setvalue(yN-log(omega+sqrt(omega*omega-1)));
else
if(rlu_(&zeroi)>.5) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
variable[process_vars(current_process,i)]
.setvalue(-yN+log(omega+sqrt(omega*omega-1)));
else
variable[process_vars(current_process,i)]
.setvalue(-yN-log(omega+sqrt(omega*omega-1)));
if (omega<lhs_limit || omega>rhs_limit) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
cout << "ERROR: " << omega << " not in " << lhs_limit << " "
<< rhs_limit << endl;
if (omega<(1+cut_size)//check to see if cut around singularity //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?); if: parameter.cc(?), parameter.cc(?)
//is violated
&& lhs_limit>(1-cut_size)
&& lhs_limit<(1+cut_size)){
cut_count++;
cut_active=true;
}
break;
default:
cerr << "invalid h function" << endl;
}
temp=cross_section[current_process];
if (temp>g_function_max(current_process,g_1)){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
numover++;
cout << "maximum =" << temp << " events " << events[current_process]
<< " current process = " << current_process
<< " current Maximums = " << g_function_max(current_process,g_1)
<< endl;
for (i=0;i<num_params;i++)
variable[i].print(); //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
redo=true;
}
else if(temp<=g_function_max(current_process,g_1)) { //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
g_sum[current_process]+=temp;
gg_sum[current_process]+=temp*temp;
}
else
cerr << "NaN detected" << endl;
double temp2;
if (temp>(temp2=rlu_(&zeroi)*g_function_max(current_process,g_1)) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
&& !cut_active){
accepted[current_process]++;
redo=false;
}
else
redo=true;
for (j=0;j<process_var_count[current_process];j++)
if (variable[process_vars(current_process,j)] //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
<cut_min[process_vars(current_process,j)])
redo=true;//check to see if we have broken any cuts. //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
}
while(redo);
generation_efficiency[current_process]
=(double)events[current_process]/(double)attempted_events[current_process];
process_cross_section[current_process]
=g_sum[current_process]/(double)attempted_events[current_process]
*(double)events[current_process]/(double)accepted[current_process];
temp=gg_sum[current_process]-pow(g_sum[current_process],2)/
(double)attempted_events[current_process];
if (temp<0)//This is really an error, but it can occur due to rounding //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
// when you have a very nice process
temp=0;
process_cross_section_std_dev[current_process]
=sqrt(temp)/(double)attempted_events[current_process];
}
void event::occur()
{//currently, only does e- e+ -> e+/- + W-/+ + neutrino(processes nr1 and nr2)
event *temp=global;
global=this;
if (!initial) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
reinitialize();
if ((process_number==mixed123 && rlu_(&zeroi)>1-branching_ratio[nr3]) || //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
process_number==nr3)
current_process.setvalue(nr3);
else if ((process_number==mixed123 && rlu_(&zeroi)>.5) || process_number==nr2) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
current_process.setvalue(nr2);
else
current_process.setvalue(nr1);
if (current_process == nr1 || current_process==nr2){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
getpoint();
}
else{
current_process.setvalue(production);
getpoint();
current_process.setvalue(decay);
getpoint();
current_process.setvalue(nr3);
}
reset_pythia();
lujets_1.k[2]=21;
lujets_1.k[3]=1;
lujets_1.k[4]=11;
lujets_1.k[5]=1;
lujets_1.k[6]=1;
lujets_1.k[4002] = 24; //exchanged W+
lujets_1.k[4003] = 12; //nu_e
lujets_1.k[4004]=79; //majorana neutrino
if (current_process==nr1){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
lujets_1.k[4005]=11; //electron
lujets_1.k[4006] = 24; //W+
}
else if (current_process==nr2){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
lujets_1.k[4005]=-11; //positron
lujets_1.k[4006] =-24; //W-
}
else {
lujets_1.k[4005]= 12; //nu_e
lujets_1.k[4006] =23; //Z
}
lujets_1.k[8002]=1;
lujets_1.k[8003]=1;
lujets_1.k[8004]=1;
lujets_1.k[8005]=5;
lujets_1.k[8006]=5;
lujets_1.k[12000]=3;
lujets_1.k[16000]=5;
lujets_1.k[12004]=6;
lujets_1.k[16004]=7;
//state of exchanged W
if (current_process==nr1 || current_process==nr2){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
lujets_1.p[12002]= in_lepton_energy[value]-En;
lujets_1.p[16002] = -sqrt(QQ);
}
else{
lujets_1.p[12002]= in_lepton_energy[value]-En3;
lujets_1.p[16002] = -sqrt(QQ3);
}
float tfloat=2*PI*rlu_(&zeroi);
//state of neutrino
if (current_process==nr1 || current_process==nr2){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
lujets_1.p[3]=-kNt*sin(tfloat);
lujets_1.p[4003]=-kNt*cos(tfloat);
lujets_1.p[8003] = -kNl;
lujets_1.p[12003]= E_nu;
}
else{
lujets_1.p[3]=-Pnt3*sin(tfloat);
lujets_1.p[4003]=-Pnt3*cos(tfloat);
lujets_1.p[8003] = Enu3*cos(direction);
lujets_1.p[12003]= Enu3;
}
lujets_1.p[16003] = 0;
//state of Heavy Majorana Neutrino
lujets_1.p[4]=-lujets_1.p[3];
lujets_1.p[4004]=-lujets_1.p[4003];
lujets_1.p[8004]=-lujets_1.p[8003];
if (current_process==nr1 || current_process==nr2) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
lujets_1.p[12004]= En;
else
lujets_1.p[12004]= En3;
lujets_1.p[16004] = majorana_mass[value];
//state of electron or neutrino
float phi=2*PI*rlu_(&zeroi);
float costheta=(rlu_(&zeroi)-.5)*2;
float theta=acos(costheta);
if (current_process==nr1 || current_process==nr2){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
lujets_1.p[5]=-sin(phi)*variable[plt];
lujets_1.p[4005]=-cos(phi)*variable[plt];
lujets_1.p[8005]=out_lepton_pl;
lujets_1.p[12005]=out_lepton_energy;
lujets_1.p[16005]=lepton_mass[value];
}
else {
lujets_1.p[5]=cos(theta)*cos(phi)*E_nu_out;
lujets_1.p[4005]=cos(theta)*sin(phi)*E_nu_out;
lujets_1.p[8005]=sin(theta)*E_nu_out;
lujets_1.p[12005]=E_nu_out;
lujets_1.p[16005]=0;
}
//state of W+/- or Z.
if (current_process==nr1 || current_process==nr2){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
lujets_1.p[6]=lujets_1.p[4]-lujets_1.p[5];
lujets_1.p[4006]=lujets_1.p[4004]-lujets_1.p[4005];
lujets_1.p[8006]=kNl-lujets_1.p[8005];
lujets_1.p[12006]=En-out_lepton_energy;
lujets_1.p[16006]=w_mass[value];
}
else {
lujets_1.p[6]=-lujets_1.p[5];
lujets_1.p[4006]=-lujets_1.p[4005];
lujets_1.p[8006]=-lujets_1.p[8005];
lujets_1.p[12006]=majorana_mass[value]-E_nu_out;
lujets_1.p[16006]=z_mass[value];
}
lujets_1.n=7;
if (current_process==nr3){ //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
ludat1_1.mstu[0]=6;
ludat1_1.mstu[1]=7;
float the=0.;
float phi=0.;
float bex=lujets_1.p[4]/lujets_1.p[12004];
float bey=lujets_1.p[4004]/lujets_1.p[12004];
float bez=lujets_1.p[8004]/lujets_1.p[12004];
lurobo_(&the,&phi,&bex,&bey,&bez);//lorentz boost the Z and nu
ludat1_1.mstu[0]=0;
ludat1_1.mstu[1]=0;
}
if (hadronisation) //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
luexec_();
global=temp;
}
C++ to HTML Conversion by ctoohtml