00001
00006 #include <stdlib.h>
00007 #include <math.h>
00008
00009 #include <bpm/bpm_messages.h>
00010 #include <bpm/bpm_process.h>
00011 #include <bpm/bpm_nr.h>
00012
00013 void _find_t0_startfit( double *wf, double ped,
00014 int peak_sample, double peak_value, double peak_fraction,
00015 int *start_sample ) {
00016
00017
00018 int S = 0;
00019 for ( S = peak_sample; S > 0; S-- ) {
00020 if ( ( ABS( wf[S] - ped ) > peak_fraction * peak_value ) &&
00021 ( ABS( wf[S] - ped ) > 6.5 ) ) {
00022 *start_sample = S;
00023 }
00024 }
00025 return;
00026 }
00027
00028
00029
00030 void _find_t0_endfit( double *wf, double ped,
00031 int peak_sample, double peak_value, double peak_fraction,
00032 int *end_sample ) {
00033
00034
00035 int S = 0;
00036 for ( S = 0; S < peak_sample; S++ ) {
00037 if ( ( ABS( wf[S] - ped ) ) < peak_fraction*peak_value ) {
00038 *end_sample = S;
00039 }
00040 }
00041 return;
00042 }
00043
00044
00045
00046 int get_t0( doublewf_t *signal, double *t0 ) {
00047
00048
00049 double adc_ped, adc_err, peak_value, sample_value, start_value;
00050 double gradient, intercept, *xval, *yval;
00051 int peak_sample, start_sample, counter, end_sample, i, retcode;
00052 double siga, sigb, chi2, q;
00053
00054 *t0 = -DBL_MAX;
00055
00056 if ( ! signal || ! t0 ) {
00057 bpm_error( "Invalid pointer arguments in get_t0(...)",
00058 __FILE__, __LINE__ );
00059 return BPM_FAILURE;
00060 }
00061
00062
00063 if( get_pedestal( signal, 20, &adc_ped, &adc_err ) == BPM_FAILURE ) {
00064 bpm_error( "Unable to retreive pedestal in get_t0(...)",
00065 __FILE__, __LINE__ );
00066 return BPM_FAILURE;
00067 }
00068
00069
00070
00071
00072 peak_value = adc_err * 4.;
00073 peak_sample = 0;
00074
00075 for ( i = 0; i < signal->ns - 10; ++i ) {
00076 sample_value = ABS( signal->wf[i] - adc_ped );
00077
00078 if ( sample_value > peak_value ) {
00079 peak_value = sample_value;
00080 peak_sample = i;
00081 }
00082
00083 }
00084
00085
00086 if ( peak_sample == 0 ) {
00087 bpm_error( "Could not find a peak in get_t0(...)",
00088 __FILE__, __LINE__ );
00089 return BPM_FAILURE;
00090 }
00091
00092
00093 _find_t0_endfit( signal->wf, adc_ped, peak_sample, peak_value, 0.9, &end_sample );
00094 _find_t0_startfit( signal->wf, adc_ped, peak_sample, peak_value, 0.1, &start_sample );
00095
00096
00097 if ( start_sample == end_sample ) {
00098 if ( bpm_verbose ) bpm_warning( "First fit initialisation failed, trying second",
00099 __FILE__, __LINE__ );
00100 _find_t0_endfit( signal->wf, adc_ped, peak_sample, peak_value, 0.95, &end_sample );
00101 _find_t0_startfit( signal->wf, adc_ped, peak_sample, peak_value, 0.05, &start_sample );
00102 }
00103
00104
00105 if ( start_sample == end_sample ) {
00106 if ( bpm_verbose ) bpm_warning( "Second fit initialisation failed, trying third",
00107 __FILE__, __LINE__ );
00108 _find_t0_endfit( signal->wf, adc_ped, peak_sample, peak_value, 0.975, &end_sample );
00109 _find_t0_startfit( signal->wf, adc_ped, peak_sample, peak_value, 0.025, &start_sample );
00110 }
00111
00112
00113 if ( ( end_sample == start_sample ) ||
00114 ( start_sample > end_sample ) ) {
00115 bpm_warning( "Cannot initialise fit, returning end_sample in get_t0(...)",
00116 __FILE__, __LINE__ );
00117 *t0 = (double) end_sample / signal->fs;
00118 return BPM_SUCCESS;
00119 }
00120
00121
00122
00123 xval = (double*) calloc( end_sample - start_sample + 1, sizeof(double) );
00124 yval = (double*) calloc( end_sample - start_sample + 1, sizeof(double) );
00125
00126 if ( ! ( xval && yval ) ) {
00127 bpm_error( "Coudn't allocate memory in get_t0(...)",
00128 __FILE__, __LINE__ );
00129 return BPM_FAILURE;
00130 }
00131
00132 counter = 0;
00133 for ( i = start_sample; i < end_sample + 1; i++ ) {
00134 yval[counter] = signal->wf[i] - adc_ped;
00135 xval[counter++] = (double) i;
00136 }
00137
00138
00139 if ( nr_fit( xval, yval, end_sample - start_sample + 1, NULL, 0, &intercept, &gradient,
00140 &siga, &sigb, &chi2, &q ) == BPM_FAILURE ) {
00141 bpm_error( "T0 straight line fit failed in get_t0(...)", __FILE__, __LINE__ );
00142 *t0 = -DBL_MAX;
00143 retcode = BPM_FAILURE;
00144 } else {
00145
00146 if( ABS( gradient ) > 0. ) {
00147 *t0 = ( peak_value * 0.5 - intercept ) / gradient;
00148 if ( *t0 < 0 || *t0 > signal->ns ) {
00149 bpm_error( "Fitted t0 value out of range in get_t0(...)",
00150 __FILE__, __LINE__ );
00151 *t0 = -DBL_MAX;
00152 retcode = BPM_FAILURE;
00153 } else {
00154
00155 *t0 /= signal->fs;
00156 retcode = BPM_SUCCESS;
00157 }
00158 } else {
00159 bpm_error( "Gradient in t0 fit equals 0 in get_t0(...)",
00160 __FILE__, __LINE__ );
00161 *t0 = -DBL_MAX;
00162 retcode = BPM_FAILURE;
00163 }
00164
00165 }
00166
00167
00168 free( xval );
00169 free( yval );
00170
00171 return retcode;
00172 }