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 20 of file TSFit.cc.

References plshdim, sdim, and findQualityFiles::size.

20  {
21  sdim = size;
22  plshdim = size_sh;
23  // sample_flag = new int[size];
24 // t = new double[sdim];
25 // z = new double[sdim];
26 // f = new double[sdim];
27 // acc = new double[sdim];
28 // adfmx = new double[sdim];
29 // adcp = new double[sdim];
30 // maskp3 = new double[sdim];
31 // corel = new double[sdim];
32 // nbcor = new double[sdim];
33 // int k;
34 // ff = new double *[sdim];
35 // for(k=0;k<sdim;k++)ff[k] = new double[4];
36 // der = new double *[sdim];
37 // for(k=0;k<sdim;k++)der[k] = new double[5];
38 
39 
40 }
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 257 of file TSFit.cc.

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

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

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

Referenced by fit_third_degree_polynomial().

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

References errmat, i, j, and sdim.

86  {
87  // Input: noise_initial value noise (in adc channels) as read in the
88  // data base.
89 /*------------------------------------------------------------------------*/
90 
91  int i, j;
92  //double one_over_noisesq;
93 
94  //one_over_noisesq = 1. / ( noise_initialvalue * noise_initialvalue );
95  for(i=0;i<sdim;i++ ){
96  for(j=0;j<sdim;j++){
97  errmat[i][j] = noise_initialvalue;
98  //errmat[i][j] = 0.;
99  }
100  //errmat[i][i] = one_over_noisesq;
101  }
102 }
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 189 of file TSFit.cc.

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

Referenced by fpol3dg().

189  {
190  // inversion of a positive definite symetric matrix of size n
191 
192  int i, j, k, jj;
193  double r, s;
194  double deter=0;
195 
196  /* initialisation */
197 
198  if( n > matdim ){
199  printf(
200  "ERROR : trying to use TSFit::inverms with size %d( max allowed %d\n",
201  n, matdim );
202  return -999.;
203  }
204 
205  int zero = 0;
206  memset( (char *)al, zero, 8*n*n );
207  memset( (char *)be, zero, 8*n*n );
208  /*
209  for(i=0;i<n;i++){
210  for(j=0;j<n;j++){
211  al[i][j] = 0.;
212  be[i][j] = 0.;
213  }
214  }
215  */
216  /* decomposition en vecteurs sur une base orthonormee */
217  al[0][0] = sqrt( g[0][0] );
218  for(i=1;i<n;i++){
219  al[i][0] = g[0][i] / al[0][0];
220  for(j=1;j<=i;j++){
221  s=0.;
222  for(k=0;k<=j-1;k++){
223  s = s + al[i][k] * al[j][k];
224  }
225  r= g[i][j] - s;
226  if( j < i ) al[i][j] = r / al[j][j];
227  if( j == i ) al[i][j] = sqrt( r );
228  }
229  }
230  /* inversion de la matrice al */
231  be[0][0] = 1./al[0][0];
232  for(i=1;i<n;i++){
233  be[i][i] = 1. / al[i][i];
234  for(j=0;j<i;j++){
235  jj=i-j-1;
236  s=0.;
237  for(k=jj+1;k<=i;k++){
238  s=s+ be[i][k] * al[k][jj];
239  }
240  be[i][jj]=-s/al[jj][jj];
241  }
242  }
243  /* calcul de la matrice ginv */
244  for(i=0;i<n;i++){
245  for(j=0;j<n;j++){
246  s=0.;
247  for(k=0;k<n;k++){
248  s=s+ be[k][i] * be[k][j];
249  }
250  ginv[i][j]=s;
251  }
252  }
253 
254  return deter;
255 }
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:46
int j
Definition: DBlmapReader.cc:9
int k[5][pyjets_maxn]
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 42 of file TSFit.cc.

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

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