Page MenuHomeHEPForge
Files F19250771

No OneTemporary
Actions

Size
22 KB
Referenced Files
None
Subscribers
None

View Options

Index: branches/DMNLO/micrOMEGAs-2.4.1/MSSM/PatrickcMSSM.cpp
===================================================================
--- branches/DMNLO/micrOMEGAs-2.4.1/MSSM/PatrickcMSSM.cpp (revision 1031)
+++ branches/DMNLO/micrOMEGAs-2.4.1/MSSM/PatrickcMSSM.cpp (revision 1032)
@@ -1,702 +1,668 @@
+
//============================= My test program for DM@NLO ========================================
//----------------------------------------------- Header
// Includes etc.
#define MASSES_INFO
#define CONSTRAINTS
#define OMEGA
#define INDIRECT_DETECTION
#define CDM_NUCLEON
#define RGE spheno
#include"../sources/micromegas.h"
#include"../sources/micromegas_aux.h"
#include"lib/pmodel.h"
#include <string>
using namespace std;
#define SUGRAMODEL_(A) A ## SUGRA
#define SUGRAMODEL(A) SUGRAMODEL_(A)
#define PRINTRGE_(A) printf(" Spectrum calculator is %s\n", #A)
#define PRINTRGE(A) PRINTRGE_(A)
// Global Variables
// This defines the RD calculation option (0=MO result (CalcHep), 1=tree (DM@NLO), 2=full one-loop (DM@NLO))
static int RDoption, DMNLOAufrufe, CalcHepAufrufe, IntErrors;
static bool DMNLOReady;
static double ChannelFilter;
// This array contains the 30 most relevant processes for a given relic calculation
typedef struct { double weight; string prtcl1; string prtcl2; string prtcl3; string prtcl4;} bChannel;
static bChannel PatChannel[30];
// External functions are included
extern "C"
{
void chichi2qq_(double &Pcm, int *mflags,double &cutoff, double *myresult);
void ntnt2bb_(double &Pcm,double &mbinp,double &mtinp,int *iflag,double *myresult);
- void directdetection_(double &E_R, double &pcm, double &costheta, int &Supp_Mode, int &Loop_Mode, int &Denom_Mode, double &Qout, double *DDresults);
- void looptest_(double &Danger, double *myresult);
- void looptest2_(double &Danger, double *myresult);
+ void directdetection_(double &pcm, double &costheta, int &Supp_Mode, int &Loop_Mode, int &Denom_Mode, int &Box_Mode, double &Qout, double *DDresults);
}
//----------------------------------------------- ImproveCrossSection
// Here it is checked if DM@NLO is able to improve the cross section from Calchep
// Interface between DM@NLO and micrOMEGAS
// Additional functions for improveCrossSection etc.
bool CanImprove(long n1, long n2, long n3, long n4);
bool ShallImprove(long n1, long n2, long n3, long n4);
void IntError(double Tree, double RealCorrections);
void SetFlags(int *mflags, long n1, long n2, long n3, long n4);
double DMNLOContribution();
double EstimateChannel(string s1, string s2, string s3, string s4);
long Name2PDG(string s);
int DMNLOChannels();
-
+
void improveCrossSection(long n1,long n2,long n3,long n4,double Pcm, double *addr)
{
double test[9], cutoff = 0.0;
int mflags[11];
CalcHepAufrufe++;
if(RDoption>0 && DMNLOReady && CanImprove(n1, n2, n3, n4)){
DMNLOAufrufe++;
SetFlags(mflags, n1, n2, n3, n4);
chichi2qq_(Pcm,mflags,cutoff,test);
*addr=test[0];
if(RDoption==2 && ShallImprove(n1, n2, n3, n4)){
IntError(test[0],test[6]);
*addr = test[8];
}
}
}
bool CanImprove(long n1, long n2, long n3, long n4)
{
//Neu_i + Neu_j -> q qbar
if((n1 == 1000022 || n1 == 1000023 || n1 == 1000025 || n1 == 1000035) && (n2 == 1000022 || n2 == 1000023 || n2 == 1000025 || n2 == 1000035)){
if ((n3==1 && n4==-1) || (n3==2 && n4==-2) || (n3==3 && n4==-3) || (n3==4 && n4==-4) || (n3==5 && n4==-5) || (n3==6 && n4==-6)
|| (n4==1 && n3==-1) || (n4==2 && n3==-2) || (n4==3 && n3==-3) || (n4==4 && n3==-4) || (n4==5 && n3==-5) || (n4==6 && n3==-6)) {
// if ((n3==1 && n4==-1) || (n3==2 && n4==-2) || (n3==3 && n4==-3) || (n3==4 && n4==-4)
// || (n4==1 && n3==-1) || (n4==2 && n3==-2) || (n4==3 && n3==-3) || (n4==4 && n3==-4)) {
// printf("CanImprove = True \n");
return true;
}
}
//Neu_i + Cha_j -> q q'bar
if(((n1 == 1000024 || n1 == 1000037 || n1 == -1000024 || n1 == -1000037) && (n2 == 1000022 || n2 == 1000023 || n2 == 1000025 || n2 == 1000035))
|| ((n2 == 1000024 || n2 == 1000037 || n2 == -1000024 || n2 == -1000037) && (n1 == 1000022 || n1 == 1000023 || n1 == 1000025 || n1 == 1000035))){
if ((n3==2 && n4==-1) || (n3==4 && n4==-3) || (n3==6 && n4==-5) || (n4==2 && n3==-1) || (n4==4 && n3==-3) || (n4==6 && n3==-5)) {
// if ((n3==2 && n4==-1) || (n3==4 && n4==-3) || (n4==2 && n3==-1) || (n4==4 && n3==-3)) {
// printf("CanImprove = True \n");
return true;
}
}
//Cha_i + Cha_j -> q qbar
if(((n1 == 1000024 || n1 == 1000037) && (n2 == -1000024 || n2 == -1000037)) || ((n1 == -1000024 || n1 == -1000037) && (n2 == 1000024 || n2 == 1000037))){
if ((n3==1 && n4==-1) || (n3==2 && n4==-2) || (n3==3 && n4==-3) || (n3==4 && n4==-4) || (n3==5 && n4==-5) || (n3==6 && n4==-6)
|| (n4==1 && n3==-1) || (n4==2 && n3==-2) || (n4==3 && n3==-3) || (n4==4 && n3==-4) || (n4==5 && n3==-5) || (n4==6 && n3==-6)) {
// if ((n3==1 && n4==-1) || (n3==2 && n4==-2) || (n3==3 && n4==-3) || (n3==4 && n4==-4)
// || (n4==1 && n3==-1) || (n4==2 && n3==-2) || (n4==3 && n3==-3) || (n4==4 && n3==-4)) {
// printf("CanImprove = True \n");
return true;
}
}
return false;
}
bool ShallImprove(long n1, long n2, long n3, long n4)
{
int i;
for(i=0; PatChannel[i].weight>ChannelFilter ;i++){
if (strcmp(PatChannel[i].prtcl1.c_str(), pdg2name(n1))==0 && strcmp(PatChannel[i].prtcl2.c_str(), pdg2name(n2))==0
&& strcmp(PatChannel[i].prtcl3.c_str(), pdg2name(n3))==0 && strcmp(PatChannel[i].prtcl4.c_str(), pdg2name(n4))==0) {
return true;
}
}
return false;
}
void IntError(double Tree, double RealCorrections)
{
if(abs(RealCorrections)>0.5*abs(Tree)){
printf("\n\n NUMERICAL INTEGRATION OF 2->3 PROCESSES FAILED!!! \n\n");
IntErrors++;
}
}
void SetFlags(int *mflags, long n1, long n2, long n3, long n4)
{
//Neu_i + Neu_j -> q qbar
if((n1 == 1000022 || n1 == 1000023 || n1 == 1000025 || n1 == 1000035) && (n2 == 1000022 || n2 == 1000023 || n2 == 1000025 || n2 == 1000035)){
mflags[0] = 0;
if(n1 == 1000022) mflags[1] = 1;
if(n1 == 1000023) mflags[1] = 2;
if(n1 == 1000025) mflags[1] = 3;
if(n1 == 1000035) mflags[1] = 4;
if(n2 == 1000022) mflags[2] = 1;
if(n2 == 1000023) mflags[2] = 2;
if(n2 == 1000025) mflags[2] = 3;
if(n2 == 1000035) mflags[2] = 4;
if(n3==1 || n4==1){
mflags[4] = 1;
mflags[3] = 4;
}
if(n3==2 || n4==2){
mflags[4] = 1;
mflags[3] = 3;
}
if(n3==3 || n4==3){
mflags[4] = 2;
mflags[3] = 4;
}
if(n3==4 || n4==4){
mflags[4] = 2;
mflags[3] = 3;
}
if(n3==5 || n4==5){
mflags[4] = 3;
mflags[3] = 4;
}
if(n3==6 || n4==6){
mflags[4] = 3;
mflags[3] = 3;
}
}
//Neu_i + Cha_j -> q q'bar
if(((n1 == 1000024 || n1 == 1000037 || n1 == -1000024 || n1 == -1000037) && (n2 == 1000022 || n2 == 1000023 || n2 == 1000025 || n2 == 1000035))
|| ((n2 == 1000024 || n2 == 1000037 || n2 == -1000024 || n2 == -1000037) && (n1 == 1000022 || n1 == 1000023 || n1 == 1000025 || n1 == 1000035))){
mflags[0] = 1;
mflags[3] = 3;
if(n1 == 1000024 || n1 == -1000024) mflags[1] = 1;
if(n1 == 1000037 || n1 == -1000037) mflags[1] = 2;
if(n2 == 1000022) mflags[2] = 1;
if(n2 == 1000023) mflags[2] = 2;
if(n2 == 1000025) mflags[2] = 3;
if(n2 == 1000035) mflags[2] = 4;
if(n2 == 1000024 || n2 == -1000024) mflags[1] = 1;
if(n2 == 1000037 || n2 == -1000037) mflags[1] = 2;
if(n1 == 1000022) mflags[2] = 1;
if(n1 == 1000023) mflags[2] = 2;
if(n1 == 1000025) mflags[2] = 3;
if(n1 == 1000035) mflags[2] = 4;
if(n3==2 || n4 ==2){
mflags[4] = 1;
}
if(n3==4 || n4 ==4){
mflags[4] = 2;
}
if(n3==6 || n4 ==6){
mflags[4] = 3;
}
}
//Cha_i + Cha_j -> q qbar
if(((n1 == 1000024 || n1 == 1000037) && (n2 == -1000024 || n2 == -1000037)) || ((n1 == -1000024 || n1 == -1000037) && (n2 == 1000024 || n2 == 1000037))){
mflags[0] = 2;
if(n1 == 1000024) mflags[1] = 1;
if(n1 == 1000037) mflags[1] = 2;
if(n2 == -1000024) mflags[2] = 1;
if(n2 == -1000037) mflags[2] = 2;
if(n2 == 1000024) mflags[1] = 1;
if(n2 == 1000037) mflags[1] = 2;
if(n1 == -1000024) mflags[2] = 1;
if(n1 == -1000037) mflags[2] = 2;
if(n3==1 || n4==1){
mflags[4] = 1;
mflags[3] = 4;
}
if(n3==2 || n4==2){
mflags[4] = 1;
mflags[3] = 3;
}
if(n3==3 || n4==3){
mflags[4] = 2;
mflags[3] = 4;
}
if(n3==4 || n4==4){
mflags[4] = 2;
mflags[3] = 3;
}
if(n3==5 || n4==5){
mflags[4] = 3;
mflags[3] = 4;
}
if(n3==6 || n4==6){
mflags[4] = 3;
mflags[3] = 3;
}
}
//REMAINING CHICHI2QQ DEFAULT SETTINGS
// diagonalize (0 = Masses from MO, 1 = Diagonalize)
mflags[5] = 0;
// iflux (1 = sigma, 0 = v.sigma)
mflags[6] = 1;
// only tree (1 = only tree, 0 = Loops)
mflags[7] = 1;
if(RDoption==2 && ShallImprove(n1, n2, n3, n4)) mflags[7] = 0;
// Kinematics in MO Cross Section (0 = Ours, 1 = MO)
mflags[8] = 0;
// Resummation Style (0 Old DMNLO Style, 1 Spira/Patrick Style)
mflags[9] = 1;
// NNLO Resummation (1 = on, 0 = off)
mflags[10] = 1;
// Note: The use NNLO Res. = 1 sets Res. Style to 1 automatically.
//printf("Flags 0 1 2 3 4 = %i %i %i %i %i \n", mflags[0], mflags[1], mflags[2], mflags[3], mflags[4]);
}
double DMNLOContribution()
{
double result = 0.;
long n1, n2, n3, n4;
for(int i=0; PatChannel[i].weight>ChannelFilter ;i++){
n1 = Name2PDG(PatChannel[i].prtcl1);
n2 = Name2PDG(PatChannel[i].prtcl2);
n3 = Name2PDG(PatChannel[i].prtcl3);
n4 = Name2PDG(PatChannel[i].prtcl4);
if (CanImprove(n1, n2, n3, n4)) {
result += PatChannel[i].weight;
}
}
return result;
}
int DMNLOChannels()
{
int result = 0;
long n1, n2, n3, n4;
for(int i=0; PatChannel[i].weight>ChannelFilter ;i++){
n1 = Name2PDG(PatChannel[i].prtcl1);
n2 = Name2PDG(PatChannel[i].prtcl2);
n3 = Name2PDG(PatChannel[i].prtcl3);
n4 = Name2PDG(PatChannel[i].prtcl4);
if (CanImprove(n1, n2, n3, n4)) {
result++;
}
}
return result;
}
long Name2PDG(string s)
{
long result = 0;
if (s=="d") return 1;
if (s=="u") return 2;
if (s=="s") return 3;
if (s=="c") return 4;
if (s=="b") return 5;
if (s=="t") return 6;
if (s=="D") return -1;
if (s=="U") return -2;
if (s=="S") return -3;
if (s=="C") return -4;
if (s=="B") return -5;
if (s=="T") return -6;
if (s=="~o1") return 1000022;
if (s=="~o2") return 1000023;
if (s=="~o3") return 1000025;
if (s=="~o4") return 1000035;
if (s=="~1+") return 1000024;
if (s=="~2+") return 1000037;
if (s=="~1-") return -1000024;
if (s=="~2-") return -1000037;
return result;
}
double EstimateChannel(string s1, string s2, string s3, string s4)
{
double result = 0.;
for(int i=0; i<30 ;i++){
if (s1==PatChannel[i].prtcl1 && s2==PatChannel[i].prtcl2 && s3==PatChannel[i].prtcl3 && s4==PatChannel[i].prtcl4) {
result += PatChannel[i].weight;
}
}
return result;
}
//====================================== Begin Main Function ======================================
int main(int argc, char** argv)
{
//----------------------------------------------- Declare variables
// General variables
int fast=1;
double Beps=1.E-5, cutoff=0.0;
double OmegaMO,OmegaTree,OmegaFull,Xf;
int err, err2;
char mess[10];
char lsp[5];
double mychannel[7];
double m0,mhf,a0, tb;
double gMG1, gMG2, gMG3, gAl, gAt, gAb, sgn, gMHu, gMHd, gMl2;
double gMl3, gMr2, gMr3, gMq2, gMq3, gMu2, gMu3, gMd2, gMd3;
// My parameters
double mhiggs, mneut, mbottom, mCH;
int Zaehler = 1;
int indi, indj, itype, i;
// Parameters for my cross section and propagators
double t1,test[9];
int mflags[11], mflags2[8];
//----------------------------------------------- Create and open files
FILE * Ausgabe1;
FILE * Ausgabe2;
FILE * Ausgabe3;
FILE * Ausgabe4;
FILE * Ausgabe5;
Ausgabe1 = fopen("output1.txt", "w");
Ausgabe2 = fopen("output2.txt", "w");
Ausgabe4 = fopen("output4.txt", "w");
Ausgabe5 = fopen("output5.txt", "w");
//----------------------------------------------- Read cMSSM parameters and calculate spectrum
// Standard DD Scenario m0 = 500, mhf = 500, a0 = 0, tb = 12, sgn = +1
if(argc<5)
{
printf("\nStandard GUT-scale parameters (to be given): \n");
printf(" m0, mhf, A0, tanb, sgn(mu) \n\n");
exit(1);
}
else // mSUGRA parameters
{
sscanf(argv[1],"%lf",&m0);
sscanf(argv[2],"%lf",&mhf);
sscanf(argv[3],"%lf",&a0);
sscanf(argv[4],"%lf",&tb);
sscanf(argv[5],"%lf",&sgn);
}
if(sgn!=-1) sgn=1;
// Top-Masse einlesen
assignValW("Mtp",172.4);
//------------ Massen und Kopplungen setzen
// Gaugino-Massen werden gesetzt
gMG1 = mhf, gMG2 = mhf, gMG3 = mhf;
// Skalare Massen und Kopplungen an der GUT-Skala werden festgelegt.
gAl = a0, gAt = a0, gAb = a0;
gMl2 = m0, gMl3 = m0, gMr2 = m0, gMr3 = m0;
gMq2 = m0, gMq3 = m0, gMu2 = m0, gMd2 = m0, gMu3 = m0, gMd3 = m0;
gMHu = m0, gMHd = m0;
// Weitergabe an SPheno
err2 = SUGRAMODEL(RGE) (tb, gMG1, gMG2, gMG3, gAl, gAt, gAb, sgn, gMHu, gMHd, gMl2, gMl3, gMr2, gMr3, gMq2, gMq3, gMu2, gMu3, gMd2, gMd3);
//----------------------------------------------- Tests
// Test if model is ok
if(err2<0){
OmegaMO = (double)err;
OmegaTree = (double)err;
OmegaFull = (double)err;
}
// Test if LSP is neutralino
else
{
err = sortOddParticles(lsp);
if(strcmp(lsp,"~o1")!=0)
{
OmegaMO = -99.0;
OmegaTree = -99.0;
OmegaFull = -99.0;
}
//================================ Compute relic density ==========================================
else
{
if(err) { printf("Can't calculate %s ",mess); return 1;}
// Prerun
Ausgabe3 = fopen("output3.txt", "w+");
RDoption = 0; // SET RD OPTION HERE!
ChannelFilter = 0.02; // SET CHANNEL FILTER HERE!
DMNLOReady = false;
darkOmegaFO(&Xf,1,1.E-5);
printChannels(Xf,0.,Beps,1,Ausgabe3);
rewind(Ausgabe3);
DMNLOReady = true;
DMNLOAufrufe = 0;
CalcHepAufrufe = 0;
IntErrors = 0;
printf("Prerun finished \n");
// printf("Freeze Out XF = %8.5e \n",Xf);
// Set PatChannel for ShallImprove
char percentsign[3];
char prtcl1[3];
char prtcl2[3];
char prtcl3[3];
char prtcl4[3];
// printf("\nThe 30 most relevant channels are: \n");
for(i=0; i<30;i++){
fscanf(Ausgabe3,"%lf %s %s %s %s %s", &PatChannel[i].weight, percentsign, prtcl1, prtcl2, prtcl3, prtcl4);
PatChannel[i].weight = PatChannel[i].weight/100.;
PatChannel[i].prtcl1 = prtcl1;
PatChannel[i].prtcl2 = prtcl2;
PatChannel[i].prtcl3 = prtcl3;
PatChannel[i].prtcl4 = prtcl4;
// printf("%lf%% %s + %s -> %s + %s \n", 100.*PatChannel[i].weight, PatChannel[i].prtcl1.c_str(),
// PatChannel[i].prtcl2.c_str(), PatChannel[i].prtcl3.c_str(), PatChannel[i].prtcl4.c_str());
}
fclose(Ausgabe3);
printf("PatChannel set \n");
// Actual Relic Calculation
OmegaMO=darkOmega(&Xf,fast,Beps);
// RDoption=1;
// OmegaTree=darkOmega(&Xf,fast,Beps);
// RDoption=2;
// OmegaFull=darkOmega(&Xf,fast,Beps);
printf("Actual run finished \n");
// Find Masses
mbottom = findValW("MbMb");
mneut = findValW("MNE1");
mhiggs = findValW("Mh");
mCH = findValW("MHc");
// Control output
printf("\n");
printf(" ==== Recent Progress ==== \n");
printf("omega*h2 = %8.5e \n",OmegaMO);
printf("m0 = %8.5e \n",m0);
printf("mhf = %8.5e \n",mhf);
printf("RD = %i \n", RDoption);
printf("DMNLOAufrufe = %i \n", DMNLOAufrufe);
printf("CalcHepAufrufe = %i \n", CalcHepAufrufe);
printf("IntErrors = %i \n", IntErrors);
// Output in file
fprintf(Ausgabe1, " %i %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %i \n", Zaehler,
m0, mhf, OmegaMO, OmegaTree, OmegaFull, DMNLOContribution(), mneut, findValW("MNE2"),findValW("MC1"), mhiggs, mCH, findValW("MSt1"), IntErrors);
fprintf(Ausgabe2, " %i %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %i \n",
Zaehler, m0, mhf, DMNLOContribution(),
EstimateChannel("~o1", "~o1", "b", "B"), EstimateChannel("~o1", "~o1", "t", "T"),
EstimateChannel("~o1", "~o2", "b", "B"), EstimateChannel("~o1", "~o2", "t", "T"),
EstimateChannel("~o1", "~o3", "b", "B"), EstimateChannel("~o1", "~o3", "t", "T"),
EstimateChannel("~o2", "~o2", "b", "B"), EstimateChannel("~o2", "~o2", "t", "T"),
EstimateChannel("~o1", "~1+", "t", "B"), EstimateChannel("~o1", "~1+", "S", "c"),
EstimateChannel("~o1", "~1+", "u", "D"),
EstimateChannel("~o2", "~1+", "t", "B"), EstimateChannel("~o2", "~1+", "S", "c"),
EstimateChannel("~o2", "~1+", "u", "D"),
EstimateChannel("~1+", "~1-", "b", "B"), EstimateChannel("~1+", "~1-", "t", "T"),
EstimateChannel("~t1", "~T1", "G", "G"), EstimateChannel("~t1", "~t1", "t", "t"),
EstimateChannel("~o1", "~t1", "G", "t"), EstimateChannel("~t1", "~T1", "W+", "W-"),
EstimateChannel("~o1", "~t1", "t", "h"), EstimateChannel("~t1", "~T1", "h", "h"),
DMNLOChannels());
}
}
// }
// }
//=================================== My Cross Section ============================================
/*
t1 = 100.0;
// Loop over t1 (kinetic energy)
// for(t1=1.0;t1<=500.0;t1+=1.0)
// {
//----------------------------------------------- Set flags
for(indi=1; indi<=1; indi+=1)
{
for(indj=1; indj<=1; indj+=1)
{
for(itype=3; itype<=3; itype+=1)
{
//CHICHI2QQ CASE
// chargino (0-no charginos, 1- neutralino & chargino initial state, 2-charginos)
mflags[0] = 0;
// neutralino or chargino index - indi (in mixed initial state - chargino)
mflags[1] = indi;
// neutralino or chargino index - indj
mflags[2] = indj;
// type (1=neutrino, 2=electron, 3=up, 4=down) in mixed states it is up-type quark
mflags[3] = itype;
// generation
mflags[4] = 1;
// diagonalize (0 = Masses from MO, 1 = Diagonalize)
mflags[5] = 0;
// iflux (1=sigma 0=v.sigma)
mflags[6] = 1;
// only tree (1 = only tree, 0 = Loops)
mflags[7] = 0;
// Kinematics in MO Cross Section (0=Ours, 1=MO)
mflags[8] = 0;
// Resummation Style (0 Old DMNLO Style, 1 Spira/Patrick Style)
mflags[9] = 1;
// NNLO Resummation (1 = on, 0 = off)
mflags[10] = 1;
// Note: The use NNLO Res. = 1 sets Res. Style to 1 automatically.
// CalcHEP Aufruf
double cosmin=-1.0, cosmax=1.0, calchep;
double hbar_c2 = 389379660;
numout* channel;
channel=newProcess("~o1,~o1->t,T","n1n1_tT"); // Defining process
if(channel){
int ntot, err;
procInfo1(channel,&ntot,NULL,NULL);
calchep= cs22(channel,ntot,t1,cosmin,cosmax,&err); // Calculate cross section with CalcHEP
calchep= calchep/hbar_c2; // Conversion in 1/GeV^2
}
// Calculate cross section and write it in a file
chichi2qq_(t1,mflags,cutoff,test);
fprintf(Ausgabe4, "%i %i %i %i %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e \n",
Zaehler, indi, indj, itype, t1, calchep, test[0], test[1], test[2], test[3], test[4], test[5], test[6], test[7], test[8]);
// Control output
printf("\n");
printf(" ==== Recent Progress ==== \n");
printf("Zaehler = %i \n", Zaehler);
printf("t1 = %8.5e \n",t1);
printf("tb = %8.5e \n",tb);
printf("indi = %i \n",indi);
printf("indj = %i \n",indj);
printf("itype = %i \n",itype);
Zaehler++;
// End all the parameter loops
}
}
}
// }
*/
//=================================== Direct Detection ============================================
- double E_R = 10.0; //Recoil Energy in KeV
double pcm = 5.0; //Momentum in CMS frame in GeV
double costheta = 0.5; //Scattering angle
- double Qout = 10; //Renormalization scale
+ double Qout = 1000; //Renormalization scale
int Supp_Mode = 0;
int Loop_Mode = 1;
- int Denom_Mode = 0;
+ int Denom_Mode = 1;
+ int Box_Mode = 0;
double DDresults[32];
// for (Qout = 10; Qout<=10; Qout+=500) {
// for (pcm = 5; pcm<=5; pcm+=100) {
- directdetection_(E_R, pcm, costheta, Supp_Mode, Loop_Mode, Denom_Mode, Qout, DDresults);
+ directdetection_(pcm, costheta, Supp_Mode, Loop_Mode, Denom_Mode, Box_Mode, Qout, DDresults);
printf("Qout = %f \n", Qout);
printf("pcm = %f \n", pcm);
fprintf(Ausgabe5, "%8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e \n",
DDresults[0], DDresults[1], DDresults[2], DDresults[3], DDresults[4], DDresults[5], DDresults[6], DDresults[7]);
fprintf(Ausgabe5, "%8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e \n",
DDresults[8], DDresults[9], DDresults[10], DDresults[11], DDresults[12], DDresults[13], DDresults[14], DDresults[15]);
fprintf(Ausgabe5, "%8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e \n",
DDresults[16], DDresults[17], DDresults[18], DDresults[19], DDresults[20], DDresults[21], DDresults[22], DDresults[23]);
fprintf(Ausgabe5, "%8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e \n",
DDresults[24], DDresults[25], DDresults[26], DDresults[27], DDresults[28], DDresults[29], DDresults[30], DDresults[31]);
// }
-// }
-
-
-/*
- int count;
- double Danger;
- double loopresults[20];
- int Zaehler2 = 1;
-
- for(count = 1000; count >= -1000; count -= 1){
- Danger = count*0.0001;
- looptest_(Danger, loopresults);
- fprintf(Ausgabe5, "%i %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e %8.9e \n",
- Zaehler2, Danger, loopresults[0], loopresults[1], loopresults[2], loopresults[3],
- loopresults[4], loopresults[5], loopresults[6], loopresults[7],loopresults[8],
- loopresults[9], loopresults[10], loopresults[11], loopresults[12], loopresults[13],
- loopresults[14], loopresults[15], loopresults[16], loopresults[17], loopresults[18],
- loopresults[19]);
- Zaehler2++;
- } */
-
-
-// int count;
-// double Danger;
-// double loopresults[20];
-// int Zaehler2 = 1;
-
-// for(count = 10; count >= 0; count -= 1){
-// Danger = 100.0;
-// Danger = -0.;
-// looptest2_(Danger, loopresults);
-// fprintf(Ausgabe5, "%i %8.9e %8.9e \n", Zaehler2, Danger, loopresults[0]);
-// Zaehler2++;
-// }
-
+// }
+
//==================================== End of program =============================================
// Close files
fclose(Ausgabe1);
fclose(Ausgabe2);
fclose(Ausgabe4);
fclose(Ausgabe5);
return 0;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Sep 30, 5:44 AM (7 h, 13 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6566313
Default Alt Text
(22 KB)

Event Timeline

HEPForge.org · .