intTwoHDM::whichinput=1;//which type of input you will provide (1=hadronic)
TwoHDM::TwoHDM():MassBeta(2){
//CONSTRUCTOR
alphaSM=-atan(1);
bsgamma_max=4.1;
bsgamma_min=3.;
Ecm2=2000.;
Ecm7=7000.;
Ecm8=8000.;
#ifdef HBMODE
HiggsBoundsAnalysis.rewrite_values(nHzero,nHplus,whichanalyses,whichinput);//sets the values needed to initialize HiggsBounds (will be initialized in UserInitCalcs)
#endif
//Loading numerical tables for cross-sections
LoadAll2HDMTables();
///////////////////////////////
//Initialise arrays to zero
//Higgs couplings to vector bosons
for(inti=0;i<nHzero;++i){
Rhvv[i]=0;
for(intj=0;j<nHzero;++j)
Rhhv[i][j]=0;
for(intt=0;t<4;++t)
for(intf=0;f<3;++f)
Rhf[i][f]=0;
}
for(inti=0;i!=8;++i){
mass_bounds[i]=false;
L_bounds[i]=false;
phi_bounds[i]=true;
}
phi_bounds[2]=false;
phi_bounds[6]=false;
}
voidTwoHDM::readfile(std::istream&in){
//This function is called from the overloaded operator >> when a TwoHDM input file is provided.
//On the base class, Potential::readfile simply calls openinput() on "in".
//Since TwoHDM files only contain the relevante variables, "this" readfile opens an auxiliary file in the same format has Potential::readfile
//and then replaces the relevante variables bty the values in "in".
//check of all independent parameters were defined in the input file
for(inti=0;i!=5;++i)
if(!mass_bounds[i]){
cerr<<"ScannerS ERROR: Unable to find a necessary parameter in the 2HDM input file. This model ALWAYS requires you define the range for mH, mh, mA, mH+, tanbeta and m12^2. These are independent parameters that will be randomly generated."<<endl;
exit(EXIT_FAILURE);
}
if((!L_bounds[7])||(!tanbeta_bounds)){
cerr<<"ScannerS ERROR: Unable to find a necessary parameter in the 2HDM input file. This model ALWAYS requires you define the range for mH, mh, mA, mH+, tanbeta and m12^2. These are independent parameters that will be randomly generated."<<endl;
exit(EXIT_FAILURE);
}
//inert model
if(PhiPar.Vecmin[0]==0&&PhiPar.Vecmax[0]==0){
if(alpha_bounds){
cerr<<"ScannerS ERROR: While using the inert model (tanbeta==0) there are no mixing directions so alpha is not defined."<<endl;
exit(EXIT_FAILURE);
}
if(Cp.L.Vecmin[7]!=0||Cp.L.Vecmax[7]!=0){
cerr<<"ScannerS ERROR: In the inert model (tanbeta==0), m12^2 must also be equal to 0."<<endl;
exit(EXIT_FAILURE);
}
if((!L_bounds[4])||(!L_bounds[5])){
cerr<<"ScannerS ERROR: Unable to find a necessary parameter in the 2HDM input file. In the inert model (tanbeta==0) you must define lamdba2 and lamdba3. These are independent parameters that will be randomly generated."<<endl;
exit(EXIT_FAILURE);
}
//In the inert model there is no mixing so if tanbeta==0, there is no alpha
MixingPars.Vec.resize(0);
MixingPars.Vecmin.resize(0);
MixingPars.Vecmax.resize(0);
}
//non-inert models
else{
if(!alpha_bounds){
cerr<<"ScannerS ERROR: Unable to find a necessary parameter in the 2HDM input file. Non-inert models require you to define the mixing angle alpha. This is an independent parameter that will be randomly generated."<<endl;
exit(EXIT_FAILURE);
}
if(Cp.Mass.Vecmax[1]>Cp.Mass.Vecmin[0]){
std::clog<<endl<<"ScannerS Warning: the mass range for the light Higgs (mh) is superior or overlaps that of the heavy Higgs (mH)."<<endl;
std::clog<<"-> Resetting internal alpha to generate physical alpha in [-Pi/2:Pi/2] to cover all possible cases and so that ordering is consistent with input file. "<<endl;
std::clog<<"-> Do you wish to continue? [Y/N]"<<endl;;
stringchoice;
do{
cin>>choice;
if(choice=="N"){
std::clog<<"Scan aborted by the user!"<<endl;
exit(EXIT_FAILURE);
}
elseif(choice!="Y")
std::clog<<"Please type exacly Y or N."<<endl;
}while(choice!="Y");
//resetting alpha ranges to [0;2Pi] if the user wished to continue the run
std::clog<<endl<<"The ranges of parameters to be scanned are:"<<endl;
std::clog<<endl<<"Ranges for vacuum expectation values:"<<endl;
for(size_ti=0;i!=Phi.size();++i){
if(phi_bounds[i])
std::clog<<"phi"<<i<<" in ["<<Phi.vgetmin(i)<<":"<<Phi.vgetmax(i)<<"]"<<endl;
}
std::clog<<endl<<"Ranges for couplings:"<<endl;
for(size_ti=0;i!=Cp.L.size();++i){
if(L_bounds[i])
std::clog<<"L"<<i<<" in ["<<Cp.L.vgetmin(i)<<":"<<Cp.L.vgetmax(i)<<"]"<<endl;
}
std::clog<<endl<<"Ranges for physical state masses:"<<endl;
for(size_ti=0;i!=Cp.Mass.size();++i){
if(mass_bounds[i])
std::clog<<"m"<<i<<" in ["<<Cp.Mass.vgetmin(i)<<":"<<Cp.Mass.vgetmax(i)<<"]"<<endl;
}
std::clog<<endl<<"Range for alpha (mixing parameterization):"<<endl;
std::clog<<"alpha in ["<<MixingPars.Vecmin[0]<<":"<<MixingPars.Vecmax[0]<<"]"<<endl;
std::clog<<endl<<"Range for tanbeta (v2/v6 parameterization):"<<endl;
std::clog<<"tanbeta in ["<<PhiPar.Vecmin[0]<<":"<<PhiPar.Vecmax[0]<<"]"<<endl;
std::clog<<endl;
std::clog<<"The seed you are using for the random number generator is: "<<r.getseed()<<endl<<endl;
}
boolTwoHDM::CheckBounds(){
//unlike Potential::CheckBounds(), this only checks dependent parameters for which the range was explicity declared in the 2HDM input file. if they were not declared in input file (and were loaded of 2HDM.potential by default) L_bounds will be false and the parameters will not be checked
//The only dependend parameters (noticing that in this class we use paramerizations for both the mixing and the VeV angle)
//WARNING: This routine assumes the two scalars that get vevs in the decompositions of the two doublets are phi[2] and phi[6]!
voidTwoHDM::PhysicalStates2HDM(){
//This routines returns the masses of the scalars in the 2HDM as well as the alpha mixing angle for the 2 CP even Higgs bosons and the beta angle between the vevs.
// INPUT:
// Mass - internal (vector like) array of masses
// Mixing - internal (matrix like) mixing matrix
// Phi - internal (vector like) array of field values
// OUTPUT:
// mHheavy - Mass of the heavy CP even Higgs
// mHlight - Mass of the light CP even Higgs
// mHc - Mass of the charged Higgs
// mA - Mass of the CP odd Higgs
// alpha - mixing angle for the 2 CP even Higgs bosons
// beta - angle between the vevs
mA=Cp.Mass[2];//Mass of the pseudo scalar (numbering checked in VERBOSE mode)
mHcharged=Cp.Mass[3];//Mass of the charged Higgs (numbering checked in VERBOSE mode. Note Mass[4] is the same because it corresponds to the imaginary part).
//This routines returns the masses of the scalars in the 2HDM as well as the alpha mixing angle for the 2 CP even Higgs bosons and the beta angle between the vevs.
/// INPUT:
// Mass - internal (vector like) array of masses
// Mixing - internal (matrix like) mixing matrix
// Phi - internal (vector like) array of field values
/// OUTPUT:
// mHheavy - Mass of the heavy CP even Higgs
// mhlight - Mass of the light CP even Higgs
// mHc - Mass of the charged Higgs
// mA - Mass of the CP odd Higgs
// alpha - mixing angle for the 2 CP even Higgs bosons
// beta - angle between the vevs
mA=Mass[2];//Mass of the pseudo scalar (numbering checked in VERBOSE mode)
mHc=Mass[3];//Mass of the charged Higgs (numbering checked in VERBOSE mode. Note Mass[4] is the same because it corresponds to the imaginary part).
std::cerr<<"Error in double InternalTableSusHiPP(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Tables for Ecm = "<<Ecm<<" not available. Ecm must be in [7000:14000] GeV"<<endl;
std::cerr<<"Error in double InternalTableSusHiPP(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Tables for mHparticle = "<<mHparticle<<" not available. Allowed range is [30:1000] GeV."<<endl;
std::cerr<<"Error in double InternalTableSusHiPP(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Tables for mHparticle = "<<mHparticle<<" not available. Allowed range is [30:1000] GeV."<<endl;
exit(EXIT_FAILURE);
}
}
doublegT,gB;
doublecosbeta=1/sqrt(1+tanbeta*tanbeta);
doublesinbeta=tanbeta*cosbeta;
if(Type2HDM==1||Type2HDM==3){
if(Hparticle==0){
gT=cos(alpha)/sinbeta;
gB=gT;
}
elseif(Hparticle==1){
gT=1./tanbeta;
gB=-gT;
}
elseif(Hparticle==2){
gT=sin(alpha)/sinbeta;
gB=gT;
}
else{
std::cerr<<"Error in double InternalTableSusHiPP(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Bad Hparticle = "<<Hparticle<<". Allowed values are 0,1,2 for h, A or H."<<endl;
exit(EXIT_FAILURE);
}
}
elseif(Type2HDM==2||Type2HDM==4){
if(Hparticle==0){
gT=cos(alpha)/sinbeta;
gB=-sin(alpha)/cosbeta;
}
elseif(Hparticle==1){
gT=1./tanbeta;
gB=tanbeta;
}
elseif(Hparticle==2){
gT=sin(alpha)/sinbeta;
gB=cos(alpha)/cosbeta;
}
else{
std::cerr<<"Error in double InternalTableSusHiPP(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Bad Hparticle = "<<Hparticle<<". Allowed values are 0,1,2 for h, A or H."<<endl;
exit(EXIT_FAILURE);
}
}
else{
std::cerr<<"Error in double InternalTableSusHiPP(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Bad Type2HDM = "<<Type2HDM<<". Allowed values are 1,2,3,4."<<endl;
std::cerr<<"Error in double InternalTableSusHiPPbar(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Tables for Ecm = "<<Ecm<<" not available. Ecm must be in [500:2000] GeV"<<endl;
std::cerr<<"Error in double InternalTableSusHiPPbar(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Tables for mHparticle = "<<mHparticle<<" not available. Allowed range is [30:1000] GeV."<<endl;
std::cerr<<"Error in double InternalTableSusHiPPbar(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Tables for mHparticle = "<<mHparticle<<" not available. Allowed range is [30:1000] GeV."<<endl;
exit(EXIT_FAILURE);
}
}
doublegT,gB;
doublecosbeta=1/sqrt(1+tanbeta*tanbeta);
doublesinbeta=tanbeta*cosbeta;
if(Type2HDM==1||Type2HDM==3){
if(Hparticle==0){
gT=cos(alpha)/sinbeta;
gB=gT;
}
elseif(Hparticle==1){
gT=1./tanbeta;
gB=-gT;
}
elseif(Hparticle==2){
gT=sin(alpha)/sinbeta;
gB=gT;
}
else{
std::cerr<<"Error in double InternalTableSusHiPPbar(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Bad Hparticle = "<<Hparticle<<". Allowed values are 0,1,2 for h, A or H."<<endl;
exit(EXIT_FAILURE);
}
}
elseif(Type2HDM==2||Type2HDM==4){
if(Hparticle==0){
gT=cos(alpha)/sinbeta;
gB=-sin(alpha)/cosbeta;
}
elseif(Hparticle==1){
gT=1./tanbeta;
gB=tanbeta;
}
elseif(Hparticle==2){
gT=sin(alpha)/sinbeta;
gB=cos(alpha)/cosbeta;
}
else{
std::cerr<<"Error in double InternalTableSusHiPPbar(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Bad Hparticle = "<<Hparticle<<". Allowed values are 0,1,2 for h, A or H."<<endl;
exit(EXIT_FAILURE);
}
}
else{
std::cerr<<"Error in double InternalTableSusHiPPbar(double Ecm, double mHparticle, double alpha,double tanbeta,int Hparticle,int Type2HDM): Bad Type2HDM = "<<Type2HDM<<". Allowed values are 1,2,3,4."<<endl;