GammaWare  Head Version for release 0.9
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FitSpek.cpp
Go to the documentation of this file.
1 // FitSpek.cpp: implementation of the CFitSpek class.
2 //
4 
5 #include <memory.h>
6 #include <stdlib.h>
7 #include <cmath>
8 #include <algorithm>
9 
10 #include "FitSpek.h"
11 
12 //#ifndef max
13 //#define max(a,b) (((a) > (b)) ? (a) : (b))
14 //#endif
15 
16 //#ifndef min
17 //#define min(a,b) (((a) < (b)) ? (a) : (b))
18 //#endif
19 
20 //#ifdef _DEBUG
21 //#undef THIS_FILE
22 //static char THIS_FILE[]=__FILE__;
23 //#define new DEBUG_NEW
24 //#endif
25 
27 // Construction/Destruction //
29 
30 CFitSpek::CFitSpek()
31 {
32  m_bValid = m_bFunct = m_bBckgr = false;
33  m_bLine = m_bFlat = m_bSlope = false;
34  m_bGauss = m_bExp = m_bSin = false;
35  m_bStep = m_bTailL = m_bTailR = false;
36  m_pRes = NULL;
37  m_pErr = NULL;
38  m_pPar = NULL;
39  m_pDatVal = NULL;
40  m_pDatErr = NULL;
41 
42  // the default settings
43  // [Bit1 enabled] [Bit2 all equal]
44  m_pGFitPar[0] = 1; // 0 Bg0
45  m_pGFitPar[1] = 0; // 1 Bg1 not enabled
46  m_pGFitPar[2] = 1; // 2+0 Amplitude
47  m_pGFitPar[3] = 1; // 2+1 Position
48  m_pGFitPar[4] = 3; // 2+2 FWMH all equal
49  m_pGFitPar[5] = 3; // 2+3 Step all equal
50  m_pGFitPar[6] = 2; // 2+4 TailL all equal but not enabled
51  m_pGFitPar[7] = 2; // 2+5 TailR all equal but not enabled
52 }
53 
54 CFitSpek::~CFitSpek()
55 {
56  Reset();
57 }
58 
59 void CFitSpek::Reset()
60 {
61  m_bValid = m_bFunct = m_bBckgr = false;
62  m_bLine = m_bFlat = m_bSlope = false;
63  m_bGauss = m_bExp = m_bSin = false;
64  m_bStep = m_bTailL = m_bTailR = false;
65  if(m_pRes) {delete [] m_pRes; m_pRes = NULL;}
66  if(m_pErr) {delete [] m_pErr; m_pErr = NULL;}
67  if(m_pPar) {delete [] m_pPar; m_pPar = NULL;}
68  if(m_pDatVal) {delete [] m_pDatVal; m_pDatVal = NULL;}
69  if(m_pDatErr) {delete [] m_pDatErr; m_pDatErr = NULL;}
70 }
71 
72 void CFitSpek::GFitParDlg()
73 {
74 #if 0
75  CFitSpekDlg setUp;
76 
77  setUp.m_bB0 = m_pGFitPar[0]&1;
78  setUp.m_bB1 = m_pGFitPar[1]&1;
79  setUp.m_bAmpl = m_pGFitPar[2]&1;
80  setUp.m_bPosi = m_pGFitPar[3]&1;
81  setUp.m_bFWHM = m_pGFitPar[4]&1;
82  setUp.m_bStep = m_pGFitPar[5]&1;
83  setUp.m_bTailL = m_pGFitPar[6]&1;
84  setUp.m_bTailR = m_pGFitPar[7]&1;
85  setUp.m_bFWHMeq = m_pGFitPar[4]&2;
86  setUp.m_bTailLeq = m_pGFitPar[6]&2;
87  setUp.m_bTailReq = m_pGFitPar[7]&2;
88 
89  if(setUp.DoModal() != IDOK)
90  return;
91 
92  m_pGFitPar[0] = setUp.m_bB0 ? 1 : 0;
93  m_pGFitPar[1] = setUp.m_bB1 ? 1 : 0;
94  m_pGFitPar[2] = setUp.m_bAmpl ? 1 : 0;
95  m_pGFitPar[3] = setUp.m_bPosi ? 1 : 0;
96  m_pGFitPar[4] = setUp.m_bFWHM ? 1 : 0;
97  m_pGFitPar[5] = setUp.m_bStep ? 1 : 0;
98  m_pGFitPar[6] = setUp.m_bTailL ? 1 : 0;
99  m_pGFitPar[7] = setUp.m_bTailR ? 1 : 0;
100  m_pGFitPar[4] |= setUp.m_bFWHMeq ? 2 : 0;
101  m_pGFitPar[5] |= 2 ;
102  m_pGFitPar[6] |= setUp.m_bTailLeq ? 2 : 0;
103  m_pGFitPar[7] |= setUp.m_bTailReq ? 2 : 0;
104 #endif
105 
106 }
107 
109 // Fit of a stright line //
111 
112 bool CFitSpek::CalcStrightLine(float *pSpek, int chanA, int chanB, int doSlope)
113 {
114  // nothing valid yet
115  m_bValid = false;
116 
117  float *pd;
118 
119  if(chanA < chanB) {
120  m_nChanA = chanA;
121  m_nChanB = chanB;
122  }
123  else {
124  m_nChanA = chanB;
125  m_nChanB = chanA;
126  }
127  m_nBgChanA = m_nChanA;
128  m_nBgChanB = m_nChanB;
129  m_nNdat = m_nChanB - m_nChanA + 1;
130 
131  // tilted line (needs at least 2 channels)
132  if(doSlope&1) {
133  if(m_nNdat < 2) return false;
134 
135  double s1 = 0., sx = 0., sxx = 0., sy = 0., syx = 0.;
136  double x, y, e;
137  pd = pSpek + m_nChanA;
138  if(doSlope&2) {
139  // use only the two extreme channels
140  for(int ii = 0; ii < m_nNdat; ii+=m_nNdat-1) {
141  y = pd[ii];
142  x = ii+0.5; // considera il centro del canale
143  e = 1.; // andrebbe sostituito con l'errore statistico
144  s1 += e;
145  sx += e*x;
146  sxx += e*x*x;
147  sy += e*y;
148  syx += e*y*x;
149  }
150  }
151  else {
152  // use all channels in region
153  x = 0.5; // considera il centro del canale
154  for(int ii = 0; ii < m_nNdat; ii++, x += 1.) {
155  y = *pd++;
156  e = 1.; // andrebbe sostituito con l'errore statistico
157  s1 += e;
158  sx += e*x;
159  sxx += e*x*x;
160  sy += e*y;
161  syx += e*y*x;
162  }
163  }
164  double deter = s1*sxx-sx*sx;
165  if(deter >= 0.) {
166  m_dBack0 = (sxx*sy - sx*syx)/deter;
167  m_dBack1 = (-sx*sy + s1*syx)/deter;
168  m_dBack0Err = sqrt(sxx/deter);
169  m_dBack1Err = sqrt(s1 /deter);
170 // double eback01 = -sx /deter;
171  m_bFlat = true;
172  m_bSlope = true;
173  double sy1 = 0.;
174  double sy2 = 0.;
175  x = 0.5; // considera il centro del canale
176  pd = pSpek + m_nChanA;
177  for(int ii = 0; ii < m_nNdat; ii++, x += 1.) {
178  double dy = *pd++ - (m_dBack0 + m_dBack1*x);
179  sy1 += fabs(dy);
180  sy2 += dy*dy;
181  }
182  m_dBackL1 = sy1 / m_nNdat;
183  m_dBackL2 = sqrt(sy2 / m_nNdat);
184  }
185  }
186 
187  // Flat line
188  else {
189  double sy0 = 0.;
190  pd = pSpek + m_nChanA;
191  if(doSlope&2) {
192  // use only the two extreme channels
193  for(int ii = 0; ii < m_nNdat; ii+=m_nNdat-1)
194  sy0 += pd[ii];
195  m_dBack0 = sy0 / 2;
196  }
197  else {
198  // use all channels in region
199  for(int ii = 0; ii < m_nNdat; ii++)
200  sy0 += *pd++;
201  m_dBack0 = sy0 / m_nNdat;
202  }
203  double sy1 = 0.;
204  double sy2 = 0.;
205  pd = pSpek + m_nChanA;
206  for(int ii = 0; ii < m_nNdat; ii++) {
207  double dy = *pd++ - m_dBack0;
208  sy1 += fabs(dy);
209  sy2 += dy*dy;
210  }
211  m_dBackL1 = sy1 / m_nNdat;
212  m_dBackL2 = sqrt(sy2 / m_nNdat);
213  m_dBack1 = m_dBack1Err = 0.;
214  m_bFlat = true;
215  m_bSlope = false;
216  }
217 
218  m_bValid = true;
219  m_bLine = true;
220  m_bBckgr = true;
221  m_bFunct = false;
222 
223  // puntatori alle funzioni risultato del calcolo della retta
224  m_pfFitFunc = &CFitSpek::LineFunc;
225  m_pfFitBack = &CFitSpek::LineFunc;
226 
227  return m_bValid;
228 }
229 
230 // metodo pubblico (attraverso membro puntatore)
231 // che lavora sui risultati del fit
232 double CFitSpek::LineFunc(double xx)
233 {
234  if(!m_bLine) return 0.;
235 
236  double val = 0;
237  if(m_bFlat ) val += m_dBack0;
238  if(m_bSlope) val += m_dBack1*(xx-m_nChanA);
239 
240  return val;
241 }
242 
244 // Calcola i primi momenti (Area, Posizione, Sigma) e il loro errore //
245 // sottrae un fondo lineare calcolato dai primi e ultimi nChanBack canali //
247 
248 bool CFitSpek::CalcPeakMoments(float *pSpek, int chanA, int chanB, int nChanBack)
249 {
250  // nothing valid yet
251  Reset();
252 
253  if(chanA < chanB) {
254  m_nChanA = chanA;
255  m_nChanB = chanB;
256  }
257  else {
258  m_nChanA = chanB;
259  m_nChanB = chanA;
260  }
261  m_nBgChanA = m_nChanA;
262  m_nBgChanB = m_nChanB;
263  m_nNdat = m_nChanB - m_nChanA + 1;
264 
265  int ndat = m_nNdat - 2 * nChanBack;
266  if(ndat < 3)
267  return false;
268 
269 // i primi e gli ultimi nChanBack canali vengono usati per calcolare un fondo inclinato
270 
271  m_dBack0 = 0.;
272  m_dBack1 = 0.;
273  if(nChanBack > 0) {
274  double s1 = 0., sx = 0., sxx = 0., sy = 0., syx = 0.;
275  double x, y, e;
276  for(int ii = m_nChanA; ii <= m_nChanB; ii++) {
277  if(ii < m_nChanA+nChanBack || ii > m_nChanB-nChanBack) {
278  y = pSpek[ii];
279  x = ii - m_nChanA + 0.5;
280  e = 1.; // andrebbe sostituito con l'errore statistico ???
281  s1 += e;
282  sx += e*x;
283  sxx += e*x*x;
284  sy += e*y;
285  syx += e*y*x;
286  }
287  }
288  double deter = s1*sxx-sx*sx;
289  if(deter >= 0.) {
290  m_dBack0 = (sxx*sy - sx*syx)/deter;
291  m_dBack1 = (-sx*sy + s1*syx)/deter;
292  m_dBack0Err = sxx/deter;
293  m_dBack1Err = s1 /deter;
294 // double er01 = -sx /deter;
295  }
296  }
297  m_bFlat = m_dBack0 != 0;
298  m_bSlope = m_dBack1 != 0;
299  m_bBckgr = m_bFlat || m_bSlope;
300  m_bLine = true;
301 
302  double dMV0, dMV1;
303  double dKV2, dKV3, dKV4;
304 
305  // Area e Posizione calcolati sottraendo il fondo (senza propagazione dell'errore)
306  // Posizioni riferite all'inizio di m_nChanA, con riferimento al centro del canale
307  // tenere la zona di lavoro sufficientemente stretta per evitare che le deviazioni
308  // dal fondo in zone lontano dal picco pesino troppo
309  // (sarebbe meglio fare due passaggi escludendo i canali piu' lontani di ~3 FWHM)
310 
311  dMV0 = dMV1 = 0.;
312  double xii, xyy, back;
313  for(int ii = m_nChanA+nChanBack; ii <= m_nChanB-nChanBack; ii++) {
314  xii = ii + 0.5 - m_nChanA;
315  back = m_dBack0 + xii*m_dBack1;
316  xyy = pSpek[ii];
317  xyy -= back; // sottrae il fondo
318  dMV0 += xyy; // M0 ===> Area
319  dMV1 += xyy * xii; // M1 ===> Posizione
320  }
321  if(dMV0 == 0.)
322  return false;
323  dMV1 /= dMV0;
324 
325  // i prossimi 3 momenti (calcolati rispetto al baricentro dMV1)
326 
327  dKV2 = dKV3 = dKV4 = 0;
328  double xny;
329  for(int ii = m_nChanA+nChanBack; ii <= m_nChanB-nChanBack; ii++) {
330  xii = ii + 0.5 - m_nChanA;
331  back = m_dBack0 + xii*m_dBack1;
332  xyy = pSpek[ii];
333  xyy -= back; // sottrae il fondo
334  xyy = max(0., xyy); // lo sviluppo in momenti vale solo per dati positivi!!!
335  xii -= dMV1; // momenti rispetto al baricentro
336  dKV2 += (xny = xyy*xii*xii); // K2 ==> K2 = Varianza = sigma^2
337  dKV3 += (xny *= xii); // K3 ==> Skewness = K3 / sigma^3
338  dKV4 += (xny *= xii); // K4 ==> Kurtosis = K4 / sigma^4 - 3
339  }
340  dKV2 = dKV2 / dMV0;
341  dKV3 = dKV3 / dMV0;
342  dKV4 = dKV4 / dMV0;
343 
344  m_nNpeaks = 1; // risultati riportati come fit gaussiano a 1 picco
345  m_pRes = new CFitRes [m_nNpeaks];
346  m_pErr = new CFitRes [m_nNpeaks];
347 
348  m_pRes[0].Area = dMV0;
349  m_pRes[0].Posi = dMV1 + m_nChanA;
350  m_pRes[0].Sigma = m_pRes[0].Area ? sqrt(dKV2) : 0.;
351  m_pRes[0].Fwhm = m_pRes[0].Sigma * s2fwhm;
352  m_pRes[0].Fw05 = m_pRes[0].Sigma * s2fwhm;
353  m_pRes[0].Fw01 = m_pRes[0].Sigma * s2fwtm;
354  m_pRes[0].W01L = sqrt(log(10.)/log(2.));
355  m_pRes[0].W01R = sqrt(log(10.)/log(2.));
356  m_pRes[0].Ampli = m_pRes[0].Sigma ? m_pRes[0].Area/(m_pRes[0].Sigma*sqrt2pi) : 0.;
357  m_bStep = false;
358  m_bTailL = false;
359  m_bTailR = false;
360 
361 
363  // Err(Kn) = sqrt( (K2n - Kn*Kn)/Area ) //
365 
366  m_pErr[0].Area = sqrt(dMV0); // ??
367  m_pErr[0].Posi = sqrt((dKV2 )/ dMV0); // K1 = 0 per definizione;
368  m_pErr[0].Sigma = sqrt((dKV4 - dKV2*dKV2)/ dMV0)/(2*sqrt(dKV2));
369  m_pErr[0].Fwhm = m_pErr[0].Sigma * s2fwhm;
370 
371  m_bValid = true;
372  m_bFunct = true;
373  m_bGauss = true;
374 
375  // puntatori alle funzioni risultato del calcolo dei momenti
376  // coincidono con quella del fit gaussiano
377  m_pfFitFunc = &CFitSpek::GFitFunc; // background + all peaks
378  m_pfFitPeak = &CFitSpek::GFitPeak; // one particular peak
379  m_pfFitBack = &CFitSpek::LineFunc; // background only
380 
381  // ora chalcola un chi2
382  m_dChi2 = 0.;
383  for(int ii = m_nChanA+nChanBack; ii <= m_nChanB-nChanBack; ii++) {
384  double dy = pSpek[ii] - FitFunc(ii+0.5);
385  dy *= dy;
386  dy /= pSpek[ii] ? pSpek[ii] : 1;
387  m_dChi2 += dy;
388  }
389  m_dChi2 /= m_nChanB-nChanBack - (m_nChanA+nChanBack);
390  return m_bValid;
391 }
392 
394 // Fit of a gaussian peak with tails over a stright+step background //
396 
397 /*
398  X dove si vuole calcolare (riferito all'inizio dello spettro)
399  chanA Z canale di riferimento
400  back0 B0 fondo costante
401  back1 B1 slope del fondo riferito a chanA
402  ampli A ampiezza del picco
403  posi P posizione del picco (riferita all'inizio dello spettro)
404  sigma W varianza del picco
405  step S ampiezza step (frazione dell'unita',positiva) ToDo
406  tailL L attacco coda sinistra (in unita' di sigma, negativo) ToDo
407  tailR R attacco coda destra (in unita' di sigma, positivo) ToDo
408 
409  yy = (X-P)/W
410 
411  Funct(yy) = B + SUMpeaks(A * ( F1 + F2 ))
412 
413  B = B0 + B1 * (X-Z)
414 
415  = exp(-0.5*L*(2*yy-L)) per yy < L
416  F1 = exp(-0.5*yy^2) per L <= yy <= R
417  = exp(-0.5*R*(2*yy-R)) per R < yy
418 
419  F2 = S/(1+exp(yy))^2
420 
421 */
422 
423 // Questo obbrobrio e' una traduzione di un programma Fortran (recal_cob ...)
424 // Il metodo e' complicato dalla possibilita' di imporre limiti ai parametri
425 // I parametri sono controllati da m_pGFitPar che va caricata dall'esterno con GFitParSet()
426 int CFitSpek::CalcGaussianFit(float *pSpek, int chanA, int chanB, std::vector<double>&vPeaks)
427 {
428  // nothing valid yet
429  Reset();
430 
431  // puntatori alle funzioni per Curfit()
432  m_pfFunct = &CFitSpek::GFfunF; // funzione
433  m_pfDeriv = &CFitSpek::GFderF; // derivata
434 
435  m_nOffB = 2;
436  m_nOffN = 6;
437 
438  int ii;
439  int ierr=0;
440  double chi2, chi2n;
441 
442  // Verify fit region
443  if(chanA < chanB) {
444  m_nChanA = chanA;
445  m_nChanB = chanB;
446  }
447  else {
448  m_nChanA = chanB;
449  m_nChanB = chanA;
450  }
451  m_nBgChanA = m_nChanA;
452  m_nBgChanB = m_nChanB;
453 
454  // nothing to fit if region empty
455  bool empty = true;
456  for(ii = m_nChanA; ii <= m_nChanB; ii++) {
457  if(pSpek[ii]) {
458  empty = false;
459  break;
460  }
461  }
462  if(empty) {
463  ierr = -1;
464  return ierr;
465  }
466 
467  // keep only peaks internal to fit region
468  double *lPeaks = new double[std::max(1,(int)vPeaks.size())];
469  int nPeaks = 0;
470  for(size_t np = 0; np < vPeaks.size(); np++) {
471  if(vPeaks[np] > m_nChanA && vPeaks[np] < m_nChanB) {
472  lPeaks[nPeaks++] = vPeaks[np];
473  }
474  }
475 
476  // if no peak present, set one at max amplitude channel
477  if(nPeaks == 0) {
478  int pm = (m_nChanA+m_nChanB+1)/2;
479  float vm = pSpek[pm];
480  for(ii = m_nChanA+1; ii < m_nChanB-1; ii++) {
481  if(pSpek[ii] > vm) {
482  pm = ii;
483  vm = pSpek[pm];
484  }
485  }
486  lPeaks[0] = pm + 0.5;
487  nPeaks = 1;
488  }
489 
490  m_nNdat = m_nChanB - m_nChanA + 1; // margins included
491  m_nNpeaks = nPeaks;
492  m_nNpar = m_nOffB + m_nOffN * m_nNpeaks;
493 
494  // create parameters structure
495  m_pPar = new CFitPar [m_nNpar];
496 
497  // the data to fit
498  m_pDatVal = new double [m_nNdat];
499  m_pDatErr = new double [m_nNdat];
500 
501  // prepare vector of data and 1/(sigma^2)
502  float *ps = pSpek + m_nChanA;
503  for(ii = 0; ii < m_nNdat; ii++) {
504  double yy = *ps++;
505  m_pDatVal[ii] = yy;
506  m_pDatErr[ii] = (yy) ? (1./fabs(yy)) : (1.);
507  }
508 
510  // estimates for the first fit round //
511  // limiting number of free parameters //
513 
514  double bgL = (m_pDatVal[0 ] + m_pDatVal[1 ])/2.;
515  double bgR = (m_pDatVal[m_nNdat-1] + m_pDatVal[m_nNdat-2])/2.;
516 
517  // background level [0]
518  // background slope [1]
519  int bflag = 0;
520  if(m_pGFitPar[0]&1) bflag |= 1; // BG0 enabled ?
521  if(m_pGFitPar[1]&1) bflag |= 2; // BG1 enabled ?
522  switch (bflag) {
523  case 0:
524  m_bFlat = false;
525  m_bSlope = false;
526  m_pPar[0].Value = 0.;
527  m_pPar[1].Value = 0.;
528  break;
529  case 1:
530  m_bFlat = true;
531  m_bSlope = false;
532  m_pPar[0].Value = (bgR+bgL)/2;
533  m_pPar[1].Value = 0.;
534  break;
535  case 2:
536  m_bFlat = false;
537  m_bSlope = true;
538  m_pPar[0].Value = 0.;
539  m_pPar[1].Value = (bgR-bgL)/m_nNdat;
540  break;
541  case 3:
542  m_bFlat = true;
543  m_bSlope = true;
544  m_pPar[0].Value = bgL;
545  m_pPar[1].Value = (bgR-bgL)/m_nNdat;
546  break;
547  }
548 
549  m_pPar[0].Status = (m_bFlat ) ? p_Free : p_Fixed;
550  m_pPar[1].Status = (m_bSlope) ? p_Free : p_Fixed;
551 
552  m_pPar[0].Bounds = p_BNone;
553  m_pPar[1].Bounds = p_BNone;
554 
555  bgL = m_pPar[0].Value;
556  bgR = m_pPar[1].Value;
557 
558  // offs+0 amplitude
559  // offs+1 positions
560  // offs+2 sigma
561  // offs+3 step
562  // offs+4 tailL
563  // offs+5 tailR
564 
565  // set-up the starting values
566  for(int np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
567  double ppos = lPeaks[np] - m_nChanA;
568  int ipos = int(ppos);
569  double ampl = m_pDatVal[ipos];
570  double step = 0.;
571  double back = bgL + ppos*bgR;
572  double ampl2 = (ampl+back)/2.;
573  int width = 1;
574  for(int iL = ipos-1; iL > 0; iL--) { // estimate of peak width
575  if(m_pDatVal[iL] < ampl2) break;
576  width++;
577  }
578  for(int iR = ipos+1; iR < m_nNdat-1; iR++) {
579  if(m_pDatVal[iR] < ampl2) break;
580  width++;
581  }
582  m_pPar[offs+0].Value = ampl - back;
583  m_pPar[offs+1].Value = ppos;
584  m_pPar[offs+2].Value = width/s2fwhm;
585  m_pPar[offs+3].Value = step;
586  m_pPar[offs+4].Value = -100; // non existing
587  m_pPar[offs+5].Value = 100; // non existing
588  }
589 
590  // replace sigmas with weighted average
591  double sigma = 0;
592  double total = 0;
593  for(int np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
594  sigma += fabs(m_pPar[offs+0].Value)*m_pPar[offs+2].Value;
595  total += fabs(m_pPar[offs+0].Value);
596  }
597  sigma /= total;
598  sigma = max(sigma, 1./sqrt(12.)); // 1 channel uniformly distributed
599  for(int np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
600  m_pPar[offs+2].Value = sigma;
601  }
602 
603  // now set the controls for the first fit
604  for(int np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
605  // amplitude --> free but positive
606  m_pPar[offs+0].Status = p_Free;
607  m_pPar[offs+0].LinkTo = offs+0;
608  m_pPar[offs+0].Bounds = p_BLow;
609  m_pPar[offs+0].Valmin = 0.;
610 
611  // position --> free, but inside fit region
612  m_pPar[offs+1].Status = p_Free;
613  m_pPar[offs+1].LinkTo = offs+1;
614  m_pPar[offs+1].Bounds = p_BBoth;
615  m_pPar[offs+1].Valmin = 0.;
616  m_pPar[offs+1].Valmax = m_nNdat-1.;
617 
618  // sigma --> all equal to the first
619  if(np == 0) {
620  m_pPar[offs+2].Status = p_Free;
621  m_pPar[offs+2].LinkTo = offs+2;
622  m_pPar[offs+2].Bounds = p_BBoth;
623  m_pPar[offs+2].Valmin = 1./sqrt(12.); // 1 channel uniformly distributed
624  m_pPar[offs+2].Valmax = m_nNdat/(s2fwhm*2.);
625  }
626  else {
627  m_pPar[offs+2].Status = p_Linked;
628  m_pPar[offs+2].LinkTo = m_nOffB+2;
629  }
630 
631  // step --> fixed
632  m_bStep = false;
633  m_pPar[offs+3].Status = p_Fixed;
634 
635  // tail left --> fixed
636  m_bTailL = false;
637  m_pPar[offs+4].Status = p_Fixed;
638 
639  // tail right --> fixed
640  m_bTailR = false;
641  m_pPar[offs+5].Status = p_Fixed;
642  }
643 
644  // count number of free parameters
645  m_nFpar = 0;
646  for(int np = 0; np < m_nNpar; np++) {
647  if(m_pPar[np].Status == p_Free) m_nFpar++;
648  }
649 
650  if(m_nFpar < 1) {
651  ierr = -3; // nothing to fit
652  goto cleanup;
653  }
654 
655  if(m_nNdat - m_nFpar < 1) {
656  ierr = -4; // insufficient degrees of freedom
657  goto cleanup;
658  }
659 
663 
664  chi2 = Curfit(m_nNdat, m_nNpar, m_nFpar);
665  if(chi2 < 0) {
666  ierr = -5;
667  goto cleanup;
668  }
669 
673 
674  for(int np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
675 
676  // amplitude --> free but positive
677  m_pPar[offs+0].Status = p_Free;
678  m_pPar[offs+0].LinkTo = offs+0;
679  m_pPar[offs+0].Bounds = p_BLow;
680  m_pPar[offs+0].Valmin = 0.;
681 
682  // position --> free, but inside fit region
683  m_pPar[offs+1].Status = p_Free;
684  m_pPar[offs+1].LinkTo = offs+1;
685  m_pPar[offs+1].Bounds = p_BBoth;
686  m_pPar[offs+1].Valmin = 1.;
687  m_pPar[offs+1].Valmax = m_nNdat-2.;
688 
689  // sigma --> free; possibly all equal to first
690  if(np == 0 || !(m_pGFitPar[4]&2)) {
691  m_pPar[offs+2].Value = fabs(m_pPar[offs+2].Value);
692  m_pPar[offs+2].Status = p_Free;
693  m_pPar[offs+2].LinkTo = offs+2;
694  m_pPar[offs+2].Bounds = p_BBoth;
695  m_pPar[offs+2].Valmin = 1./sqrt(12.); // 1 channel uniformly distributed
696  m_pPar[offs+2].Valmax = m_nNdat/(s2fwhm*2.);
697  }
698  else {
699  m_pPar[offs+2].Value = fabs(m_pPar[offs+2].Value);
700  m_pPar[offs+2].Status = p_Linked;
701  m_pPar[offs+2].LinkTo = m_nOffB+2;
702  }
703 
704  // step --> if free can be all equal to that of first peak
705  if(m_pGFitPar[5]&1) {
706  m_bStep = true;
707  if(np == 0 || !(m_pGFitPar[5]&2)) {
708  m_pPar[offs+3].Status = p_Free;
709  m_pPar[offs+3].LinkTo = offs+3;
710  m_pPar[offs+3].Bounds = p_BBoth;
711  m_pPar[offs+3].Valmin = 0;
712  m_pPar[offs+3].Valmax = 1;
713  }
714  else {
715  m_pPar[offs+3].Status = p_Linked;
716  m_pPar[offs+3].LinkTo = m_nOffB+3;
717  }
718  }
719 
720  // tail left --> if free can be all equal to that of first peak
721  if(m_pGFitPar[6]&1) {
722  m_bTailL = true;
723  m_pPar[offs+4].Value = -2.; // attacco iniziale a -2 sigma (~FWHM)
724  if(np == 0 || !(m_pGFitPar[6]&2)) {
725  m_pPar[offs+4].Status = p_Free;
726  m_pPar[offs+4].LinkTo = offs+4;
727  m_pPar[offs+4].Bounds = p_BHig|p_BLow;
728  m_pPar[offs+4].Valmax = -0.1; // upper limit of left tail
729  m_pPar[offs+4].Valmin = -50.; // lower limit of left tail
730  }
731  else {
732  m_pPar[offs+4].Status = p_Linked;
733  m_pPar[offs+4].LinkTo = m_nOffB+4;
734  }
735  }
736 
737  // tail right --> if free can be all equal to that of first peak
738  if(m_pGFitPar[7]&1) {
739  m_bTailR = true;
740  m_pPar[offs+5].Value = 2.; // attacco iniziale a +2 sigma (~FWHM)
741  if(np == 0 || !(m_pGFitPar[7]&2)) {
742  m_pPar[offs+5].Status = p_Free;
743  m_pPar[offs+5].LinkTo = offs+5;
744  m_pPar[offs+5].Bounds = p_BLow|p_BHig;
745  m_pPar[offs+5].Valmin = 0.1; // lower limit of right tail
746  m_pPar[offs+5].Valmax = 50.; // upper limit of right tail
747  }
748  else {
749  m_pPar[offs+5].Status = p_Linked;
750  m_pPar[offs+5].LinkTo = m_nOffB+5;
751  }
752  }
753  }
754 
755  // count number of free parameters
756  m_nFpar = 0;
757  for(int np = 0; np < m_nNpar; np++) {
758  if(m_pPar[np].Status == p_Free) m_nFpar++;
759  }
760 
761  if(m_nFpar < 1) {
762  ierr = -6; // nothing to fit
763  goto cleanup;
764  }
765 
766  if(m_nNdat - m_nFpar < 1) {
767  ierr = -7; // insufficient degrees of freedom
768  goto cleanup;
769  }
770 
771  chi2n = Curfit(m_nNdat, m_nNpar, m_nFpar);
772  if(chi2n < 0) {
773  if(chi2n != -2) { // -2 vuol dire che non e' riuscito a fare meglio
774  ierr = -8;
775  goto cleanup;
776  }
777  chi2n = chi2;
778  }
779 
780  m_dChi2 = chi2n/(m_nNdat-m_nFpar);
781 
782  m_bValid = true;
783  m_bGauss = true;
784  m_bLine = m_bFlat || m_bSlope;
785  m_bFunct = true;
786  m_bBckgr = m_bLine;
787 
788  // save result in data_members
789  GF_res();
790 
791  // puntatori alle funzioni risultato del fit
792  m_pfFitFunc = &CFitSpek::GFitFunc; // background + all peaks
793  m_pfFitBack = &CFitSpek::GFitBack; // background only
794  m_pfFitPeak = &CFitSpek::GFitPeak; // one particular peak
795 
796 cleanup:
797  delete [] m_pPar; m_pPar = NULL;
798  delete [] m_pDatVal; m_pDatVal = NULL;
799  delete [] m_pDatErr; m_pDatErr = NULL;
800  delete [] lPeaks;
801 
802  return (m_bValid) ? (m_nNpeaks) : (ierr);
803 }
804 
805 int CFitSpek::CalcGaussianFit(float *pSpek, int chanA, int chanB)
806 {
807  std::vector<double> vPeaks;
808  if(chanA > chanB)
809  std::swap(chanA, chanB);
810  vPeaks.push_back(chanA-1); // outside fit region so that CalcGaussianFit will assumeand auto locate 1 peak
811  return CalcGaussianFit(pSpek, chanA, chanB, vPeaks);
812 }
813 
814 // the private methods for the gaussian fit
815 // works in offset-zero
816 double CFitSpek::GFfunF(double xx, double *xpar)
817 {
818 /*
819  X dove si vuole calcolare (riferito all'inizio della regione di fit)
820  Par(0) = B0 fondo costante
821  Par(1) = B1 fondo slope
822  Par(2) = A ampiezza
823  Par(3) = P posizione
824  Par(4) = W sigma
825  Par(5) = S ampiezza step
826  Par(6) = L attacco coda sinistra (distanza da P, misurata in unit� di sigma; L < 0)
827  Par(7) = R attacco coda destra (distanza da P, misurata in unit� di sigma; R > 0)
828 
829  yy = (X-P)/W distanza da P in unit� di sigma
830 
831  Funct(yy) = B0 + B1*X + A * ( F1 + F2 )
832 
833  = exp(-0.5*L*(2*yy-L)) yy <= L
834  F1 = exp(-0.5*yy^2) L <= yy <= R
835  = exp(-0.5*R*(2*yy-R)) R <= yy
836 
837  F2 = S/(1+exp(yy))^2
838 */
839 
840  double amp, pos, sig, yy, yl, yr, ex, f1, f2, ey, eyp1i;
841 
842  int np, offs;
843  double fsum = 0.;
844  for(np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
845  amp = xpar[offs+0];
846  pos = xpar[offs+1];
847  sig = xpar[offs+2];
848  yy = (xx-pos)/sig;
849  yl = xpar[offs+4];
850  yr = xpar[offs+5];
851 
852  if(m_bTailL && yy < yl) {
853  ex = -0.5*yl*(2.*yy-yl);
854  }
855  else if(m_bTailR && yy > yr) {
856  ex = -0.5*yr*(2.*yy-yr);
857  }
858  else {
859  ex = -0.5*yy*yy;
860  }
861  f1 = exp(ex);
862 
863  f2 = 0;
864  if (m_bStep) {
865  ey = exp(yy);
866  eyp1i = 1./(1.+ey);
867  f2 = xpar[offs+3]*eyp1i*eyp1i;
868  }
869 
870  fsum += amp*(f1+f2);
871  }
872  return xpar[0] + xpar[1]*xx + fsum;
873 
874 }
875 
876 // calcola le derivate parziali della funzione nel punto i
877 void CFitSpek::GFderF(int id, double *xpar, double *deriv, double &deltay, double &weight)
878 {
879  double amp, pos, sig, stp, f0, f1, f2, f3;
880  double yy, yl, yr, ex, ey, dd;
881 
882  double xx = id + 0.5; // derivatives calculated at mid channel
883 
884  deltay = m_pDatVal[id] - GFfunF(xx, xpar);
885  weight = m_pDatErr[id];
886 
887  memset(deriv, 0, m_nNpar*sizeof(double));
888 
889  deriv[0] = 1.; // background level
890  deriv[1] = xx; // background slope
891 
892  int np, offs;
893  for(np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
894  amp = xpar[offs+0];
895  pos = xpar[offs+1];
896  sig = xpar[offs+2];
897  stp = xpar[offs+3];
898  yl = xpar[offs+4];
899  yr = xpar[offs+5];
900  yy = (xx-pos)/sig;
901 
902  if(m_bTailL && yy < yl) {
903  ex = -0.5*yl*(2.*yy-yl);
904  f1 = exp(ex);
905  f0 = yl/sig;
906  deriv[offs+4] = amp*f1*(yl-yy);
907  }
908  else if(m_bTailR && yy > yr) {
909  ex = -0.5*yr*(2.*yy-yr);
910  f1 = exp(ex);
911  f0 = yr/sig;
912  deriv[offs+5] = amp*f1*(yr-yy);
913  }
914  else {
915  ex = -0.5*yy*yy;
916  f1 = exp(ex);
917  f0 = yy/sig;
918  }
919 
920  f2 = f3 = 0;
921  if(m_bStep) {
922  ey = exp(yy);
923  f3 = 1./(1.+ey);
924  f2 = f3*f3;
925  f3 *= 2*ey/sig;
926  deriv[offs+3] = amp*f2;
927  }
928 
929  dd = amp*(f0*f1 + stp*f2*f3);
930 
931  deriv[offs+0] = f1+stp*f2; // amplitude
932  deriv[offs+1] = dd; // position
933  deriv[offs+2] = dd*yy; // sigma
934  }
935 
936  return;
937 }
938 
939 // Salva i risultati in m_pRes e m_pErr
940 void CFitSpek::GF_res()
941 {
942  m_dBack0 = m_pPar[0].Value; m_dBack0Err = m_pPar[0].Error;
943  m_dBack1 = m_pPar[1].Value; m_dBack1Err = m_pPar[1].Error;
944 
945  m_pRes = new CFitRes [m_nNpeaks];
946  m_pErr = new CFitRes [m_nNpeaks];
947 
948  int np, offs;
949  for(np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
950  m_pRes[np].Ampli = m_pPar[offs+0].Value;
951  m_pRes[np].Area = m_pPar[offs+0].Value * m_pPar[offs+2].Value * sqrt2pi;
952  m_pRes[np].Posi = m_pPar[offs+1].Value+m_nChanA;
953  m_pRes[np].Sigma = m_pPar[offs+2].Value;
954  m_pRes[np].Fwhm = m_pPar[offs+2].Value * s2fwhm;
955  m_pRes[np].Fw05 = m_pPar[offs+2].Value * s2fwhm;
956  m_pRes[np].Fw01 = m_pPar[offs+2].Value * s2fwtm;
957  m_pRes[np].W01L = sqrt(log(10.)/log(2.));
958  m_pRes[np].W01R = sqrt(log(10.)/log(2.));
959  m_pRes[np].Step = m_pPar[offs+3].Value;
960  m_pRes[np].TailL = m_pPar[offs+4].Value;
961  m_pRes[np].TailR = m_pPar[offs+5].Value;
962 
963  m_pErr[np].Ampli = m_pPar[offs+0].Error;
964  m_pErr[np].Area = m_pRes[np].Area*sqrt(pow(m_pPar[offs+0].Error/m_pPar[offs+0].Value,2.) +
965  pow(m_pPar[offs+2].Error/m_pPar[offs+2].Value,2.));
966  m_pErr[np].Posi = m_pPar[offs+1].Error;
967  m_pErr[np].Sigma = m_pPar[offs+2].Error;
968  m_pErr[np].Fwhm = m_pPar[offs+2].Error * s2fwhm;
969  m_pErr[np].Step = m_pPar[offs+3].Error;
970  m_pErr[np].TailL = m_pPar[offs+4].Error;
971  m_pErr[np].TailR = m_pPar[offs+5].Error;
972  }
973 
975  //for(np = 0; np < m_nNpeaks; np++) {
976  // double xxl = m_pRes[np].Posi - 10.*m_pRes[np].Fwhm;
977  // double xxr = m_pRes[np].Posi + 10.*m_pRes[np].Fwhm;
978  // double xxx;
979  // double area = 0., yval = 0.;
980  // for(xxx = xxl; xxx <= xxr; xxx += 0.1) {
981  // yval = GFitPeak(xxx, np);
982  // area += yval;
983  // }
984  // area = area/10;
985  // continue;
986  //}
987 
988  // aree e larghezze gia' calcolate analiticamente se gaussiane senza code
989  if(!m_bTailL && !m_bTailR)
990  return;
991 
992  // area dei picchi calcolata usando erf() e l'integrale analitico delle code esponenziali
993  for(np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
994  double area = 0.;
995  if(m_bTailL) {
996  double L = fabs(m_pPar[offs+4].Value);
997  double a = exp(-0.5*L*L)/L;
998  double b = sqrt2pi/2*GF_erf(L/m_sqrt2);
999  area += a + b;
1000  }
1001  else {
1002  area += sqrt2pi/2;
1003  }
1004  if(m_bTailR) {
1005  double R = fabs(m_pPar[offs+5].Value);
1006  double a = exp(-0.5*R*R)/R;
1007  double b = sqrt2pi/2*GF_erf(R/m_sqrt2);
1008  area += a + b;
1009  }
1010  else {
1011  area += sqrt2pi/2;
1012  }
1013  area *= m_pPar[offs+0].Value*m_pPar[offs+2].Value;
1014  m_pErr[np].Area *= area/m_pRes[np].Area; // errore solo ri-scalato ???
1015  m_pRes[np].Area = area;
1016  }
1017 
1018  for(np = 0; np < m_nNpeaks; np++) {
1019  m_pRes[np].Fw05 = GF_y2w(np, 0.5);
1020  m_pRes[np].Fw01 = GF_y2w(np, 0.1);
1021  m_pRes[np].W01L = GF_y2w(np, 0.1, -1)/(m_pRes[np].Fwhm/2);
1022  m_pRes[np].W01R = GF_y2w(np, 0.1, +1)/(m_pRes[np].Fwhm/2);
1023  }
1024 
1025  return;
1026 
1027 }
1028 
1029 // NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
1030 // Pag 221.
1031 // erfcc returns the complementary error function erfc(x)
1032 // with fractional error everywhere less than 1.2 10-7
1033 // Here it has been transformed to erf() and double
1034 double CFitSpek::GF_erf(double x)
1035 {
1036  double t,z,ans;
1037  z=fabs(x);
1038  t=1.0/(1.0+0.5*z);
1039  ans=t*exp(-z*z-1.26551223+
1040  t*( 1.00002368+
1041  t*( 0.37409196+
1042  t*( 0.09678418+
1043  t*(-0.18628806+
1044  t*( 0.27886807+
1045  t*(-1.13520398+
1046  t*( 1.48851587+
1047  t*(-0.82215223+
1048  t* 0.17087277 )))))))) );
1049 //return x >= 0.0 ? ans : 2.0-ans; // erfc
1050  return x >= 0.0 ? 1-ans : ans-1; // erf
1051 }
1052 
1053 // find actual width of the peak at level yval considering the tails
1054 double CFitSpek::GF_y2w(int np, double yval, int side)
1055 {
1056  const double precision = 1.e-6;
1057 
1058  if(!m_bGauss || np < 0 || np > m_nNpeaks) return 0.;
1059  if(yval <=0 || yval >= 1.) return 0.;
1060 
1061  CFitRes *pRes = m_pRes+ np;
1062  if(pRes->Ampli == 0) return 0.;
1063 
1064  double ff = 1./pRes->Ampli;
1065  double x0 = pRes->Posi;
1066 
1067  double dl = pRes->Sigma*sqrt(2.*log(1./yval));
1068  double xl = x0 - dl;
1069  double yl = GFitPeak(xl, np)*ff;
1070 
1071  while(fabs(yl/yval-1) > precision) {
1072  dl /= 4.;
1073  if(yl > yval) {
1074  while(yl > yval) {
1075  xl -= dl;
1076  yl = GFitPeak(xl, np)*ff;
1077  }
1078  }
1079  else {
1080  while(yl < yval) {
1081  xl += dl;
1082  yl = GFitPeak(xl, np)*ff;
1083  }
1084  }
1085  }
1086 
1087  double dr = pRes->Sigma*sqrt(2.*log(1./yval));
1088  double xr = x0 + dr;
1089  double yr = GFitPeak(xr, np)*ff;
1090 
1091  while(fabs(yr/yval-1) > precision) {
1092  dr /= 4.;
1093  if(yr > yval) {
1094  while(yr > yval) {
1095  xr += dr;
1096  yr = GFitPeak(xr, np)*ff;
1097  }
1098  }
1099  else {
1100  while(yr < yval) {
1101  xr -= dr;
1102  yr = GFitPeak(xr, np)*ff;
1103  }
1104  }
1105  }
1106 
1107  if(side < 0)
1108  return x0 - xl;
1109  else if (side > 0)
1110  return xr - x0;
1111  else
1112  return xr - xl;
1113 }
1114 
1115 // metodo pubblico (attraverso membro puntatore)
1116 // che lavora sui risultati del fit e calcola la funzione+fondo
1117 double CFitSpek::GFitFunc(double xx)
1118 {
1119  double yy, f1, f2, fsum, back, val;
1120 
1121  if(!m_bGauss) return 0.;
1122 
1123  fsum = 0.;
1124  for(int np = 0; np < m_nNpeaks; np++) {
1125  yy = (xx - m_pRes[np].Posi) / m_pRes[np].Sigma;
1126 
1127  f1 = 0.;
1128  if(m_bTailL && yy < m_pRes[np].TailL) {
1129  double yl = m_pRes[np].TailL;
1130  f1 = exp(-0.5*yl*(2.*yy-yl));
1131  }
1132  else if(m_bTailR && yy > m_pRes[np].TailR) {
1133  double yr = m_pRes[np].TailR;
1134  f1 = exp(-0.5*yr*(2.*yy-yr));
1135  }
1136  else {
1137  f1 = exp(-0.5*yy*yy);
1138  }
1139 
1140  f2 = 0.;
1141  if(m_bStep) {
1142  double ey = exp(yy);
1143  double ez = 1. / (1. + ey);
1144  f2 = m_pRes[np].Step * ez * ez;
1145  }
1146 
1147  fsum += m_pRes[np].Ampli * (f1 + f2);
1148  }
1149 
1150  back = 0;
1151  if(m_bFlat ) back += m_dBack0;
1152  if(m_bSlope) back += m_dBack1*(xx-m_nChanA);
1153 
1154  val = back + fsum;
1155  return val;
1156 
1157 }
1158 
1159 // metodo pubblico (attraverso membro puntatore)
1160 // che lavora sui risultati del fit e calcola la funzione per un picco
1161 double CFitSpek::GFitPeak(double xx, int np)
1162 {
1163  double yy, f1;
1164 
1165  if(!m_bGauss || np < 0 || np > m_nNpeaks) return 0.;
1166 
1167  yy = (xx - m_pRes[np].Posi) / m_pRes[np].Sigma;
1168 
1169  f1 = 0.;
1170  if(m_bTailL && yy < m_pRes[np].TailL) {
1171  double yl = m_pRes[np].TailL;
1172  f1 = exp(-0.5*yl*(2.*yy-yl));
1173  }
1174  else if(m_bTailR && yy > m_pRes[np].TailR) {
1175  double yr = m_pRes[np].TailR;
1176  f1 = exp(-0.5*yr*(2.*yy-yr));
1177  }
1178  else {
1179  f1 = exp(-0.5*yy*yy);
1180  }
1181 
1182  f1 *= m_pRes[np].Ampli;
1183 
1184  return f1;
1185 
1186 }
1187 
1188 // metodo pubblico (attraverso membro puntatore)
1189 // che lavora sui risultati del fit e calcola il solo fondo
1190 double CFitSpek::GFitBack(double xx)
1191 {
1192  double yy, fsum, back, val;
1193 
1194  if(!m_bGauss) return 0.;
1195 
1196  fsum = 0.;
1197  if(m_bStep) {
1198  for(int np = 0; np < m_nNpeaks; np++) {
1199  yy = (xx - m_pRes[np].Posi) / m_pRes[np].Sigma;
1200  double ey = exp(yy);
1201  double ez = 1. / (1. + ey);
1202  double f2 = m_pRes[np].Step * ez * ez;
1203  fsum += m_pRes[np].Ampli * f2;
1204  }
1205  }
1206 
1207  back = 0;
1208  if(m_bFlat ) back += m_dBack0;
1209  if(m_bSlope) back += m_dBack1*(xx-m_nChanA);
1210 
1211  val = back + fsum;
1212  return val;
1213 
1214 }
1215 
1216 // Minimizzazione del chi2 con il metodo Marquardt, secondo Bevington-1 pag. 237
1217 // Puntatori alle funzioni m_pfFunct, m_pfDeriv sono membri privati della classe CFitSpek
1218 // Matrici [][] usate come 1D-array indirizzate per righe (C-style) con la macro M12(i1,i2,l2)
1219 
1220 double CFitSpek::Curfit(int ndat, int npar, int mpar)
1221 {
1222  double *vvalue = new double [npar]; // vettore dei parametri di fit
1223  double *verror = new double [npar]; // vettore degli errori dei parametri di fit
1224  double *deriv = new double [npar]; // vettore delle derivate calcolato in m_pfDeriv()
1225 
1226  double *beta = new double [mpar]; // Vettore del gradiente
1227  double *alpha = new double [mpar*mpar]; // Matrice di curvatura
1228  double *sqrp = new double [mpar]; // sqrt(alpha[i][i])
1229  double *array = new double [mpar*mpar]; // matrice (modificata) di curvatura
1230 
1231  int i, j, k;
1232  double dyi,wwi, wwiderj;
1233  double dpar;
1234 
1235  // riporta i parametri di fit
1236  for(i = 0; i < npar; i++) vvalue[i] = m_pPar[i].Value;
1237 
1238  double chisqr = Chisqr(vvalue); // chi2 iniziale
1239  double chisqt, chisq1, deter;
1240  double flamda = 0.001; // valore iniziale variabile controllo
1241  int nsteps = 0; // conta il numero di passi effettuati
1242 
1244  // max 20 iterazione se chi2 aumenta dopo la minimizzazione
1245 
1246  for(int iter = 0; iter < 20; iter++) {
1247  // il vecchio valore di chi2
1248  chisqt = chisqr;
1249 
1250  // accumula alpha e beta
1251  memset(beta, 0, sizeof(double)*mpar);
1252  memset(alpha, 0, sizeof(double)*mpar*mpar);
1253  for(i = 0; i < ndat; i++) {
1254  // calcola le derivate
1255  (this->*m_pfDeriv)(i, vvalue, deriv, dyi, wwi);
1256 
1257  // impacchetta le derivate
1258  if(mpar < npar) Manage(deriv, 1);
1259 
1260  // accumula alpha e beta sugli mpar parametri liberi
1261  for(j = 0; j < mpar; j++) {
1262  if(deriv[j]) {
1263  wwiderj = wwi*deriv[j];
1264  beta[j] += dyi*wwiderj;
1265  for(k = 0; k <= j; k++) {
1266  alpha[M12(j,k,mpar)] += deriv[k]*wwiderj;
1267  }
1268  }
1269  }
1270  }
1271  // riempie il triangolo superiore
1272  for(j=0; j < mpar; j++) {
1273  for(k=0; k < j; k++) {
1274  alpha[M12(k,j,mpar)] = alpha[M12(j,k,mpar)];
1275  }
1276  }
1277 
1278  // con questa modifica degli elementi diagonali nulli, e' possibile continuare
1279  // a usare tutta la matrice e cioe' anche righe e colonne di parametri fissi
1280  // e' ancora necessaria ???
1281  for(j = 0; j < mpar; j++) {
1282  if(alpha[M12(j,j,mpar)] < 1.e-15)
1283  alpha[M12(j,j,mpar)] = 1.e-15;
1284  }
1285 
1286  // salva le normalizzazioni per array
1287  for(j = 0; j < mpar; j++) {
1288  sqrp[j] = sqrt(alpha[M12(j,j,mpar)]);
1289  }
1290 
1291  // loop su flambda
1292  while(true) {
1293 
1294  // calcola la matrice di curvatura modificata
1295  for(j = 0; j < mpar; j++) {
1296  for(k = 0; k < mpar; k++) {
1297  array[M12(j,k,mpar)] = alpha[M12(j,k,mpar)]/(sqrp[j]*sqrp[k]) ;
1298  }
1299  array[M12(j,j,mpar)] = 1. + flamda ;
1300  }
1301 
1302  // inversione con gauss-jordan
1303  deter = Matinv(array, mpar, mpar);
1304  if(deter <= 0.) { // errore
1305  chisqr = -1.;
1306  goto cleanup;
1307  }
1308 
1309  // impacchetta le variabili per poter calcolare gli incrementi
1310  if(mpar < npar) Manage(vvalue, 2);
1311 
1312  // calcola gli incrementi
1313  for(j = 0; j < mpar; j++) {
1314  dpar = 0.;
1315  for(k = 0; k < mpar; k++) {
1316  dpar += array[M12(j,k,mpar)]*beta[k]/(sqrp[j]*sqrp[k]);
1317  }
1318  vvalue[j] += dpar;
1319  }
1320 
1321  // spacchettamento delle variabili
1322  if(mpar < npar) Manage(vvalue, 3);
1323 
1324  // verifica dei parametri legati
1325  Manage(vvalue, 4);
1326 
1327  // il nuovo valore di chi2
1328  chisq1 = Chisqr(vvalue);
1329 
1330  // se diminuito, accetta i valori ed eventualmente re-itera
1331  if(chisq1 < chisqr) {
1332  chisqr = chisq1;
1333  break;
1334  }
1335 
1336  // scarta gli incrementi appena calcolati
1337  for(i = 0; i < npar; i++)
1338  vvalue[i] = m_pPar[i].Value;
1339  // ri-prova con flambda piu' grande
1340  flamda *= 10.;
1341  if(flamda > 1.e20) // ma solo fino a questo valore
1342  goto cleanup;
1343 
1344  } // fine loop su flambda
1345 
1346  // accetta i nuovi parametri e decide se re-iterare
1347  for(i = 0; i < npar; i++) m_pPar[i].Value = vvalue[i];
1348  nsteps++;
1349 
1350  // calcola gli errori (impacchettati!!)
1351  for(j = 0; j < mpar; j++) {
1352  if(alpha[M12(j,j,mpar)] == 0)
1353  verror[j] = 0.;
1354  else
1355  verror[j] = sqrt(array[M12(j,j,mpar)]/alpha[M12(j,j,mpar)]);
1356  }
1357  //spacchetta gli errori
1358  if(mpar < npar) Manage(verror, 5);
1359 
1360  if((chisqt-chisqr)/chisqr < 0.001) // chi2 stabilizzato --> fine delle iterazioni
1361  break;
1362 
1363  // diminuisce flambda (fino a 10^-7) e re-itera
1364  if(flamda > 1.e-6)
1365  flamda *= 0.1;
1366 
1367  } // fine loop sulle iterazioni
1368 
1369 cleanup:
1370  delete [] deriv;
1371  delete [] beta;
1372  delete [] alpha;
1373  delete [] sqrp;
1374  delete [] array;
1375  delete [] vvalue;
1376  delete [] verror;
1377 
1378  return (nsteps) ? (chisqr) : (-2);
1379 }
1380 
1381 // calcolo del chiquadro
1382 double CFitSpek::Chisqr(double *xpar)
1383 {
1384  double xx = 0.5;
1385  double chi=0;
1386  for(int ii = 0; ii < m_nNdat; ii++, xx += 1.) {
1387  double xval = (this->*m_pfFunct)(xx, xpar);
1388  double dval = m_pDatVal[ii] - xval;
1389  chi += m_pDatErr[ii] * dval*dval;
1390  }
1391  return chi;
1392 }
1393 
1394 // gestione dell'impacchettamento e dei limiti
1395 void CFitSpek::Manage(double *xpar, int flag)
1396 {
1397  int j, k;
1398 
1399  switch (flag) {
1400  case 1:
1401  // somma le derivate dei parametri legati
1402  for(j = 0; j < m_nNpar; j++) {
1403  if(m_pPar[j].Status == p_Linked) {
1404  xpar[m_pPar[j].LinkTo] += xpar[j];
1405  }
1406  }
1407  // impacchettamento delle derivate
1408  k = 0;
1409  for(j = 0; j < m_nNpar; j++) {
1410  if(m_pPar[j].Status == p_Free) {
1411  xpar[k] = xpar[j];
1412  k++;
1413  }
1414  }
1415  break;
1416  case 2:
1417  // impacchettamento dei parametri liberi
1418  k = 0;
1419  for(j = 0; j < m_nNpar; j++) {
1420  if(m_pPar[j].Status == p_Free) {
1421  xpar[k] = xpar[j];
1422  m_pPar[j].PackId = k;
1423  k++;
1424  }
1425  }
1426  break;
1427  case 3:
1428  // spacchettamento dei parametri liberi
1429  for(j = m_nNpar-1; j >=0; j--) {
1430  if(m_pPar[j].Status == p_Free)
1431  xpar[j] = xpar[m_pPar[j].PackId];
1432  }
1433  // riprende quelli legati e quelli fissi
1434  for(j = m_nNpar-1; j >=0; j--) {
1435  if(m_pPar[j].Status == p_Fixed)
1436  xpar[j] = m_pPar[j].Value;
1437  else if(m_pPar[j].Status == p_Linked)
1438  xpar[j] = xpar[m_pPar[j].LinkTo];
1439  }
1440  break;
1441  case 4:
1442  // verifica i limiti imposti ai parametri liberi
1443  for(j = 0; j < m_nNpar; j++) {
1444  if(m_pPar[j].Status == p_Free && m_pPar[j].Bounds != p_BNone) {
1445  if(m_pPar[j].Bounds & p_BLow)
1446  xpar[j] = max(m_pPar[j].Valmin, xpar[j]);
1447  if(m_pPar[j].Bounds & p_BHig)
1448  xpar[j] = min(m_pPar[j].Valmax, xpar[j]);
1449  }
1450  }
1451  // aggiorna i valori dei parametri legati ad altri parametri
1452  // che possono essere stati modificati dai limiti
1453  for(j = 0; j < m_nNpar; j++) {
1454  if(m_pPar[j].Status == p_Linked) {
1455  xpar[j] = xpar[m_pPar[j].LinkTo];
1456  }
1457  }
1458  break;
1459  case 5:
1460  // spacchettamento degli errori
1461  for(j = m_nNpar-1; j >=0; j--) {
1462  if(m_pPar[j].Status == p_Free)
1463  m_pPar[j].Error = xpar[m_pPar[j].PackId];
1464  else
1465  m_pPar[j].Error = 0;
1466  }
1467  for(j = 0; j < m_nNpar; j++) {
1468  if(m_pPar[j].Status == p_Fixed)
1469  m_pPar[j].Error = 0.;
1470  else if(m_pPar[j].Status == p_Linked)
1471  m_pPar[j].Error = m_pPar[m_pPar[j].LinkTo].Error;
1472  }
1473  break;
1474  }
1475 }
1476 
1477 // inversione di matrice simmetrica definita positiva
1478 // con il metodo di eliminazione Gauss-Jordan e pivot
1479 // Bevington-1 pag.302
1480 // Matrice [mpar][mpar] usa un unico blocco di memoria (vettore) con elementi
1481 // della matrice indirizzati per righe (C-style) con la macro M12(i1,i2,mpar)
1482 // Viene invertita solo la sottomatrice [npar][npar]
1483 
1484 double CFitSpek::Matinv(double *array, int mpar, int npar)
1485 {
1486  int ii, jj, kk;
1487  double det, amax, damax, dd;
1488 
1489  if(npar > mpar) return 0.;
1490  if(npar < 1) return 0.;
1491  if(npar == 1) {
1492  amax = array[0];
1493  if(amax == 0) return 0.;
1494  array[0] = 1./amax;
1495  return amax;
1496  }
1497 
1498  int *iik = new int[npar];
1499  int *jjk = new int[npar];
1500 
1501  det = 1.;
1502  for(kk = 0; kk < npar; kk++) {
1503  // find largest element in rest of array
1504  amax = 0;
1505  for(ii = kk; ii < npar; ii++) {
1506  for(jj = kk; jj < npar; jj++) {
1507  dd = array[M12(ii,jj,mpar)];
1508  if(fabs(dd) > fabs(amax)) {
1509  amax = dd;
1510  iik[kk] = ii;
1511  jjk[kk] = jj;
1512  }
1513  }
1514  }
1515 
1516  if(amax == 0.) {
1517  det = 0.;
1518  goto cleanup;
1519  }
1520  // interchange rows and columns to bring amax into array[kk][kk]
1521  damax = 1./amax;
1522  ii = iik[kk];
1523  if(ii > kk) {
1524  for(jj = 0; jj < npar; jj++) {
1525  dd = array[M12(kk,jj,mpar)];
1526  array[M12(kk,jj,mpar)] = array[M12(ii,jj,mpar)];
1527  array[M12(ii,jj,mpar)] = -dd;
1528  }
1529  }
1530  jj = jjk[kk];
1531  if(jj > kk) {
1532  for(ii = 0; ii < npar; ii++) {
1533  dd = array[M12(ii,kk,mpar)];
1534  array[M12(ii,kk,mpar)] = array[M12(ii,jj,mpar)];
1535  array[M12(ii,jj,mpar)] = -dd;
1536  }
1537  }
1538 
1539  // accumulate elements of inverse matrix
1540  for(ii = 0; ii < npar; ii++) {
1541  if(ii != kk)
1542  array[M12(ii,kk,mpar)] *= -damax;
1543  }
1544  for(ii = 0; ii < npar; ii++) {
1545  if(ii == kk) continue;
1546  for(jj = 0; jj < npar; jj++) {
1547 // if(ii != kk && jj != kk)
1548  if(jj != kk)
1549  array[M12(ii,jj,mpar)] += array[M12(ii,kk,mpar)] * array[M12(kk,jj,mpar)];
1550  }
1551  }
1552  for(jj = 0; jj < npar; jj++) {
1553  if(jj != kk)
1554  array[M12(kk,jj,mpar)] *= damax;
1555  }
1556  array[M12(kk,kk,mpar)] = damax;
1557  det *= amax;
1558  }
1559 
1560  // restore ordering of matrix;
1561  for(kk = npar-1; kk >= 0; kk--) {
1562  jj = iik[kk];
1563  if(jj > kk) {
1564  for(ii = 0; ii < npar; ii++) {
1565  dd = array[M12(ii,kk,mpar)];
1566  array[M12(ii,kk,mpar)] = -array[M12(ii,jj,mpar)];
1567  array[M12(ii,jj,mpar)] = dd;
1568  }
1569  }
1570  ii = jjk[kk];
1571  if(ii > kk) {
1572  for(jj = 0; jj < npar; jj++) {
1573  dd = array[M12(kk,jj,mpar)];
1574  array[M12(kk,jj,mpar)] = -array[M12(ii,jj,mpar)];
1575  array[M12(ii,jj,mpar)] = dd;
1576  }
1577  }
1578  }
1579 
1580 cleanup:
1581  delete [] iik;
1582  delete [] jjk;
1583  return det;
1584 }
1585 
1587 // Calcola la costante di decadimento esponenziale //
1588 // su un fondo piatto calcolato in [bgA--bgB] //
1590 
1591 bool CFitSpek::CalcExponential(float *pSpek, int chA, int chB, int bgA, int bgB)
1592 {
1593  // nothing valid yet
1594  Reset();
1595 
1596  m_dBack0 = 0.;
1597  m_dBack1 = 0.;
1598  m_dBack0Err = 0.;
1599  m_dBack1Err = 0.;
1600 
1601  // fondo piatto
1602  int ndat;
1603  if(bgA >-1 && bgB > -1) {
1604  if(bgA > bgB) {int tmp = bgA; bgA = bgB; bgB = tmp;}
1605  ndat = bgB-bgA+1;
1606  for(int nn = bgA; nn <= bgB; nn++)
1607  m_dBack0 += pSpek[nn];
1608  m_dBack0 /= ndat;
1609  m_bFlat = true;
1610  m_bSlope = false;
1611  }
1612 
1613  if(chA > chB) {int tmp = chA; chA = chB; chB = tmp;}
1614  ndat = chB-chA+1;
1615  if(ndat < 2)
1616  return false;
1617 
1618  double s1 = 0., sx = 0., sxx = 0., sy = 0., syx = 0.;
1619  double x, y, e;
1620  x = 0.5; // considera il centro del canale
1621  float * pd = pSpek + chA;
1622  for(int ii = chA; ii <= chB; ii++, x += 1.) {
1623  y = *pd++ - m_dBack0;
1624  if(y > 0) {
1625  y = log(y);
1626  e = 1.; // come trattare l'errore ??
1627  s1 += e;
1628  sx += e*x;
1629  sxx += e*x*x;
1630  sy += e*y;
1631  syx += e*y*x;
1632  }
1633  }
1634  double deter = s1*sxx-sx*sx;
1635  if(deter < 0.)
1636  return false;
1637 
1638  m_bValid = true;
1639  m_bExp = true;
1640  m_bBckgr = m_bFlat || m_bSlope;
1641  m_bFunct = true;
1642  m_bLine = true;
1643  m_nNpeaks = 0;
1644 
1645  m_nChanA = chA;
1646  m_nChanB = chB;
1647  m_nBgChanA = bgA;
1648  m_nBgChanB = bgB;
1649 
1650  m_dExpAmpli = (sxx*sy - sx*syx)/deter;
1651  m_dExpDecay = (-sx*sy + s1*syx)/deter;
1652  m_dExpAmpli = exp(m_dExpAmpli);
1653  m_dExpDecay = -1./m_dExpDecay;
1654 
1655  // puntatori alle funzioni risultato del fit
1656  m_pfFitFunc = &CFitSpek::ExpFunc; // exp + back
1657  m_pfFitPeak = &CFitSpek::ExpPeak; // exp only
1658  m_pfFitBack = &CFitSpek::LineFunc; // background only
1659 
1660 
1661  return true;
1662 }
1663 
1664 // metodo pubblico (attraverso membro puntatore)
1665 // che lavora sui risultati del fit e ritorna esponenziale + fondo
1666 double CFitSpek::ExpFunc(double xx)
1667 {
1668  if(!m_bExp) return 0.;
1669 
1670  double val;
1671  xx -= m_nChanA;
1672  val = m_dExpAmpli * exp(-xx/m_dExpDecay);
1673 
1674  if(m_bFlat)
1675  val += m_dBack0;
1676 
1677  return val;
1678 }
1679 
1680 // metodo pubblico (attraverso membro puntatore)
1681 // che lavora sui risultati del fit e ritorna esponenziale + fondo
1682 double CFitSpek::ExpPeak(double xx, int /*nn*/)
1683 {
1684  if(!m_bExp) return 0.;
1685 
1686  double val;
1687  xx -= m_nChanA;
1688  val = m_dExpAmpli * exp(-xx/m_dExpDecay);
1689 
1690  return val;
1691 }
1692 
1694 // Fit of sine wave over a stright background //
1696 
1697 /*
1698  X dove si vuole calcolare (riferito all'inizio dello spettro)
1699  chanA Z canale di riferimento
1700  back0 B0 fondo costante
1701  back1 B1 slope del fondo riferito a chanA
1702  ampli A ampiezza della sinusoide
1703  omega B pulsazione della sinusoide
1704  phase C fase della sinusoide
1705 
1706 
1707  Funct(X) = B0 + B1 * (X-Z) + A * sin(B*X + C)
1708 
1709 */
1710 
1711 // adattato da CalcGaussianFit
1712 int CFitSpek::CalcSinusoid(float *pSpek, int chanA, int chanB, double period, bool bFlat, bool bSlope)
1713 {
1714  // nothing valid yet
1715  Reset();
1716 
1717  // puntatori alle funzione per Curfit()
1718  m_pfFunct = &CFitSpek::SFfunF; // funzione
1719  m_pfDeriv = &CFitSpek::SFderF; // derivata
1720 
1721  m_nOffB = 2;
1722  m_nOffN = 3;
1723 
1724  int ii, np;
1725  int ierr=0;
1726  double chi2, chi2n;
1727 
1728  // Verify fit region
1729  if(chanA < chanB) {
1730  m_nChanA = chanA;
1731  m_nChanB = chanB;
1732  }
1733  else {
1734  m_nChanA = chanB;
1735  m_nChanB = chanA;
1736  }
1737  m_nBgChanA = m_nChanA;
1738  m_nBgChanB = m_nChanB;
1739 
1740  // nothing to fit if region empty
1741  bool empty = true;
1742  for(ii = m_nChanA; ii <= m_nChanB; ii++) {
1743  if(pSpek[ii]) {
1744  empty = false;
1745  break;
1746  }
1747  }
1748  if(empty) {
1749  ierr = -1;
1750  return ierr;
1751  }
1752 
1753  // nothing to fit if frequency too high
1754  if(period < 2) {
1755  ierr = -2;
1756  return ierr;
1757  }
1758 
1759  m_nNdat = m_nChanB - m_nChanA + 1; // margins included
1760  m_nNpeaks = 1;
1761  m_nNpar = m_nOffB + m_nOffN * m_nNpeaks;
1762 
1763  // create parameters structure
1764  m_pPar = new CFitPar [m_nNpar];
1765 
1766  // the data to fit
1767  m_pDatErr = new double [m_nNdat];
1768  m_pDatVal = new double [m_nNdat];
1769 
1770  // prepare vector of data and 1/(sigma^2)
1771  float *ps = pSpek + m_nChanA;
1772  for(ii = 0; ii < m_nNdat; ii++) {
1773  double yy = *ps++;
1774  m_pDatVal[ii] = yy;
1775  m_pDatErr[ii] = 1.;
1776  }
1777 
1779  // estimates for the first fit round //
1780  // limiting number of free parameters //
1782 
1783  double vsum = 0.;
1784  double vmin = m_pDatVal[0];
1785  double vmax = m_pDatVal[0];
1786  for(ii = 0; ii < m_nNdat; ii++) {
1787  double dd = m_pDatVal[ii];
1788  vsum += dd;
1789  vmin = min(vmin, dd);
1790  vmax = max(vmax, dd);
1791  }
1792  double bg0 = vsum / m_nNdat;
1793  double bg1 = 0.;
1794 
1795  m_bFlat = bFlat;
1796  m_bSlope = bSlope;
1797 
1798  // background level [0]
1799  if(m_bFlat) {
1800  m_pPar[0].Value = bg0;
1801  m_pPar[0].Status = p_Free;
1802  }
1803  else {
1804  m_pPar[0].Value = 0;
1805  m_pPar[0].Status = p_Fixed;
1806  }
1807 
1808  // backgrounr slope [1] --> no slope
1809  bSlope = m_bSlope;
1810  m_bSlope = false;
1811  m_pPar[1].Value = 0.;
1812  m_pPar[1].Status = p_Fixed;
1813 
1814  bg0 = m_pPar[0].Value;
1815  bg1 = m_pPar[1].Value;
1816 
1817  int offs; // keep this mechanism to add other components
1818  for(np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
1819  // amplitude --> limited;
1820  m_pPar[offs+0].Value = (vmax - vmin)/ 2;;
1821  m_pPar[offs+0].Status = p_Free;
1822  // frequency --> free
1823  m_pPar[offs+1].Value = 2*m_pi/period;
1824  m_pPar[offs+1].Status = p_Free;
1825  // phase --> free
1826  m_pPar[offs+2].Value = 0;
1827  m_pPar[offs+2].Status = p_Free;
1828  }
1829 
1830  // count number of free parameters
1831  m_nFpar = 0;
1832  for(np = 0; np < m_nNpar; np++) {
1833  if(m_pPar[np].Status == p_Free) m_nFpar++;
1834  }
1835 
1836  if(m_nFpar < 1) {
1837  ierr = -3; // nothing to fit
1838  goto cleanup;
1839  }
1840 
1841  if(m_nNdat - m_nFpar < 1) {
1842  ierr = -4; // insufficient degrees of freedom
1843  goto cleanup;
1844  }
1845 
1849 
1850  chi2 = Curfit(m_nNdat, m_nNpar, m_nFpar);
1851  if(chi2 < 0) {
1852  ierr = -5;
1853  goto cleanup;
1854  }
1855 
1859 
1860  // background level [0]
1861  m_pPar[0].Status = (m_bFlat) ? p_Free : p_Fixed;
1862 
1863  // backgrounr slope [1]
1864  m_bSlope = bSlope;
1865  m_pPar[1].Status = (m_bSlope) ? p_Free : p_Fixed;
1866 
1867  for(np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
1868  // amplitude --> free;
1869  m_pPar[offs+0].Status = p_Free;
1870  // frequency --> free
1871  m_pPar[offs+1].Status = p_Free;
1872  // phase --> free
1873  m_pPar[offs+2].Status = p_Free;
1874  }
1875 
1876  // count number of free parameters
1877  m_nFpar = 0;
1878  for(np = 0; np < m_nNpar; np++) {
1879  if(m_pPar[np].Status == p_Free) m_nFpar++;
1880  }
1881 
1882  if(m_nFpar < 1) {
1883  ierr = -6; // nothing to fit
1884  goto cleanup;
1885  }
1886 
1887  if(m_nNdat - m_nFpar < 1) {
1888  ierr = -7; // insufficient degrees of freedom
1889  goto cleanup;
1890  }
1891 
1892  chi2n = Curfit(m_nNdat, m_nNpar, m_nFpar);
1893  if(chi2n < 0) {
1894  if(chi2n != -2) { // -2 vuol dire che non e' riuscito a fare meglio
1895  ierr = -8;
1896  goto cleanup;
1897  }
1898  chi2n = chi2;
1899  }
1900 
1901  m_dChi2 = chi2n/(m_nNdat-m_nFpar);
1902  m_bValid = true;
1903  m_bSin = true;
1904  m_bBckgr = m_bFlat || m_bSlope;
1905  m_bLine = true;
1906  m_bFunct = true;
1907 
1908  // save result in data_menbers
1909  m_dBack0 = m_pPar[0].Value; m_dBack0Err = m_pPar[0].Error;
1910  m_dBack1 = m_pPar[1].Value; m_dBack1Err = m_pPar[1].Error;
1911 
1912  m_pRes = new CFitRes [m_nNpeaks];
1913  m_pErr = new CFitRes [m_nNpeaks];
1914 
1915  for(np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
1916  m_pRes[np].Ampli = m_pPar[offs+0].Value;
1917  m_pRes[np].Omega = m_pPar[offs+1].Value;
1918  m_pRes[np].Phase = m_pPar[offs+2].Value;
1919 
1920  m_pErr[np].Ampli = m_pPar[offs+0].Error;
1921  m_pErr[np].Omega = m_pPar[offs+1].Error;
1922  m_pErr[np].Phase = m_pPar[offs+2].Error;
1923  }
1924 
1925  // puntatori alle funzioni risultato del fit
1926  m_pfFitFunc = &CFitSpek::SFitFunc; // sinusoid + background
1927  m_pfFitPeak = &CFitSpek::SFitWave; // the sinusoid
1928  m_pfFitBack = &CFitSpek::LineFunc; // background only
1929 
1930 cleanup:
1931  delete [] m_pPar; m_pPar = NULL;
1932  delete [] m_pDatVal; m_pDatVal = NULL;
1933  delete [] m_pDatErr; m_pDatErr = NULL;
1934 
1935  return (m_bValid) ? (m_nNpeaks) : (ierr);
1936 }
1937 
1938 // the private methods for the sinusoidal fit
1939 // works in offset-zero
1940 double CFitSpek::SFfunF(double xx, double *xpar)
1941 {
1942 /*
1943  X dove si vuole calcolare (riferito all'inizio della regione di fit)
1944  Par(0) = B0 fondo costante
1945  Par(1) = B1 fondo slope
1946  Par(2) = A ampiezza
1947  Par(3) = B frequenza
1948  Par(4) = C fase
1949 
1950  Funct(X) = B0 + B1*X + A * sin(B*X + C)
1951 */
1952 
1953  double ampli, omega, phase, f1;
1954 
1955  int np, offs;
1956  double fsum = 0.;
1957  for(np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
1958  ampli = xpar[offs+0];
1959  omega = xpar[offs+1];
1960  phase = xpar[offs+2];
1961  f1 = ampli*sin(omega*xx+phase);
1962  fsum += f1;
1963  }
1964  return xpar[0] + xpar[1]*xx + fsum;
1965 }
1966 
1967 // calcola le derivate parziali della funzione nel punto i
1968 void CFitSpek::SFderF(int id, double *xpar, double *deriv, double &deltay, double &weight)
1969 {
1970  double ampli, omega, phase, angle, cosvv;
1971 
1972  double xx = id + 0.5; // derivatives calculated at mid channel
1973 
1974  deltay = m_pDatVal[id] - SFfunF(xx, xpar);
1975  weight = 1.;
1976 
1977  memset(deriv, 0, m_nNpar*sizeof(double));
1978 
1979  deriv[0] = 1.; // background level
1980  deriv[1] = xx; // background slope
1981 
1982  int np, offs;
1983  for(np = 0, offs = m_nOffB; np < m_nNpeaks; np++, offs += m_nOffN) {
1984  ampli = xpar[offs+0];
1985  omega = xpar[offs+1];
1986  phase = xpar[offs+2];
1987  angle = omega*xx+phase;
1988  cosvv = cos(angle) * ampli;
1989 
1990  deriv[offs+0] = sin(angle);
1991  deriv[offs+1] = cosvv * xx ;
1992  deriv[offs+2] = cosvv;
1993  }
1994 
1995  return;
1996 }
1997 
1998 // metodo pubblico (attraverso membro puntatore)
1999 // che lavora sui risultati del fit e calcola la funzione+fondo
2000 double CFitSpek::SFitFunc(double xx)
2001 {
2002  double yy, f1, fsum, back, val;
2003 
2004  if(!m_bSin) return 0.;
2005 
2006  fsum = 0.;
2007  yy = xx - m_nChanA;
2008  for(int np = 0; np < m_nNpeaks; np++) {
2009  f1 = m_pRes[np].Ampli * sin( m_pRes[np].Omega * yy + m_pRes[np].Phase);
2010  fsum += f1;
2011  }
2012 
2013  back = 0;
2014  if(m_bFlat ) back += m_dBack0;
2015  if(m_bSlope) back += m_dBack1*(xx-m_nChanA);
2016 
2017  val = back + fsum;
2018  return val;
2019 
2020 }
2021 
2022 // metodo pubblico (attraverso membro puntatore)
2023 // che lavora sui risultati del fit e calcola la funzione per un picco
2024 double CFitSpek::SFitWave(double xx, int np)
2025 {
2026  double yy, f1;
2027 
2028  if(!m_bValid || np < 0 || np > m_nNpeaks) return 0.;
2029 
2030  yy = xx - m_nChanA;
2031  f1 = m_pRes[np].Ampli * sin( m_pRes[np].Omega * yy + m_pRes[np].Phase);
2032 
2033  return f1;
2034 
2035 }
b
TBrowser * b
Definition: adf/Macros/Deprecated/Load.C:128
CFitSpek::m_pfFunct
double(CFitSpek::* m_pfFunct)(double, double *)
Definition: FitSpek.h:208
CFitSpek::m_pfFitFunc
double(CFitSpek::* m_pfFitFunc)(double)
Definition: FitSpek.h:203
p_BHig
const int p_BHig
Definition: FitSpek.h:29
__LoadWatchers__::tmp
TString tmp
Definition: GANILLoadWatchers.C:47
CFitSpek::FitFunc
double FitFunc(double x)
Definition: FitSpek.h:123
CFitSpek::m_pfFitBack
double(CFitSpek::* m_pfFitBack)(double)
Definition: FitSpek.h:204
m_pi
const double m_pi
Definition: FitSpek.h:41
p_Fixed
const int p_Fixed
Definition: FitSpek.h:25
m_sqrt2
const double m_sqrt2
Definition: FitSpek.h:42
CFitSpek::CalcGaussianFit
int CalcGaussianFit(float *pSpek, int chanA, int chanB, std::vector< double > &vPeaks)
Definition: FitSpek.cpp:426
CFitSpek::CalcStrightLine
bool CalcStrightLine(float *pSpek, int chanA, int chanB, int doSlope)
Definition: FitSpek.cpp:112
CFitSpek::CalcSinusoid
int CalcSinusoid(float *pSpek, int chanA, int chanB, double period, bool bFlat, bool bSlope)
Definition: FitSpek.cpp:1712
s2fwtm
const double s2fwtm
Definition: FitSpek.h:45
p_BBoth
const int p_BBoth
Definition: FitSpek.h:30
CFitSpek::m_pfDeriv
void(CFitSpek::* m_pfDeriv)(int, double *, double *, double &, double &)
Definition: FitSpek.h:209
CFitSpek::Sigma
double Sigma(int np)
Definition: FitSpek.h:90
CFitSpek::CalcPeakMoments
bool CalcPeakMoments(float *pSpek, int chanA, int chanB, int nChanBack)
Definition: FitSpek.cpp:248
CFitSpek::GFitParDlg
void GFitParDlg()
Definition: FitSpek.cpp:72
p_Linked
const int p_Linked
Definition: FitSpek.h:24
CFitSpek::CFitSpek
CFitSpek()
Definition: FitSpek.cpp:30
M12
#define M12(i1, i2, l2)
Definition: FitSpek.h:21
CFitSpek::m_pfFitPeak
double(CFitSpek::* m_pfFitPeak)(double, int)
Definition: FitSpek.h:205
p_BLow
const int p_BLow
Definition: FitSpek.h:28
CFitSpek::~CFitSpek
virtual ~CFitSpek()
Definition: FitSpek.cpp:54
p_BNone
const int p_BNone
Definition: FitSpek.h:27
s2fwhm
const double s2fwhm
Definition: FitSpek.h:44
CFitSpek::CalcExponential
bool CalcExponential(float *pSpek, int chanA, int chanB, int bgA=-1, int bgB=-1)
Definition: FitSpek.cpp:1591
p_Free
const int p_Free
Definition: FitSpek.h:23
CFitSpek::Reset
void Reset()
Definition: FitSpek.cpp:59
sqrt2pi
const double sqrt2pi
Definition: FitSpek.h:43