// Author: Ricardo Eusebi, FNAL
// phone: (630) 840 4326 
// email: eusebi@fnal.gov
// Created: sometime in 2004
// This code is in the process of being PUBLISHED.
// DO NOT PUBLISH OR DISTRIBUTE WITHOUT THE CONSENT OF RICARDO EUSEBI

#ifndef TOPANDHIGGSCOMMON_DEF
#define TOPANDHIGGSCOMMON_DEF
#include "iostream.h"

// -------------------------------------------------------------------
//      Some stuff common to Carena's and Brendan papers
// -------------------------------------------------------------------


Double_t 
Inte(Double_t a, Double_t b, Double_t c){

  //Check troubling  cases first
  if (fabs(a)==fabs(b) && fabs(c) ==0) { /* cout<<"a=b; c=0"<<endl;*/ return 1.0/(a*a);}
  if (fabs(a)==fabs(c) && fabs(b) ==0) { /* cout<<"a=c; b=0"<<endl;*/ return 1.0/(a*a);}
  if (fabs(b)==fabs(c) && fabs(a) ==0) { /* cout<<"b=c; a=0"<<endl;*/ return 1.0/(b*b);}
  if (fabs(a)==fabs(b) && fabs(b) ==fabs(c)) { /* cout<<"a=b=c"<<endl;   */ return 0.5/(a*a);}

  if (fabs(a)==fabs(b) && c!=0){
    //cout<<"a=b; c!=0"<<endl;
    Double_t b2mc2 = b*b-c*c;
    Double_t b2mc22 = pow(b2mc2,2);
    return (c*c * log(pow(c/b,2))/b2mc22 + 1.0/b2mc2);
  }else if (fabs(a)==fabs(c) && b!=0){
    //cout<<"a=b; c!=0"<<endl;
    Double_t c2mb2 = c*c-b*b;
    Double_t c2mb22 = pow(c2mb2,2);
    return (b*b * log(pow(b/c,2))/c2mb22 + 1.0/c2mb2);
  }else if (fabs(b)==fabs(c) && a!=0){
    //cout<<"a=b; c!=0"<<endl;
    Double_t b2ma2 = b*b-a*a;
    Double_t b2ma22 = pow(b2ma2,2);
    return (a*a * log(pow(a/b,2))/b2ma22 + 1.0/b2ma2);
  }


  Double_t Num,Den;
  Num  = 2*pow(a*b,2)*log(fabs(a/b));
  Num += 2*pow(b*c,2)*log(fabs(b/c));
  Num += 2*pow(c*a,2)*log(fabs(c/a));
  Den = (a*a-b*b)*(b*b-c*c)*(a*a-c*c);
 
  return Num/Den;

}//Inte

Double_t  
Lambda(Double_t a, Double_t b, Double_t c){

  Double_t res = a*a + b*b +  c*c;
  res -= 2*(a*b + a*c+ c*b);
  return res;

}// Lambda

//Calculate the beta0 from page 62, Eq 22. CDFNote 4295
Double_t  beta0(Double_t nf){  return (33.0 - 2.0*nf)/3.0; /* 3.0 */ } //beta0

//Calculate the Delta from page 62, Eq 23. CDFNote 4295
Double_t  DeltaFromalphaZ (Double_t nf){

  Double_t alpha_s = 0.119; //0.119
  Double_t mu2   = pow(91.19,2);
  Double_t expo  = 4* PI /(alpha_s*beta0(nf));
  
  return sqrt(mu2*exp(-expo));

}//Delta

//Calculate the alpha strong according to the usual formulae
Double_t alpha_strong(Int_t nf, Double_t Q){

  Double_t res1 = 2. * beta0(nf) * log(Q/DeltaFromalphaZ(nf));

  return 4 * TMath::Pi() /res1;

}//alpha_strong




//Taken from Carena's paper hep-ph 9912516 , page 5
// Note that it is ill defined when the parameters have the same magnitude
//Calculate the Evolution Factor. See Carena's paper page 9, Eq 19.
Double_t EvolutionFactor(Int_t nf, Double_t Q2, Double_t Q1){

  Double_t alphaQ1 = alpha_strong(nf,Q1);
  Double_t alphaQ2 = alpha_strong(nf,Q2);
  Double_t df = 4./beta0(nf); // == 12/(33-2*nf)

  Double_t JfNum = - (8982. - 504*nf  + 40 * nf * nf);
  Double_t JfDen = 3 * pow ( 33-2*nf  ,2);
  Double_t Jf = JfNum/JfDen;

  Double_t res1 = pow(alphaQ2/alphaQ1,df);
  Double_t res2 = 1 + (alphaQ1-alphaQ2)* Jf /(4*TMath::Pi());

  return res1 * res2;

}//EvolutionFactor

//Make the running b mass
Double_t 
Running_bquark_Mass(Double_t Q){

  if (Q>=4.25)
    return EvolutionFactor(5,Q,4.25) * 4.25;

  //Not sure how to define it for low Q.
  return 0.0;

}//Running_bquark_Mass



//Make the running b mass
Double_t 
Running_bquark_Mass_test(Double_t Q){

  if (Q>=4.25)
    return EvolutionFactor(6,Q*Q,175) * EvolutionFactor(5,175,4.25) * 4.25;

  //Not sure how to define it for low Q.
  return 0.0;

}//Running_bquark_Mass





#endif