Page MenuHomeHEPForge

No OneTemporary

Index: branches/sHdecay_CxSM_and_RxSM/ScannerSUserRxSM.cpp
===================================================================
--- branches/sHdecay_CxSM_and_RxSM/ScannerSUserRxSM.cpp (revision 109)
+++ branches/sHdecay_CxSM_and_RxSM/ScannerSUserRxSM.cpp (revision 110)
@@ -1,255 +1,304 @@
#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 RxSM ////////////////////
//////////////////////////////////////////////////////////
// These functions will only be called when you use a //
// TwoHDM input file. //
//////////////////////////////////////////////////////////
void RxSM::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
//// Setting up the output print options
//// leave the 1st argument at 1 if you want print points 1 by 1
if(DM_candidate)
setprint(50,22);
else
setprint(50,34);
//// 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);
LoadAll2HDMTables();
#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 RxSM::UserAnalysis(PhiRef & Phi,LambdaRef & L,MassRef & Mass, MmixingRef & Mixing){
//////////////////////////////////////////////////////
// ENTER CODE FOR YOUR TESTS/OUTPUT DURING THE SCAN //
//////////////////////////////////////////////////////
//--------------------------------------------------//
if(not DM_candidate)//Test not needed for DM case
//Test Electroweak precision observables
if(not TestSTU_2sigma_xSM2_mixall(Mass,Mixing))
return false;
//Colider tests
#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
return true;
}
void RxSM::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;
for(size_t i=0;i!=2;++i){
points2print[store_vector_index][data]=Cp.Mass[h[i]];
++data;
}
points2print[store_vector_index][data]=Phi[2];
++data;
points2print[store_vector_index][data]=Phi[4];
++data;
for(size_t i=0;i!=5;++i){
points2print[store_vector_index][data]=Cp.L[i];
++data;
}
points2print[store_vector_index][data]=alpha;
++data;
/*for(int i=0;i!=nHzero;++i){
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;
}
//branching ratios for higgs to 2 higgs decay
if(DM_candidate){
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjinvisible[0];
++data;
}
else{
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjhihi[0][1];
++data;
points2print[store_vector_index][data]=HiggsBoundsAnalysis.BR_hjhihi[1][0];
}*/
}
void RxSM::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 RxSM::CheckStability(LambdaRef & L)const{
///////////////////////////////////////////////////////////////////
///// RxSM 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. /////
///////////////////////////////////////////////////////////////////
- // WARNING!!!! THIS IS ONLY A SUFICIENT CONDITION!!!!!!!!!1
- if(L[1]<0 || L[2]<0 || L[4]<0){
+ if(L[1]<0 || L[4]<0 || L[2]<-sqrt(L[1]*L[4]/6.)){
#ifdef VERBOSE
std::clog<<endl<<"---> ScannerS: Point rejected by RxSM::CheckStability(LambdaRef & L)."<<endl;
#endif
return false;
}
return true;
}
bool RxSM::CheckGlobal(PhiRef & Phi0,LambdaRef & L0,Potential & V)const{
///////////////////////////////////////////////////////////////////
///// RxSM 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;
+ //Case where no EWSB occurs H=0 S=0
+ Phi[2]=0e0;
+ Phi[4]=0e0;
+ Vtest=V();
+ if(Vtest<V0){
+ return false;
+ }
+
+ //clog<<"pass1"<<endl;
+ //Case without Higgs mechanism (H=0)
+ double Rsq=-6*L[3]/L[4];
+ if(Rsq>0){
+ Phi[2]=0e0;
+ Phi[4]=sqrt(Rsq);
+ Vtest=V();
+ if(Vtest<V0)
+ return false;
+ }
+ //Dark matter phase (tested when broken phase)
+ if(not DM_candidate){
+ double Hsq=-L[0]/L[1];
+ Phi[2]=sqrt(Hsq);
+ Phi[4]=0;
+ Vtest=V();
+ if(Vtest<V0)
+ return false;
+ }
+
+ //Broken phase (tested when DM case)
+ if(DM_candidate){
+ Rsq=(L[0]*L[2]/L[1]-L[3])/(L[4]/6.-L[2]*L[2]/L[1]);
+ double Hsq=-L[0]/L[1]-L[2]/L[1]*Rsq;
+ Phi[2]=sqrt(Hsq);
+ Phi[4]=sqrt(Rsq);
+ Vtest=V();
+ if(Vtest<V0)
+ return false;
+ }
+
return true;
}
//////////////////////////////////////////////////////////
//////// RxSM default parameterizations //////////
//////////////////////////////////////////////////////////
// These functions will be called when you use the RxSM //
// class. RxSM.potential does not have parametrizations //
// on by default, so they will not be called. //
//////////////////////////////////////////////////////////
void RxSM::MyPhiParametrization(const PhiParamVec & PhiPar,PhiVec & Phi){
///////////////////////////////////////////////////////////
/////// WARNING: THIS FUNCTION WILL NOT BE CALLED BY /////
/////// DEFAULT SINCE THE RxSM.potential FILE /////
/////// DOESN'T HAVE THIS PARAMETRIZATION ON /////
///////////////////////////////////////////////////////////
}
void RxSM::MyInternalMixing(const PhiParamVec & PhiPar,const PhiVec & Phi,MixingparamVec & MixPar,vector< vector<double> > & MixInternal){
///////////////////////////////////////////////////////////
/////// WARNING: THIS FUNCTION WILL NOT BE CALLED BY /////
/////// DEFAULT SINCE THE RxSM.potential FILE /////
/////// DOESN'T HAVE THIS PARAMETRIZATION ON /////
///////////////////////////////////////////////////////////
}
void RxSM::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 RxSM.potential FILE /////
/////// DOESN'T HAVE THIS PARAMETRIZATION ON /////
///////////////////////////////////////////////////////////
IndependentMs.push_back(0);
IndependentMs.push_back(1);
Mass[0]=MLPar[0];
Mass[1]=MLPar[1];
h[0]=0;
h[1]=1;
if(DM_candidate && Mixing[0][2]==0){
//clog<<"Warning: Flipped case"<<endl;
Mass[0]=MLPar[1];
Mass[1]=MLPar[0];
h[0]=1;
h[1]=0;
}
else if(not DM_candidate &&MLPar[0]>MLPar[1]){
h[1]=0;
h[0]=1;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:38 PM (9 h, 51 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023669
Default Alt Text
(11 KB)

Event Timeline