#define DATUM "2010_11_09" /**************************************************** This is the file '2Dfit.c.2010_11_09', genuine name: '2Dfit.c'. The aim of this software is the analysis of the experimental results obtained at the beamtimes in July 2009 December 2009 July 2010 at the surface diffraction beamline of the Max-Planck-Institut fuer Metallforschung at the ANKA synchrotron facility in Karlsruhe. The program '2Dfit' was written at the Max-Planck-Institut fuer Metallforschung (Stuttgart). The model function in this program consists of three peaks, all three are the sum of two 2D-Gaussians with the same peak coordinats. The backgrund is constant. The last version with the single 2D-Gaussian plus constant background is the version 2Dfit.c.2010_11_04. Authors: Agnes Szoekefalvi-Nagy, Janos Major (major@mf.mpg.de) The newest version of the source code '2Dfit.c' can freely be obtained from the web page of the Max Planck Institute for Metals Research (www.mf.mpg.de), page of the Department for Low Dimensional and Metastable Materials, page "software". The program '2Dfit' is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. The included Fortran fit routine minuit is a GPL (GNU General Public Licence, http://www.gnu.org/copyleft/gpl.html) software of CERN (Geneva, Switzerland) (http://consult.cern.ch/writeup/minuit/). The c-minuit coupling of the Fortran minuit routine with the C language program, the c-minuit system was developed by Andras Major (http://c-minuit.sourceforge.net). For installation of the minuit C bibliothek see the man page of the superfit program and the file c-minuit-CVS-20031201.tar.gz (web page of the Max Planck Institute for Metals Research (www.mf.mpg.de), page of the Department for Low Dimensional and Metastable Materials, page "software"). This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You can obtain a copy of the GNU General Public License from http://www.gnu.org/copyleft/gpl.html; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. ************************************************** First version: 2Dfit.c.2010.09.04 which is based on the program: file_read.c.2010_08_30 The standard site for this file is: $HOME/c-minuit// To compile and link please use: make install *************************************************/ /* for monitoring all errors: int ERROR = 0; for stopping at the first error: int ERROR = 1; */ int ERROR = 0; #define __NO_MATH_INLINES /* makes slower but, may be, more precise */ #include #include #include #include #include #include #include "cfortran.h" #include "minuitfcn.h" #include "minuit.h" #define EICHUNG 0.00804 /* mythen detector: degree/channel */ #define MAX_FILE 5100 /* max. number of files */ #define MAX_PIXEL 5000 /* max. number a pixels along the ppm axes */ #define MAX_CHANNEL 1400 /* max. number of detector channels*/ #define MAX_LEN 256 /* max. number of input characters */ #define EFFZERO 1e-40 /* my definition of zero */ #define NEARZERO 1e-10 /* zero for the matrix functions */ #define LIM_EXP 70. /* limit for the function exp */ #define SCREEN_WRITE 50 /* screen write out period during interpolation */ #define PLOTZERO 1.e-6 /* smaller values are not plotted in gnuplot */ #define RAD 0.0174532925199433 /* 1 degree in radian */ #define PI 3.1415926535897932384626433832795028841971693993751 #define MAX_MESSPUNKT 10000000 /* max. number of measuring points */ #define MAX_FIT 50000 /* max. number of fitted plots */ #define NM 1000 /* max. number of fit parameters */ void gnuplot_read_dat_files_t (void); void gnuplot_read_dat_files_q (void); /* function prototypes: */ void message (char *text, char *text1); void error (char *text); void terminate (char *text); void wait (char *text); void read_dat_files(void); void read_fit_dat_files(void); double gauss2D (double A, double x, double y, double x0, double y0, double sx, double sy, double theta); void fcn(int npar,double* grad,double* fcnval,double* x,int iflag,void* futil); int tr(int); double power(double, double); void header (int); /* global declarations: */ int mode; double start[NM], step[NM], min[NM], max[NM]; /* do not change this line */ char parname[NM][MAX_LEN]; /* names of the MINUIT pars */ double fcn_end, fcnred_end; /* minfit results */ int n_par; /* minfit results */ double HH,KK,LL,MON[50000],DET[50000], MMON[50000], MMONITOR; /* fit data */ double A1, H1, K1, SH1, SK1, theta1, B; double A2, H2, K2, SH2, SK2, theta2; double A3, H3, K3, SH3, SK3, theta3; double A4, SH4, SK4, theta4; double A5, SH5, SK5, theta5; double A6, SH6, SK6, theta6; double A1_fit[MAX_FIT], H1_fit[MAX_FIT], K1_fit[MAX_FIT]; double SH1_fit[MAX_FIT], SK1_fit[MAX_FIT], B_fit[MAX_FIT]; double A1_err[MAX_FIT], H1_err[MAX_FIT], K1_err[MAX_FIT]; double SH1_err[MAX_FIT], SK1_err[MAX_FIT], theta1_fit[MAX_FIT]; double theta1_err[MAX_FIT], B_err[MAX_FIT]; double A2_fit[MAX_FIT], H2_fit[MAX_FIT], K2_fit[MAX_FIT]; double SH2_fit[MAX_FIT], SK2_fit[MAX_FIT]; double A2_err[MAX_FIT], H2_err[MAX_FIT], K2_err[MAX_FIT]; double SH2_err[MAX_FIT], SK2_err[MAX_FIT], theta2_fit[MAX_FIT]; double theta2_err[MAX_FIT]; double A3_fit[MAX_FIT], H3_fit[MAX_FIT], K3_fit[MAX_FIT]; double SH3_fit[MAX_FIT], SK3_fit[MAX_FIT]; double A3_err[MAX_FIT], H3_err[MAX_FIT], K3_err[MAX_FIT]; double SH3_err[MAX_FIT], SK3_err[MAX_FIT], theta3_fit[MAX_FIT]; double theta3_err[MAX_FIT]; double A4_fit[MAX_FIT]; double A4_err[MAX_FIT]; double SH4_fit[MAX_FIT], SK4_fit[MAX_FIT]; double SH4_err[MAX_FIT], SK4_err[MAX_FIT], theta4_fit[MAX_FIT]; double theta4_err[MAX_FIT]; double A5_fit[MAX_FIT]; double A5_err[MAX_FIT]; double SH5_fit[MAX_FIT], SK5_fit[MAX_FIT]; double SH5_err[MAX_FIT], SK5_err[MAX_FIT], theta5_fit[MAX_FIT]; double theta5_err[MAX_FIT]; double A6_fit[MAX_FIT]; double A6_err[MAX_FIT]; double SH6_fit[MAX_FIT], SK6_fit[MAX_FIT]; double SH6_err[MAX_FIT], SK6_err[MAX_FIT], theta6_fit[MAX_FIT]; double theta6_err[MAX_FIT]; double chired_fit[MAX_FIT]; double Hf [MAX_FIT], Kf [MAX_FIT], Lf [MAX_FIT]; double FCNRED[MAX_FIT]; int counts [MAX_FILE][MAX_CHANNEL]; /* experimental results */ double q_x[MAX_FILE][MAX_CHANNEL]; double q_y[MAX_FILE][MAX_CHANNEL]; double eichung=EICHUNG; /* PSD degree/channel */ double channel_m; /* PSD reference channel tth=0->direct beam*/ double PPM[MAX_PIXEL][MAX_PIXEL]; /* array to be transfered to a ppm file */ double GNU[MAX_PIXEL][MAX_PIXEL]; /* array to be transferred to gnuplot */ double a_i[MAX_FILE][MAX_CHANNEL]; /* the corresponding q_x values */ double a_f[MAX_FILE][MAX_CHANNEL]; /* the corresponding alpha_j values */ double H[MAX_MESSPUNKT], K[MAX_MESSPUNKT], L[MAX_MESSPUNKT]; double maxdm=0; /* MAX(det/mon) */ int Time[MAX_MESSPUNKT], mon[MAX_MESSPUNKT], det[MAX_MESSPUNKT]; char *datfile; /* first (common) part of the file names */ char *datpath; /* path of the file names */ char ppmfile [MAX_LEN]; char psfile [MAX_LEN]; char pngfile [MAX_LEN]; char gnuplotfile [MAX_LEN]; char gnuplot2Dfile [MAX_LEN]; int number, t, tmax, tfile; double COUNTSperMMON [5000][50000]; /* exp. counts[shot][time]*/ double threshold; /* between 0 and 1: only counts > threshold*maxdm are included in the fit */ int errorcount = 0; /* for the function error */ int gnuplotversion; /* 0 if >= 4.0; >0 if <=3.7 */ int fcncount=0; char headerchar[2000]; /********************** function: main ***************/ int main(int argc, char *argv[]){ char *kontr; if (argc != 7) {message("\n You have made a mistake in the input,", "please try again using the format:\n"); message (" ", "2Dfit \n"); message (" =","1: reads data files and makes nothing else"); message (" ","2: makes 1D graphic (png) files by gnuplot "); message (" "," (axes: t, counts/mean_monitor) "); message (" "," (for autoscaling please edit the source code)"); message (" ","3: makes 2D graphic (png) files by gnuplot "); message (" "," (axes: H, K; colour: counts/mean_monitor) "); message (" "," (for autoscaling please edit the source code)"); message (" ","4: fit + 2D_DATA png files "); message (" ","5: fit"); message (" =","path of the data files (e.g., . or ../ or DAT)"); message (" =","first (common) part of the file names"); message (" =","number of mesh points = number of files"); message (" ","0: all available files"); message (" =","number of time points = number of lines in the files"); message (" ","0: all available time points (lines)"); message ("=","between 0 and 1: only COUNTSperMMON>threshold*maxdm"); message (" ","will be included in the fit (0: no threshold)"); message ( "Attention: all files produced by the present program will be located in\n", " the directory where the program is started. Existing files with"); message ( " the same name will be overwritten without warning.",""); fprintf (stderr, "Program version: 2Dfit.c." DATUM "\n"); message (" ", "2Dfit \n"); terminate (" input error ");} mode = strtod(argv[1], &kontr); if(kontr == argv[1]) terminate(" wrong first parameter\n"); datpath = argv[2]; message("path of the data files= ", datpath); datfile = argv[3]; message("common part of the file names=", datfile); number = strtod(argv[4], &kontr); if(kontr == argv[4]) terminate(" wrong 4th parameter\n"); tmax = strtod(argv[5], &kontr); if(kontr == argv[5]) terminate(" wrong 5th parameter\n"); threshold = strtod(argv[6], &kontr); if(kontr == argv[6]) terminate(" wrong 5th parameter\n"); if ((threshold<0)||(threshold>1)) terminate(" threshold parameter is out of range\n"); /* the gnuplot version will be asked: */ {FILE *pipe; pipe = popen("gnuplot", "w"); fprintf(pipe, "set terminal png \n"); gnuplotversion = pclose(pipe); if (gnuplotversion == 0) message("A new gnuplot version is present.","\n"); else message ("An old gnuplot version is present: no png files will be available.", "\n");} /* ready: gnuplot version */ switch (mode){ case 1: /* this case only for a general test */ read_dat_files (); break; case 2: read_dat_files (); /* maxdm (range of the counts/mon) and MMON[ii] = mean_monitor will be defined */ gnuplot_read_dat_files_t (); break; case 3: read_dat_files (); /* maxdm (range of the counts) and MMON[ii] = mean_monitor will be defined */ gnuplot_read_dat_files_q (); break; case 4: /* fit the 2D Gaussian function to the constant t data*/ read_dat_files (); /* maxdm (range of the counts/mon) and writes ***t.png 2D files and MMON[ii] = mean_monitor will be defined */ read_fit_dat_files (); break; case 5: /* fit the 2D Gaussian function to the constant t data*/ read_dat_files (); /* maxdm (range of the counts/mon) and MMON[ii] = mean_monitor will be defined */ read_fit_dat_files (); break; default: terminate(" Illegal mode value. "); break;} message ("Program version: file_read.c." DATUM, "\n" ); error ("terminal"); return 0;} /********************** function: main ***************/ /*************************** header ****************************/ void header (int i) {/* prepares the string headerfile to be copied as a header for the out files i = 1 for fitting i = 2 for not fitting */ char c[2000]; time_t curtime; curtime = time (NULL); strcpy(headerchar, "# "); strcpy(c, asctime (localtime (&curtime))); strcat(headerchar,c); strcpy(c, "# First (common) part of the data file names: \n"); strcat(headerchar,c); strcpy(c, "# "); strcat(headerchar,c); strcpy(c, datpath); strcat(headerchar,c); strcpy(c, datfile); strcat(headerchar,c); strcpy(c, "\n"); strcat(headerchar,c); sprintf (c, "# The initial time values in experimental units: \n"); strcat(headerchar,c); sprintf (c, "# Further initial parameters: \n"); strcat(headerchar,c); sprintf (c, "# mode=%1i number=%4i 1=%4i number=%4i threshold=%.2f \n", mode, number, 1, number, threshold); strcat(headerchar,c); strcpy(c, "# This file was produced by the software '' version " DATUM "\n"); strcat(headerchar,c); // if (i==1) { //if (strcmp("MODEL","model1")); //strcpy(c, "# The used model function for the fit is the 'model1',\n" // "# which is the stretched exponential fuction \n"); //if (strcmp("MODEL","model2")); //strcpy(c, "The used model function for the fit is the 'model2',\n" // "# which is the sum of two exponentials \n"); //if (strcmp("MODEL","model3")); //strcpy(c, "# The used model function for the fit is the 'model3',\n" // "# which is the sum of three exponentials;\n" // "# the third one can be delayed and gives zero for negative \n" // "# arguments (Heaviside function) \n"); // strcat(headerchar,c); } strcpy(c, "# \n"); strcat(headerchar,c);} /*************************** header ****************************/ /*************************** gauss2D ****************************/ double gauss2D(double A, double x , double y, double x0, double y0, double sx, double sy, double theta) { /* model function: 2D Gaussian function A = amplitude; sx = sigma_x; sy = sigma_y theta = rotation angle of the ellipse in degrees */ double dummy, a, b, c; theta = theta * PI / 180.; /* degree -> radian */ a = cos(theta)*cos(theta)/(2.*sx*sx) + sin(theta)*sin(theta)/(2.*sy*sy); b = -sin(2.*theta)/(4.*sx*sx) + sin(2.*theta)/(4.*sy*sy); c = sin(theta)*sin(theta)/(2.*sx*sx) + cos(theta)*cos(theta)/(2.*sy*sy); if (sx < EFFZERO) terminate ("sigma_x is too small in model function"); if (sy < EFFZERO) terminate ("sigma_y is too small in model function"); // dummy = -(x-x0)*(x-x0)/(2.*sx*sx) -(y-y0)*(y-y0)/(2.*sy*sy); dummy = -a*(x-x0)*(x-x0)+2.*b*(x-x0)*(y-y0)-c*(y-y0)*(y-y0); if (dummy < -100) dummy =0; else dummy = exp(dummy)/(2.*PI*sx*sy); return (A*dummy);} /*************************** gauss2D ****************************/ /**************************** fcn *****************************/ void fcn(int npar,double* grad,double* fcnval,double* x,int iflag,void* futil) { int i; double Dev; A1 = x[0]; H1 = x[1]; K1 = x[2]; SH1 = x[3]; SK1 = x[4]; theta1 = x[5]; A2 = x[6]; H2 = x[7]; K2 = x[8]; SH2 = x[9]; SK2 = x[10]; theta2 = x[11]; A3 = x[12]; H3 = x[13]; K3 = x[14]; SH3 = x[15]; SK3 = x[16]; theta3 = x[17]; A4 = x[18]; SH4 = x[19]; SK4 = x[20]; theta4 = x[21]; A5 = x[22]; SH5 = x[23]; SK5 = x[24]; theta5 = x[25]; A6 = x[26]; SH6 = x[27]; SK6 = x[28]; theta6 = x[29]; B = x[30]; /* here begins the function to be minimized (*fcnval=) : */ *fcnval=0.; /* Maximum Likelihood fit with Poisson statistics for (i=1; i<=number; i++) { Mod = MODEL(A,Hf[i],Kf[i],H0,K0,SH,SK,theta)+B; Exp = COUNTSperMMON[i][t] * MMON[i]; *fcnval+= Mod - Exp; if (Exp>1.e-6) *fcnval+=Exp*log(Exp/Mod); } *fcnval = *fcnval*2.; */ /* least squares fit for non-constant monitor intensities: */ for (i=1; i<=number; i++) { if (COUNTSperMMON[i][t]>maxdm*threshold){ Dev = COUNTSperMMON[i][t] - (gauss2D(A1,Hf[i],Kf[i],H1,K1,SH1,SK1,theta1) + gauss2D(A2,Hf[i],Kf[i],H2,K2,SH2,SK2,theta2) + gauss2D(A3,Hf[i],Kf[i],H3,K3,SH3,SK3,theta3) + gauss2D(A4,Hf[i],Kf[i],H1,K1,SH4,SK4,theta4) + gauss2D(A5,Hf[i],Kf[i],H2,K2,SH5,SK5,theta5) + gauss2D(A6,Hf[i],Kf[i],H3,K3,SH6,SK6,theta6) + B); *fcnval+=MMON[i]*MMON[i]*Dev*Dev/(COUNTSperMMON[i][t]*MMON[i]+1.0);}} if (*fcnval == 0) terminate("*fcnval == 0 --> the fit does not work properly! \n"); /* print to screen: fit progress: fcncount = fcncount + 1; if (fcncount%SCREEN_WRITE==0) printf("fcn status: (%4i) fcn = %12.4f fcn_red = %12.4f \n", fcncount, *fcnval, *fcnval/(tmax-npar)); */ if (iflag==3){fcn_end = *fcnval; fcnred_end = *fcnval/(number-npar); n_par=npar;}} /**************************** fcn *****************************/ /********************** function: read_fit_dat_files ***************/ void read_fit_dat_files () { /* experimental counts will be read in: */ int i, ii, lineindex, j; char filename [MAX_LEN]; char line [MAX_LEN+1]; float fH, fK, fL; int imon, idet; char Dummy[MAX_LEN]; int dummy; double ddummy, monitor; int idummy, maxcall, parameter; char strategy[20], sdummy[MAX_LEN+1]; FILE *file; FILE *pipe; char gfilename [MAX_LEN]; for (ii=1; ii<=number; ii++) { sprintf(filename, "%s/%s%i%s", datpath, datfile, ii, ".dat"); file=fopen(filename,"r"); lineindex=0; if (file==NULL) terminate(" The desired file does not exist. "); t = 0; while (fgets (line, MAX_LEN, file) != NULL) { lineindex=lineindex+1; if (lineindex==9) { sscanf (line, "%s%f%f%f", Dummy, &fH, &fK, &fL); Hf[ii] = fH; Kf[ii] = fK; Lf[ii] = fL; /* printf("%i H= %f K=%f L=%f time=%i mon=%i det=%i \n ", index,H[index],K[index],L[index],time[index],mon[index],det[index]); */} else if (lineindex>12) { sscanf (line, "%i%i", &imon, &idet); t=t+1; COUNTSperMMON [ii][t] = idet/MMON [ii]; } } fclose(file);} tfile = t; printf("########## tfile=%i tmax=%i ##########\n", tfile, tmax); for (t=1; t<=tmax; t++) { fcncount = 0; MNINIT(5,6,7); /* do not change this line, please */ /* minuit parameter definitions : */ i = 1; sprintf(parname[i],"A1"); start[i]=0.00029; step[i]=0.0001; min[i]=1.e-8; max[i]=1; i = 2; sprintf(parname[i],"H1"); start[i]=1.87; step[i]=.1; min[i]=1.865; max[i]=1.895; i = 3; sprintf(parname[i],"K1"); start[i]=.0005; step[i]=.000001; min[i]=-0.001; max[i]=0.01; i = 4; sprintf(parname[i],"SH1"); start[i]=.01; step[i]=.001; min[i]=.0001; max[i]=0.08; i = 5; sprintf(parname[i],"SK1"); start[i]=.0001; step[i]=.00001; min[i]=.0001; max[i]=0.001; i = 6; sprintf(parname[i],"theta1"); start[i]=.0; step[i]=0.1; min[i]=-40.; max[i]=40.; i = 7; sprintf(parname[i],"A2"); start[i]=0.000; step[i]=0.000; min[i]=1.e-8; max[i]=1; i = 8; sprintf(parname[i],"H2"); start[i]=1.87; step[i]=.0; min[i]=1.865; max[i]=1.895; i = 9; sprintf(parname[i],"K2"); start[i]=.0005; step[i]=.00000; min[i]=-0.001; max[i]=0.01; i =10; sprintf(parname[i],"SH2"); start[i]=.01; step[i]=.00; min[i]=.0001; max[i]=0.08; i =11; sprintf(parname[i],"SK2"); start[i]=.0001; step[i]=.0000; min[i]=.0001; max[i]=0.001; i =12; sprintf(parname[i],"theta2"); start[i]=.0; step[i]=.0; min[i]=-40.; max[i]=40.; i =13; sprintf(parname[i],"A3"); start[i]=0.000; step[i]=0.000; min[i]=1.e-8; max[i]=1; i =14; sprintf(parname[i],"H3"); start[i]=1.87; step[i]=.0; min[i]=1.865; max[i]=1.895; i =15; sprintf(parname[i],"K3"); start[i]=.0005; step[i]=.00000; min[i]=-0.001; max[i]=0.01; i =16; sprintf(parname[i],"SH3"); start[i]=.01; step[i]=.001; min[i]=.0001; max[i]=0.08; i =17; sprintf(parname[i],"SK3"); start[i]=.0001; step[i]=.0000; min[i]=.0001; max[i]=0.001; i =18; sprintf(parname[i],"theta3"); start[i]=.0; step[i]=.0; min[i]=-40.; max[i]=40.; i =19; sprintf(parname[i],"A4"); start[i]=0.000; step[i]=0.0001; min[i]=1.e-8; max[i]=1; i =20; sprintf(parname[i],"SH4"); start[i]=.01; step[i]=.001; min[i]=.0001; max[i]=0.08; i =21; sprintf(parname[i],"SK4"); start[i]=.0001; step[i]=.00001; min[i]=.0001; max[i]=0.001; i =22; sprintf(parname[i],"theta4"); start[i]=.0; step[i]=.01; min[i]=-40.; max[i]=40.; i =23; sprintf(parname[i],"A5"); start[i]=0.000; step[i]=0.000; min[i]=1.e-8; max[i]=1; i =24; sprintf(parname[i],"SH5"); start[i]=.01; step[i]=.00; min[i]=.0001; max[i]=0.08; i =25; sprintf(parname[i],"SK5"); start[i]=.0001; step[i]=.0000; min[i]=.0001; max[i]=0.001; i =26; sprintf(parname[i],"theta5"); start[i]=.0; step[i]=.0; min[i]=-40.; max[i]=40.; i =27; sprintf(parname[i],"A6"); start[i]=0.000; step[i]=0.000; min[i]=1.e-8; max[i]=1; i =28; sprintf(parname[i],"SH6"); start[i]=.01; step[i]=.00; min[i]=.0001; max[i]=0.08; i =29; sprintf(parname[i],"SK6"); start[i]=.0001; step[i]=.0000; min[i]=.0001; max[i]=0.001; i =30; sprintf(parname[i],"theta6"); start[i]=.0; step[i]=.0; min[i]=-40.; max[i]=40.; i =31; sprintf(parname[i], "B"); start[i]=.685; step[i]=0; min[i]=1e-3; max[i]=100.; for (j=1; j<=31; j++) MNPARM(j, parname[j], start[j], step[j], min[j], max[j], dummy); parameter = 2; /* 0: minimum number of fcn calls 1: (minuit default) 2: most precise results */ sprintf(strategy, "%s %i", "SET STRATEGY", parameter); MNCOMD(minuitfcn, strategy, dummy, NULL); /* default: SET WARNINGS SET NOWARNINGS */ sprintf(strategy, "%s", "SET NOWARNINGS"); MNCOMD(minuitfcn, strategy, dummy, NULL); parameter = 0; /* 0: minimum output 1: default value, normal output 2: additional output 3: maximum output*/ sprintf(strategy, "%s %i", "SET PRINTOUT", parameter); MNCOMD(minuitfcn, strategy, dummy, NULL); parameter = 1; /* 0.5: -log Likelihood (1 sigma) 1 : (minuit default) chisquare function (1 sigma) 4 : chisquare function (2 sigma) */ sprintf(strategy, "%s %i", "SET ERRORDEF", parameter); MNCOMD(minuitfcn, strategy, dummy, NULL); maxcall = 5000; sprintf(strategy, "%s %i", "SIMPLEX", maxcall); MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ maxcall = 5000; sprintf(strategy, "%s %i", "MINIMIZE", maxcall); MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ sprintf(strategy, "%s %i", "CALL", 3); MNCOMD(minuitfcn, strategy, dummy, NULL); /* results of the fit */ for (j=1; j<=31; j++) /* fit results */ MNPOUT(j, sdummy, start[j], step[j], ddummy, ddummy, idummy); A1_fit [t] = start[ 1]; H1_fit[t] = start[ 2]; K1_fit[t] = start[ 3]; SH1_fit[t] = start[ 4]; SK1_fit[t] = start[ 5]; theta1_fit[t] = start[ 6]; A2_fit [t] = start[ 7]; H2_fit[t] = start[ 8]; K2_fit[t] = start[ 9]; SH2_fit[t] = start[10]; SK2_fit[t] = start[11]; theta2_fit[t] = start[12]; A3_fit [t] = start[13]; H3_fit[t] = start[14]; K3_fit[t] = start[15]; SH3_fit[t] = start[16]; SK3_fit[t] = start[17]; theta3_fit[t] = start[18]; A4_fit [t] = start[19]; SH4_fit[t] = start[20]; SK4_fit[t] = start[21]; theta4_fit[t] = start[22]; A5_fit [t] = start[23]; SH5_fit[t] = start[24]; SK5_fit[t] = start[25]; theta5_fit[t] = start[26]; A6_fit [t] = start[27]; SH6_fit[t] = start[28]; SK6_fit[t] = start[29]; theta6_fit[t] = start[30]; B_fit [t] = start[31]; A1_err [t] = step [ 1]; H1_err[t] = step [ 2]; K1_err[t] = step [ 3]; SH1_err[t] = step [ 4]; SK1_err[t] = step [ 5]; theta1_err[t] = step [ 6]; A2_err [t] = step [ 7]; H2_err[t] = step [ 8]; K2_err[t] = step [ 9]; SH2_err[t] = step [10]; SK2_err[t] = step [11]; theta2_err[t] = step [12]; A3_err [t] = step [13]; H3_err[t] = step [14]; K3_err[t] = step [15]; SH3_err[t] = step [16]; SK3_err[t] = step [17]; theta3_err[t] = step [18]; A4_err [t] = step [19]; SH4_err[t] = step [20]; SK4_err[t] = step [21]; theta4_err[t] = step [22]; A5_err [t] = step [23]; SH5_err[t] = step [24]; SK5_err[t] = step [25]; theta5_err[t] = step [26]; A6_err [t] = step [27]; SH6_err[t] = step [28]; SK6_err[t] = step [29]; theta6_err[t] = step [30]; B_err [t] = step [31]; FCNRED[t] = fcnred_end; printf("########## tfile=%i tmax=%i n_par=%i fcnred=%g ##########\n", tfile, tmax, n_par, fcnred_end);} /* writes the fit result file: */ header(1); file = fopen ("fitresult", "w"); // fprintf (file, "%s", headerchar); fprintf (file, "# 1:A1 2:A1err 3:H1 4:H1err 5:K1 6:K1err " "7:SH1 8:SH1err 9:SK1 10:SK1err 11:theta1 12:theta1_err 13:max\n" "#14:A2 15:A2err 16:H2 17:H2err 18:K2 19:K2err " "20:SH2 21:SH2err 22:SK2 23:SK2err 24:theta2 25:theta2_err 26:max\n" "#27:A3 28:A3err 29:H3 30:H3err 31:K3 32:K3err " "33:SH3 34:SH3err 35:SK3 36:SK3err 37:theta3 38:theta3_err 39:max\n" "#40:A4 41:A4err " "42:SH4 43:SH4err 44:SK4 45:SK4err 46:theta4 47:theta4_err 48:max\n" "#49:A5 50:A5err " "51:SH5 52:SH5err 53:SK5 54:SK5err 55:theta5 56:theta5_err 57:max\n" "#58:A5 59:A5err " "60:SH6 61:SH6err 62:SK6 63:SK6err 64:theta6 65:theta6_err 66:max\n" "#67:B 68:Berr 69:maxdm 70:fcnred \n"); for (t=1; t<=tmax; t++){ if (MMON[t]>EFFZERO) monitor = MMON[t]; else monitor = 1.; fprintf (file, " %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e " " %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e " " %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e " " %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e " " %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e " " %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e " " %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e \n", A1_fit [t], A1_err [t], H1_fit[t], H1_err[t], K1_fit[t], K1_err[t], SH1_fit[t], SH1_err[t], SK1_fit[t], SK1_err[t], theta1_fit [t], theta1_err [t], A1_fit [t]/(2.*PI*SH1_fit[t]*SK1_fit[t]), A2_fit [t], A2_err [t], H2_fit[t], H2_err[t], K2_fit[t], K2_err[t], SH2_fit[t], SH2_err[t], SK2_fit[t], SK2_err[t], theta2_fit [t], theta2_err [t], A2_fit [t]/(2.*PI*SH2_fit[t]*SK2_fit[t]), A3_fit [t], A3_err [t], H3_fit[t], H3_err[t], K3_fit[t], K3_err[t], SH3_fit[t], SH3_err[t], SK3_fit[t], SK3_err[t], theta3_fit [t], theta3_err [t], A3_fit [t]/(2.*PI*SH3_fit[t]*SK3_fit[t]), A4_fit [t], A4_err [t], SH4_fit[t], SH4_err[t], SK4_fit[t], SK4_err[t], theta4_fit [t], theta4_err [t], A4_fit [t]/(2.*PI*SH4_fit[t]*SK4_fit[t]), A5_fit [t], A5_err [t], SH5_fit[t], SH5_err[t], SK5_fit[t], SK5_err[t], theta5_fit [t], theta5_err [t], A5_fit [t]/(2.*PI*SH5_fit[t]*SK5_fit[t]), A6_fit [t], A6_err [t], SH6_fit[t], SH6_err[t], SK6_fit[t], SK6_err[t], theta6_fit [t], theta6_err [t], A6_fit [t]/(2.*PI*SH6_fit[t]*SK6_fit[t]), B_fit [t], B_err [t], maxdm, FCNRED[t]);} fprintf (file, "# 1:A1 2:A1err 3:H1 4:H1err 5:K1 6:K1err " "7:SH1 8:SH1err 9:SK1 10:SK1err 11:theta1 12:theta1_err 13:max\n" "#14:A2 15:A2err 16:H2 17:H2err 18:K2 19:K2err " "20:SH2 21:SH2err 22:SK2 23:SK2err 24:theta2 25:theta2_err 26:max\n" "#27:A3 28:A3err 29:H3 30:H3err 31:K3 32:K3err " "33:SH3 34:SH3err 35:SK3 36:SK3err 37:theta3 38:theta3_err 39:max\n" "#40:A4 41:A4err " "42:SH4 43:SH4err 44:SK4 45:SK4err 46:theta4 47:theta4_err 48:max\n" "#49:A5 50:A5err " "51:SH5 52:SH5err 53:SK5 54:SK5err 55:theta5 56:theta5_err 57:max\n" "#58:A5 59:A5err " "60:SH6 61:SH6err 62:SK6 63:SK6err 64:theta6 65:theta6_err 66:max\n" "#67:B 68:Berr 69:maxdm 70:fcnred \n"); fclose(file); for (t=1; t<=tmax; t++){ sprintf(gfilename, "%s%04i%s", datfile,t,"t.gnu"); header(2); file = fopen (gfilename, "w"); fprintf (file, "%s", headerchar); for (ii=1; ii<=number; ii++){ fprintf (file, " %.8f \t %.8f \t %f \n", Hf[tr(ii)], Kf[tr(ii)], COUNTSperMMON [tr(ii)][t]); if (ii%31==0) fprintf (file, " \n");} fclose(file);} if (mode==4) for (t=1; t<=tmax; t++) { sprintf( filename, "%s%04i%s", datfile,t,"t.png"); sprintf(gfilename, "%s%04i%s", datfile,t,"t.gnu"); pipe = popen("gnuplot -geometry 900x900+10+10", "w"); /* fprintf (pipe, "set logscale zcb \n"); */ fprintf (pipe, "set title '<%s t=%i micro_s>' \n", filename, t); /* fprintf(pipe, "set terminal postscript color solid lw 2\n"); */ fprintf(pipe, "set terminal png \n"); /* fprintf(pipe, "set terminal gif \n"); */ fflush(pipe); fprintf(pipe, "set output '%s' \n", filename); fprintf(pipe, "set view 60,45 \n"); /* fprintf(pipe, "set view 0,0 \n"); */ fflush(pipe); /* version constant scaling: */ fprintf (pipe, "splot [][][0:%f] '%s' t '' with pm3d \n", ceil(maxdm), gfilename); /* version auto scaling: */ /* fprintf (pipe, "splot '%s' t '' with pm3d \n", gfilename); */ fclose (pipe);} for (t=1; t<=tmax; t++){ sprintf(gfilename, "%s%04i%s", datfile,t,"f.gnu"); header(2); file = fopen (gfilename, "w"); fprintf (file, "%s", headerchar); for (ii=1; ii<=number; ii++){ fprintf (file, " %.8f \t %.8f \t %f \n", Hf[tr(ii)], Kf[tr(ii)], (gauss2D(A1_fit[t],Hf[tr(ii)],Kf[tr(ii)],H1_fit[t],K1_fit[t], SH1_fit[t],SK1_fit[t],theta1_fit[t]) + gauss2D(A2_fit[t],Hf[tr(ii)],Kf[tr(ii)],H2_fit[t],K2_fit[t], SH2_fit[t],SK2_fit[t],theta2_fit[t]) + gauss2D(A3_fit[t],Hf[tr(ii)],Kf[tr(ii)],H3_fit[t],K3_fit[t], SH3_fit[t],SK3_fit[t],theta3_fit[t]) + gauss2D(A4_fit[t],Hf[tr(ii)],Kf[tr(ii)],H1_fit[t],K1_fit[t], SH4_fit[t],SK4_fit[t],theta4_fit[t]) + gauss2D(A5_fit[t],Hf[tr(ii)],Kf[tr(ii)],H2_fit[t],K2_fit[t], SH5_fit[t],SK5_fit[t],theta5_fit[t]) + gauss2D(A6_fit[t],Hf[tr(ii)],Kf[tr(ii)],H3_fit[t],K3_fit[t], SH6_fit[t],SK6_fit[t],theta6_fit[t]) + B_fit[t])); if (ii%31==0) fprintf (file, " \n");} fclose(file);} pipe = popen("gnuplot", "w"); fprintf (pipe, "set term png \n"); fprintf (pipe, "set view 90,0 \n"); fprintf (pipe, "set xyplane 0.05 \n"); for (t=1; t<=tmax; t++){ fprintf (pipe, "set output '%s%04i%s' \n", datfile,t,"tf.png"); fprintf (pipe, "splot '%s%04i%s' w l, '%s%04i%s' w l \n", datfile,t,"t.gnu", datfile,t,"f.gnu"); fflush(pipe);} fclose (pipe); return;} /********************** function: read_fit_dat_files ***************/ /**************** function: gnuplot_read_dat_files_t ***************/ void gnuplot_read_dat_files_t () { /* constructs 1D graphic files by gnuplot: file parameter: q; plot axis: time */ int ii, mtime, lineindex; char filename [MAX_LEN]; char gfilename [MAX_LEN]; char line [MAX_LEN+1]; float fH, fK, fL; char dummy[MAX_LEN]; FILE *file; FILE *pipe; /* H, K, L values will be read in: */ for (ii=1; ii<=number; ii++){ sprintf(filename, "%s/%s%i%s", datpath, datfile,ii,".dat"); /* here will be chosen: which graphics format will be used: */ /* sprintf(gfilename, "%s%03i%s", datfile,ii,".ps"); */ sprintf(gfilename, "%s%03i%s", datfile,ii,".png"); /* sprintf(gfilename, "%s%03i%s", datfile,ii,".gif"); */ printf("%i ", ii); fflush(stdout); file=fopen(filename,"r"); mtime=0; lineindex=0; if (file==NULL) terminate(" The desired file does not exist. "); while (fgets (line, MAX_LEN, file) != NULL) { lineindex=lineindex+1; if (lineindex==9) sscanf (line, "%s%f%f%f", dummy, &fH, &fK, &fL); /* printf("%i H= %f K=%f L=%f Time=%i mon=%i det=%i \n ", index,H[index],K[index],L[index],Time[index],mon[index],det[index]); */} fclose(file); /**** ez itt meg rossz for (t=1; t<=tmax; t++){ sprintf(gfilename, "%s%04i%s", datfile,t,"q.gnu"); header(2); file = fopen (gfilename, "w"); fprintf (file, "%s", headerchar); for (ii=first; ii<=last; ii++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(ii)], Kf[tr(ii)], COUNTSperMMON [tr(ii)][t]); if (ii%31==0) fprintf (file, " \n");} fclose(file);} ****/ /* gnuplot settings: */ pipe = popen("gnuplot -geometry 900x900+10+10", "w"); /* fprintf (pipe, "set logscale z \n"); */ fprintf (pipe, "set title '<%s H=%f K=%f L=%f>' \n", filename,fH,fK,fL); /* fprintf(pipe, "set terminal postscript color solid lw 2\n"); */ fprintf(pipe, "set terminal png \n"); /* fprintf(pipe, "set terminal gif \n"); */ fflush(pipe); fprintf(pipe, "set output '%s' \n", gfilename); fflush(pipe); /* gnuplot command: */ /* version constant scaling: */ fprintf (pipe, "plot [1:%d][0:%.0f] '%s' using ($0):($2)/'%f' t '' w l lw 2 lt 1 \n", tmax, ceil(maxdm), filename, MMON [ii]); /* version auto scaling: */ /* fprintf (pipe, "plot [1:%d][0:] '%s' using ($0):($2)/'%f' t '' w l lw 2 lt 1 \n", tmax, filename, MMON [ii]); */ fclose (pipe); } printf("\n"); return;} /**************** function: gnuplot_read_dat_files_t ***************/ /**************** function: gnuplot_read_dat_files_q ***************/ void gnuplot_read_dat_files_q () { /* constructs 2D graphic and gnuplot files by gnuplot: file parameter: time; plot axes: q */ int t, ii, lineindex, imon, idet; char filename [MAX_LEN]; char gfilename [MAX_LEN]; char line [MAX_LEN+1]; float fH, fK, fL; char dummy[MAX_LEN]; FILE *file; FILE *pipe; for (ii=1; ii<=number; ii++) { sprintf(filename, "%s/%s%i%s", datpath, datfile,ii,".dat"); file=fopen(filename,"r"); t=0; lineindex=0; if (file==NULL) terminate(" The desired file does not exist. "); while (fgets (line, MAX_LEN, file) != NULL) { lineindex=lineindex+1; if (lineindex==9) { sscanf (line, "%s%f%f%f", dummy, &fH, &fK, &fL); Hf[ii] = fH; Kf[ii] = fK; Lf[ii] = fL; /* printf("%i H= %f K=%f L=%f time=%i mon=%i det=%i \n ", index,H[index],K[index],L[index],time[index],mon[index],det[index]); */} else if (lineindex>12) { sscanf (line, "%i%i", &imon, &idet); t=t+1; COUNTSperMMON [ii][t] = idet/MMON [ii]; } } fclose(file);} for (t=1; t<=tmax; t++){ sprintf(gfilename, "%s%04i%s", datfile,t,"t.gnu"); header(2); file = fopen (gfilename, "w"); fprintf (file, "%s", headerchar); for (ii=1; ii<=number; ii++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(ii)], Kf[tr(ii)], COUNTSperMMON [tr(ii)][t]); if (ii%31==0) fprintf (file, " \n");} fclose(file);} for (t=1; t<=tmax; t++) { printf("%i ", t); fflush(stdout); /* here will be chosen: which graphics format will be used: */ /* sprintf(gfilename, "%s%04i%s", datfile,t,"t.ps"); */ sprintf( filename, "%s%04i%s", datfile,t,"t.png"); sprintf(gfilename, "%s%04i%s", datfile,t,"t.gnu"); /* sprintf(gfilename, "%s%04i%s", datfile,ii,"t.gif"); */ /* gnuplot settings: */ pipe = popen("gnuplot -geometry 900x900+10+10", "w"); /* fprintf (pipe, "set logscale zcb \n"); */ fprintf (pipe, "set title '<%s t=%i * time unit>' \n", filename, t); /* fprintf(pipe, "set terminal postscript color solid lw 2\n"); */ fprintf(pipe, "set terminal png \n"); /* fprintf(pipe, "set terminal gif \n"); */ fflush(pipe); fprintf(pipe, "set output '%s' \n", filename); /*fprintf(pipe, "set view 60,45 \n");*/ fprintf(pipe, "set view 0,0 \n"); fflush(pipe); /* gnuplot command: */ /* version constant scaling: */ fprintf (pipe, "set zrange [0:%f] \n", maxdm); /***correct ???***/ fprintf (pipe, "set cbrange [0:%f] \n", ceil(maxdm)); fprintf (pipe, "splot [][][0:%f] '%s' t '' with pm3d \n", ceil(maxdm), gfilename); /* version fixed scaling: */ /*fprintf (pipe, "set zrange [0:80] \n"); fprintf (pipe, "set cbrange [0:80] \n"); fprintf (pipe, "splot [:][:][:] '%s' t '' with pm3d \n", gfilename);*/ /* version auto scaling: */ /* fprintf (pipe, "splot '%s' t '' with pm3d \n", gfilename); */ fclose (pipe);} printf("\n"); return;} /**************** function: gnuplot_read_dat_files_q ***************/ /****************************** function: tr ***********************/ int tr (int i) {/*reorganising of the data file order: K will be monotonous*/ int j; if (i <= 310) j = i; else if (i <= 651) j = i + 310; else j = i - 341; return i; /* nem rendez at */ /*return j; atrendez */ } /****************************** function: tr ***********************/ /***************************** function: power *********************/ double power(double x, double p) {/* calculates x^p */ return exp(p * log(x));} /***************************** function: power *********************/ /********************** function: read_dat_files ***************/ void read_dat_files () { /* experimental counts will be read in: */ int ii, index=0, mtime=0, lineindex; char filename [MAX_LEN]; char line [MAX_LEN+1]; float fH, fK, fL; /* segedvaltozo */ int imon, idet; /* segedvaltozo egesztipusu */ char dummy[MAX_LEN]; FILE *file; double maxdmt[10000]; double maxH[10000], maxK[10000], maxL[10000]; for (ii=0; ii<10000; ii++) maxdmt[ii]=0; if (number == 0) {ii = 0; do {ii = ii + 1; sprintf(filename, "%s/%s%i%s", datpath, datfile,ii,".dat"); file = fopen(filename, "r"); if (file != NULL) fclose(file);} while (file != NULL); number = ii - 1; if (number == 0) terminate ("no data files are available"); else fprintf (stderr, "The number of data files: %d \n", number); } for (ii=1; ii<=number; ii++){ sprintf(filename, "%s/%s%i%s", datpath, datfile,ii,".dat"); /* printf("%s \n", filename); */ printf("%i ", ii); fflush(stdout); file=fopen(filename,"r"); mtime=0; lineindex=0; MMONITOR = 0; if (file==NULL) terminate(" The desired file does not exist. "); while (fgets (line, MAX_LEN, file) != NULL) { lineindex=lineindex+1; if (lineindex==9) sscanf (line, "%s%f%f%f", dummy, &fH, &fK, &fL); else if (lineindex>12) {sscanf (line, "%i%i", &imon, &idet); mtime=mtime+1; index=index+1; Time[index]=mtime; H[index]=fH; K[index]=fK; L[index]=fL; mon[index]=imon; MMONITOR = MMONITOR + imon; det[index]=idet; if (maxdm<(double)idet/(double)imon) maxdm=(double)idet/(double)imon; /* printf("%i H= %f K=%f L=%f Time=%i mon=%i det=%i \n ", index,H[index],K[index],L[index],Time[index],mon[index],det[index]); */ if (maxdmt[mtime]<(double)idet/(double)imon){ maxdmt[mtime]=(double)idet/(double)imon; maxH[mtime]=fH; maxK[mtime]=fK; maxL[mtime]=fL; } } } MMON[ii] = MMONITOR/(double) mtime; fclose(file); } if (tmax == 0) {tmax = mtime; fprintf (stderr, "\nThe number of time points: %d \n", tmax);} sprintf(filename, "maxtest.asc"); file=fopen(filename,"w"); for (ii=1; ii<=mtime; ii++) fprintf(file,"%.4f \t %.4f \t %.4f \t %.4f \n", maxdmt[ii], maxH[ii], maxK[ii], maxL[ii]); fclose(file); printf ("\nMAX(det/mon)= %f \n",maxdm); return;} /********************** function: read_dat_files ***************/ /********************** function: terminate ***************/ void terminate (char *text){ /* Text bei Fehler */ message("terminate:", text); exit (1);} /********************** function: terminate ***************/ /********************** function: message ***************/ void message (char *text1, char *text2){ /*Text beim Suchen*/ fprintf(stderr, "%s %s \n", text1, text2);} /********************** function: message ***************/ /********************** function: error ***************/ void error (char *ref){ /* reports on runtime errors Usage: error ("label"); #include and #include and int errorcount = 0; (global) are necessary.*/ #ifdef _FENV_H int raised; #endif /* fenv.h*/ char text[MAX_LEN]; strcpy(text, "errno error message ("); strcat(text, ref); strcat(text, ")"); perror(text); #ifdef _FENV_H raised = fetestexcept (FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID); /*(FE_INEXACT | FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);*/ if (raised) {errorcount = errorcount + 1; sprintf(text, "%-3i", errorcount); strcat(text, "fp error message ("); strcat(text, ref); strcat(text, "):"); /*if (raised & FE_INEXACT) strcat(text, " INEXACT ");*/ if (raised & FE_DIVBYZERO) strcat(text, " DIVBYZERO "); if (raised & FE_UNDERFLOW) strcat(text, " UNDERFLOW "); if (raised & FE_OVERFLOW) strcat(text, " OVERFLOW "); if (raised & FE_INVALID) strcat(text, " INVALID "); /* putchar('\a'); fflush(stdout); */ message (text, " "); whistle(); if (ERROR) terminate (" Floating point error detected.");} feclearexcept (FE_ALL_EXCEPT); #endif /* fenv.h*/ errno = 0;} /********************** function: error ***************/ /********************** wait ********************/ void wait (char *text){ char ch[MAX_LEN]; fprintf(stderr, "Further with RETURN: %s", text); fgets(ch, MAX_LEN-1, stdin);} /********************** wait ********************/ /********************** whistle ********************/ void whistle (void){putchar('\a'); fflush(stdout);} /********************** whistle ********************/