This section contains model files for MicroMegas that have been created with FeynRules and then exported in CalcHEP format with subsequent slight corrections to make the model work in CalcHEP (sometimes CalcHEP is picky on the order of function definitions etc). The model files listed here can be used in a newProject created project workspace inside MicroMegas 2.2 and higher without modifications. The separate subsections for separate models describe the publication according to which they are implemented as well as variable explanations.
These tips are based on the micrOMEGAs 2.0 manual which can be found here.
Default model parameters are defined in work/model/vars1.mdl file. To change these variables you can use the assignVal(name,val) function. If you want to print out parameters, then you can either print out all parameters with for example printVar(stdout) or find the variable first with findVal(name,&val) and then use the printf facility in C.
The findVal function can also be used to find the value of calculated variables e.g. vacuum expectation value or some new particle mass etc, however to make this work the respective parameter needs to have a * prefixed in front of the variable name in the work/models/func1.mdl or it will just not be found.
Before any other calculations can happen the model parameters have to be set to their needed values and the function sortOddParticles(message) has to be called where message is a string variable containing the error message in case something goes wrong. Remember, that this function has to be called always after you change the model parameters as otherwise you might get inconsistent results. Once the odd particles are sorted you can find out which particle was defined as the dark matter candidate from the message variable and the mass of the particle from the lopmass_() function e.g.
printf("Dark Matter candidate is %s with mass %1.2f GeV\n",mess,lopmass_());
The routine which is why MicroMegas was first conceived is the calculation of the relic density, which you can do with the following function:
darkOmega(&Xf, fast, Beps)
where Xf is a pointer to a double precision variable, fast is either 1 or 0 depending if you want a fast or precise calculation. To the first order fast gives you the same result within 1% so if that's ok, then it's the preferred method. Beps defines a cutoff for which channels will be added to the calculation. If you just want the LOP itself, then set Beps to 1. If you want to consider most channels then a recommended value is around 10^-4 - 10^-6. The function returns the Omega h^2 value, Xf will be M_LOP/Tf, where Tf is the freezout temperature.
If you want to know which channels contributed to the calculation, then the function printChannels can be used, the arguments the function takes are:
printChannels(Xf,cut,Beps,prcnt,file) Xf - as returned by darkOmega cut - lowest value to be printed Beps is the same as in darkOmega prcnt - if non-zero, then the output will be in percentages file - file where to print this information or stdout to print to the output screen
Calculating the spectrums of photons, positrons, anti-protons and neutrinos is easily possible with the calcSpectrum command. However note that the command uses pre-built tables from pythia hadronizations to get the final spectra based on relative contributions from different channels. AFAIK a newer version for spectrum calculation is on the way, but as of the writing of this document was not yet available in MicroMegas.
calcSpectrum(v,outP,tab,&err) v - relative velocity in natural units (e.g. 0.001 in our galaxy ~300km/s) outP - particle for which the spectrum is calculated (0,1,2,3,4,5 e.g. gamma,e+,p-,nu_e,nu_mu,nu_tau) tab - array of doubles which contains 250 elements, essentially defined as double tab[250]; err - pointer to an int which will be 0 if everything was ok The function itself returns the relative cross section <sigma v>
The spectrum you can then use for what ever reasons you need. The most likely use is to propagate it with GalProp and for that you first need to save the spectrum to a file, which can be done using the following function:
spectrTable(tab,fname,mess,Xmin,N) tab - the array with the spectrum from calcSpectrum function fname - filename, where to write mess - string that describes the spectrum Xmin - minimal energy to consider in the units of E/M_LOP N - number of points (<300)
If you want to use the spectrum directly during code execution, then you can use zInterp(x,tab) function where x=log(E/M) that you want and tab is the spectrum. The result is dN/dx where N is the number of particles.
Although MicroMegas supports also cross section calculations they aren't trivial. More comprehensive information is availabel in MicroMegas documentation as referenced above. The one situation that I was interested in is calculation of the relative cross sections of the final annihilation channels based on the current model parameters inside the program. One can do that based on an example given in the sample MSSM code. I'll give an altered example with descriptions below:
numout* cc; // variable pointing to the compiled library
double Pcm; // Center of mass energy
cc=newProcess("","omg__t_si__t_si"); // Start using the pre-compiled library, you have to look
// up the name in work/so_generated/ folder e.g. work/so_generated/omg__t_si__t_si.so
if(cc) // If the library is available
{ int ntot,l;
char * name[4];
double mass[4];
procInfo1(cc,&ntot,NULL,NULL); // ntot is filled with number of sub-processes available
for(l=1;l<=ntot; l++) // iterate over all subprocesses
{ int err;
double cs;
procInfo2(cc,l,name,mass); // Fills an array of particle names and masses for that particular sub-process
if(l==1) { Pcm=mass[0]*v/2; printf("(Pcm=%.2E)\n",Pcm);}
// Just some printout for the first sub-process. The particle in slot 0 is always the LOP
printf("%3s,%3s -> %3s %3s ",name[0],name[1],name[2],name[3]); // Print the actual process
cs= cs22(cc,l,Pcm,-1.,1.,&err);
// Calculate the cross section using the library in cc, subprocess l, c.m.s energy Pcm.
// -1 .. 1 in cos(theta) phase space e.g. all
if(err) printf("Error\n"); // If there is an error just say so
else if(cs==0.) printf("Zero\n"); // If the cross setion comes out zero print Zero
else printf("%.2E [pb] ( sigma*v=%.2E [cm^3/sec] ) \n",cs, cs*v*2.9979E-26);
// Print out the cross section in natural units as well as units of cm^3/s
}
}
an example output of that code snippet:
(Pcm=1.81E-01) ~Si,~Si -> H H 2.85E+02 [pb] ( sigma*v=1.71E-26 [cm^3/sec] ) ~Si,~Si -> w W 2.57E+01 [pb] ( sigma*v=1.54E-27 [cm^3/sec] ) ~Si,~Si -> Z Z 2.11E+01 [pb] ( sigma*v=1.26E-27 [cm^3/sec] ) ~Si,~Si -> b B 7.62E-01 [pb] ( sigma*v=4.57E-29 [cm^3/sec] ) ~Si,~Si -> t T 1.87E+01 [pb] ( sigma*v=1.12E-27 [cm^3/sec] ) ~Si,~Si -> c C 6.96E-02 [pb] ( sigma*v=4.17E-30 [cm^3/sec] ) ~Si,~Si -> ta TA 3.63E-02 [pb] ( sigma*v=2.18E-30 [cm^3/sec] )