file: ../SURF/event.h


#include "parameter.h"
#include "f2c.h"
#include "LEDA/array2.h"
#include "assert.h"

//created by John Langford (jl@crush.caltech.edu)
//This class's role is producing an event involving Heavy Majorana Neutrinos
//using monte carlo methodology.
//It has many modifiable parameters.

//Because this class depends on the fortran pythia libraries it is nonreintrant

//these are defined to make it impossible to set an invalid parameter
typedef enum {active,inactive} final_state_type;
typedef enum {nr1,nr2,production,decay,nr3,mixed123} process_type;

typedef enum {cnst,v_inv,v_inv_v_inv,exp_neg_v_Av,v_exp_neg_v_v_Av,
omega_to_yl} h_type;
//used to increase efficiency of calculations.  part of the general monte carlo
//method

//defined for code readability.
typedef enum {plt,t,yl,full_t,time} var_type;

const int value=0;//Several parameters are stored in arrays with the first 
const int squared=1;//value = val and second value = val^2

const int maj2_minus_val2=2;//more constants defined to access regular systems
const int maj2_plus_2val2=3;//of equations
const int gamma=4;

const int A=0;//A,B,C stand for moderately large equations
const int B=1;
const int C=2;

const int num_params=6;//plt,t,yl,phi,time (to decay)

const int num_procs=5;//nr1,nr2,nr3.  Or,in more readable terms: N->e- + W+
const int num_real_procs=4;
//N->e+ + W- and N production then decay N->v+Z

const double fermi=.0000116637;

class event {
private:
  //internal things
  bool initial;
  parameter QQ;         //in GeV (virtual mass ^ 2)
  parameter QQ3;         //in GeV (virtual mass ^ 2)
  int kc;               //jetset's internal code
  parameter g_1;        //used to tell which side of the singularity you are on
  array2<var_type> process_vars;
  array<int> process_var_count;
  array<double> g_sum;  //sum of all g calls
  array<double> gg_sum; //sum of all g calls squared
  
  //parameters modifiable but which cause a recalculation of g_max
  array<parameter> in_lepton_momentum;          //in Gev, and GeV^2
  array<parameter> in_lepton_energy;
  array<parameter> lepton_mass;
  array<parameter> majorana_mass; //in GeV
  process_type process_number;
  double safety_factor;//multiplied to maximum of G found
  array<h_type> h_function; //plt,t,yl
  double cut_size; //distance from singularity
  //the following have 3rd and 4th entry = maj2_minus_<2nd>, maj2_plus_2<2nd>
  //5th entry=gamma<1st>
  array<parameter> z_mass;
  array<parameter> w_mass;
  parameter full_gamma;
  
  array<parameter> A_value;  //{Ayl,At,Apl} constants used to tune monte carlo
  parameter theta_w;         //weak angle
  parameter mixing;          //zhi_nu_N in B+G  assumed to be .1 for now
  parameter Vzhi;            //Vzhi_e_N * V_e_v in B+G assumed to be .1 for now

  void reset_pythia();//set pythia variables to 0's
  void reinitialize();//called when an event is about to occur if anything has         //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
                      // changed
  void maximizer();   //finds maximum of g function's
  friend void gwrap(integer *npar,doublereal *deriv,doublereal *difsig,
		    doublereal *xf,integer *nflag,doublereal *futil);
                      //required by the minuit routines
  void getpoint();
  //gets a phase space point of the variables

  void event::set_singular_var();//sets yl and plt to be along the singularity

  //function which calculate dependent variables
  friend double hmult();
  friend double square(); 
  friend double makeu();
  friend double makeu3();
  friend double make_out_lepton_pl();
  friend double make_out_lepton_energy();
  friend double makeQQ();
  friend double makeQQ3();
  friend double make_mix_theta( );
  friend double makekNt();
  friend double make_lhs_limit();
  friend double make_rhs_limit();
  friend double makeA();
  friend double makeB();
  friend double makeC();
  friend double make_first_term();
  friend double make_nr1_cross();
  friend double make_nr2_cross();
  friend double make_direction();
  friend double make_decay_time_cross();
  friend double make_majorana_production_cross();
  friend double makemt();
  friend double makeEn();
  friend double make_En3();
  friend double make_Enu3();
  friend double make_Pnt3();
  friend double makeyN();
  friend double makekNl();
  friend double make_E_nu();
  friend double make_E_nu_out();
  friend double make_kNt_min();
  friend double make_W_mass();
  friend double times();
  friend double make_br12();
  friend double plus2();
  friend double minus();
  friend double makeAt();
  friend double makeAyl();
  friend double makeAplt();
  friend double make_A_time();
  friend double make_kinematic_min_yl(); 
  friend double make_kinematic_max_yl(); 
  friend double make_kinematic_max_plt(); 
  friend double make_kinematic_max_t(); 
  friend double make_kinematic_min_full_t(); 
  friend double make_kinematic_min_t(); 
  friend double make_br3(); 
  friend double make_s();
  friend double make_sqrt_s();
  friend double	make_gamma_z();
  friend double	make_full_gamma();
  friend double make_gamma_w();
  friend double make_w_square();
  friend double make_z_square();
  friend double make_major_square();
  friend double make_ilm();
  friend double make_ilms();
  friend double make_iles();
  friend double make_ilmas();
  friend double make_m2p2w2();
  friend double make_m2p2z2();
  friend double make_m2mw2();
  friend double make_m2mz2();
public:
  //Variables that you may read but probably shouldn't set (used by system)
  array<double> process_cross_section; //{nr1,nr2,nr3} in pb
  array<double> process_cross_section_std_dev; //{nr1,nr2,nr3} in pb
  array<double> generation_efficiency; //{nr1,nr2,nr3}
  array2<double> g_function_max;//{nr1,nr2,nr3} maximum of g functions found
  int cut_count;//the number of times the cut about the singularity is 
  // violated.
  array<int> accepted;//number of accepted events
  array<int> events;//number of events
  array<int> attempted_events;//The number of attempts to generate events 
  //(many are invalid)
  parameter out_lepton_pl;//out lepton variables
  parameter out_lepton_energy;
  int numover;
  //lots of dependent parameters.  For their definition it is probably 
  //best to look in friends.cc.  Names correspond to that in B + G
  array<parameter> sqrt_s;
  array<parameter> letters;//A,B,C in B+G
  parameter first_term;
  parameter u;
  parameter u3;
  parameter lhs_limit;
  parameter rhs_limit;
  parameter mix_theta;
  parameter kNl;
  parameter mt;
  parameter kNt;
  parameter kNt_min;
  parameter En;
  parameter En3;
  parameter Enu3;
  parameter E_nu_out;
  parameter Pnt3;
  parameter E_nu;
  parameter yN;
  array<parameter> kinematic_min;//{plt,t,yl}
  array<parameter> kinematic_max;
  array<parameter> branching_ratio;
  parameter direction;
  array<parameter> cross_section;
  parameter current_process;
  
  //Variables you may set to change functionality
  //cuts
  array<double> cut_min;//{plt,t,yl}  General cuts on the cross-sections
  array<double> cut_max;
  //process control variables
  final_state_type final_state;
  array<parameter> variable;//plt,t,yl
  bool hadronisation;
  void in_lepton_E_set(double in)
       {in_lepton_energy[value].setvalue(in);initial=false;};
  double in_lepton_E_read(){return in_lepton_energy[value];};
  void lepton_mass_set(double in)
       {lepton_mass[value].setvalue(in);initial=false;};
  double lepton_mass_read(){return lepton_mass[value];};
  void major_mass_set(double in)
       {majorana_mass[value].setvalue(in);initial=false;};
  double major_mass_read(){return majorana_mass[value];};
  void process_number_set(process_type in)
       {process_number=in;current_process.setvalue(in);initial=false;};
  int process_number_read(){return process_number;};
  void safety_factor_set(double in)
  {
    int i; 
    assert(in>=1);//a lower safety factor simply doesn't make sense.
    for (i=0;i<num_procs;i++){
      if (g_function_max(i,true)!=0)        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
	g_function_max(i,true)*=in/safety_factor;
      if (g_function_max(i,false)!=0)        //OVERLOAD CALL: if: parameter.cc(?), parameter.cc(?)
	g_function_max(i,false)*=in/safety_factor;
    }
    safety_factor=in;
  };
  void h_function_plt_set(h_type in){h_function[plt]=in;initial=false;};
  h_type h_function_plt_read(){return h_function[plt];};
  void h_function_t_set(h_type in){h_function[t]=in;initial=false;};
  h_type h_function_t_read(){return h_function[t];};
  double safety_factor_read(){return safety_factor;};
  void cut_size_set(double in){cut_size=in;initial=false;};
  double cut_size_read(){return cut_size;};
  void A_value_plt_set(double in){A_value[plt].setvalue(in);initial=false;};
  double A_value_plt_read(){return A_value[plt];};
  void A_value_t_set(double in){A_value[t].setvalue(in);initial=false;};
  double A_value_t_read(){return A_value[t];};
  event();
  void print();        //OVERLOAD CALL: print: event.cc(event), parameter.cc(parameter)
  void occur();
};

//used by all of the friend functions.  Actual location is in friend.cc
extern event* global;



Back to Source File Index


C++ to HTML Conversion by ctoohtml