ATLAS Offline Software
Public Member Functions | Private Member Functions | List of all members
LArWaveHelper Class Reference

#include <LArWaveHelper.h>

Collaboration diagram for LArWaveHelper:

Public Member Functions

LArWave translate (const LArWave &theWave, int nShift, double baseline=0.) const
 
LArWave Dtranslate (const LArWave &theWave, double tShift, double baseline=0.) const
 
unsigned int getMin (const LArWave &theWave) const
 return index of minimum sample More...
 
unsigned int getMax (const LArWave &theWave) const
 return index of maximum sample More...
 
double getDMax (const LArWave &theWave, double &tMax) const
 return amplitude aproximation from poly2 fit around maxima, and it's time More...
 
LArWave derive (const LArWave &theWave) const
 crude derivative More...
 
LArWave derive_smooth (const LArWave &theWave) const
 smoothed derivative More...
 
LArWave autocorr (const LArWave &theWave) const
 autocorrelation function (not normalized) More...
 
LArWave subSample (const LArWave &theWave, unsigned Nfirst, unsigned deltaN) const
 
double getBaseline (const LArWave &theWave, unsigned nBase) const
 
double getSumTail (const LArWave &theWave, unsigned iFirst) const
 
double getSumRegion (const LArWave &theWave, unsigned iFirst, unsigned iLast) const
 
double getSumSquareTail (const LArWave &theWave, unsigned iFirst) const
 
double getSumSquareRegion (const LArWave &theWave, unsigned iFirst, unsigned iLast) const
 
double getArea (const LArWave &theWave) const
 
double getWidth (const LArWave &theWave) const
 
double getMaxAmp (const LArWave &theWave) const
 
double getT0 (const LArWave &theWave) const
 
double getT5 (const LArWave &theWave) const
 
unsigned getStart (const LArWave &theWave) const
 
unsigned getZeroCross (const LArWave &theWave) const
 
double getJitter (const LArWaveCumul &theWave) const
 
LArWaveDerivedQuantitiesP getDerivedQuantities (const LArWaveCumul &theWave) const
 
std::vector< double > linfit (const LArWave &theWave, unsigned iFirst, unsigned iLast, double &rho) const
 
std::vector< double > expfit (const LArWave &theWave, unsigned iFirst, unsigned iLast, double &rho) const
 
std::vector< double > polyfit (const LArWave &theWave, unsigned iFirst, unsigned iLast, unsigned Ndeg) const
 
std::vector< LArWave > linearMasterWave (const std::vector< const LArWave * > &vWaves, const std::vector< double > &vAmpli) const
 

Private Member Functions

unsigned fact (unsigned N) const
 
std::vector< double > linfit (const std::vector< double > &X, const std::vector< double > &Y, double &rho) const
 
std::vector< double > expfit (const std::vector< double > &X, const std::vector< double > &Y, double &rho) const
 
std::vector< double > polyfit (const std::vector< double > &X, const std::vector< double > &Y, unsigned Ndeg) const
 
unsigned get_fit_vectors (const LArWave &theWave, unsigned iFirst, unsigned iLast, std::vector< double > &X, std::vector< double > &Y) const
 

Detailed Description

Definition at line 14 of file LArWaveHelper.h.

Member Function Documentation

◆ autocorr()

LArWave LArWaveHelper::autocorr ( const LArWave &  theWave) const

autocorrelation function (not normalized)

Definition at line 322 of file LArWaveHelper.cxx.

323 {
324  std::vector<double> wave = theWave.getWave() ;
325  double dt = theWave.getDt() ;
326  unsigned int nsamples = theWave.getSize() ;
327  std::vector<double> ac ;
328  ac.resize(nsamples,0.0) ;
329  for ( unsigned int i=0 ; i<nsamples ; i++ ) {
330  for ( unsigned int k=0 ; k<nsamples-i ; k++ ) {
331  ac[i] += wave[i]*wave[i+k] ;
332  }
333  ac[i] *= dt ;
334  }
335  return LArWave( ac , dt ) ;
336 }

◆ derive()

LArWave LArWaveHelper::derive ( const LArWave &  theWave) const

crude derivative

Definition at line 256 of file LArWaveHelper.cxx.

257 {
258  std::vector<double> wave = theWave.getWave() ;
259  const double dt = theWave.getDt() ;
260  unsigned int nsamples = theWave.getSize() ;
261  std::vector<double> dwave ;
262  dwave.resize(nsamples,0.0) ;
263  if ( nsamples > 1 ) {
264  const double inv_dt = 1. / dt;
265  dwave[0] = ( wave[1] - wave[0] ) * inv_dt ;
266  for ( unsigned int i=1 ; i<nsamples-1 ; i++ ) {
267  dwave[i] = ( wave[i+1] - wave[i-1] ) * inv_dt * 0.5;
268  }
269  dwave[nsamples-1] = ( wave[nsamples-1] - wave[nsamples-2] ) * inv_dt ;
270  }
271  return LArWave(dwave,dt) ;
272 }

◆ derive_smooth()

LArWave LArWaveHelper::derive_smooth ( const LArWave &  theWave) const

smoothed derivative

Definition at line 274 of file LArWaveHelper.cxx.

275 {
276  /* ================================================================
277  * 7 points first derivative
278  * using Savitsky-Golay smoothing filtering (m = polynomial order)
279  * ================================================================ */
280  static const double c[3][4][7] =
281  { { { -0.107, -0.071, -0.036, 0.000, 0.036, 0.071, 0.107 } ,
282  { 0.012, -0.071, -0.107, -0.095, -0.036, 0.071, 0.226 } ,
283  { 0.131, -0.071, -0.179, -0.190, -0.107, 0.071, 0.345 } ,
284  { 0.250, -0.071, -0.250, -0.286, -0.179, 0.071, 0.464 } } , // m = 2
285  { { 0.087, -0.266, -0.230, 0.000, 0.230, 0.266, -0.087 } ,
286  { 0.123, -0.183, -0.218, -0.095, 0.075, 0.183, 0.115 } ,
287  { -0.008, 0.067, -0.040, -0.190, -0.246, -0.067, 0.484 } ,
288  { -0.306, 0.484, 0.306, -0.286, -0.734, -0.484, 1.020 } } , // m = 3
289  { { 0.087, -0.266, -0.230, 0.000, 0.230, 0.266, -0.087 } ,
290  { -0.049, 0.219, -0.276, -0.439, 0.018, 0.584, -0.057 } ,
291  { -0.079, 0.234, -0.063, -0.333, -0.270, 0.099, 0.413 } ,
292  { 0.269, -0.856, 0.497, 0.863, -0.543, -1.824, 1.594 } } } ; // m = 4
293  const int m=3 ; // choose between 2 and 4
294  const int ipoly = m-2 ;
295  //
296  std::vector<double> wave = theWave.getWave() ;
297  const double dt = theWave.getDt() ;
298  const double inv_dt = 1. / dt;
299  unsigned int nsamples = theWave.getSize() ;
300  std::vector<double> dwave ;
301  dwave.resize(nsamples,0.0) ;
302  //
303  for ( unsigned int i=0 ; i<nsamples ; i++ ) {
304  if ( i<=2 ) { // first 3 points
305  int ifilt = 3 - i ;
306  for ( unsigned int j=0 ; j<=6 ; j++ )
307  dwave[i] += -c[ipoly][ifilt][6-j] * wave[j] ;
308  } else if ( i>=nsamples-3 ) { // last 3 points
309  int ifilt = nsamples - i + 2;
310  for ( unsigned int j=nsamples-7 ; j<=nsamples-1 ; j++ )
311  dwave[i] += c[ipoly][ifilt][j-(nsamples-7)] * wave[j] ;
312  } else { // intermediate points
313  int ifilt = 0 ;
314  for ( unsigned int j=i-3 ; j<=i+3 ; j++ )
315  dwave[i] += c[ipoly][ifilt][j-(i-3)] * wave[j] ;
316  }
317  dwave[i] *= inv_dt ;
318  }
319  return LArWave( dwave , dt ) ;
320 }

◆ Dtranslate()

LArWave LArWaveHelper::Dtranslate ( const LArWave &  theWave,
double  tShift,
double  baseline = 0. 
) const

Definition at line 28 of file LArWaveHelper.cxx.

30 {
31  unsigned int nsamples = theWave.getSize() ;
32  std::vector<double> trWave ;
33  trWave.resize( nsamples );
34  int ibin = 0;
35  const double Dt = theWave.getDt();
36  const double inv_Dt = 1. / Dt;
37  if(fabs(tShift) > Dt) { // move first a number of bins
38  ibin = int(tShift * inv_Dt);
39  for ( unsigned int i=0 ; i<theWave.getSize() ; ++i ) {
40  if ( (int)i + ibin >= 0 && (int)i + ibin < (int)nsamples) trWave[i+ibin] = theWave.getSample(i) ;
41  }
42  if(ibin > 0) {
43  for ( unsigned int i=0 ; i<(unsigned int)ibin ; ++i ) {
44  trWave[i] = baseline ;
45  }
46  } else {
47  for ( unsigned int i=nsamples+ibin; i<nsamples; ++i) {
48  trWave[i] = baseline;
49  }
50  }
51  } else {
52  for ( unsigned int i=0 ; i<theWave.getSize() ; ++i ) trWave[i] = theWave.getSample(i);
53  }
54 
55  // now move trWave for part of bin
56  std::vector<double> tranWave ;
57  tranWave.resize(trWave.size());
58  double tmpshift = tShift - ibin * Dt;
59  if(tmpshift > 0.) {
60  tranWave[0] = baseline;
61  for(unsigned int i=1; i<trWave.size(); ++i) {
62  tranWave[i] = trWave[i] - tmpshift * (trWave[i] - trWave[i-1]) * inv_Dt;
63  }
64  } else {
65  tranWave[trWave.size() - 1] = baseline;
66  for(int i=trWave.size()-2; i>=0; --i) {
67  tranWave[i] = trWave[i] - tmpshift * (trWave[i] - trWave[i+1]) * inv_Dt;
68  }
69  }
70 
71  return LArWave( tranWave , Dt ) ;
72 }

◆ expfit() [1/2]

std::vector< double > LArWaveHelper::expfit ( const LArWave &  theWave,
unsigned  iFirst,
unsigned  iLast,
double &  rho 
) const

Definition at line 479 of file LArWaveHelper.cxx.

481 {
482  std::vector<double> X,Y ;
483  get_fit_vectors(theWave,iFirst,iLast,X,Y) ;
484  return expfit(X,Y,rho) ;
485 }

◆ expfit() [2/2]

std::vector< double > LArWaveHelper::expfit ( const std::vector< double > &  X,
const std::vector< double > &  Y,
double &  rho 
) const
private

Definition at line 576 of file LArWaveHelper.cxx.

578 {
579  std::vector<double> C ;
580  C.resize(0) ;
581  rho = 0 ;
582  unsigned N = X.size() ;
583  if ( N>=2 && Y.size()==N ) {
584  std::vector<double> logY ;
585  logY.resize(N) ;
586  for ( unsigned i=0;i<N;i++ ) {
587  if ( Y[i]*Y[0] <= 0.0 ) {
588  fprintf(stderr,"ERROR: exp. fit of zero-crossing wave\n") ;
589  return C ;
590  }
591  logY[i] = log(fabs(Y[i])) ;
592  }
593  C = linfit(X,logY,rho) ;
594  C[0] = exp(C[0]) ;
595  if ( Y[0] < 0.0 ) C[0] = - C[0] ;
596  }
597  return C ;
598 }

◆ fact()

unsigned LArWaveHelper::fact ( unsigned  N) const
inlineprivate

Definition at line 63 of file LArWaveHelper.h.

63 { return ( N<2 ) ? 1 : N*fact(N-1) ; };

◆ get_fit_vectors()

unsigned LArWaveHelper::get_fit_vectors ( const LArWave &  theWave,
unsigned  iFirst,
unsigned  iLast,
std::vector< double > &  X,
std::vector< double > &  Y 
) const
private

Definition at line 661 of file LArWaveHelper.cxx.

663 {
664  unsigned length = theWave.getSize() ;
665  if ( length == 0 ) {
666  X.resize(0) ;
667  Y.resize(0) ;
668  return 0 ;
669  } else {
670  if ( iLast >= length ) iLast = length - 1 ;
671  if ( iFirst > iLast ) iFirst = iLast ;
672  unsigned vlength = iLast - iFirst + 1 ;
673  X.resize(vlength) ;
674  Y.resize(vlength) ;
675  double dt = theWave.getDt() ;
676  for ( unsigned i=iFirst ; i<=iLast ; i++ ) {
677  unsigned k = i - iFirst ;
678  X[k] = i * dt ;
679  Y[k] = theWave.getSample(i) ;
680  }
681  return vlength ;
682  }
683 }

◆ getArea()

double LArWaveHelper::getArea ( const LArWave &  theWave) const

Definition at line 141 of file LArWaveHelper.cxx.

142 {
143  const unsigned int N=theWave.getSize();
144  double surf=0.;
145  const unsigned Nbase=5;
146  double amax=-9.999E+99;
147  double base = getBaseline(theWave,Nbase);
148  for (unsigned int i=1; i<N ; i++) {
149  const double asample=theWave.getSample(i);
150  if((asample-base)>amax)amax = asample-base;
151  }
152 
153  const double inv_amax = 1. / amax;
154  for (unsigned int i=0; i<N ; i++){
155  const double asample=theWave.getSample(i)-base;
156  if(asample>0)surf+=asample*inv_amax;
157  }
158 
159  return surf;
160 }

◆ getBaseline()

double LArWaveHelper::getBaseline ( const LArWave &  theWave,
unsigned  nBase 
) const

Definition at line 351 of file LArWaveHelper.cxx.

352 {
353  if (nBase == 0) return 0;
354  double base = 0.0 ;
355  for ( unsigned i=0 ; i<nBase ; i++ ) base += theWave.getSample(i) ;
356  return base/nBase ;
357 }

◆ getDerivedQuantities()

LArWaveDerivedQuantitiesP LArWaveHelper::getDerivedQuantities ( const LArWaveCumul &  theWave) const

Definition at line 456 of file LArWaveHelper.cxx.

457 {
458  float baseline = getBaseline(theWave,10);
459  float maxAmp = getMaxAmp(theWave);
460  // please note the 'dt' normalization in the following!
461  float dt = theWave.getDt();
462  float tmaxAmp = dt * getMax(theWave);
463  float width = dt * getWidth(theWave);
464  float rT0 = dt * getT0(theWave);
465  float posLobe = dt * getSumRegion(theWave,0,getZeroCross(theWave))/maxAmp;
466  float jitter = getJitter(theWave);
467  return LArWaveDerivedQuantitiesP(baseline, maxAmp, tmaxAmp, width, rT0, posLobe, jitter, theWave.getFlag() );
468 }

◆ getDMax()

double LArWaveHelper::getDMax ( const LArWave &  theWave,
double &  tMax 
) const

return amplitude aproximation from poly2 fit around maxima, and it's time

Definition at line 103 of file LArWaveHelper.cxx.

104 {
105 // const unsigned int N = theWave.getSize() ;
106  unsigned int imax = getMax(theWave) ;
107  double amax = -9.999E+99 ;
108  std::vector<double> pfit = polyfit(theWave, imax - 4, imax + 4, 2);
109  if(pfit.size() < 3 || pfit[2] == 0.) { tMax = -9.999E+99; return amax; }
110  tMax = - pfit[1] / (2 * pfit[2]);
111  amax = pfit[0];
112  return amax ;
113 }

◆ getJitter()

double LArWaveHelper::getJitter ( const LArWaveCumul &  theWave) const

Definition at line 433 of file LArWaveHelper.cxx.

434 {
435  // Calculate the squared error of the CaliWave
436  std::vector<double> errorsSquared(theWave.getErrors());
437  for (unsigned int iError = 0; iError < errorsSquared.size(); ++iError) {
438  if (errorsSquared[iError] == 0.0)
439  return 0.; // no error in the Wave, cannot compute jitter
440  errorsSquared[iError] = errorsSquared[iError] * errorsSquared[iError];
441  }
442  // Get the wave derivative
443  LArWave derivedWave = derive(theWave);
444  LArWave derivativeSquared = derivedWave * derivedWave;
445  // Fit jitter
446  double rho = 0.;
447  std::vector<double> theLinearFit = linfit(derivativeSquared.getWave(),errorsSquared,rho);
448  // Output
449  if ( theLinearFit[1] > 0. && theLinearFit[1]==theLinearFit[1] ) // skip negative value and (hopefully) NAN
450  return sqrt(theLinearFit[1]);
451  else
452  return 0.;
453 }

◆ getMax()

unsigned int LArWaveHelper::getMax ( const LArWave &  theWave) const

return index of maximum sample

Definition at line 89 of file LArWaveHelper.cxx.

90 {
91  const unsigned int N = theWave.getSize() ;
92  unsigned int imax = N + 1 ;
93  double amax = -9.999E+99 , asample ;
94  for ( unsigned int i=1 ; i<N ; i++ ) {
95  if ( (asample=theWave.getSample(i)) > amax ) {
96  imax = i ;
97  amax = asample ;
98  }
99  }
100  return imax ;
101 }

◆ getMaxAmp()

double LArWaveHelper::getMaxAmp ( const LArWave &  theWave) const

Definition at line 129 of file LArWaveHelper.cxx.

130 {
131  const unsigned int N = theWave.getSize() ;
132  double amax = -9.999E+99;
133  for ( unsigned int i=1; i<N ; ++i ) {
134  double asample = theWave.getSample(i);
135  if ( asample>amax )
136  amax = asample;
137  }
138  return amax;
139 }

◆ getMin()

unsigned int LArWaveHelper::getMin ( const LArWave &  theWave) const

return index of minimum sample

Definition at line 74 of file LArWaveHelper.cxx.

75 {
76  const unsigned int N = theWave.getSize() ;
77  unsigned int imin = N + 1 ;
78  double amin = +9.999E+99 , asample ;
79  for ( unsigned int i=1 ; i<N ; i++ ) {
80  if ( (asample=theWave.getSample(i)) < amin ) {
81  imin = i ;
82  amin = asample ;
83  }
84  }
85  return imin ;
86 }

◆ getStart()

unsigned int LArWaveHelper::getStart ( const LArWave &  theWave) const

Definition at line 413 of file LArWaveHelper.cxx.

414 {
415  const unsigned Nbase = 5 ;
416  const double Threshold = 0.001 ;
417  double base = getBaseline(theWave,Nbase) ;
418  unsigned iMax = getMax(theWave) ;
419  double thr = base + Threshold * ( theWave.getSample( iMax ) - base ) ;
420  for ( unsigned i=iMax-1 ; i>0 ; i-- )
421  if (theWave.getSample(i)<thr) return i-1 ;
422  return 0 ;
423 }

◆ getSumRegion()

double LArWaveHelper::getSumRegion ( const LArWave &  theWave,
unsigned  iFirst,
unsigned  iLast 
) const

Definition at line 374 of file LArWaveHelper.cxx.

374  {
375 
376  if (iLast>theWave.getSize())
377  iLast=theWave.getSize();
378  if (iFirst >= iLast) return 0.0;
379 
380  double sum = 0.0;
381  std::vector<double>::const_iterator it=theWave.getWave().begin()+iFirst;
382  std::vector<double>::const_iterator it_e=theWave.getWave().begin()+iLast;
383  for(;it!=it_e;++it) {
384  sum+=(*it);
385  }
386 
387  return sum;
388 }

◆ getSumSquareRegion()

double LArWaveHelper::getSumSquareRegion ( const LArWave &  theWave,
unsigned  iFirst,
unsigned  iLast 
) const

Definition at line 396 of file LArWaveHelper.cxx.

397 {
398  if (iLast>theWave.getSize())
399  iLast=theWave.getSize();
400  if (iFirst >= iLast) return 0.0;
401 
402  double sum = 0.0;
403  std::vector<double>::const_iterator it=theWave.getWave().begin()+iFirst;
404  std::vector<double>::const_iterator it_e=theWave.getWave().begin()+iLast;
405  for(;it!=it_e;++it) {
406  sum+=((*it)*(*it));
407  }
408 
409  return sum;
410 }

◆ getSumSquareTail()

double LArWaveHelper::getSumSquareTail ( const LArWave &  theWave,
unsigned  iFirst 
) const

Definition at line 391 of file LArWaveHelper.cxx.

392 {
393  return getSumSquareRegion(theWave,iFirst,theWave.getSize());
394 }

◆ getSumTail()

double LArWaveHelper::getSumTail ( const LArWave &  theWave,
unsigned  iFirst 
) const

Definition at line 359 of file LArWaveHelper.cxx.

360 {
361  return getSumRegion(theWave,iFirst,theWave.getSize());
362 }

◆ getT0()

double LArWaveHelper::getT0 ( const LArWave &  theWave) const

Definition at line 162 of file LArWaveHelper.cxx.

163 {
164  const unsigned Nbase=5;
165  double rT0=0.;
166  double asamp1=0.;
167  double asamp2=0.;
168  double base = getBaseline(theWave,Nbase);
169  double amax = getMaxAmp(theWave);
170  unsigned int imax=getMax(theWave);
171 
172  for (unsigned int i=1; i<imax ; i++) {
173  asamp1=theWave.getSample(i);
174  asamp2=theWave.getSample(i-1);
175  if((asamp1-base)>amax*0.10 && (asamp2-base)<amax*0.10)
176  rT0 = i-1+(0.10*amax-(asamp2-base))/(asamp1-asamp2);
177  }
178 
179 
180  return rT0;
181 }

◆ getT5()

double LArWaveHelper::getT5 ( const LArWave &  theWave) const

Definition at line 183 of file LArWaveHelper.cxx.

184 {
185 // const unsigned Nbase=5;
186  double rT5=0.;
187 // double base = getBaseline(theWave,Nbase);
188  double tmax;
189  double amax = getDMax(theWave, tmax);
190  unsigned int imax=getMax(theWave);
191 
192  for (unsigned int i=imax - 20; i>0 ; --i) {
193  if(theWave.getSample(i) < 0.05*amax) {
194  std::vector<double> pfit = polyfit(theWave, i-3, i+3, 3);
195  for(double t=(i+2)*theWave.getDt(); t>=(i-2)*theWave.getDt(); t-= 0.05) {
196  if(pfit[0]+pfit[1]*t+pfit[2]*t*t+pfit[3]*t*t*t < 0.05*amax) {
197  rT5 = t+0.05;
198  break;
199  }
200  }
201  break;
202  }
203  }
204 
205  return rT5;
206 }

◆ getWidth()

double LArWaveHelper::getWidth ( const LArWave &  theWave) const

Definition at line 209 of file LArWaveHelper.cxx.

210 {
211  const unsigned int N = theWave.getSize() ;
212  unsigned int imax= 1;
213  const unsigned Nbase=5;
214  double amax=-9.999E+99;
215  double base = getBaseline(theWave,Nbase);
216  for (unsigned int i=1; i<N ; i++) {
217  const double& asample=theWave.getSample(i);
218  if((asample-base)>amax){
219  imax = i;
220  amax = asample-base ;
221  }
222  }
223  if (amax == 0.) {
224  // the width is ill-defined in that case.
225  return 0.;
226  }
227 
228  double wleft=0.;
229  double wright=0.;
230  double width=0.;
231  double asamp1=0.;
232  double asamp2=0.;
233 
234  for (unsigned int i=1; i<imax; i++){
235  asamp1=theWave.getSample(i);
236  asamp2=theWave.getSample(i-1);
237  if((asamp1-base)>=amax/2. && (asamp2-base)<=amax/2.)
238  {
239  wleft= i-1+(amax/2.-(asamp2-base))/(asamp1-asamp2);
240 
241  }
242  }
243  for (unsigned int i=imax; i<N; i++){
244  asamp1=theWave.getSample(i);
245  asamp2=theWave.getSample(i-1);
246  if((asamp1-base)<=amax/2. && (asamp2-base)>=amax/2.) {
247  wright= i-1+((asamp2-base)-amax/2.)/(asamp2-asamp1);
248 
249  }
250  }
251  width=wright-wleft;
252 
253  return width;
254 }

◆ getZeroCross()

unsigned int LArWaveHelper::getZeroCross ( const LArWave &  theWave) const

Definition at line 425 of file LArWaveHelper.cxx.

426 {
427  for ( unsigned i=getMax(theWave)+1 ; i<theWave.getSize() ; i++)
428  if (theWave.getSample(i)<0) return i ;
429  return 0 ;
430 }

◆ linearMasterWave()

std::vector< LArWave > LArWaveHelper::linearMasterWave ( const std::vector< const LArWave * > &  vWaves,
const std::vector< double > &  vAmpli 
) const

Definition at line 498 of file LArWaveHelper.cxx.

500 {
501  std::vector<LArWave> MWandDAC0 ;
502  unsigned nWaves = vWaves.size() ;
503  if ( nWaves != vAmpli.size() ) {
504  std::cout << " linearMasterWave ERROR: wrong number of amplitudes!" << std::endl ;
505  return MWandDAC0 ;
506  }
507  double dt = vWaves[0]->getDt() ;
508  unsigned len = vWaves[0]->getSize() ;
509  for ( unsigned k=1 ; k<nWaves ; k++ ) {
510  if ( vWaves[k]->getDt() != dt ) {
511  std::cout << " linearMasterWave ERROR: waves don't have same dt" << std::endl ;
512  return MWandDAC0 ;
513  }
514  if ( vWaves[k]->getSize() != len ) {
515  std::cout << " linearMasterWave ERROR: waves don't have same length" << std::endl ;
516  return MWandDAC0 ;
517  }
518  }
519  LArWave mw(len,dt) ;
520  LArWave dac0(len,dt) ;
521  for ( unsigned i=0 ; i<len ; i++ ) {
522  std::vector<double> ramp ;
523  for ( unsigned k=0 ; k<nWaves ; k++ ) {
524  ramp.push_back(vWaves[k]->getSample(i)) ;
525  }
526  double rho ;
527  std::vector<double> par = linfit(vAmpli,ramp,rho) ;
528  if (par.size() != 2) {
529  std::cout << " linearMasterWave ERROR: linear fit failed" << std::endl ;
530  return MWandDAC0 ;
531  }
532  mw[i] = par[1] ;
533  dac0[i] = par[0] ;
534  }
535  MWandDAC0.push_back(dac0) ;
536  MWandDAC0.push_back(mw) ;
537  return MWandDAC0 ;
538 }

◆ linfit() [1/2]

std::vector< double > LArWaveHelper::linfit ( const LArWave &  theWave,
unsigned  iFirst,
unsigned  iLast,
double &  rho 
) const

Definition at line 471 of file LArWaveHelper.cxx.

473 {
474  std::vector<double> X,Y ;
475  get_fit_vectors(theWave,iFirst,iLast,X,Y) ;
476  return linfit(X,Y,rho) ;
477 }

◆ linfit() [2/2]

std::vector< double > LArWaveHelper::linfit ( const std::vector< double > &  X,
const std::vector< double > &  Y,
double &  rho 
) const
private

Definition at line 545 of file LArWaveHelper.cxx.

547 {
548  std::vector<double> C ;
549  C.resize(0) ;
550  rho = 0 ;
551  unsigned N = X.size() ;
552  if ( N<2 || Y.size()!=N ) {
553  // fprintf(stderr,"ERROR: lin. fit on less than 2 points\n") ;
554  } else {
555  double Sx=0. , Sy=0. , Sx2=0. , Sxy=0. , S=0. , Sy2=0. ;
556  for (unsigned i=0;i<N;i++) {
557  S += 1 ;
558  Sx += X[i] ;
559  Sy += Y[i] ;
560  Sx2 += X[i] * X[i] ;
561  Sxy += X[i] * Y[i] ;
562  Sy2 += Y[i] * Y[i] ;
563  }
564  C.resize(2) ;
565  C[1] = ( S*Sxy - Sx*Sy ) / ( S*Sx2 - Sx*Sx ) ;
566  C[0] = ( Sy*Sx2 - Sxy*Sx ) / ( S*Sx2 - Sx*Sx ) ;
567  if ( C[1] != C[1] || C[0] != C[0] ) { // protection agains C[i] being nan
568  C.resize(0) ;
569  return C ;
570  }
571  rho = ( S*Sxy - Sx*Sy ) / sqrt( (S*Sx2-Sx*Sx) * (S*Sy2-Sy*Sy) ) ;
572  }
573  return C ;
574 }

◆ polyfit() [1/2]

std::vector< double > LArWaveHelper::polyfit ( const LArWave &  theWave,
unsigned  iFirst,
unsigned  iLast,
unsigned  Ndeg 
) const

Definition at line 487 of file LArWaveHelper.cxx.

489 {
490  std::vector<double> X,Y ;
491  get_fit_vectors(theWave,iFirst,iLast,X,Y) ;
492  return polyfit(X,Y,Ndeg) ;
493 }

◆ polyfit() [2/2]

std::vector< double > LArWaveHelper::polyfit ( const std::vector< double > &  X,
const std::vector< double > &  Y,
unsigned  Ndeg 
) const
private

Definition at line 600 of file LArWaveHelper.cxx.

603 {
604  // to be implemented - return to avoid compiler warning RDS
605  std::vector<double> dummy;
606  if(Ndeg > 2 || (Y.empty()) || (X.size() != Y.size())) {
607  return dummy;
608  }
609  switch (Ndeg) {
610  case 0: { dummy.resize(1);
611  double sum = 0.;
612  for(unsigned int i=0; i<Y.size(); ++i) sum += Y[i];
613  dummy[0] = sum / Y.size();
614  return (dummy); }
615  case 1: { double rho; return linfit(X,Y,rho);}
616  case 2: { // fill the matrix XT*X
617  double a11,a12,a13,a21,a22,a23,a31,a32,a33;
618  a11=1.;
619  double sum =0., sum2=0., sum3=0., sum4=0.;
620  for(unsigned int i=0; i<X.size(); ++i) {
621  sum += X[i];
622  sum2 += X[i]*X[i];
623  sum3 += X[i]*X[i]*X[i];
624  sum4 += X[i]*X[i]*X[i]*X[i];
625  }
626  a12 = a21 = sum;
627  a13 = a22 = a31 = sum2;
628  a23 = a32 = sum3;
629  a33 = sum4;
630  // Inverse
631  double ai11,ai12,ai13,ai21,ai22,ai23,ai31,ai32,ai33;
632  const double det =
633  a11*(a33*a22-a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12-a22*a13);
634  if(fabs(det) < 1.e-6) return (dummy); // Not invertable
635  ai11 = (a33*a22 - a32*a23);
636  ai12 = ai21 = -(a33*a12-a32*a13);
637  ai22 = a33*a11 - a31*a13;
638  ai13 = ai31 = a23*a12 - a22*a13;
639  ai23 = ai32 = -(a23*a11-a21*a13);
640  ai33 = a22*a11-a21*a12;
641  sum = 0.; sum2=0.; sum3=0.;
642  for(unsigned int i=0; i<X.size(); ++i) {
643  sum += Y[i];
644  sum2 += X[i]*Y[i];
645  sum3 += X[i]*X[i]*Y[i];
646  }
647  dummy.resize(3);
648  const double inv_det = 1. / det;
649  dummy[0] = (ai11*sum + ai12*sum2 + ai13*sum3) * inv_det;
650  dummy[1] = (ai21*sum + ai22*sum2 + ai23*sum3) * inv_det;
651  dummy[2] = (ai31*sum + ai32*sum2 + ai33*sum3) * inv_det;
652 
653  return (dummy);
654  }
655  default : { // will not happen
656  return (dummy);
657  }
658  }
659 }

◆ subSample()

LArWave LArWaveHelper::subSample ( const LArWave &  theWave,
unsigned  Nfirst,
unsigned  deltaN 
) const

Definition at line 338 of file LArWaveHelper.cxx.

339 {
340  unsigned length = theWave.getSize() ;
341  if ( length <= 0 ) return LArWave() ;
342  unsigned Nsamp = (length-Nfirst + (deltaN-1)) / deltaN ;
343  LArWave w( Nsamp , deltaN*theWave.getDt() ) ;
344  for (unsigned i=0;i<Nsamp;i++) {
345  unsigned j = Nfirst + i * deltaN ;
346  w[i] = theWave.getSample(j) ;
347  }
348  return w ;
349 }

◆ translate()

LArWave LArWaveHelper::translate ( const LArWave &  theWave,
int  nShift,
double  baseline = 0. 
) const

Definition at line 11 of file LArWaveHelper.cxx.

13 {
14  unsigned int nsamples = theWave.getSize() ;
15  std::vector<double> trWave ;
16  trWave.resize( nsamples+nShift ) ;
17  for ( unsigned int i=0 ; i<theWave.getSize() ; i++ ) {
18  if ( (int)i + nShift >= 0 ) trWave[i+nShift] = theWave.getSample(i) ;
19  }
20  if ( nShift > 0 ) {
21  for ( unsigned int i=0 ; i<(unsigned int)nShift ; i++ ) {
22  trWave[i] = baseline ;
23  }
24  }
25  return LArWave( trWave , theWave.getDt() ) ;
26 }
The documentation for this class was generated from the following files:
LArWave
Definition: LArWave.h:31
base
std::string base
Definition: hcg.cxx:78
LArWaveHelper::getSumRegion
double getSumRegion(const LArWave &theWave, unsigned iFirst, unsigned iLast) const
Definition: LArWaveHelper.cxx:374
LArWave::getSize
size_t getSize() const
number of time samples
Definition: LArWave.h:62
LArWave::getFlag
unsigned getFlag() const
flag: ...
Definition: LArWave.h:178
TRTCalib_Extractor.det
det
Definition: TRTCalib_Extractor.py:36
LArWaveHelper::getWidth
double getWidth(const LArWave &theWave) const
Definition: LArWaveHelper.cxx:209
LArWave::getWave
const std::vector< double > & getWave() const
Wave parameters.
Definition: LArWave.h:167
LArWaveHelper::getBaseline
double getBaseline(const LArWave &theWave, unsigned nBase) const
Definition: LArWaveHelper.cxx:351
DMTest::C
C_v1 C
Definition: C.h:26
skel.it
it
Definition: skel.GENtoEVGEN.py:407
LArWave::getDt
const double & getDt() const
delta time
Definition: LArWave.h:50
get_generator_info.stderr
stderr
Definition: get_generator_info.py:40
LArWaveHelper::getJitter
double getJitter(const LArWaveCumul &theWave) const
Definition: LArWaveHelper.cxx:433
LArWaveHelper::getMaxAmp
double getMaxAmp(const LArWave &theWave) const
Definition: LArWaveHelper.cxx:129
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
LArWaveHelper::getMax
unsigned int getMax(const LArWave &theWave) const
return index of maximum sample
Definition: LArWaveHelper.cxx:89
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
JetTiledMap::S
@ S
Definition: TiledEtaPhiMap.h:44
LArWaveHelper::getDMax
double getDMax(const LArWave &theWave, double &tMax) const
return amplitude aproximation from poly2 fit around maxima, and it's time
Definition: LArWaveHelper.cxx:103
Monitored::X
@ X
Definition: HistogramFillerUtils.h:24
LArWaveHelper::get_fit_vectors
unsigned get_fit_vectors(const LArWave &theWave, unsigned iFirst, unsigned iLast, std::vector< double > &X, std::vector< double > &Y) const
Definition: LArWaveHelper.cxx:661
LArWave::getSample
const double & getSample(const unsigned int i) const
Amplitude per time bin.
Definition: LArWave.h:53
IBLCalibrationConfig.thr
thr
Definition: IBLCalibrationConfig.py:39
getSize
int getSize(std::map< std::string, std::vector< std::string >> &collection, const std::string &object)
Definition: SUSYToolsAlg.cxx:1577
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
lumiFormat.i
int i
Definition: lumiFormat.py:85
CaloNoise_fillDB.dt
dt
Definition: CaloNoise_fillDB.py:56
LArWaveHelper::getSumSquareRegion
double getSumSquareRegion(const LArWave &theWave, unsigned iFirst, unsigned iLast) const
Definition: LArWaveHelper.cxx:396
baseline
@ baseline
Definition: SUSYToolsTester.cxx:94
CalibDbCompareRT.dummy
dummy
Definition: CalibDbCompareRT.py:59
LArWaveHelper::fact
unsigned fact(unsigned N) const
Definition: LArWaveHelper.h:63
LArWaveDerivedQuantitiesP
Definition: LArWaveDerivedQuantitiesP.h:10
imax
int imax(int i, int j)
Definition: TileLaserTimingTool.cxx:33
LArWaveHelper::linfit
std::vector< double > linfit(const LArWave &theWave, unsigned iFirst, unsigned iLast, double &rho) const
Definition: LArWaveHelper.cxx:471
Monitored::Y
@ Y
Definition: HistogramFillerUtils.h:24
LArWaveHelper::polyfit
std::vector< double > polyfit(const LArWave &theWave, unsigned iFirst, unsigned iLast, unsigned Ndeg) const
Definition: LArWaveHelper.cxx:487
LArWaveHelper::getZeroCross
unsigned getZeroCross(const LArWave &theWave) const
Definition: LArWaveHelper.cxx:425
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:28
ReadOfcFromCool.nsamples
nsamples
Definition: ReadOfcFromCool.py:115
python.CaloAddPedShiftConfig.default
default
Definition: CaloAddPedShiftConfig.py:43
Threshold
Threshold
Definition: TRIGGERidentity.h:18
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
Base_Fragment.width
width
Definition: Sherpa_i/share/common/Base_Fragment.py:59
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
LArWaveHelper::derive
LArWave derive(const LArWave &theWave) const
crude derivative
Definition: LArWaveHelper.cxx:256
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:198
python.compressB64.c
def c
Definition: compressB64.py:93
length
double length(const pvec &v)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:26
LArWaveHelper::expfit
std::vector< double > expfit(const LArWave &theWave, unsigned iFirst, unsigned iLast, double &rho) const
Definition: LArWaveHelper.cxx:479
LArWaveHelper::getT0
double getT0(const LArWave &theWave) const
Definition: LArWaveHelper.cxx:162
LArWaveCumul::getErrors
const std::vector< double > & getErrors() const
error vector
Definition: LArWaveCumul.h:138
fitman.rho
rho
Definition: fitman.py:532
fitman.k
k
Definition: fitman.py:528
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
Generated on Mon Sep 29 2025 21:40:53 for ATLAS Offline Software by   doxygen 1.8.18