#ifndef PANDPYTHCONSTS_H
#define PANDPYTHCONSTS_H

//for pandorarun beam type...
const int epluseminus  = 1;
const int gammagamma   = 2;
const int eminuseminus = 3;

//copied from hepdata.h from pandora....

const double PI = 3.14159265358979323846;
const double TINY =  1.0e-30;


template<class T>
inline
T SQR(T a){return a*a;}

//double ABS(double a); 
//int  ABS(int a);

//double MIN(double a,double b);
//double  MAX(double a, double b);
//double SIGN(double a, double b);

//int MIN(int a, int b);
//int  MAX(int a, int b);
//int SIGN(int a, int b);

//int INT(double a);

//void therror(char error_text[]);

const double alpha0 = 1.0/137.036;
const double alpha = 1/128.14;
const double alphaw = 1/29.60;
const double sstw = alpha/alphaw;
const double cstw = 1.0 - sstw;
const double mw = 80.4;
const double mz = 91.187;
     /*  kinematic mass values  */
const double mt = 175.0;
const double mb = 4.2;
const double mc = 1.2;
const double mtau = 1.777;
const double me = 0.510999e-3;
      /*  m_MSbar at Q = 2 mt   */
const double mtMSb = 159.0;
const double mbMSb = 2.68;
const double mcMSb = 0.54;
/*  for renormalization group running, switch b0 at 2 mtMSb:  */
inline double bzero(double Q){ return  (Q < 2.0*mtMSb) ?  7.333 : 7.0 ; }
const double alphasZ = 0.119;
const double alphas2mt =
	alphasZ/(1.0 + (7.33/(2.0*PI))*alphasZ*log(2.0*mtMSb/mz));
inline double alphas(double Q){
   return   alphas2mt/(1.0 + (bzero(Q)/(2.0*PI))*alphas2mt*log(Q/(2.0*mtMSb)));
}

// conversion from (GeV)^-2 to pb,fb:
const double pb  = 0.389380e9;
const double fb  = 0.389380e12;

// conversion from GeV^-1 to fm:
const double fm = 0.197327053;
const double cm = 0.197327053e-13;

// widths of t, W, Z:

const double WGqfactor = 3.0 *  (1.0 + alphas(sqrt(mw))/PI
                                  + 1.41 * SQR(alphas(sqrt(mw))/PI) );
const double ZGqfactor = 3.0 *  (1.0 + alphas(sqrt(mz))/PI
                                  + 1.41 * SQR(alphas(sqrt(mz))/PI) );
const double tGqfactor =  (1.0  - (0.866 * alphas(mt) 
                                         + 1.69 * SQR(alphas(mt))));
                         /* Czarnecki and Melnikov  hep-ph/9806244 */
const double NfW = (3.0 + 2.0 * WGqfactor);
const double NfZ = 3.0 *  SQR(0.5) + 3.0 *  SQR(0.5 - sstw) +  3.0 * SQR(sstw)
    + 2.0 *  (SQR(0.5 - (2.0/3.0)*sstw) + SQR((2.0/3.0)*sstw))* ZGqfactor 
  +(3.0 - 0.017)*(SQR(0.5 - (1.0/3.0)*sstw) + SQR((1.0/3.0)*sstw))* ZGqfactor;
const double GammaW = 2.08;
const double GammaZ = 2.496;
const double Gammat =  (alphaw* mt/16.0) * SQR(1.0 - SQR(mw/mt))
                                    *( 2.0 + SQR(mt/mw)) * tGqfactor;

/*  use the following function to compute the Breit-Wigner distributions
       for particles with widths :    0 < x < 1   */
inline double XTOMS(double x, double M, double G){
                return M*M + M*G * tan((PI- 2.0*G/M)*(x - 0.5)); }
/*  the following minimum values are also useful :  */
const double MINmw = sqrt(XTOMS(0.0,mw,GammaW));
const double MINmz = sqrt(XTOMS(0.0,mz,GammaZ));
const double MINmt = sqrt(XTOMS(0.0,mt,Gammat));
const double SMALLWIDTH = 1.0e-9;


/*  ID variables for beam particles   */

const int electron = 11;
const int positron = -11;
const int photon = 22;
 
/*  numerical constants   */


const double third = 1.0/3.0;
const double tthird = 2.0/3.0;


/* `Only' variables used to restrict decay patterns    */

const int allStates = 546;
const int leptonsOnly = 547;
const int quarksOnly =  548;
const int lightquarksOnly = 549;
const int heavyquarksOnly = 550;
const int bcOnly = 551;
const int cOnly = 552;
const int bOnly = 553;
const int eOnly = 554;
const int muOnly = 555;
const int tauOnly = 556;
const int invisibleOnly = 557;

const int glueOnly = 660;;
const int cglueOnly = 661;
const int tOnly = 662;
const int bcglueOnly = 663;
const int WWOnly = 664;
const int ZZOnly = 665;
const int WZOnly = 666;
const int ggOnly = 667;
const int gZOnly = 668;

// copied from ebeams.h from pandoraV2.1...
/*    note polarization = +1.0 denotes complete right-handed polarization */
 
const int NLC500 = 500;
const int NLC1000 = 1000;
const int NLC1500 = 1500;

const int NLC500A = 501;
const int NLC500B = 500;
const int NLC500C = 502;
const int NLC500H = 503;
const int NLC1000A = 1001;
const int NLC1000B = 1000;
const int NLC1000C = 1002;
const int NLC1000H = 1003;

const int JLC500 = 500;
const int JLC1000 = 1000;
const int JLC1500 = 1500;

const int JLC500A = 501;
const int JLC500B = 500;
const int JLC500C = 502;
const int JLC500H = 503;
const int JLC1000A = 1001;
const int JLC1000B = 1000;
const int JLC1000C = 1002;
const int JLC1000H = 1003;

const int TESLA500 = 510;
const int TESLA800 = 810;

const int CLIC500 = 520;
const int CLIC1000 = 1020;
const int CLIC3000 = 3020;
const int CLIC5000 = 5020;

const int EE500 = 499;

#endif