Page MenuHomeHEPForge
Files F19250920

No OneTemporary
Actions

Size
15 KB
Referenced Files
None
Subscribers
None

View Options

Index: tags/ScannerSv1_1_1/ScannerSUserxSM2.cpp
===================================================================
--- tags/ScannerSv1_1_1/ScannerSUserxSM2.cpp (revision 89)
+++ tags/ScannerSv1_1_1/ScannerSUserxSM2.cpp (revision 90)
@@ -1,432 +1,433 @@
#include "ScannerScore/ScannerSUser.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 xSM ////////////////////
//////////////////////////////////////////////////////////
// These functions will only be called when you use a //
// TwoHDM input file. //
//////////////////////////////////////////////////////////
void xSM2::UserInitCalcs(void){
//////////////////////////////////////////////////////////////////////////
// ENTER HERE CODE FOR YOUR INITIAL CALCULATIONS BEFORE THE SCAN STARTS //
//////////////////////////////////////////////////////////////////////////
//// Setting up the output print options
//// leave the 1st argument at 1 if you want print points 1 by 1
setprint(1,15);
//// Setting up the output print options
//// leave the 1st argument at 1 if you want print points 1 by 1
//setprint(1,49);
//// 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(0);
//MIX_ALL example
//setprint(50,49);
//DM example
//setprint(50,39);
#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 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 xSM2::UserAnalysis(PhiRef & Phi,LambdaRef & L,MassRef & Mass, MmixingRef & Mixing){
//////////////////////////////////////////////////////
// ENTER CODE FOR YOUR TESTS/OUTPUT DURING THE SCAN //
//////////////////////////////////////////////////////
//Test Electroweak precision observables
if(not TestSTU_2sigma_xSM2_mixall(Mass,Mixing))
return false;
//computes scalar widths for this point
compute_scalar_widths();
#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
if(not PrepareHiggsBoundsInput())
return false;
int HBresult,chan,ncombined;
double obsratio;
HiggsBoundsAnalysis.run(HBresult,chan,obsratio,ncombined);
#ifdef HSMODE
//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(HBresult==0 || Pvalue<0.0027)
return false; //Exclude point if excluded from non-observation in the HB search channels or if probability of model (calculated with higgs signals) given the observed signal rates is too low (3sigma). Note that 2sigma is more standard, however, we can recover the 1 and 2 sigma regions by post-processing using the provided pvalue on output.
#endif
#endif
#endif
#ifdef MICROMEGASMODE
//if micromegas is on AND we have DM candidate
if(DM_candidate){
assignMOvalues();
// Compute dark matter contributions from dark matter particle 1 (lightest)
double Xf,cut=0.01,Beps=1e-5;
// First the relic density contribution
Omega=darkOmega(&Xf,1,Beps);
if(Omega>0.13 || Omega<0)
return false;
// spin_amplitudes: 0: pA0; 1: pAspin; 2: nA0; 3:nAspin
double spin_amplitudes[4][2];
if(not test_nucleon_amplitudes(Omega,sigmaSI3,Cp.Mass[2],spin_amplitudes,sigmaup)){
return false;
}
}
#endif
}
void xSM2::what2print(){
//////////////////////////////////////////////////////////////////////////
// CHOSE WHAT YOU WANT TO BE PRINTED TO THE OUTPUT FILE HERE //
//////////////////////////////////////////////////////////////////////////
// YOU NEED TO USE setprint() IN UserInitCalcs() TO DEFINE HOW MUCH //
// DATA YOU WANT TO PRINT PER POINT (AND HOW MANY POINTS TO STORE //
// BEFORE PRINTING). //
//////////////////////////////////////////////////////////////////////////
size_t data=0;
//Masses
for(size_t i=0;i!=3;++i){
points2print[store_vector_index][data]=Cp.Mass[i];
++data;
}
//VEVs
for(size_t i=4;i!=6;++i){
points2print[store_vector_index][data]=Phi[i];
++data;
}
//couplings
for(size_t i=0;i!=7;++i){
points2print[store_vector_index][data]=Cp.L[i];
++data;
}
//Mixing angles
for(size_t i=0;i!=3;++i){
points2print[store_vector_index][data]=Mmixing[i][2];
++data;
}
/*//////////////////////// MIX_ALL example /////////////////////
size_t data=0;
for(size_t i=0;i!=3;++i){
points2print[store_vector_index][data]=Cp.Mass[i];
++data;
}
for(size_t i=4;i!=6;++i){
points2print[store_vector_index][data]=Phi[i];
++data;
}
for(size_t i=0;i!=7;++i){
points2print[store_vector_index][data]=Cp.L[i];
++data;
}
#ifdef HBMODE
//Branching Ratios & decay width
for(int i=0;i!=3;++i){
points2print[store_vector_index][data]=Mmixing[i][2];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjss[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjcc[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjbb[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjmumu[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjtautau[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjWW[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjZZ[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjZga[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjgaga[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjgg[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.MhGammaTot[i];
++data;
}
#ifdef HSMODE
points2print[store_vector_index][data]=Pvalue;
++data;
#endif
#endif
////////////////////////---------------------////////////////////*/
/*////////////////////////// DM example /////////////////////////
size_t data=0;
for(size_t i=0;i!=3;++i){
points2print[store_vector_index][data]=Cp.Mass[i];
++data;
}
for(size_t i=4;i!=5;++i){
points2print[store_vector_index][data]=Phi[i];
++data;
}
for(size_t i=0;i!=7;++i){
points2print[store_vector_index][data]=Cp.L[i];
++data;
}
#ifdef HBMODE
//Branching Ratios & decay width
for(int i=0;i!=2;++i){
points2print[store_vector_index][data]=Mmixing[i][2];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjss[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjcc[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjbb[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjmumu[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjtautau[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjWW[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjZZ[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjZga[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjgaga[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjgg[i];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.MhGammaTot[i];
++data;
}
points2print[store_vector_index][data]=Mmixing[2][2];
++data;
#ifdef HSMODE
points2print[store_vector_index][data]=Pvalue;
++data;
#endif
#endif
#ifdef MICROMEGASMODE
points2print[store_vector_index][data]=Omega;
++data;
#endif
points2print[store_vector_index][data]=sigmaSI3;
++data;
points2print[store_vector_index][data]=sigmaup;
++data;
////////////////////////---------------------////////////////////*/
}
void xSM2::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 xSM2::CheckStability(LambdaRef & L)const{
///////////////////////////////////////////////////////////////////
///// xSM2 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[1]<0 || L[4]<0 || (L[2]<0 && pow(L[2],2)>L[1]*L[4])){
#ifdef VERBOSE
std::clog<<endl<<"---> ScannerS: Point rejected by xSM2::CheckStability(LambdaRef & L)."<<endl;
#endif
return false;
}
return true;
}
bool xSM2::CheckGlobal(PhiRef & Phi0,LambdaRef & L0,Potential & V)const{
///////////////////////////////////////////////////////////////////
///// xSM2 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. /////
///////////////////////////////////////////////////////////////////
//Bind references to auxiliary potential Vaux
PhiRef Phi(V);
LambdaRef L(V);
//Copy values from current parameter space point
Phi.copy(Phi0);
L.copy(L0);
//// Compute potential at the current minimum
double V0=V();
double Vtest;
double * atemp;
atemp = new double [3+1];
double * ztemp;
ztemp = new double [2*3];
//Possible minimum with H=0 and A=0
atemp[0]=L[6]*sqrt(2);
atemp[1]=L[3]/2+L[5]/2;
atemp[2]=0;
atemp[3]=L[4]/4;
polyroots(atemp,ztemp,3);
for(int i=0;i<3;++i){
if(abs(ztemp[2*i+1])<1e-8*abs(ztemp[2*i])){
Phi[2]=0e0;
Phi[4]=ztemp[2*i];
+ Phi[5]=0e0;
Vtest=V();
if(Vtest<V0)
return false;
}
Phi.copy(Phi0);
}
//Possible minimum with H=0 and A and S non-zero
double vx = -L[6]*sqrt(2)/L[5];
double squarevy= 2/L[4]*( L[5] - L[3] - L[4]*L[6]*L[6]/L[5]/L[5] );
if(squarevy>0){
Phi[2]=0;
Phi[4]=vx;
Phi[5]=sqrt(squarevy);
Vtest=V();
if(Vtest<V0)
return false;
}
Phi.copy(Phi0);
//Possible minimum with H ans S nonzero and A=0
atemp[0]=2*L[6]*sqrt(2);
atemp[1]=L[3]+L[5]-L[2]*L[0]/L[1];
atemp[2]=0;
atemp[3]=L[4]/2e0-L[2]*L[2]/2e0/L[1];
polyroots(atemp,ztemp,3);
for(int i=0;i<3;i++){
if(abs(ztemp[2*i+1])<1e-8*abs(ztemp[2*i])){
Phi[4]=ztemp[2*i];
double Xtemp=-(2*L[0]+L[2]*Phi[4]*Phi[4])/L[1];
if(Xtemp>0){
Phi[2]=sqrt(Xtemp);
Phi[5]=0e0;
Vtest=V();
if(Vtest<V0)
return false;
}
}
Phi.copy(Phi0);
}
//Possible minimum with H nonzero and A nonzero
vx=-L[6]*sqrt(2)/L[5];
squarevy= 2*(L[5]*L[5]*(L[1]*(L[5]-L[3])+L[0]*L[2])-L[4]*L[6]*L[6]*L[1]+L[2]*L[2]*L[6]*L[6])/(L[4]*L[5]*L[5]*L[1]-L[2]*L[2]*L[5]*L[5]);
//2*( L[1]*L[5]- L[1]*L[3]+ ( -L[4]*L[1]+L[2]*L[2] )*L[6]*L[6]/L[5/L[5] +L[0]*L[2])/(L[1]*L[4]-L[2]*L[2]);
double vsquare=-(2*L[0] + L[2] *vx*vx + L[2]*squarevy )/L[1];
if(squarevy>0){
Phi[2]=sqrt(vsquare);
Phi[4]=vx;
Phi[5]=sqrt(squarevy);
Vtest=V();
if(Vtest<V0){
#ifdef VERBOSE
std::clog<<endl<<"---> ScannerS: Point rejected by xSM2::CheckGlobal(PhiRef & Phi0,LambdaRef & L0,Potential & V)."<<endl;
#endif
return false;
}
}
delete [] atemp;
delete [] ztemp;
return true;
}
bool RGE_Stop_criteria(const std::vector<double> & y,void *pvoid){
///////////////////////////////////////////////////////////
/////// This function is can be used with the /////
/////// RGETools to evolve the potential couplings /////
///////////////////////////////////////////////////////////
return true;
}
//////////////////////////////////////////////////////////
//////// xSM2 default parameterizations //////////
//////////////////////////////////////////////////////////
// These functions will be called when you use the xSM2 //
// class. xSM2.potential does not have parametrizations //
// on by default, so they will not be called. //
//////////////////////////////////////////////////////////
void xSM2::MyPhiParametrization(const PhiParamVec & PhiPar,PhiVec & Phi){
///////////////////////////////////////////////////////////
/////// WARNING: THIS FUNCTION WILL NOT BE CALLED BY /////
/////// DEFAULT SINCE THE xSM2.potential FILE /////
/////// DOESN'T HAVE THIS PARAMETRIZATION ON /////
///////////////////////////////////////////////////////////
}
void xSM2::MyInternalMixing(const PhiParamVec & PhiPar,const PhiVec & Phi,MixingparamVec & MixPar,vector< vector<double> > & MixInternal){
///////////////////////////////////////////////////////////
/////// WARNING: THIS FUNCTION WILL NOT BE CALLED BY /////
/////// DEFAULT SINCE THE xSM2.potential FILE /////
/////// DOESN'T HAVE THIS PARAMETRIZATION ON /////
///////////////////////////////////////////////////////////
}
void xSM2::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 xSM2.potential FILE /////
/////// DOESN'T HAVE THIS PARAMETRIZATION ON /////
///////////////////////////////////////////////////////////
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Sep 30, 5:45 AM (11 h, 46 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6566343
Default Alt Text
(15 KB)

Event Timeline

HEPForge.org · .