#include #include #include #include #include #include "util.h" int main(void) { double gain[6], c[50], dc[50], area[50], darea[50]; double ac, dac, eg, deg, aeg, daeg, eff, w; float deff = 0.0f, pars[10]; int meff, i, j[50], j1, j2; int nc, iorder, nlines, nterms; char title[80], nf[80], nam[50][30], ans[80], line[128]; FILE *file1, *file2; #define IFMT "%d%lf%lf%lf%lf%*24c%28c" #define OFMT "%3d %10.4f %8.4f %10.0f %8.0f %10.4f %8.4f %.28s\n" int effeval(float *, double, double *); int calc_energy(double, double, double *, double *, double *, int); /* ENERGY - a program to calculate energies for gf3.sto files and divide by detector efficiencies */ printf( "\n\n" "WELCOME TO ENERGY.\n" "This program uses the information in .sto files to calculate energies\n" " and intensities for fitted peaks. The running of this program is a\n" " necessary step in creating input (.eng) files for the programs legft,\n" " rdm and rdmfit. The .sto file is created using the SA command in gf3.\n" " If you have fitted one spectrum there will be a line in your .sto file\n" " for each fitted peak. If you have fitted a series of n spectra\n" " (e.g. a set of angular distribution spectra) and used the SA command\n" " properly, there will be a set of n lines in your .sto file for each\n" " peak. For more than 20 spectra, you will have to do some editing.\n\n" "The energy calibration will be taken from a .aca or .cal file, made\n" " by running the program encal. The energy will be calculated for each\n" " line in the .sto file. In addition, if there is more than one line per\n" " peak (i.e. more than one fitted spectrum), an average energy will be\n" " calculated. You may also divide the peak areas by the detector\n" " efficiency to obtain the gamma-ray intensity. If you do so, the\n" " efficiency parameters will be taken from a .aef or .eff file,\n" " generated by running the program effit.\n" " C version D. C. Radford Sept 1999\n"); /* ask for data file name */ while (1) { cask("gf2/gf3 .sto file = ? (default .ext = .sto)", nf, 80); setext(nf, ".sto", 80); if (!(file1 = open_readonly(nf))) continue; /* ask for output file name */ cask("Output file name = ? (default .ext = .eng)", nf, 80); setext(nf, ".eng", 80); file2 = open_new_file(nf, 1); /* ask for energy calibration file name and read parameters */ while (!cask("Encal parameter file = ? (default .ext = .aca, .cal)", nf, 80) || read_cal_file(nf, title, &iorder, gain)); nterms = iorder + 1; while (1) { nc = cask("Number of lines per peak = ?", ans, 80); if (inin(ans, nc, &nlines, &j1, &j2)) continue; if (nlines <= 0) nlines = 1; if (nlines <= 50) break; printf("Maximum is 50.\n"); } meff = -1; if (caskyn("Do you want to divide the areas by the efficiency? (Y/N)")) { meff = 1; /* ask for effit output file name and read parameters */ while (!cask("Effit parameter file = ? (default .ext = .aef, .eff)", nf, 80) || read_eff_file(nf, title, pars)); deff = exp(pars[9]) - 1.f; printf(" Percentage error on efficiency = %.4f\n", deff * 100.f); fprintf(file2, " No. Centroid +- error Intensity +- error" " Energy +- error Sp.name\n"); } else { fprintf(file2, " No. Centroid +- error Area +- error" " Energy +- error Sp.name\n"); } if (!fgets(line, 120, file1)) goto NEXTFIL; while (1) { ac = 0.0f; dac = 0.0f; i = 0; while (i < nlines) { if (!fgets(line, 120, file1)) goto NEXTFIL; if (strstr(line, "Centroid") || sscanf(line, IFMT, &j[i], &c[i], &dc[i], &area[i], &darea[i], nam[i]) != 6 || c[i] == 0.0) continue; /* if centroid == 0.0 ignore line */ if (dc[i] <= 0.0) { printf("Errors must be nonzero; zero values set to 1.0 chs.\n"); dc[i] = 1.0; } if (nlines > 1) { w = 1.0f / (dc[i] * dc[i]); ac += w * c[i]; dac += w; } i++; } if (nlines > 1) { if (ac == 0.0f) continue; /* if centroid == 0.0 ignore lines */ ac /= dac; dac = 1.0f / sqrt(dac); eff = 1.0f; calc_energy(ac, dac, &aeg, &daeg, gain, nterms); if (meff > 0) { if (effeval(pars, aeg, &eff)) goto NEXTFIL; for (i = 0; i < nlines; ++i) { area[i] = area[i] * 100.0f / eff; darea[i] = darea[i] * 100.0f / eff; darea[i] = sqrt(darea[i]*darea[i] + area[i]*deff * area[i]*deff); } } for (i = 0; i < nlines; ++i) { calc_energy(c[i], dc[i], &eg, °, gain, nterms); fprintf(file2, OFMT, j[i], c[i], dc[i], area[i], darea[i], eg, deg, nam[i]); } if (meff > 0) { fprintf(file2, " Mean energy = %10.4f +- %8.4f Effic. = %10.4f\n\n", aeg, daeg, eff); } else { fprintf(file2, " Mean energy = %10.4f +- %8.4f\n\n", aeg, daeg); } } else { i = 0; calc_energy(c[i], dc[i], &eg, °, gain, nterms); if (meff >= 0) { if (effeval(pars, eg, &eff)) goto NEXTFIL; area[i] = area[i] * 100.0f / eff; darea[i] = darea[i] * 100.0f / eff; darea[i] = sqrt(darea[i]*darea[i] + area[i]*deff * area[i]*deff); } fprintf(file2, OFMT, j[i], c[i], dc[i], area[i], darea[i], eg, deg, nam[i]); } } NEXTFIL: fclose(file1); fclose(file2); if (!caskyn("Would you like to process more files? (Y/N)")) break; } return 0; } /* energy */ /* ====================================================================== */ int effeval(float *pars, double x, double *fit) { double g, f1, f2, u3, x1, x2, u4; /* calculate the efficiency */ x1 = log(x / pars[7]); x2 = log(x / pars[8]); g = pars[6]; f1 = pars[0] + pars[1] * x1 + pars[2] * x1 * x1; f2 = pars[3] + pars[4] * x2 + pars[5] * x2 * x2; if (f1 <= 0.f || f2 <= 0.f) { printf("f1 or f2 <= 0.0 - Cannot calculate eff.\n"); return 1; } u3 = exp(-g * log(f1)); u4 = exp(-g * log(f2)); if (u3 + u4 <= 0.f) { printf("u3 + u4 <= 0.0 - Cannot calculate eff.\n"); return 1; } *fit = exp(exp(-log(u3 + u4) / g)); return 0; } /* effeval */ /* ====================================================================== */ int calc_energy(double c, double dc, double *eg, double *deg, double *gain, int nterms) { int i; *eg = gain[nterms-1]; *deg = 0.0f; for (i = nterms-1; i > 0; --i) { *eg = gain[i-1] + *eg * c; *deg = (float) i * gain[i] + *deg * c; } *deg *= dc; return 0; } /* calc_energy */