CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Public Attributes | Private Attributes
TSFit Class Reference

#include <TSFit.h>

Inheritance diagram for TSFit:

Public Member Functions

double fit_third_degree_polynomial (double *, double *)
 
double fpol3dg (int, double *, double *, double *)
 
void init_errmat (double)
 
double inverms (int, double xx[5][5], double yy[5][5])
 
void set_params (int, int, int, int, int, double, double, int, int)
 
 TSFit (int size=14, int size_sh=650)
 
virtual ~TSFit ()
 

Public Attributes

int plshdim
 
int sdim
 

Private Attributes

double acc [14]
 
double adcp [14]
 
double adfmx [14]
 
double al [5][5]
 
double alpha_th
 
double avtm
 
double be [5][5]
 
double beta_th
 
double corel [14]
 
double cov [5][5]
 
double der [14][5]
 
double errmat [14][14]
 
double f [14]
 
double ff [14][4]
 
int iinf
 
double invcov [5][5]
 
int isup
 
double maskp3 [14]
 
int n_presamples
 
int n_samples_aft_max
 
int n_samples_bef_max
 
double nbcor [14]
 
int nbr_iter_fit
 
int nbs
 
int nmxu_sto
 
double norme
 
double parfp3 [10]
 
int sample_flag [14]
 
double t [14]
 
double xki2_max
 
double z [14]
 

Detailed Description

Definition at line 16 of file TSFit.h.

Constructor & Destructor Documentation

TSFit::TSFit ( int  size = 14,
int  size_sh = 650 
)

Definition at line 19 of file TSFit.cc.

References plshdim, sdim, and findQualityFiles::size.

19  {
20  sdim = size;
21  plshdim = size_sh;
22  // sample_flag = new int[size];
23 // t = new double[sdim];
24 // z = new double[sdim];
25 // f = new double[sdim];
26 // acc = new double[sdim];
27 // adfmx = new double[sdim];
28 // adcp = new double[sdim];
29 // maskp3 = new double[sdim];
30 // corel = new double[sdim];
31 // nbcor = new double[sdim];
32 // int k;
33 // ff = new double *[sdim];
34 // for(k=0;k<sdim;k++)ff[k] = new double[4];
35 // der = new double *[sdim];
36 // for(k=0;k<sdim;k++)der[k] = new double[5];
37 
38 
39 }
int plshdim
Definition: TSFit.h:68
int sdim
Definition: TSFit.h:67
tuple size
Write out results.
virtual TSFit::~TSFit ( )
inlinevirtual

Definition at line 72 of file TSFit.h.

72 {}

Member Function Documentation

double TSFit::fit_third_degree_polynomial ( double *  bdc,
double *  ret_dat 
)

Definition at line 256 of file TSFit.cc.

References adfmx, errmat, fpol3dg(), i, iinf, isup, relval_steps::k, maskp3, nbs, norme, parfp3, sample_flag, tmp, and xki2_max.

258  {
259  // third degree polynomial fit of the pulse summit.
260  // samples are contained in array bdc and must be pedestal
261  // substracted.
262  // only samples having sample_flag >= 1 are used.
263  // the unit of time is one clock unit, that is to say 25 ns.
264  // output: ret_dat[0] = pulse height
265  // ret_dat[1] position of maximum in the sample frame in clock units
266  // ret_dat[2] adc value of the highest sample
267  // ret_dat[3] number of the highest sample
268  // ret_dat[4] lower sample number used for fitting
269  // ret_dat[5] upper sample number used for fitting
270  // errmat inverse of the error matrix
271 
272  int i;
273  int nus;
274  double xki2;
275  double tm, tmp, amp;
276 
277 
278  ret_dat[0] = -999.;
279  ret_dat[1] = -999.;
280 
281  // search the maximum
282  double val_max = 0.;
283  int imax = 0;
284  for(i=0;i<nbs;i++){
285  if( sample_flag[i] == 0 )continue;
286  if( bdc[i] > val_max ){
287  val_max = bdc[i];
288  imax = i;
289  }
290  }
291 
292  if( (val_max*val_max) * errmat[imax][imax] < 16. )return -118;
293 
294  // if( imax != 9 )printf( "imax : %d !!!!!!!!!!!!!!!!!!!!!!!!!!!\n", imax );
295 
296  if( norme == 0. )norme = val_max;
297 
298  // look for samples above 1/3 of maximum before and 1/2 after
299  double val2 = val_max / 2.;
300  double val3 = val_max / 2.;
301  int ilow = iinf;
302  int ihig = 0;
303 
304  for(i=iinf;i<=isup;i++){
305  if( sample_flag[i] >= 1 ){
306  if( ( bdc[i] < val3 ) && ( i < imax ) )ilow = i;
307  if( bdc[i] > val2 )ihig = i;
308  }
309  }
310 
311  ilow++;
312 
313  //ilow = imax - 1;
314 
315  /* le test suivant, apparemment idiot, est mis a cause des sequences 0. 2048. qui apparaissent dans certains mauvais evts JPP 11/09/00 */
316 
317  if( ihig == ilow)return -105;
318  if( ilow == imax )ilow = ilow-1;
319  //if( ihig - ilow < 3 )ihig = ilow + 3;
320  ihig = ilow + 3;
321 
322  /* printf(" third degree: ilow %d ihig %d \n",ilow,ihig); */
323  nus=0;
324  int number_of_good_samples = 0;
325  for(i=ilow;i<=ihig;i++){
326  maskp3[nus] = 0;
327  adfmx[nus] = 0.;
328 /* printf(" adc %f sample_flag %d number_of good_samples %d \n",bdc[i],sample_flag[i],number_of_good_samples); */
329  if( sample_flag[i] >= 1 ){
330  adfmx[nus] = bdc[i];
331  maskp3[nus] = 1.;
332  number_of_good_samples++;
333  }
334  nus++;
335  }
336 
337  if( number_of_good_samples < 4 ){
338  return( -106 );
339  }
340 
341  xki2 = fpol3dg( nus, &parfp3[0], &maskp3[0], &adfmx[0]);
342 
343  /* printf( "fpol3dg-----------------------------------> %f %f %f %f %f\n",
344  parfp3[0], parfp3[1], parfp3[2], parfp3[3], parfp3[4] ); */
345 
346  tm = parfp3[4] + (float)ilow;
347  amp = parfp3[5];
348 
349  if( amp * amp * errmat[0][0] < 2. )return -101.;
350  tmp = parfp3[6] + (float)ilow;
351 
352  /*
353  validation of fit quality. Most of the time the fit is done with
354  four samples, therefore there is no possible ki2 check. When more than
355  4 samples are used the ki2 is often bad. So, in order to suppress some
356  events with bad samples, a consistency check on the position of the
357  maximum and minimum of the 3rd degree polynomial is used.
358  */
359 
360  if( xki2 > xki2_max ){
361  return -102.;
362  }
363  if( (tm < (double)ilow ) || (tm > (double)ihig)){
364  return -103.;
365  }
366 
367  if( (tmp > (double)ilow ) && (tmp < (double)ihig - 1.) ){
368  return -104.;
369  }
370 
371 
372  ret_dat[0] = amp;
373  ret_dat[1] = tm;
374  ret_dat[2] = val_max;
375  ret_dat[3] = (double)imax;
376  ret_dat[4] = (double)ilow;
377  ret_dat[5] = (double)ihig;
378  ret_dat[6] = (double)tmp;
379  int k;
380  for(i=0;i<4;i++){
381  k=i+7;
382  ret_dat[k] = parfp3[i];
383  }
384 
385  return xki2;
386 }
int i
Definition: DBlmapReader.cc:9
int sample_flag[14]
Definition: TSFit.h:51
double maskp3[14]
Definition: TSFit.h:58
int nbs
Definition: TSFit.h:23
double xki2_max
Definition: TSFit.h:32
double fpol3dg(int, double *, double *, double *)
Definition: TSFit.cc:103
double norme
Definition: TSFit.h:32
double adfmx[14]
Definition: TSFit.h:56
double errmat[14][14]
Definition: TSFit.h:50
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
int isup
Definition: TSFit.h:25
int iinf
Definition: TSFit.h:25
double parfp3[10]
Definition: TSFit.h:42
double TSFit::fpol3dg ( int  nmxul,
double *  parom,
double *  mask,
double *  adc 
)

Definition at line 103 of file TSFit.cc.

References cov, delta, errmat, ff, h, i, invcov, inverms(), relval_steps::k, prof2calltree::l, alignCSCRings::s, mathSSE::sqrt(), t, and tmp.

Referenced by fit_third_degree_polynomial().

106  {
107  // fit third degree polynomial
108  // nmxul array adc[] length
109  // parom return parameters (a0,a1,a2,a3,pos max,height)
110  // fplo3dg uses only the diagonal terms of errmat[][]
111  // errmat inverse of the error matrix
112 
113  int i, k, l;
114  double h, t2, tm, delta, tmp;
115  double xki2, dif, difmx, deglib;
116  double bv[4], s;
117 
118  deglib=(double)nmxul - 4.;
119  for(i=0;i<nmxul;i++){
120  t[i]=i;
121  ff[i][0]=1.;
122  ff[i][1]=t[i];
123  ff[i][2]=t[i]*t[i];
124  ff[i][3]=ff[i][2]*t[i];
125  }
126  /* computation of covariance matrix */
127  for(k=0;k<4;k++){
128  for(l=0;l<4;l++){
129  s=0.;
130  for(i=0;i<nmxul;i++){
131  s=s+ff[i][k]*ff[i][l]*errmat[i][i]*mask[i];
132  }
133  cov[k][l]=s;
134  }
135  s=0.;
136  for(i=0;i<nmxul;i++){
137  s=s+ff[i][k]*adc[i]*errmat[i][i]*mask[i];
138  }
139  bv[k]=s;
140  }
141  /* parameters */
142  inverms( 4, cov, invcov );
143  for(k=0;k<4;k++){
144  s=0.;
145  for(l=0;l<4;l++){
146  s=s+bv[l]*invcov[l][k];
147  }
148  parom[k]=s;
149  }
150 
151  if( parom[3] == 0. ){
152  parom[4] = -1000.;
153  parom[5] = -1000.;
154  parom[6] = -1000.;
155  return 1000000.;
156  }
157  /* worst hit and ki2 */
158  xki2=0.;
159  difmx=0.;
160  for(i=0;i<nmxul;i++ ){
161  t2=t[i]*t[i];
162  h= parom[0]+parom[1]*t[i]+parom[2]*t2+parom[3]*t2*t[i];
163  dif=(adc[i]-h)*mask[i];
164  xki2=xki2+dif*dif*errmat[i][i];
165  if(dif > difmx) {
166  difmx=dif;
167  }
168  }
169  if(deglib > 0.5) xki2=xki2/deglib;
170  /* amplitude and maximum position */
171  delta=parom[2]*parom[2]-3.*parom[3]*parom[1];
172  if(delta > 0.){
173  delta=sqrt(delta);
174  tm=-(delta+parom[2])/(3.*parom[3]);
175  tmp=(delta-parom[2])/(3.*parom[3]);
176  }
177  else{
178  parom[4] = -1000.;
179  parom[5] = -1000.;
180  parom[6] = -1000.;
181  return xki2;
182  }
183  parom[4]=tm;
184  parom[5]= parom[0]+parom[1]*tm+parom[2]*tm*tm+parom[3]*tm*tm*tm;
185  parom[6]=tmp;
186  return xki2;
187 }
int adc(sample_type sample)
get the ADC sample (12 bits)
dbl * delta
Definition: mlp_gen.cc:36
double invcov[5][5]
Definition: TSFit.h:37
int i
Definition: DBlmapReader.cc:9
T sqrt(T t)
Definition: SSEVec.h:48
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
double inverms(int, double xx[5][5], double yy[5][5])
Definition: TSFit.cc:188
double errmat[14][14]
Definition: TSFit.h:50
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
double cov[5][5]
Definition: TSFit.h:37
double ff[14][4]
Definition: TSFit.h:62
double t[14]
Definition: TSFit.h:52
void TSFit::init_errmat ( double  noise_initialvalue)

Definition at line 85 of file TSFit.cc.

References errmat, i, j, and sdim.

85  {
86  // Input: noise_initial value noise (in adc channels) as read in the
87  // data base.
88 /*------------------------------------------------------------------------*/
89 
90  int i, j;
91  //double one_over_noisesq;
92 
93  //one_over_noisesq = 1. / ( noise_initialvalue * noise_initialvalue );
94  for(i=0;i<sdim;i++ ){
95  for(j=0;j<sdim;j++){
96  errmat[i][j] = noise_initialvalue;
97  //errmat[i][j] = 0.;
98  }
99  //errmat[i][i] = one_over_noisesq;
100  }
101 }
int i
Definition: DBlmapReader.cc:9
int j
Definition: DBlmapReader.cc:9
int sdim
Definition: TSFit.h:67
double errmat[14][14]
Definition: TSFit.h:50
double TSFit::inverms ( int  n,
double  xx[5][5],
double  yy[5][5] 
)

Definition at line 188 of file TSFit.cc.

References al, be, g, i, j, findQualityFiles::jj, relval_steps::k, gen::n, alignCSCRings::r, alignCSCRings::s, and mathSSE::sqrt().

Referenced by fpol3dg().

188  {
189  // inversion of a positive definite symetric matrix of size n
190 
191  int i, j, k, jj;
192  double r, s;
193  double deter=0;
194 
195  /* initialisation */
196 
197  if( n > matdim ){
198  printf(
199  "ERROR : trying to use TSFit::inverms with size %d( max allowed %d\n",
200  n, matdim );
201  return -999.;
202  }
203 
204  int zero = 0;
205  memset( (char *)al, zero, 8*n*n );
206  memset( (char *)be, zero, 8*n*n );
207  /*
208  for(i=0;i<n;i++){
209  for(j=0;j<n;j++){
210  al[i][j] = 0.;
211  be[i][j] = 0.;
212  }
213  }
214  */
215  /* decomposition en vecteurs sur une base orthonormee */
216  al[0][0] = sqrt( g[0][0] );
217  for(i=1;i<n;i++){
218  al[i][0] = g[0][i] / al[0][0];
219  for(j=1;j<=i;j++){
220  s=0.;
221  for(k=0;k<=j-1;k++){
222  s = s + al[i][k] * al[j][k];
223  }
224  r= g[i][j] - s;
225  if( j < i ) al[i][j] = r / al[j][j];
226  if( j == i ) al[i][j] = sqrt( r );
227  }
228  }
229  /* inversion de la matrice al */
230  be[0][0] = 1./al[0][0];
231  for(i=1;i<n;i++){
232  be[i][i] = 1. / al[i][i];
233  for(j=0;j<i;j++){
234  jj=i-j-1;
235  s=0.;
236  for(k=jj+1;k<=i;k++){
237  s=s+ be[i][k] * al[k][jj];
238  }
239  be[i][jj]=-s/al[jj][jj];
240  }
241  }
242  /* calcul de la matrice ginv */
243  for(i=0;i<n;i++){
244  for(j=0;j<n;j++){
245  s=0.;
246  for(k=0;k<n;k++){
247  s=s+ be[k][i] * be[k][j];
248  }
249  ginv[i][j]=s;
250  }
251  }
252 
253  return deter;
254 }
int i
Definition: DBlmapReader.cc:9
#define matdim
Definition: TSFit.h:9
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
double be[5][5]
Definition: TSFit.h:38
T sqrt(T t)
Definition: SSEVec.h:48
int j
Definition: DBlmapReader.cc:9
double al[5][5]
Definition: TSFit.h:38
void TSFit::set_params ( int  n_samples,
int  niter,
int  n_presmpl,
int  sample_min,
int  sample_max,
double  time_of_max,
double  chi2_max,
int  nsbm,
int  nsam 
)

Definition at line 41 of file TSFit.cc.

References alpha_th, avtm, beta_th, iinf, isup, relval_steps::k, n_presamples, n_samples_aft_max, n_samples_bef_max, nbr_iter_fit, nbs, norme, sample_flag, and xki2_max.

44  {
45  // parameters initialisation of the fit package
46  // nsbm n_samples_bef_max, nsam n_samples_aft_max
47 
48  nbs = n_samples;
49  nbr_iter_fit = niter;
50  n_presamples = n_presmpl;
51  iinf = sample_min;
52  isup = sample_max;
53  avtm = time_of_max;
54  xki2_max = chi2_max;
55  n_samples_bef_max = nsbm;
56  n_samples_aft_max = nsam;
57 
58  norme = 0.;
59  alpha_th = 2.20;
60  beta_th = 1.11;
61 
62  int k;
63 
64  for(k=0;k<=nbs;k++){
65  sample_flag[k] = 0;
66  }
67 
68  for(k=sample_min;k<=sample_max;k++){
69  sample_flag[k] = 1;
70  }
71  /*
72  int lim1 = ( iinf > n_presamples ) ? n_presamples : iinf;
73  for(int k=lim1;k<=sample_max;k++){
74  sample_flag[k] = 2;
75  }
76  */
77  // printf( "sample_fag : " );
78 // for(k=0;k<=nbs;k++){
79 // printf( "%1d ", sample_flag[k] );
80 // }
81 // printf( "\n" );
82 
83 }
double avtm
Definition: TSFit.h:26
int sample_flag[14]
Definition: TSFit.h:51
int n_samples_bef_max
Definition: TSFit.h:27
double alpha_th
Definition: TSFit.h:35
int nbs
Definition: TSFit.h:23
double xki2_max
Definition: TSFit.h:32
int n_presamples
Definition: TSFit.h:24
double norme
Definition: TSFit.h:32
double beta_th
Definition: TSFit.h:35
int nbr_iter_fit
Definition: TSFit.h:36
int n_samples_aft_max
Definition: TSFit.h:28
int isup
Definition: TSFit.h:25
int iinf
Definition: TSFit.h:25

Member Data Documentation

double TSFit::acc[14]
private

Definition at line 55 of file TSFit.h.

double TSFit::adcp[14]
private

Definition at line 57 of file TSFit.h.

double TSFit::adfmx[14]
private

Definition at line 56 of file TSFit.h.

Referenced by fit_third_degree_polynomial().

double TSFit::al[5][5]
private

Definition at line 38 of file TSFit.h.

Referenced by inverms().

double TSFit::alpha_th
private

Definition at line 35 of file TSFit.h.

Referenced by set_params().

double TSFit::avtm
private

Definition at line 26 of file TSFit.h.

Referenced by set_params().

double TSFit::be[5][5]
private

Definition at line 38 of file TSFit.h.

Referenced by inverms().

double TSFit::beta_th
private

Definition at line 35 of file TSFit.h.

Referenced by set_params().

double TSFit::corel[14]
private

Definition at line 59 of file TSFit.h.

double TSFit::cov[5][5]
private

Definition at line 37 of file TSFit.h.

Referenced by fpol3dg().

double TSFit::der[14][5]
private

Definition at line 63 of file TSFit.h.

double TSFit::errmat[14][14]
private

Definition at line 50 of file TSFit.h.

Referenced by fit_third_degree_polynomial(), fpol3dg(), and init_errmat().

double TSFit::f[14]
private
double TSFit::ff[14][4]
private

Definition at line 62 of file TSFit.h.

Referenced by fpol3dg().

int TSFit::iinf
private

Definition at line 25 of file TSFit.h.

Referenced by fit_third_degree_polynomial(), and set_params().

double TSFit::invcov[5][5]
private

Definition at line 37 of file TSFit.h.

Referenced by fpol3dg().

int TSFit::isup
private

Definition at line 25 of file TSFit.h.

Referenced by fit_third_degree_polynomial(), and set_params().

double TSFit::maskp3[14]
private

Definition at line 58 of file TSFit.h.

Referenced by fit_third_degree_polynomial().

int TSFit::n_presamples
private

Definition at line 24 of file TSFit.h.

Referenced by set_params().

int TSFit::n_samples_aft_max
private

Definition at line 28 of file TSFit.h.

Referenced by set_params().

int TSFit::n_samples_bef_max
private

Definition at line 27 of file TSFit.h.

Referenced by set_params().

double TSFit::nbcor[14]
private

Definition at line 60 of file TSFit.h.

int TSFit::nbr_iter_fit
private

Definition at line 36 of file TSFit.h.

Referenced by set_params().

int TSFit::nbs
private

Definition at line 23 of file TSFit.h.

Referenced by fit_third_degree_polynomial(), and set_params().

int TSFit::nmxu_sto
private

Definition at line 34 of file TSFit.h.

double TSFit::norme
private

Definition at line 32 of file TSFit.h.

Referenced by fit_third_degree_polynomial(), and set_params().

double TSFit::parfp3[10]
private

Definition at line 42 of file TSFit.h.

Referenced by fit_third_degree_polynomial().

int TSFit::plshdim

Definition at line 68 of file TSFit.h.

Referenced by TSFit().

int TSFit::sample_flag[14]
private

Definition at line 51 of file TSFit.h.

Referenced by fit_third_degree_polynomial(), and set_params().

int TSFit::sdim

Definition at line 67 of file TSFit.h.

Referenced by init_errmat(), and TSFit().

double TSFit::t[14]
private

Definition at line 52 of file TSFit.h.

Referenced by fpol3dg().

double TSFit::xki2_max
private

Definition at line 32 of file TSFit.h.

Referenced by fit_third_degree_polynomial(), and set_params().

double TSFit::z[14]
private

Definition at line 53 of file TSFit.h.

Referenced by geometryXMLparser.Alignable::pos().