Page MenuHomeHEPForge

ScannerSUser.cpp
No OneTemporary

ScannerSUser.cpp

#include "ScannerScore/ScannerSUser.h"
#include "ScannerSUserAux.h"
using namespace std;
//////////////////////////////////////////////////////////
////// Function contents to be defined by the user ///////
//////////////////////////////////////////////////////////
/// WARNING: Do not change the prototypes of these ///
/// template functions since they are used by the core ///
/// routines. Only change their contents! ///
//////////////////////////////////////////////////////////
void 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*abs(sin(PhiPar[0]));
Phi[6]=246*abs(cos(PhiPar[0]));
}
void MyInternalMixing(const PhiParamVec & PhiPar,const PhiVec & Phi,MixingparamVec & MixPars,vector< vector<double> > & MixInternal,RandGen & r){
//////////////////////////
// Variables:
// PhiPar[] Vector of VEV parameters
// Phi[] Vector of VEVs
// MixPars[] 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 MixPars[] 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 (MixPars) 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 MixPars (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 MixPars[i]=function(Phi[0],...,PhiPars[0],...)
MixPars[0]=PhiPar[0]-acos(-1)/2e0;
// ii) Expressions for the mixing matrix elements in the form MixInternal[i][j]=Function(MixPars[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(MixPars[0]);
MixInternal[0][1]=-sin(MixPars[0]);
MixInternal[1][0]=sin(MixPars[0]);
MixInternal[1][1]=cos(MixPars[0]);
}
void MyCoupMassRelations(const PhiParamVec & PhiPar,const PhiVec & Phi,const vector< vector<double> > & Mixing,const MixingparamVec & MixPars,vector<double> & L,vector<double> & Mass, MLparamVec & MLparams, vector<double> & IndependentMs,vector<double> & IndependentLs,RandGen & r){
//////////////////////////
// Variables:
// PhiPar[] Vector of VEV parameters
// Phi[] Vector of VEVs
// MixPars[] 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 MixPars[] 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 (MLparams) whose ranges have been specified in the Mathematica notebook.
//Again, we allow here the use the option to actually change the values of MLparams (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]=MLparams[0];
Mass[2]=MLparams[0];
}
bool CheckStability(LambdaRef & L){
if(L[3]>0 && L[4]>0 && L[5]+ sqrt(L[3]*L[4])>0 && L[5]+L[6]-abs(L[7]) + sqrt(L[3]*L[4]) >0 ){
return false;
}
return true;
}
bool CheckGlobal(PhiRef & Phi0,LambdaRef & L0,Potential & V){
//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);
///////////////////////////////////////////
///// Your calculations here! /////////////
///////////////////////////////////////////
//// Compute discriminant D
double kd = pow((L[3]/L[4]),0.25);
double Disc=(L[0]-kd*kd*L[1])/(Phi[6]/Phi[2]-kd);
if(Disc < 0 || Disc == 0)
return false;
Phi.copy(Phi0);
return true;
}
bool UserAnalysis(PhiRef & Phi,LambdaRef & L,MassRef & Mass, MmixingRef & Mixing){
/// decide first which one is the heavy mixing Higgs
double mHlight=Mass[0],mHheavy=Mass[1],alpha,mA=Mass[2],mHcharged=Mass[3];
if(mHlight>mHheavy){
mHlight=Mass[1];
mHheavy=Mass[0];
alpha=acos(Mixing[0][2]);
}
else
alpha=acos(Mixing[1][2]);
///////////////////////////////////
//////Example of Superiso call ////
///////////////////////////////////
//Create Input File for superiso
int Type2HDM=2;
// 2HDM sub-type
//Type I = 1
//Type II = 2
//Type III () = 3
//Type IV () = 4
double tanbeta=Phi[6]/Phi[2];
CreateInputFileSuperiso2HDM(mHlight,mHheavy,mA,mHcharged,alpha,tanbeta,Type2HDM);
//Example of na observable calculated with superiso
cout<< "bsgamma="<<bsgamma_calculator(superisofile)<<endl;
///////////////////////////////////
///// Example of SuShi call ///////
//////////////////////////////////
int Hparticle=0/*m,A,H*/,pp_ppbar=0,order=1;
double CM_energy=8000;
CreateInputFileSuShi2HDM(Hparticle,pp_ppbar,order,CM_energy,mHlight, mHheavy,mA,mHcharged,alpha,Type2HDM);
double xsecggh_out,errxsecggh_out,xsecbbh_out;
sushixsection_(order,xsecggh_out,errxsecggh_out,xsecbbh_out);
cout <<"----> " <<xsecggh_out <<"\t" << errxsecggh_out << "\t" << xsecbbh_out<< endl;
////////////////////////////////////////
///// Example of Micromegas call ///////
///////////////////////////////////////
////// Dark matter
/// Relic density
//ForceUG=1; /* to Force Unitary Gauge assign 1 */
/*assignVal("Mh",mHlight);
assignVal("MHH",mHheavy);
assignVal("MAs",mA);
assignVal("MHc",mHcharged);
assignVal("m12sqr",L[0]);
assignVal("sb",tanbeta/sqrt(1+tanbeta*tanbeta));
assignVal("sa",sin(alpha));
char cdmName[10];
int err=sortOddParticles(cdmName);
if(err) { cerr<<"Can't calculate "<<cdmName<< endl; exit(-1);}
///// Compute dark matter contributions from dark matter particle 1 (lightest)
double Omega,Xf,cut=0.01,Beps=1e-5;
/// First the relic density contribution
Omega=darkOmega(&Xf,1,1e-5);
/////////Messages for verbose mode: Comment!
cout << endl<<"------------------------------------------------------------"<< endl;
cout << "**** Omega1 = " << Omega<<"\t Mcdm = "<< Mcdm<< endl;
printChannels(Xf,0.01,Beps,1,stdout);
if(Omega>0.112 || Omega<0)
return false;
double pA0[2],pA5[2],nA0[2],nA5[2];
double Nmass=0.939; //nucleon mass
double SCcoeff;
double err1=nucleonAmplitudes(NULL, pA0,pA5,nA0,nA5);
SCcoeff=4/M_PI*3.8937966E8*pow(Nmass*Mcdm/(Nmass+ Mcdm),2.); // Value given by micromegas
double sigmaSI3=SCcoeff*pA0[0]*pA0[0];
cout <<"SigmaSI = "<<sigmaSI3<< endl;
if(mA>=exp(1.9)&& mA<=exp(6.9)){
double sigmaup=Interp(log(mA),Xenon100);
//clog <<"sigmaup = "<< sigmaup<< endl;
if(sigmaSI3*Omega/0.112>sigmaup)
return false;
}*/
////////////////////////////////////////
////////// Final print out//////////////
////////////////////////////////////////
cout<< "--- Masses ---"<< endl;
for(size_t i=0; i!=Phi.size();++i){
cout <<Mass[i]<< endl;
}
cout<< "--- Couplings ---"<< endl;
for(size_t i=0; i!=L.size();++i){
cout << L[i]<< endl;
}
cout<< "--- VEVs ---"<< endl;
for(size_t i=0; i!=Phi.size();++i){
cout << Phi[i]<< endl;
}
cout<<endl;
return true;
}

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 3:47 PM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4012450
Default Alt Text
ScannerSUser.cpp (11 KB)

Event Timeline