Page MenuHomeHEPForge

ScannerSUserTwoHDM.cpp
No OneTemporary

ScannerSUserTwoHDM.cpp

#include "ScannerScore/ScannerSModels/2HDM/ScannerSTwoHDM.h"
using namespace std;
//////////////////////////////////////////////////////////
////// Function contents to be defined by the user ///////
//////////////////////////////////////////////////////////
// WARNING: Don't change the prototypes of the template //
// functions in this file since they are used by the //
// core routines. Only change their contents! //
//////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////
////////////// Member Function of TwoHDM /////////////////
//////////////////////////////////////////////////////////
// These functions will only be called when you use a //
// TwoHDM input file. //
//////////////////////////////////////////////////////////
void TwoHDM::UserInitCalcs(void){
//////////////////////////////////////////////////////////////////////////
// ENTER HERE CODE FOR YOUR INITIAL CALCULATIONS BEFORE THE SCAN STARTS //
//////////////////////////////////////////////////////////////////////////
//// TEST general SUSHI TABLES
//#include "ScannerScore/test_includes/tests_sushi_tables_LHC_noEW.cpp"
//exit(0);
/////// END OF TEST ///
//// Setting up the output print options
//// leave the 1st argument at 1 if you want print points 1 by 1
setprint(50,67);
//// Use this vector to print/keep any extra data you calculate on UserAnalysis
//// Keep it at 0 if you have no use for it
extra_data.resize(6);
#ifdef HBMODE
//// Initialising HiggsBounds!
HiggsBoundsAnalysis.InitialiseHiggsBounds();
#ifdef HSMODE
//// Initialising HiggsSignals!
//int Np=0; //Number of free parameters of the model
initialize_higgssignals_latestresults_(nHzero,nHplus);
//setup_nparam_(Np); // Set number of extra free parameters for the model (no call defaults to 0)
int pdf=2; //choose probability density function of Higgs mass arround the measured central value (HiggsSignals)
setup_pdf_(pdf);
#endif
#endif
}
bool TwoHDM::UserAnalysis(PhiRef & Phi,LambdaRef & L,MassRef & Mass, MmixingRef & Mixing){
////////////////////////////////////////////////////////////////
// ENTER CODE FOR YOUR TESTS/OUTPUT DURING THE SCAN //
////////////////////////////////////////////////////////////////
// User defined analysis function to be executed once for //
// each valid parameter space point that was generated. This //
// routine serves to accept or reject the point for example //
// through testing experimental constraints. //
////////////////////////////////////////////////////////////////
PhysicalStates2HDM();
if(not TestSTU_2sigma_2HDM(mHlight,mHheavy,mA,mHcharged,alpha,beta))
return false;
// 1- Testing low energy stuff
////1.1- Zbbar at LEP
if(not TestZbbarbLEP())
return false;
////1.2 flavour Physics observables
if(not Test_bs_gamma())
return false;
// 2- Collider physics
#ifdef HBMODE
//Current implementation uses Hdecay to compute branching ratios for HiggsBounds. Alternatively you can set the values of HiggsBoundsAnalysis using your own calculations.
#ifdef HDECAYMODE
PrepareHiggsBoundsInput();
//Call to run HiggsBounds
int HBresult,chan, ncombined;
double obsratio;
HiggsBoundsAnalysis.run(HBresult,chan,obsratio,ncombined);
if(HBresult==0)
return false; //Exclude point if excluded from non-observation in the HB search channels.
#ifdef HSMODE //WARNING: This is currently OFF as the interpretation of the HiggSignals output is more difficult (see CheckSignalStrengths_gen below).
//Call to run Higgs Signals
int mode=1; //input
//output variables
double csqmu;
double csqmh;
double csqtot;
int nobs;
run_higgssignals_(mode,csqmu,csqmh,csqtot,nobs,Pvalue);
if(Pvalue<0.0027)
return false; //Exclude point if the observed signal rates is too low (3sigma). //WARNING: This is currently OFF if HSMODE is not set, as the interpretation of the HiggSignals output is more difficult (see CheckSignalStrengths_gen below).
#endif
// Test signal rates for the SM Higgs at the LHC7+8 using the ATLAS+CMS combination
//check at two sigma
constexpr int nsigma = 3;
//with a seperation of 5GeV
constexpr double separation = 5.;
std::map<std::string,double> mu_results;
if(not CheckSignalStrengths_gen(nsigma,separation,HiggsBoundsAnalysis,mu_results))
return false;
#endif
#endif
return true;
}
void TwoHDM::what2print(){
/////////////////////////// TwoHDM example ///////////////////////
// input parameters
points2print[store_vector_index][0]=mHheavy;
points2print[store_vector_index][1]=mHlight;
points2print[store_vector_index][2]=mA;
points2print[store_vector_index][3]=mHcharged;
points2print[store_vector_index][4]=alpha;
points2print[store_vector_index][5]=tanbeta;
//couplings
for(int i=0;i<Cp.L.size();++i)
points2print[store_vector_index][6+i]=Cp.L[i];
points2print[store_vector_index][14]=bsgamma;
#ifdef HBMODE
//lightHiggs - h
//Branching ratios and total width
points2print[store_vector_index][15]=HiggsBoundsAnalysis.BR_hjss[0];
points2print[store_vector_index][16]=HiggsBoundsAnalysis.BR_hjcc[0];
points2print[store_vector_index][17]=HiggsBoundsAnalysis.BR_hjbb[0];
points2print[store_vector_index][18]=HiggsBoundsAnalysis.BR_hjmumu[0];
points2print[store_vector_index][19]=HiggsBoundsAnalysis.BR_hjtautau[0];
points2print[store_vector_index][20]=HiggsBoundsAnalysis.BR_hjWW[0];
points2print[store_vector_index][21]=HiggsBoundsAnalysis.BR_hjZZ[0];
points2print[store_vector_index][22]=HiggsBoundsAnalysis.BR_hjZga[0];
points2print[store_vector_index][23]=HiggsBoundsAnalysis.BR_hjgaga[0];
points2print[store_vector_index][24]=HiggsBoundsAnalysis.BR_hjgg[0];
points2print[store_vector_index][25]=*HiggsBoundsAnalysis.BR_hjhihi[0];
points2print[store_vector_index][26]=HiggsBoundsAnalysis.MhGammaTot[0];
//Pseudo scalar - A
//Branching ratios and total width
points2print[store_vector_index][27]=HiggsBoundsAnalysis.BR_hjss[1];
points2print[store_vector_index][28]=HiggsBoundsAnalysis.BR_hjcc[1];
points2print[store_vector_index][29]=HiggsBoundsAnalysis.BR_hjbb[1];
points2print[store_vector_index][30]=HiggsBoundsAnalysis.BR_hjmumu[1];
points2print[store_vector_index][31]=HiggsBoundsAnalysis.BR_hjtautau[1];
points2print[store_vector_index][32]=HiggsBoundsAnalysis.BR_hjWW[1];
points2print[store_vector_index][33]=HiggsBoundsAnalysis.BR_hjZZ[1];
points2print[store_vector_index][34]=HiggsBoundsAnalysis.BR_hjZga[1];
points2print[store_vector_index][35]=HiggsBoundsAnalysis.BR_hjgaga[1];
points2print[store_vector_index][36]=HiggsBoundsAnalysis.BR_hjgg[1];
points2print[store_vector_index][37]=HiggsBoundsAnalysis.MhGammaTot[1];
//Heavy Higgs -
//Branching ratios and total width
points2print[store_vector_index][38]=HiggsBoundsAnalysis.BR_hjss[2];
points2print[store_vector_index][39]=HiggsBoundsAnalysis.BR_hjcc[2];
points2print[store_vector_index][40]=HiggsBoundsAnalysis.BR_hjbb[2];
points2print[store_vector_index][41]=HiggsBoundsAnalysis.BR_hjmumu[2];
points2print[store_vector_index][42]=HiggsBoundsAnalysis.BR_hjtautau[2];
points2print[store_vector_index][43]=HiggsBoundsAnalysis.BR_hjWW[2];
points2print[store_vector_index][44]=HiggsBoundsAnalysis.BR_hjZZ[2];
points2print[store_vector_index][45]=HiggsBoundsAnalysis.BR_hjZga[2];
points2print[store_vector_index][46]=HiggsBoundsAnalysis.BR_hjgaga[2];
points2print[store_vector_index][47]=HiggsBoundsAnalysis.BR_hjgg[2];
points2print[store_vector_index][48]=HiggsBoundsAnalysis.BR_hjhihi[2][0];
points2print[store_vector_index][49]=HiggsBoundsAnalysis.BR_hjhihi[2][1];
points2print[store_vector_index][50]=HiggsBoundsAnalysis.MhGammaTot[2];
//Charged Higgs - H^+
//Branching ratios and total width
points2print[store_vector_index][51]=mHcharged;
points2print[store_vector_index][52]=HiggsBoundsAnalysis.BR_tWpb[0];
points2print[store_vector_index][53]=HiggsBoundsAnalysis.BR_tHpjb[0];
points2print[store_vector_index][54]=HiggsBoundsAnalysis.BR_Hpjcs[0];
points2print[store_vector_index][55]=HiggsBoundsAnalysis.BR_Hpjcb[0];
points2print[store_vector_index][56]=HiggsBoundsAnalysis.BR_Hptaunu[0];
points2print[store_vector_index][57]=HiggsBoundsAnalysis.MHplusGammaTot[0];
//Hlight
//Tevatron
points2print[store_vector_index][58]=HiggsBoundsAnalysis.CS_tev_hj_ratio[0];
//LHC7
points2print[store_vector_index][59]=HiggsBoundsAnalysis.CS_lhc7_hj_ratio[0];
//LHC8
points2print[store_vector_index][60]=HiggsBoundsAnalysis.CS_lhc8_hj_ratio[0];
//A
//Tevatron
points2print[store_vector_index][61]=HiggsBoundsAnalysis.CS_tev_hj_ratio[1];
//LHC7
points2print[store_vector_index][62]=HiggsBoundsAnalysis.CS_lhc7_hj_ratio[1];
//LHC8
points2print[store_vector_index][63]=HiggsBoundsAnalysis.CS_lhc8_hj_ratio[1];
//Hheavy
//Tevatron
points2print[store_vector_index][64]=HiggsBoundsAnalysis.CS_tev_hj_ratio[2];
//LHC7
points2print[store_vector_index][65]=HiggsBoundsAnalysis.CS_lhc7_hj_ratio[2];
//LHC8
points2print[store_vector_index][66]=HiggsBoundsAnalysis.CS_lhc8_hj_ratio[2];
//HiggsSignals Pvalue
#ifdef HSMODE
points2print[store_vector_index][67]=Pvalue;
#endif
#endif
//-----------------------------------------------------//
int data=67;
//extra_data
#ifdef HSMODE
data=68;
#endif
for(size_t i=0;i!=extra_data.size();++i){
points2print[store_vector_index][data]=extra_data[i];
++data;
}
//---------------------------------------------------/*/
}
void TwoHDM::UserFinalCalcs(void){
////////////////////////////////////////////////////////////////////////
// ENTER HERE CODE FOR YOUR FINAL CALCULATIONS AFTER THE SCAN IS DONE //
////////////////////////////////////////////////////////////////////////
// Function for User defined final calculations to be executed after //
// the scan is finished //
////////////////////////////////////////////////////////////////////////
}
bool TwoHDM::CheckStability(LambdaRef & L)const{
///////////////////////////////////////////////////////////////////
///// 2HDM Potential Stability function /////
///////////////////////////////////////////////////////////////////
///// THIS FUNCTION IS ON A USER.CPP FILE MAINLY TO SHOW THE /////
///// USER WHAT THE SCAN IS CALCULATING. WE WOULD ADVISE YOU /////
///// NOT TO ALTER THE CONTENTS OF THIS FUNCTION. /////
///////////////////////////////////////////////////////////////////
if(L[3]>0 && L[4]>0 && L[5]+sqrt(L[3]*L[4])>0 && L[5]+L[6]-abs(L[2])+sqrt(L[3]*L[4])>0){
return true;
}
else{
#ifdef VERBOSE
std::clog<<endl<<"---> ScannerS: Point rejected by TwoHDM::CheckStability(LambdaRef & L)."<<endl;
#endif
return false;
}
}
bool TwoHDM::CheckGlobal(PhiRef & Phi,LambdaRef & L,Potential & V)const{
///////////////////////////////////////////////////////////////////
///// 2HDM Global Minimum function /////
///////////////////////////////////////////////////////////////////
///// THIS FUNCTION IS ON A USER.CPP FILE MAINLY TO SHOW THE /////
///// USER WHAT THE SCAN IS CALCULATING. WE WOULD ADVISE YOU /////
///// NOT TO ALTER THE CONTENTS OF THIS FUNCTION. /////
///////////////////////////////////////////////////////////////////
double kd = pow((L[3]/L[4]),0.25);
double Disc=L[7]*(L[0]-kd*kd*L[1])*(Phi[6]/Phi[2]-kd);
if(Disc <= 0){
#ifdef VERBOSE
std::clog<<endl<<"---> ScannerS: Point rejected by TwoHDM::CheckGlobal(PhiRef & Phi,LambdaRef & L,Potential & V)."<<endl;
#endif
return false;//If condition not met, reject point
}
else
return true;
}
//////////////////////////////////////////////////////////
//////// 2HDM default parameterizations //////////
//////////////////////////////////////////////////////////
// These functions will be called when you use the 2HDM //
// class in order to define alpha and tanbeta. //
//////////////////////////////////////////////////////////
void TwoHDM::MyPhiParametrization(const PhiParamVec & PhiPar,PhiVec & Phi){
//////////////////////////
// Variables:
// PhiPar[] : Vector of VEV parameters
// Phi[] : Vector of VEVs
///////////////////////////
// Description:
//The purpose of this function is to provide for a parametrization of the VEVs such that they can be generated in more generic ways other than just in a hypercubic box. For example in the 2HDM model we may want to generate h_1=246*cos(beta) and h_2=246*sin(beta), i.e. on a circle, using beta as a parameter.
//Here you must write down such a relation in the form Phi[i]=function(PhiPar[0],PhiPar[1],...). The parameter ranges have already been specified in the mathematica notebook.
//Phi[2]=246*cos(PhiPar[0]);
//Phi[6]=246*sin(PhiPar[0]);
Phi[2]=246/sqrt(1+PhiPar[0]*PhiPar[0]);
Phi[6]=246*PhiPar[0]/sqrt(1+PhiPar[0]*PhiPar[0]);
}
void TwoHDM::MyInternalMixing(const PhiParamVec & PhiPar,const PhiVec & Phi,MixingparamVec & MixPar,vector< vector<double> > & MixInternal){
//////////////////////////
// Variables:
// PhiPar[] Vector of VEV parameters
// Phi[] Vector of VEVs
// MixPar[] Vector of Mixing matrix parameters
// MixInternal[i][j] Used defined mixing matrix for the mixing block. i runs over the mass eigenstates and j runs from 0 to Nmixing-1 (Nmixing is the number or mixing states)
// r random number generator (called as r(mix,max)) in case the user whishes to change here the scan of MixPar[] more flexibly.
///////////////////////////
// Description:
//The purpose of this function is to provide for a parametrization of the mixing matrix when for simpler cases where an analytic parametrization is natural, just as for example in the 2HDM. Since at this point in the scan the VEVs have been generated (and possibly their parameters, if they have been re-parametrized), we assume a general structure where the mixing matrix parametrization can depend on such (constant) quantities (i.e. Phi, the VEVS and PhiPar, the parameters) and a new set of parameters (MixPar) whose ranges have been specified in the Mathematica notebook.
//Contrarily to the VEV re-parametrization, we allow here the use the option to actually change the values of MixPar (i.e. override the values that have been generated in the ranges defined in the mathematica notebook) so that relations between Mixing parameters and VEV parameters are possible (such as beta=alpha-Pi/2 in the decoupling limit of the 2HDM).
//Here you must write down:
// i) (this one is optional) a relation MixPar[i]=function(Phi[0],...,PhiPars[0],...)
//MixPar[0]=PhiPar[0]-acos(-1)/2e0;
// ii) Expressions for the mixing matrix elements in the form MixInternal[i][j]=Function(MixPar[0],...,Phi[0],...,PhiPars[0],...), where the index i runs over the scalar eigenstates and the index j runs over the original mixing scalar states according to the order defined in the mathematica file. IMPORTANT NOTE: You are advised to check (in VERBOSE mode) that the ordering is indeed as you think it is. Note for example that the mass eigenstates are only ordered in mass if you enforce this in the Mathematica file! The work around for example in the 2HDM mode is to split the scan in 2 scans: The first, where the SM Higgs is the Heavy (and this state mass[0] fixed to 126) and the other is the ligth one (mass[1] < 126) and then another scan where the SM Higgs is the light (and this state mass[1] fixed at 126) and the other is the Heavy one (mass[0] > 126)
MixInternal[0][0]=cos(MixPar[0]);
MixInternal[0][1]=sin(MixPar[0]);
MixInternal[1][0]=-sin(MixPar[0]);
MixInternal[1][1]=cos(MixPar[0]);
}
void TwoHDM::MyCoupMassRelations(const PhiParamVec & PhiPar,const PhiVec & Phi,const vector< vector<double> > & Mixing,const MixingparamVec & MixPar,vector<double> & L,vector<double> & Mass, MLparamVec & MLPar, vector<double> & IndependentMs,vector<double> & IndependentLs){
///////////////////////////////////////////////////////////
/////// WARNING: THIS FUNCTION WILL NOT BE CALLED BY /////
/////// DEFAULT SINCE THE TwoHDM.potential FILE /////
/////// DOESN'T HAVE THIS PARAMETRIZATION ON /////
///////////////////////////////////////////////////////////
//////////////////////////
// Variables:
// PhiPar[] Vector of VEV parameters
// Phi[] Vector of VEVs
// MixPar[] Vector of Mixing matrix parameters
// Mixing[i][j] Full mixing matrix (diagonal directions also present in the lines below the mixing directions) i runs over all mass eigenstates and j runs over all initial (gauge) scalar states according to the Mathematica ordering
// r random number generator (called as r(mix,max)) in case the user whishes to change here the scan of MixPar[] more flexibly.
///////////////////////////
// Description:
//The purpose of this function is to provide a way of imposing extra relations in the scan among the last masses and couplings that are left to be generated at the end, through a parametrization, and/or with all the other quantities that have already been generated in in the earlier stages of the scan step.
//IMPORTANT NOTE: it is crucial in this case, that you check in Verbose mode that the relations you impose below are among independent masses and couplings (in VERBOSE mode there are some lines at the end identifying the last independent quantities that are generated). If not, you may have to re-order your fields and couplings in the input mathematica notebook (the rule is that the prigram tries to leave for last, the last masses and couplings in the numbering adopted in the Mathematica notebook). In any case, below you will have to "declare" which couplings and masses you want to constrain so the code will crash and issue a warning if what you are trying to do is not allowed.
// Since at this point in the scan the VEVs, and Mixings have been generated (and possibly their parameters, if they have been re-parametrized), we assume a general structure where the relations you impose may depend on all such (constant) quantities and a new set of parameters (MLPar) whose ranges have been specified in the Mathematica notebook.
//Again, we allow here the use the option to actually change the values of MLPar (i.e. override the values that have been generated in the ranges defined in the mathematica notebook).
//Here you must write down:
// i) (WARNING: this one is NOT optional) specify here which masses and couplings (in the allowed set of independent couplings) you are going to use in your parametrization
//IndependentMs.push_back(1); //Mass[1] here for example will re-parametrised
//IndependentMs.push_back(2); //and Mass[2] as well
//and couplings (in this example we comment the next line and use only relations among masses)
IndependentLs.push_back(7);
// ii) Now write here your relations (in this example, this re-parametrization sets two masses to be equal)
//Mass[1]=MLPar[0];
//Mass[2]=MLPar[0];
L[7]=MLPar[0]*20000*Phi[2]/Phi[6];
}

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 1:43 PM (20 h, 13 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023006
Default Alt Text
ScannerSUserTwoHDM.cpp (18 KB)

Event Timeline