Skip to content

Commit 2f98f79

Browse files
author
dpscience
committed
saving all data available
saving all data available
1 parent d072387 commit 2f98f79

File tree

3 files changed

+66
-24
lines changed

3 files changed

+66
-24
lines changed

‎DReconvolutionInput.py‎

Lines changed: 14 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,16 @@
3333

3434
from DReconvolutionModel import ReconvolutionModel as reconvModel
3535

36+
#save output as *.txt file after success?
37+
__saveReconvolutionSpectrum = True
38+
__saveReconvolutionSpectrumPath = 'testData/recovolutionSpectrumOutput.txt'
39+
__saveReconvolutionSpectrumResidualPath = 'testData/recovolutionSpectrumResidualsOutput.txt'
40+
41+
#!note: IRF output is only saved if a model function is used --> __bUsingModel = True
42+
__saveReconvolutionIRF = True
43+
__saveReconvolutionIRFPath = 'testData/recovolutionIRFOutput.txt'
44+
__saveReconvolutionIRFResidualPath = 'testData/recovolutionIRFResidualsOutput.txt'
45+
3646
#expected number of components (number of exponential decay functions - LIMITED to MAX: 4):
3747
__numberOfExpDec = 2
3848

@@ -41,8 +51,8 @@
4151

4252
#expected lifetimes (tau) -> start values in [ps] (required for the levenberg marquardt fit)
4353
#note: only the first '__numberOfExpDec' related values are considered (e.g.: for __numberOfExpDec = 2 --> __expectedTau_1_in_ps AND __expectedTau_2_in_ps)
44-
__expectedTau_1_in_ps = 160.0;
45-
__expectedTau_2_in_ps = 455.0;
54+
__expectedTau_1_in_ps = 260.0;
55+
__expectedTau_2_in_ps = 1500.0;
4656
__expectedTau_3_in_ps = 160.0;
4757
__expectedTau_4_in_ps = 160.0;
4858

@@ -53,11 +63,11 @@
5363
#NOTE: Spectrum and IRF data vectors require equal length!!!
5464

5565
#file path which contains the SPECTRUM data:
56-
__filePathSpec = 'testData/spectrum_5ps.dat'
66+
__filePathSpec = 'testData/spectrum2_5ps.dat'
5767
__specDataDelimiter = '\t'
5868

5969
#file path which contains the IRF data:
60-
__filePathIRF = 'testData/irf_5ps.dat'
70+
__filePathIRF = 'testData/irf2_5ps.dat'
6171
__irfDataDelimiter = '\t'
6272

6373

‎DReconvolutionModel.py‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -58,8 +58,12 @@ def Lorentz_Cauchy(x, ampl, s, y0, x0, args=()):
5858
def Pseudovoigt1(x, ampl, a, sigma, s, y0, x0, args=()):
5959
G=np.zeros(x.size)
6060
L=np.zeros(x.size)
61-
G=(1.0/(sigma*np.sqrt(2*np.pi)))*np.exp(-0.5*((x-x0)/sigma)*((x-x0)/sigma))
62-
L=s/(np.pi*((x-x0)*(x-x0) + s*s))
61+
62+
#G=(2/sigma)*np.sqrt(np.log(2)/np.pi)*np.exp(-4*np.log(2)*((x-x0)/sigma)**2);
63+
#L=(2/np.pi)*sigma*(1/(4*(x-x0)**2+sigma**2));
64+
65+
G=(1.0/(sigma*np.sqrt(2*np.pi)))*np.exp(-0.5*((x-x0)/sigma)**2)
66+
L=s/(np.pi*((x-x0)**2 + s*s))
6367
return ampl*(a*G+(1-a)*L)+y0
6468

6569
def Pearson7(x, ampl, alpha, m, y0, x0, args=()):

‎DReconvolutionProc.py‎

Lines changed: 46 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,8 @@
4343
xSpec,ySpec = np.loadtxt(userInput.__filePathSpec, delimiter=userInput.__specDataDelimiter, unpack=True, dtype='float');
4444
xIRF,yIRF = np.loadtxt(userInput.__filePathIRF, delimiter=userInput.__irfDataDelimiter, unpack=True, dtype='float');
4545

46+
yIRFOrigin = yIRF;
47+
4648
print("shifting x to x = 0...")
4749

4850
#shift that: [0] = 0:
@@ -78,16 +80,19 @@
7880

7981
stddevIRF = 1.0
8082

81-
for i in range(0, len(xVal)-1):
83+
for i in range(0, len(xVal)):
8284
if yIRF[i] > yMaxIRF*0.5:
8385
stddevIRF = np.abs((xWhereYMaxIRF-xIRF[i]))/(2*np.sqrt(2*np.log(2)));
8486
break;
8587

86-
estimatedBkgrd = 0;
88+
estimatedBkgrd = 0;
89+
estimatedBkgrdIRF = 0;
8790
for i in range(userInput.__bkgrd_startIndex, userInput.__bkgrd_startIndex + userInput.__bkgrd_count):
88-
estimatedBkgrd += ySpec[i];
91+
estimatedBkgrd += ySpec[i];
92+
estimatedBkgrdIRF += yIRF[i];
8993

90-
estimatedBkgrd /= userInput.__bkgrd_count;
94+
estimatedBkgrd /= userInput.__bkgrd_count;
95+
estimatedBkgrdIRF /= userInput.__bkgrd_count;
9196

9297
#fit the IRF model function on data (xIRF, yIRF):
9398
if userInput.__bUsingModel:
@@ -99,7 +104,7 @@
99104
fitModelIRF.set_param_hint('ampl', min=0.0);
100105
fitModelIRF.set_param_hint('y0', min=0.0);
101106

102-
parameterListIRFFit = fitModelIRF.make_params(x=xVal, ampl=yMaxIRF, sigma=stddevIRF, y0=0, x0=xWhereYMaxIRF, args=yIRF);
107+
parameterListIRFFit = fitModelIRF.make_params(x=xVal, ampl=yMaxIRF, sigma=stddevIRF, y0=estimatedBkgrdIRF, x0=xWhereYMaxIRF, args=yIRF);
103108
#change here if you want to fix x0 and/or y0:
104109
parameterListIRFFit['x0'].vary = True;
105110
parameterListIRFFit['y0'].vary = True;
@@ -157,7 +162,7 @@
157162
fitModelIRF.set_param_hint('ampl', min=0.0);
158163
fitModelIRF.set_param_hint('y0', min=0.0);
159164

160-
parameterListIRFFit = fitModelIRF.make_params(x=xVal, ampl=yMaxIRF, s=stddevIRF, y0=0, x0=xWhereYMaxIRF, args=yIRF);
165+
parameterListIRFFit = fitModelIRF.make_params(x=xVal, ampl=yMaxIRF, s=stddevIRF, y0=estimatedBkgrdIRF, x0=xWhereYMaxIRF, args=yIRF);
161166
#change here if you want to fix x0 and/or y0:
162167
parameterListIRFFit['x0'].vary = True;
163168
parameterListIRFFit['y0'].vary = True;
@@ -216,7 +221,7 @@
216221
fitModelIRF.set_param_hint('ampl', min=0.0);
217222
fitModelIRF.set_param_hint('y0', min=0.0);
218223

219-
parameterListIRFFit = fitModelIRF.make_params(x=xVal, ampl=yMaxIRF, a=1, sigma=stddevIRF, s=stddevIRF, y0=0, x0=xWhereYMaxIRF, args=yIRF);
224+
parameterListIRFFit = fitModelIRF.make_params(x=xVal, ampl=yMaxIRF, a=0.8, sigma=stddevIRF, s=stddevIRF*0.5, y0=estimatedBkgrdIRF, x0=xWhereYMaxIRF, args=yIRF);
220225
#change here if you want to fix x0 and/or y0:
221226
parameterListIRFFit['x0'].vary = True;
222227
parameterListIRFFit['y0'].vary = True;
@@ -227,8 +232,8 @@
227232
chiSquare = resultsOfModelIRF.chisqr;
228233
redChiSquare = resultsOfModelIRF.redchi;
229234

230-
fract = (float)(resultsOfModelIRF.params['a'].value*userInput);
231-
fract_err = (float)(resultsOfModelIRF.params['a'].stderr*userInput);
235+
fract = (float)(resultsOfModelIRF.params['a'].value);
236+
fract_err = (float)(resultsOfModelIRF.params['a'].stderr);
232237

233238
sigma = (float)(resultsOfModelIRF.params['sigma'].value*userInput.__channelResolutionInPs);
234239
sigma_err = (float)(resultsOfModelIRF.params['sigma'].stderr*userInput.__channelResolutionInPs);
@@ -251,13 +256,13 @@
251256
print("background = {0} ({1})".format(yRes, yRes_err));
252257
print("");
253258
print("---------------------------------------------------------------");
254-
print("G - Gaussian: a = {0} ({})".format(fract, fract_err));
259+
print("G - Gaussian: a = {0} ({1})".format(fract, fract_err));
255260
print("---------------------------------------------------------------");
256261
print("stddev [ps] = {0} ({1})".format(sigma, sigma_err));
257262
print("FWHM [ps] = {0} ({1})".format(sigma*(2*np.sqrt(2*np.log(2))), sigma_err*(2*np.sqrt(2*np.log(2)))));
258263
print("");
259264
print("---------------------------------------------------------------");
260-
print("L - Lorentzian: (1 - a) = {0} ({})".format(1-fract, fract_err));
265+
print("L - Lorentzian: (1 - a) = {0} ({1})".format(1-fract, fract_err));
261266
print("---------------------------------------------------------------");
262267
print("s [ps] = {0} ({1})".format(s, s_err));
263268
print("---------------------------------------------------------------");
@@ -288,7 +293,7 @@
288293
fitModelIRF.set_param_hint('ampl', min=0.0);
289294
fitModelIRF.set_param_hint('y0', min=0.0);
290295

291-
parameterListIRFFit = fitModelIRF.make_params(x=xVal, ampl=yMaxIRF, alpha=stddevIRF, m=2, y0=0, x0=xWhereYMaxIRF, args=yIRF);
296+
parameterListIRFFit = fitModelIRF.make_params(x=xVal, ampl=yMaxIRF, alpha=stddevIRF, m=2, y0=estimatedBkgrdIRF, x0=xWhereYMaxIRF, args=yIRF);
292297
#change here if you want to fix x0 and/or y0:
293298
parameterListIRFFit['x0'].vary = True;
294299
parameterListIRFFit['y0'].vary = True;
@@ -302,8 +307,8 @@
302307
alpha = (float)(resultsOfModelIRF.params['alpha'].value*userInput.__channelResolutionInPs);
303308
alpha_err = (float)(resultsOfModelIRF.params['alpha'].stderr*userInput.__channelResolutionInPs);
304309

305-
m = (float)(resultsOfModelIRF.params['m'].value*userInput);
306-
m_err = (float)(resultsOfModelIRF.params['m'].stderr*userInput);
310+
m = (float)(resultsOfModelIRF.params['m'].value);
311+
m_err = (float)(resultsOfModelIRF.params['m'].stderr);
307312

308313
amplitude = (float)(resultsOfModelIRF.params['ampl'].value);
309314
amplitude_err = (float)(resultsOfModelIRF.params['ampl'].stderr);
@@ -483,7 +488,6 @@ def ExpDecay_4(x, ampl1, tau1, ampl2, tau2, ampl3, tau3, ampl4, tau4, y0, x0, ar
483488
ax2 = plt.subplot(2,1,2);
484489
ax2.set_title("Best fit: Residuals");
485490
plt.plot(xVal, resultsOfModelDecay.residual);
486-
plt.show();
487491

488492
if userInput.__numberOfExpDec == 2:
489493
print("\nrunning reconvolution with 2 component...\n")
@@ -568,7 +572,6 @@ def ExpDecay_4(x, ampl1, tau1, ampl2, tau2, ampl3, tau3, ampl4, tau4, y0, x0, ar
568572
ax2 = plt.subplot(2,1,2);
569573
ax2.set_title("Best fit: Residuals");
570574
plt.plot(xVal, resultsOfModelDecay.residual);
571-
plt.show();
572575

573576
if userInput.__numberOfExpDec == 3:
574577
print("\nrunning reconvolution with 3 component...\n")
@@ -667,7 +670,6 @@ def ExpDecay_4(x, ampl1, tau1, ampl2, tau2, ampl3, tau3, ampl4, tau4, y0, x0, ar
667670
ax2 = plt.subplot(2,1,2);
668671
ax2.set_title("Best fit: Residuals");
669672
plt.plot(xVal, resultsOfModelDecay.residual);
670-
plt.show();
671673

672674
if userInput.__numberOfExpDec == 4:
673675
print("\nrunning reconvolution with 4 component...\n")
@@ -787,5 +789,31 @@ def ExpDecay_4(x, ampl1, tau1, ampl2, tau2, ampl3, tau3, ampl4, tau4, y0, x0, ar
787789
ax2 = plt.subplot(2,1,2);
788790
ax2.set_title("Best fit: Residuals");
789791
plt.plot(xVal, resultsOfModelDecay.residual);
790-
plt.show();
792+
793+
#save data if required:
794+
if userInput.__saveReconvolutionSpectrum:
795+
ab = np.zeros(len(xVal), dtype=[('time [ps]', float), ('raw counts [#]', float), ('best fit [#]', float)]);
796+
ab['time [ps]'] = xVal*userInput.__channelResolutionInPs;
797+
ab['raw counts [#]'] = ySpec;
798+
ab['best fit [#]'] = resultsOfModelDecay.best_fit;
799+
np.savetxt(userInput.__saveReconvolutionSpectrumPath, ab, fmt='%10.3f\t%10.3f\t%10.3f', delimiter='\t', newline='\n', header='time [ps]\traw counts [#]\tbest fit [#]\n');
800+
801+
abRes = np.zeros(len(xVal), dtype=[('time [ps]', float), ('residuals [conv. level]', float)]);
802+
abRes['time [ps]'] = xVal*userInput.__channelResolutionInPs;
803+
abRes['residuals [conv. level]'] = resultsOfModelDecay.residual;
804+
np.savetxt(userInput.__saveReconvolutionSpectrumResidualPath, abRes, fmt='%10.3f\t%10.3f', delimiter='\t', newline='\n', header='time [ps]\traw counts [#]\tresiduals [conv. level]\n');
805+
806+
if userInput.__saveReconvolutionIRF and userInput.__bUsingModel:
807+
ab = np.zeros(len(xVal), dtype=[('time [ps]', float), ('raw counts [#]', float), ('best fit [#]', float)]);
808+
ab['time [ps]'] = xVal*userInput.__channelResolutionInPs;
809+
ab['raw counts [#]'] = yIRFOrigin;
810+
ab['best fit [#]'] = yIRF;
811+
np.savetxt(userInput.__saveReconvolutionIRFPath, ab, fmt='%10.3f\t%10.3f\t%10.3f', delimiter='\t', newline='\n', header='time [ps]\traw counts [#]\tbest fit [#]\n');
812+
813+
abRes = np.zeros(len(xVal), dtype=[('time [ps]', float), ('residuals [conv. level]', float)]);
814+
abRes['time [ps]'] = xVal*userInput.__channelResolutionInPs;
815+
abRes['residuals [conv. level]'] = resultsOfModelIRF.residual;
816+
np.savetxt(userInput.__saveReconvolutionIRFResidualPath, abRes, fmt='%10.3f\t%10.3f', delimiter='\t', newline='\n', header='time [ps]\traw counts [#]\tresiduals [conv. level]\n');
817+
818+
plt.show();
791819

0 commit comments

Comments
 (0)