file: ../SURF/friends.cc


#include "event.h"
#include "Cfort.h"
#include <math.h>
#include <assert.h>

//references used for these equations are the program MAJOR
//and Buchmuller and Greub, Nuclear Physics B

double make_direction()
{//from B+G equation 19
  double ret = acos(-2*global->variable[full_t]
    /(global->sqrt_s[squared]-global->majorana_mass[squared])-1);
  return ret;
}

double make_En3()
{
  return .5*(global->sqrt_s[squared]+global->majorana_mass[squared])
    /global->sqrt_s[value];
}

double make_Enu3()
{
  return .5*(global->sqrt_s[squared]-global->majorana_mass[squared])
    /global->sqrt_s[value];
}

double make_E_nu_out()
{
  return global->z_mass[maj2_minus_val2]*.5/global->majorana_mass[value];
}

double make_Pnt3()
{
  return global->Enu3*abs(sin(global->direction));
}

double makeu()
{//possible bug... how to handle negative sqrt?
  double ret=global->majorana_mass[squared]+2*global->lepton_mass[squared];
  ret-=global->sqrt_s[squared];
  ret-=global->variable[t];
  return ret;
}

double makeu3()
{//possible bug... how to handle negative sqrt?
  double ret=global->majorana_mass[squared];
  ret+=2*global->lepton_mass[squared];
  ret-=global->sqrt_s[squared];
  ret-=global->variable[full_t];
  return ret;
}

double make_out_lepton_pl()
{//major line: 716
  return sqrt(global->lepton_mass[squared]+
	      pow(global->variable[plt],2))
    *sinh(global->variable[yl]);
}

double limit(double *xf, var_type var)
{
  if (xf[0]<global->kinematic_min[var])        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
    return pow(xf[0]-global->kinematic_min[var],2);
  else if (xf[0]>global->kinematic_max[var])        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
    return pow(xf[0]-global->kinematic_max[var],2);
  else return 0.;
}

double make_out_lepton_energy()
{//EL in major line 717
  return 
    sqrt(global->lepton_mass[squared]+pow(global->variable[plt],2) 
	 +pow(global->out_lepton_pl,2));
}

double makeQQ()
{//Q2 in major line 1281
  return 
    2*global->in_lepton_energy[value]*(global->En+global->kNl)
    -global->majorana_mass[squared];
}

double makeQQ3()
{//this is still not correct
  return 
    2*global->in_lepton_energy[value]*(global->En3+global->kNl)
    -global->majorana_mass[squared];
}

double makekNt()
{//From B+G pg 354
  return 
    sqrt(global->u*global->variable[t]
	 /global->sqrt_s[squared]);
}

double make_E_nu()
{//magnitude is the same although direction will be 180 apart.
  return sqrt(global->kNt*global->kNt+global->kNl*global->kNl);
}

double make_kNt_min()
{//From B+G pg 355 
  //significant sig fig loss here
  double ret=4*global->majorana_mass[squared];
  ret*=pow(global->variable[plt],2);
  ret-=pow(global->w_mass[maj2_minus_val2],2);
  ret=ret/(4*global->variable[plt]*global->w_mass[maj2_minus_val2]);
  return 
    max(1e-6,ret);
}

double makekNl()
{//B+G 354
  if (global->current_process==nr1)        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
    return (global->variable[t]-global->u)
    /(2*global->sqrt_s[value]);
  else 
    return (-global->variable[t]+global->u)
    /(2*global->sqrt_s[value]);
}

double makemt()
{//B+G 354
  return sqrt(global->majorana_mass[squared]+pow(global->kNt,2));
}

double makeEn()
{//B+G 354
  return sqrt(pow(global->mt,2)+pow(global->kNl,2));
}

double makeyN()
{//B+G 354
  return .5*log((global->En+global->kNl)
		/(global->En-global->kNl));
}

double make_lhs_limit()
{//From B+G pg 354
  return 
    (global->w_mass[maj2_minus_val2]-2*global->variable[plt]
     *global->kNt)/(2*global->mt*global->variable[plt]);
}

double make_rhs_limit()
{//From B+G pg 354
  return 
    (global->w_mass[maj2_minus_val2]+2*global->variable[plt]
     *global->kNt)/(2*global->mt*global->variable[plt]);
}

double makeA()
{//From B+G pg 354
  if (global->current_process==nr1)        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
    return 8*global->variable[t]*(global->w_mass[maj2_minus_val2]
				  *(global->variable[t]
				    -global->majorana_mass[squared])
    +(global->majorana_mass[squared]
      -pow(global->majorana_mass[squared],2)
      /(2*global->w_mass[squared]))
				    *global->sqrt_s[value]
				    *global->variable[plt]
				    *exp(-global->variable[yl]));
  else
    return
      8*global->u*(global->w_mass[maj2_minus_val2]
		   *(global->u-global->majorana_mass[squared])
      +(global->majorana_mass[squared]
	-pow(global->majorana_mass[squared],2)
	/(2*global->w_mass[squared]))
		   *global->sqrt_s[value]
		   *global->variable[plt]
		   *exp(global->variable[yl]));
  
}

double makeB()
{//From B+G pg 354
  if (global->current_process==nr1)        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
    return 4*global->u*global->majorana_mass[squared]
      /global->w_mass[squared]
      *(global->w_mass[maj2_minus_val2]
	*(global->u-global->majorana_mass[squared])
      +(global->majorana_mass[squared]-2*global->w_mass[squared])
	*global->sqrt_s[value]*global->variable[plt]
	*exp(global->variable[yl]));
  else
    return 4*global->variable[t]*global->majorana_mass[squared]
      /global->w_mass[squared]
      *(global->w_mass[maj2_minus_val2]
	*(global->variable[t]-global->majorana_mass[squared])
      +(global->majorana_mass[squared]-2*global->w_mass[squared])
	*global->sqrt_s[value]*global->variable[plt]
	*exp(-global->variable[yl]));
}

double makeC()
{//C=A with t->u + B with u->t  B+G g 354
  double ret;
  if (global->current_process==nr1){        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
    ret=8*global->u
      *(global->w_mass[maj2_minus_val2]
	*(global->u-global->majorana_mass[squared])
      +(global->majorana_mass[squared]
	-pow(global->majorana_mass[squared],2)
	/(2*global->w_mass[squared]))
	*global->sqrt_s[value]
	*global->variable[plt]
	*exp(global->variable[yl]));
    ret+=4*global->variable[t]*global->majorana_mass[squared]
      /global->w_mass[squared]
      *(global->w_mass[maj2_minus_val2]
	*(global->variable[t]-global->majorana_mass[squared])
      +(global->majorana_mass[squared]-2*global->w_mass[squared])
	*global->sqrt_s[value]*global->variable[plt]
	*exp(-global->variable[yl]));
  }
  else {
    ret=8*global->variable[t]
      *(global->w_mass[maj2_minus_val2]
	*(global->variable[t]-global->majorana_mass[squared])
      +(global->majorana_mass[squared]
	-pow(global->majorana_mass[squared],2)
	/(2*global->w_mass[squared]))
	*global->sqrt_s[value]
	*global->variable[plt]
	*exp(-global->variable[yl]));
    ret+=4*global->u*global->majorana_mass[squared]
      /global->w_mass[squared]
      *(global->w_mass[maj2_minus_val2]
	*(global->u-global->majorana_mass[squared])
    +(global->majorana_mass[squared]-2*global->w_mass[squared])
	*global->sqrt_s[value]*global->variable[plt]
	*exp(global->variable[yl]));
  }
  return ret;
}

double make_mix_theta()
{
  return global->mixing*(.5-pow(sin(global->theta_w),2))
    /(pow(cos(global->theta_w),2)
      *(global->sqrt_s[squared]-global->z_mass[squared]));
}

double make_decay_time_cross()
{//from B+G equation 12
  double ret=global->z_mass[gamma]*global->majorana_mass[value]/global->En3;
  ret*=exp(-global->variable[time]* global->full_gamma
	   *global->majorana_mass[value]/global->En3);
  ret*=hmult();
  return ret;
}

double make_majorana_production_cross()
{//from B+G equation 15
  double ret=fermi*fermi*pow(global->w_mass[squared],2)
    /(2*PI*pow(global->sqrt_s[squared],2));
  double part1, part2,part3;
  part1=(pow(global->mixing,2)*pow(tan(global->theta_w),4)
     *(global->variable[full_t]
       *(global->variable[full_t]-global->majorana_mass[squared])
     +global->u3*(global->u3-global->majorana_mass[squared]))
     /pow(global->sqrt_s[squared]-global->z_mass[squared],2));
  part2=pow(global->mix_theta-global->Vzhi
	    /(global->variable[full_t]-global->w_mass[squared]),2)
    *global->u3*(global->u3-global->majorana_mass[squared]);
  part3=pow(global->mix_theta-global->Vzhi
	    /(global->u3-global->w_mass[squared]),2)*global->variable[full_t]
    *(global->variable[full_t]-global->majorana_mass[squared]);
  ret*=(part1+part2+part3);
  ret*=hmult();
  return 3.894e8*ret;//unit conversion again.
}

double hmult()
{
  double ret=1;
  int i;
  for (i=0;i<global->process_var_count[global->current_process];i++)
    switch (global->h_function[global->process_vars((int)global->current_process,i)]){
    case cnst:
      ret=ret*(global->kinematic_max[global->process_vars((int)global->current_process,i)]-global->kinematic_min[global->process_vars((int)global->current_process,i)]);
      break;
    case v_inv:
      ret=ret*log(global->kinematic_max[global->process_vars((int)global->current_process,i)]/global->kinematic_min[global->process_vars((int)global->current_process,i)]);
      ret=ret*global->variable[global->process_vars((int)global->current_process,i)];
      break;
    case v_inv_v_inv:
      ret=ret*(1/global->kinematic_min[global->process_vars((int)global->current_process,i)]-1/global->kinematic_max[global->process_vars((int)global->current_process,i)]);
      ret=ret*pow(global->variable[global->process_vars((int)global->current_process,i)],2);
      break;
    case exp_neg_v_Av:
      ret=ret*(exp(-global->A_value[global->process_vars((int)global->current_process,i)]*global->kinematic_min[global->process_vars((int)global->current_process,i)])
	       -exp(-global->A_value[global->process_vars((int)global->current_process,i)]*global->kinematic_max[global->process_vars((int)global->current_process,i)]))
	/global->A_value[global->process_vars((int)global->current_process,i)];
      ret=ret/exp(-global->A_value[global->process_vars((int)global->current_process,i)]*global->variable[global->process_vars((int)global->current_process,i)]);
      break;
    case v_exp_neg_v_v_Av:
      ret=ret*.5
	*(exp(-global->A_value[global->process_vars((int)global->current_process,i)]*pow(global->kinematic_min[global->process_vars((int)global->current_process,i)],2))
	  -exp(-global->A_value[global->process_vars((int)global->current_process,i)]*pow(global->kinematic_max[global->process_vars((int)global->current_process,i)],2)))
	/global->A_value[global->process_vars((int)global->current_process,i)];
      ret=ret/exp(-global->A_value[global->process_vars((int)global->current_process,i)]*pow(global->variable[global->process_vars((int)global->current_process,i)],2));
      break;
    case omega_to_yl:
      //Now the same for the yl parameter
      //http://www3.tsl.uu.se/thep/papers/TSL/ISV-92-0058.ps.gz
      ret=ret*2*PI;
      if (global->g_1){        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
	ret=ret/sqrt(global->rhs_limit);
	//This is not actually h, but rather h canceled with part of the cross
	//section
	if (global->current_process==nr1)        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
	  ret=ret*cosh(global->yN-global->variable[global->process_vars(global->current_process,i)])
	  /sqrt(1.+cosh(global->yN-global->variable[global->process_vars(global->current_process,i)]));
	else
	  ret=ret*cosh(global->yN+global->variable[global->process_vars(global->current_process,i)])
	  /sqrt(1.+cosh(global->yN+global->variable[global->process_vars(global->current_process,i)]));
      }
      else{
	ret=ret/sqrt(global->rhs_limit*global->lhs_limit);
	//ditto
	if (global->current_process==nr1)        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
	  ret=ret/abs(tanh(global->yN-global->variable[global->process_vars(global->current_process,i)]));
	else 
	  ret=ret/abs(tanh(global->yN+global->variable[global->process_vars(global->current_process,i)]));
      }
      break;
    default:
      assert(0);
      break;
    }
  return ret;
}

double make_first_term()
{//the first term in B+G equation 20
  return fermi*fermi*global->branching_ratio[global->current_process]
    *pow(global->w_mass[squared],3)*global->majorana_mass[squared]
    /(4*pow(PI*global->sqrt_s[squared]*global->w_mass[maj2_minus_val2],2)
      *global->w_mass[maj2_plus_2val2]*global->mt);
}

double make_nr1_cross()
{//B+G pg 353 note that V=1
  assert(global->current_process==nr1);
  double ret=global->first_term;
  if (global->g_1)        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
    ret=ret/sqrt(cosh(global->yN-global->variable[yl])
		 -global->lhs_limit);
  ret=
    ret*(pow(global->mix_theta-global->Vzhi
	     /(global->u-global->w_mass[squared]),2)
	 *global->letters[A]
	 +pow(global->mix_theta-global->Vzhi
	      /(global->variable[t]-global->w_mass[squared]),2)
	 *global->letters[B]
	 +pow(global->mixing/(global->sqrt_s[squared]-global->z_mass[squared]),2)
	 *pow(tan(global->theta_w),4)*global->letters[C]);
  
  //cross section done
  //now multiply by H=integral of h, and divide by h
  ret*=hmult();
  return 3.894e8*ret;//unit conversion
}

double make_nr2_cross()
{//B+G pg 353 note that V=1
  assert(global->current_process==nr2);
  double ret=global->first_term;
  if (global->g_1)        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
    ret=ret/sqrt(cosh(global->yN+global->variable[yl])
		 -global->lhs_limit);
  ret=
    ret*(pow(global->mix_theta-global->Vzhi
	     /(global->variable[t]-global->w_mass[squared]),2)
	 *global->letters[A]
	 +pow(global->mix_theta-global->Vzhi
	      /(global->u-global->w_mass[squared]),2)
	 *global->letters[B]
	 +pow(global->mixing/(global->sqrt_s[squared]-global->z_mass[squared]),2)
	 *pow(tan(global->theta_w),4)*global->letters[C]);
  
  //cross section done
  //now multiply by H=integral of h, and divide by h
  ret*=hmult();
  return 3.894e8*ret;//unit conversion
}

double make_kinematic_min_yl() 
{//B+G 354
  return global->yN-acosh(global->rhs_limit);
}

double make_kinematic_max_yl() 
{//B+G 354
  double temp;
  temp=global->yN;
  temp-=acosh(global->lhs_limit);
  return temp;
}

double make_kinematic_max_plt() 
{//B+G pg 354
  return 
    .5*(global->majorana_mass[squared]-global->w_mass[squared])
    *global->sqrt_s[value]/global->majorana_mass[squared];
}

double make_kinematic_min_t() 
{
  return (global->sqrt_s[squared]-global->majorana_mass[squared]
	  +sqrt(pow(global->sqrt_s[squared]-
		    global->majorana_mass[squared],2)
		-4*global->sqrt_s[squared]*pow(global->kNt_min,2)))/-2;
}

double make_kinematic_min_full_t() 
{
  return -(global->sqrt_s[squared]-global->majorana_mass[squared]);
}

double make_kinematic_max_t() 
{
  return (global->sqrt_s[squared]-global->majorana_mass[squared]
	  -sqrt(pow(global->sqrt_s[squared]-
		    global->majorana_mass[squared],2)
		-4*global->sqrt_s[squared]*pow(global->kNt_min,2)))/-2;
}

double make_br12() 
{//MAJOR line: 209
  return global->w_mass[gamma]
    /global->full_gamma;
}

double make_br3() 
{//MAJOR line: 209
  return 1-global->branching_ratio[nr1]-global->branching_ratio[nr2];
}

double makeAt()
{//Calculated
//MAJOR line: 217 valid for e+e-?
  return -6.93147e-5;
}

double makeAyl() 
{//MAJOR line: 222
  return 4.4
    +15.744*global->majorana_mass[squared]/
    global->sqrt_s[squared];
}

double makeAplt() 
{//MAJOR line:233
  return 46*global->majorana_mass[squared]
    /(sqrt(global->sqrt_s[value])*pow(global->w_mass[maj2_minus_val2],2));
}

double make_A_time()
{
  return global->full_gamma*global->majorana_mass[value]/global->En3;
}

double make_s() 
{//s = 4 E ^2  B+G 350
  return global->in_lepton_energy[squared]*4;
}

double make_sqrt_s() 
{// sqrt(s) = 2 E B+G 350
  double ret=(global->in_lepton_energy[value])*2;
  return ret;
}

//the remainder of these functions are general utilities that rely on depends
//instead of in.
double make_gamma_z()
{//equation 12 in B+G
  return fermi*pow(global->mixing*global->z_mass[maj2_minus_val2],2)
    *global->z_mass[maj2_plus_2val2]
    /(8*sqrt(2)*PI*pow(global->majorana_mass[value],3));
}

double make_full_gamma()
{
  return 2*global->w_mass[gamma]*global->z_mass[gamma];
}

double make_gamma_w()
{//equation 11 in B+G
  return fermi*pow(global->Vzhi*global->w_mass[maj2_minus_val2],2)
    *global->w_mass[maj2_plus_2val2]
    /(8*sqrt(2)*PI*pow(global->majorana_mass[value],3));
}

double make_w_square()
{
  return global->w_mass[value]*global->w_mass[value];
}

double make_z_square()
{
  return global->z_mass[value]*global->z_mass[value];
}

double make_major_square()
{
  return (global->majorana_mass)[value]*(global->majorana_mass)[value];
}

double make_ilm()
{
  return sqrt(global->in_lepton_momentum[squared]);
}

double make_ilms()
{
  return global->in_lepton_energy[squared]-global->lepton_mass[squared];
}

double make_iles()
{
  return global->in_lepton_energy[value]*global->in_lepton_energy[value];
}

double make_ilmas()
{
  return global->lepton_mass[value]*global->lepton_mass[value];
}

double make_m2p2w2()
{
  return global->majorana_mass[squared]+2*global->w_mass[squared];
}

double make_m2p2z2()
{
  return global->majorana_mass[squared]+2*global->z_mass[squared];
}

double make_m2mw2()
{
  return global->majorana_mass[squared]-global->w_mass[squared];
}
double make_m2mz2()
{
  return global->majorana_mass[squared]-global->z_mass[squared];
}



Back to Source File Index


C++ to HTML Conversion by ctoohtml