Repository navigation
Expand file tree
/
Copy pathpeak_comp_poisson_root.c
More file actions
562 lines (486 loc) · 15.9 KB
/
Copy pathpeak_comp_poisson_root.c
File metadata and controls
562 lines (486 loc) · 15.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
#include "peak_comp_poisson_root.h"
#include "read_config.c"
int main(int argc, char *argv[]){
int i,j,k;
FILE *expData, *simData, *results, *scalingFactors;
if(argc != 2){
printf("\npeak_comp_poisson parameter_file\n");
printf("Compares the .mca spectra designated in the parameter file "
"specified and generates cool statistics.\n\n");
exit(-1);
}
// initialize values
addBackground = 0;
for(i = 0; i < NSPECT; i++){
for(j = 0; j < S32K; j++){
expHist[i][j] = 0.;
for(k = 0; k < NSIMDATA; k++){
simHist[k][i][j] = 0.;
scaledSimHist[k][i][j] = 0.;
}
}
}
for(i = 0; i < NPAR; i++){
for(j = 0; j < NSPECT; j++){
aFinal[i][j] = 0.;
}
}
readConfigFile(argv[1]); // grab data from the config file
// check that the number of spectra being compared is fine
if(numSpectra >= NSPECT){
printf("ERROR: A spectrum number specified in the parameter file is larger "
"than the maximum value of %i. Reduce it or increase NSPECT in "
"peak_comp.h and recompile.\n",
NSPECT);
exit(-1);
}
// read in the .mca files
// check file extension of exp data and copy into float histogram
if((expData = fopen(expDataName, "r")) == NULL){
printf("ERROR: Cannot open the experiment data file %s!\n", expDataName);
exit(-1);
}
const char *dot = strrchr(expDataName, '.'); // get the file extension
if(strcmp(dot + 1, "mca") == 0){
readMCA(expData, expDataName, expHist);
}else if(strcmp(dot + 1, "fmca") == 0){
readFMCA(expData, expDataName, expHist);
}else{
printf("ERROR: Improper type of input file: %s\n", expDataName);
printf(
"Integer array (.mca) and float array (.fmca) files are supported.\n");
exit(-1);
}
fclose(expData);
for(i = 0; i < numSimData; i++){
if((simData = fopen(simDataName[i], "r")) == NULL){
printf("ERROR: Cannot open the simulated data file %s!\n", simDataName[i]);
exit(-1);
}
const char *dots = strrchr(simDataName[i], '.'); // get the file extension
if(strcmp(dots + 1, "mca") == 0){
readMCA(simData, simDataName[i], simHist[i]);
}else if(strcmp(dots + 1, "fmca") == 0){
readFMCA(simData, simDataName[i], simHist[i]);
}else{
printf("ERROR: Improper type of input file: %s\n", simDataName[i]);
printf("Integer array (.mca) and float array (.fmca) files are supported.\n");
exit(-1);
}
fclose(simData);
}
//set limits, taking re-binning into account
for(i = 0; i < numSpectra; i++){
fitStartCh[i] = startCh[i]/rebinFactor;
fitEndCh[i] = endCh[i]/rebinFactor;
}
/*for(i = 0; i < numSpectra; i++)
for(j = 0; j < S32K; j++)
for(k=0;k<numSimData;k++){
printf("sim hist %i sp %i ch %i: %f\n",k,spectrum[i],j,simHist[k][spectrum[i]][j]);
}*/
/*// read into ROOT
for(i = 0; i < numSpectra; i++)
for(j = 0; j < S32K; j++)
sim[i]->Fill(j, simHist[0][i][j]);*/
if(verbosity>0){
printf("Spectra read in...\n");
}
// do the fit
//no need to take re-binning into account here as
//this was already done at data import
double chisq = find_chisqMin();
//prepare spectra for saving to disk and/or plotting
//again, no need to take re-binning into account as
//this was already done at data import
// scale simulated data
for(k=0;k<numSimData;k++){
for(i = 0; i < numSpectra; i++){
for(j = 0; j < S32K; j++){
scaledSimHist[k][spectrum[i]][j] = aFinal[k+3][i] * simHist[k][spectrum[i]][j];
}
}
}
// add background to simulated data
for(i = 0; i < numSpectra; i++){
for(j = 0; j < S32K; j++){
if(addBackground == 1){
bgHist[spectrum[i]][j] = aFinal[0][i] + aFinal[1][i]*(double)j + aFinal[2][i]*(double)j*(double)j;
}else if(addBackground == 2){
bgHist[spectrum[i]][j] = aFinal[0][i] + aFinal[1][i]*(double)j;
}else if(addBackground == 3){
bgHist[spectrum[i]][j] = aFinal[0][i];
}else{
bgHist[spectrum[i]][j] = 0.;
}
}
}
// fit result histogram
for(i = 0; i < numSpectra; i++){
for(j = 0; j < S32K; j++){
resultsHist[spectrum[i]][j] = bgHist[spectrum[i]][j];
for(k=0;k<numSimData;k++){
resultsHist[spectrum[i]][j] += (float)scaledSimHist[k][spectrum[i]][j];
}
}
}
// print output
if(verbosity>0){
printf("Fit chisq: %.15f\n", chisq);
}else{
printf("%.15f\n", chisq);
}
// write results
if(saveResults == 1){
//write spectra
sprintf(str, "fit_poisson.fmca");
if((results = fopen(str, "w")) == NULL){
printf("ERROR: Cannot open the output file %s!\n", str);
exit(-1);
}
for(i = 0; i < numSpectra; i++){
fwrite(resultsHist[spectrum[i]], S32K * sizeof(float), 1, results);
}
fclose(results);
//save scaling factors
sprintf(str, "scalingFactors.dat");
if((scalingFactors = fopen(str, "w")) == NULL){
printf("ERROR: Cannot open the output file %s!\n", str);
exit(-1);
}
for(i = 0; i < numSpectra; i++){
fprintf(scalingFactors, "%d %.9f\n", i + 1, aFinal[0][i]);
}
fclose(scalingFactors);
}
/*printf("Scaling factors:\n");
for(i = 0; i < numSpectra; i++){
printf("%d %.9f\n", i + 1, aFinal[0][i]);
}*/
// plot results
if(plotMode >= 0){
theApp = new TApplication("App", &argc, argv);
plotSpectra();
}
return 0; // great success
}
double lrchisq(const double *par){
// chisq from likelihood ratio test
// for more information see Baker and Cousins pg. 439 and Appendix
double yi = 0.; // model
double ni = 0.; // experiment
double lrchisq = 0.; // likelihood ratio chisq
int i = 0;
int j = 0;
for(i = fitStartCh[spCurrent]; i <= fitEndCh[spCurrent]; i++){
ni = expCurrent[i]; // events in ith bin
// calculate model in the ith bin
yi=0;
for(j=0; j<numSimData; j++){
if((j>=1)&&(useRelIntensities)&&(relIntensityAvailable[j])){
yi += par[3]*par[j+3] * simCurrent[j][i]; //amplitude relative to the 1st ampliutude
}else{
yi += par[j+3] * simCurrent[j][i];
}
}
// add background if neccesary
if(addBackground == 1){
yi += par[0] + par[1]*(double)i + par[2]*(double)i*(double)i;
}else if(addBackground == 2){
yi += par[0] + par[1]*(double)i;
}else if(addBackground == 3){
yi += par[0];
}
//model cannot give less than 0 counts
if(yi<0.0){
return 1E10;
}
// evaluate chisq given input parameters
if((ni > 0.)&&(yi != 0.)){
lrchisq += (yi - ni + ni * log(ni / yi));
}else{
lrchisq += yi; // the log(0) case
}
}
lrchisq *= 2.;
/* printf("Parameters: %f %f %f %f, computed lrchisq: %f\n",par[0],par[1],par[2],par[3],lrchisq);
getc(stdin); */
return lrchisq;
}
double find_chisqMin(){
int i = 0;
int j = 0;
int k = 0;
char str[256];
double intExp,intSim;
double chisq = 0.;
if(verbosity>0)
printf("Fitting data...\n");
for(i = 0; i < numSpectra; i++){
if(verbosity>0){
printf("-Spectrum %i-\n",i+1);
}
// for more information see minimizer class documentation
// https://root.cern.ch/root/html/ROOT__Math__Minimizer.html
char minName[132] = "Minuit";
char algoName[132] = ""; // default (Migard for Minuit)
ROOT::Math::Minimizer *min = ROOT::Math::Factory::CreateMinimizer(minName, algoName);
// set tolerance , etc...
min->SetMaxFunctionCalls(10000000); // for Minuit
min->SetMaxIterations(100000);
min->SetTolerance(0.001);
if(verbosity>1){
min->SetPrintLevel(1);
}else{
min->SetPrintLevel(0); // set to 1 for more info
}
// communication to the likelihood ratio function
// (via global parameters)
for(j = 0; j < S32K; j++){
expCurrent[j] = (double)expHist[spectrum[i]][j];
}
for(k=0;k<numSimData;k++){
for(j = 0; j < S32K; j++){
simCurrent[k][j] = (double)simHist[k][spectrum[i]][j];
}
}
spCurrent = i; //needed for fit
//calculate integrals
intSim=0.;
intExp=0.;
for(j = fitStartCh[i]; j <= fitEndCh[i]; j++){
for(k=0;k<numSimData;k++){
intSim += simCurrent[k][j];
//printf("ch %i sim val: %f\n",j,simCurrent[k][j]);
}
intExp += (double)expHist[spectrum[i]][j];
}
//printf("intExp: %f, intSim: %f\n",intExp,intSim);
// create function wrapper for minmizer
// a IMultiGenFunction type
ROOT::Math::Functor lr(&lrchisq, 3+numSimData); // likelihood ratio chisq
// step size and starting variables
// may need to change for best performance
// under different running conditions
double ratio = intExp/intSim;
double variable[3+NSIMDATA];
double step[3+NSIMDATA];
for(j=0;j<3+numSimData;j++){
variable[j] = ratio/2.;
step[j] = ratio/100.;
//printf("Variable %i: val %f, step %f\n",j,variable[j],step[j]);
}
variable[1]/=100000.; //slope for linear background is usually very small
variable[2]/=100000.;
min->SetFunction(lr);
// Set pars for minimization
for(j=0;j<3+numSimData;j++){
sprintf(str, "a%i", j);
if((j>=4)&&(useRelIntensities)&&(ril[j-3]!=rih[j-3])){
min->SetLimitedVariable(j, str, variable[j], step[j],ril[j-3],rih[j-3]); //set relative intensity
}else if((j>=3)&&(forcePosAmp==1)){
min->SetLimitedVariable(j, str, variable[j], step[j],0.0,5.0*ratio); //j>=2 for amplitudes
}else if((j==1)&&(forceNegSlopeBG==1)){
min->SetLimitedVariable(j, str, -1.0*variable[j], step[j]/100., -1000.0, 0.0); //slope not allowed to be positive
}else{
min->SetVariable(j, str, variable[j], step[j]);
}
//printf("Variable %i initialized to %f, step size %f",j,variable[j], step[j]);
//printf(".\n");
}
// do the minimization
min->Minimize();
// grab parameters from minimum
const double *xs = min->X();
//if(verbosity>0)
// printf("Spectrum %i fit parameters: %f %f %f\n",i,xs[0],xs[1],xs[2]);
// assuming 3 parameters
// save pars
for(j = 0; j < 3+numSimData; j++){
if((j>=3)&&(useRelIntensities)&&(relIntensityAvailable[j-3])){
aFinal[j][i] = xs[j]*xs[3]; //intensity relative to first intensity
}else{
aFinal[j][i] = xs[j];
}
}
//print amplitudes
if(verbosity>0){
for(j = 3; j < 3+numSimData; j++){
printf("Spectrum %i, amplitude %i: %f\n",i+1,j-1,aFinal[j][i]);
}
}
chisq += min->MinValue();
}
return chisq;
}
void plotSpectra(){
int i, j, k;
TH1D *data[NSPECT];
TH1D *results[NSPECT];
TH1D *resultsBGData[NSPECT];
TH1D *resultsSimData[NSIMDATA][NSPECT];
TCanvas *c = new TCanvas("c1", "", 1618, 1000);
gStyle->SetErrorX(0.);
//create expt data TH1
char dataName[132];
for(i = 0; i < NSPECT; i++){
sprintf(dataName, "data_%2d", i);
data[i] = new TH1D(dataName, ";;", S32K/rebinFactor, 0, S32K - 1);
}
for(i = 0; i < NSPECT; i++){
for(j = 0; j < S32K; j++){
//first argument of Fill method is the x value, in
//the x-axis units, so need to rescale this as
//the source histogram may have been rebinned eariler
data[i]->Fill(j*rebinFactor, expHist[i][j]);
}
}
for(i = 0; i < NSPECT; i++){
for(j=0;j<data[i]->GetNbinsX();j++){
data[i]->SetBinError(j,sqrt(data[i]->GetBinContent(j)));
}
}
// initialize and fill results histo
char resultsName[132],resultsBGName[132],resultsSimDataName[132];
for(i = 0; i < NSPECT; i++){
sprintf(resultsName, "results_%2d", i);
results[i] = new TH1D(resultsName, ";;", S32K/rebinFactor, 0, S32K - 1);
if(plotMode>=1){
sprintf(resultsBGName, "resultsBG_%2d", i);
resultsBGData[i] = new TH1D(resultsBGName, ";;", S32K/rebinFactor, 0, S32K - 1);
for(k=0;k<numSimData;k++){
sprintf(resultsSimDataName, "resultsSimData_%2d_%2d", i, k);
resultsSimData[k][i] = new TH1D(resultsSimDataName, ";;", S32K/rebinFactor, 0, S32K - 1);
}
}
}
// be careful with indicies here
for(i = 0; i < numSpectra; i++){
for(j = 0; j < S32K; j++){
results[spectrum[i]]->Fill(j*rebinFactor, resultsHist[spectrum[i]][j]);
if(plotMode>=1){
resultsBGData[spectrum[i]]->Fill(j*rebinFactor, bgHist[spectrum[i]][j]);
for(k=0;k<numSimData;k++)
resultsSimData[k][spectrum[i]]->Fill(j*rebinFactor, scaledSimHist[k][spectrum[i]][j] + bgHist[spectrum[i]][j]);
}
}
}
// display limits
Double_t low[NSPECT], high[NSPECT];
// divide canvas to accomodate all spectra
// assumes number to plot <= 3*2=6
c->Divide(3, (int)(ceil(numSpectra/3.0)), 1E-6, 1E-6, 0);
// formatting the plot area
// be careful with indicies
for(i = 1; i <= numSpectra; i++){
c->GetPad(i)->SetRightMargin(0.01);
c->GetPad(i)->SetTopMargin(0.025);
Int_t ind = i - 1;
low[i] = startCh[ind];
high[i] = endCh[ind];
}
// plot
for(i = 1; i <= numSpectra; i++){
c->cd(i);
Int_t ind = i - 1;
// data in black
data[spectrum[ind]]->SetLineStyle(1);
data[spectrum[ind]]->SetLineWidth(2);
data[spectrum[ind]]->SetLineColor(12);
data[spectrum[ind]]->SetMarkerColor(12);
data[spectrum[ind]]->SetMarkerSize(0.8);
data[spectrum[ind]]->SetMarkerStyle(21);
data[spectrum[ind]]->GetXaxis()->SetRangeUser(low[i], high[i]);
data[spectrum[ind]]->SetStats(0);
data[spectrum[ind]]->Draw("PE1");
// simulation in red
results[spectrum[ind]]->SetLineStyle(1);
results[spectrum[ind]]->SetLineWidth(2);
results[spectrum[ind]]->SetLineColor(46);
results[spectrum[ind]]->Draw("HIST SAME");
if(plotMode>=1){
//plot individual simulated data
for(k=0;k<numSimData;k++){
resultsSimData[k][spectrum[ind]]->SetLineStyle(1);
resultsSimData[k][spectrum[ind]]->SetLineWidth(2);
resultsSimData[k][spectrum[ind]]->SetLineColor(799+k*20);
resultsSimData[k][spectrum[ind]]->Draw("HIST SAME");
}
if(plotMode!=2){
//plot background
resultsBGData[spectrum[ind]]->SetLineStyle(1);
resultsBGData[spectrum[ind]]->SetLineWidth(1);
resultsBGData[spectrum[ind]]->SetLineColor(920);
resultsBGData[spectrum[ind]]->Draw("HIST SAME");
}
}
}
c->cd(1);
// x1,y1,x2,y2
TLegend *leg;
if(plotMode>=1){
leg = new TLegend(0.70, 0.85 - numSimData*0.05, 0.90, 0.95);
}else{
leg = new TLegend(0.70, 0.85, 0.90, 0.95);
}
leg->AddEntry(data[0], "Experiment", "lp");
leg->AddEntry(results[spectrum[0]], "Simulation", "l");
if(plotMode>=1){
for(k=0;k<numSimData;k++){
sprintf(resultsName, "Branch %i", k+1);
leg->AddEntry(resultsSimData[k][spectrum[k+1]], resultsName, "l");
}
}
leg->SetFillColor(0);
leg->SetFillStyle(0);
leg->SetBorderSize(0);
leg->SetTextSize(0.045);
leg->Draw();
c->SetBorderMode(0);
c->Update();
theApp->Run(kTRUE);
}
int readMCA(FILE *inp, char *filename, float inpHist[NSPECT][S32K]){
if(verbosity>0){
printf("Reading %i spectra in MCA file: %s\n",endSpectrum,filename);
}
int mcaHist[NSPECT][S32K];
for(int i = 0; i <= endSpectrum; i++){
if(fread(mcaHist[i], S32K * sizeof(int), 1, inp) != 1){
printf("ERROR: Error reading file %s!\n", filename);
exit(-1);
}
}
for(int i = 0; i < NSPECT; i++){
int bin = -1;
for(int j = 0; j < S32K; j++){
if(j%rebinFactor==0){
bin++;
}
inpHist[i][bin] += (float)mcaHist[i][j];
}
}
return 1;
}
int readFMCA(FILE *inp, char *filename, float inpHist[NSPECT][S32K]){
if(verbosity>0){
printf("Reading %i spectra in FMCA file: %s\n",endSpectrum,filename);
}
float mcaHist[NSPECT][S32K];
for(int i = 0; i <= endSpectrum; i++){
if(fread(mcaHist[i], S32K * sizeof(float), 1, inp) != 1){
printf("ERROR: Error reading file %s!\n", filename);
exit(-1);
}
}
for(int i = 0; i < NSPECT; i++){
int bin = -1;
for(int j = 0; j < S32K; j++){
if(j%rebinFactor==0){
bin++;
}
inpHist[i][bin] += mcaHist[i][j];
}
}
return 1;
}