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 318 of file LArWaveHelper.cxx.

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

◆ derive()

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

crude derivative

Definition at line 252 of file LArWaveHelper.cxx.

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

◆ derive_smooth()

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

smoothed derivative

Definition at line 270 of file LArWaveHelper.cxx.

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

◆ 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 475 of file LArWaveHelper.cxx.

477 {
478  std::vector<double> X,Y ;
479  get_fit_vectors(theWave,iFirst,iLast,X,Y) ;
480  return expfit(X,Y,rho) ;
481 }

◆ 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 572 of file LArWaveHelper.cxx.

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

◆ 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 657 of file LArWaveHelper.cxx.

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

◆ 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 347 of file LArWaveHelper.cxx.

348 {
349  if (nBase == 0) return 0;
350  double base = 0.0 ;
351  for ( unsigned i=0 ; i<nBase ; i++ ) base += theWave.getSample(i) ;
352  return base/nBase ;
353 }

◆ getDerivedQuantities()

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

Definition at line 452 of file LArWaveHelper.cxx.

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

◆ 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 429 of file LArWaveHelper.cxx.

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

◆ 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 409 of file LArWaveHelper.cxx.

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

◆ getSumRegion()

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

Definition at line 370 of file LArWaveHelper.cxx.

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

◆ getSumSquareRegion()

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

Definition at line 392 of file LArWaveHelper.cxx.

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

◆ getSumSquareTail()

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

Definition at line 387 of file LArWaveHelper.cxx.

388 {
389  return getSumSquareRegion(theWave,iFirst,theWave.getSize());
390 }

◆ getSumTail()

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

Definition at line 355 of file LArWaveHelper.cxx.

356 {
357  return getSumRegion(theWave,iFirst,theWave.getSize());
358 }

◆ 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 
224  double wleft=0.;
225  double wright=0.;
226  double width=0.;
227  double asamp1=0.;
228  double asamp2=0.;
229 
230  for (unsigned int i=1; i<imax; i++){
231  asamp1=theWave.getSample(i);
232  asamp2=theWave.getSample(i-1);
233  if((asamp1-base)>=amax/2. && (asamp2-base)<=amax/2.)
234  {
235  wleft= i-1+(amax/2.-(asamp2-base))/(asamp1-asamp2);
236 
237  }
238  }
239  for (unsigned int i=imax; i<N; i++){
240  asamp1=theWave.getSample(i);
241  asamp2=theWave.getSample(i-1);
242  if((asamp1-base)<=amax/2. && (asamp2-base)>=amax/2.) {
243  wright= i-1+((asamp2-base)-amax/2.)/(asamp2-asamp1);
244 
245  }
246  }
247  width=wright-wleft;
248 
249  return width;
250 }

◆ getZeroCross()

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

Definition at line 421 of file LArWaveHelper.cxx.

422 {
423  for ( unsigned i=getMax(theWave)+1 ; i<theWave.getSize() ; i++)
424  if (theWave.getSample(i)<0) return i ;
425  return 0 ;
426 }

◆ linearMasterWave()

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

Definition at line 494 of file LArWaveHelper.cxx.

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

◆ linfit() [1/2]

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

Definition at line 467 of file LArWaveHelper.cxx.

469 {
470  std::vector<double> X,Y ;
471  get_fit_vectors(theWave,iFirst,iLast,X,Y) ;
472  return linfit(X,Y,rho) ;
473 }

◆ 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 541 of file LArWaveHelper.cxx.

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

◆ polyfit() [1/2]

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

Definition at line 483 of file LArWaveHelper.cxx.

485 {
486  std::vector<double> X,Y ;
487  get_fit_vectors(theWave,iFirst,iLast,X,Y) ;
488  return polyfit(X,Y,Ndeg) ;
489 }

◆ 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 596 of file LArWaveHelper.cxx.

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

◆ subSample()

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

Definition at line 334 of file LArWaveHelper.cxx.

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

◆ 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:370
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
LArWave::getSize
size_t getSize() const
number of time samples
Definition: LArWave.h:62
LArWave::getFlag
unsigned getFlag() const
flag: ...
Definition: LArWave.h:178
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
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:347
DMTest::C
C_v1 C
Definition: C.h:26
skel.it
it
Definition: skel.GENtoEVGEN.py:396
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:429
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:657
LArWave::getSample
const double & getSample(const unsigned int i) const
Amplitude per time bin.
Definition: LArWave.h:53
getSize
int getSize(std::map< std::string, std::vector< std::string >> &collection, const std::string &object)
Definition: SUSYToolsAlg.cxx:1578
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
lumiFormat.i
int i
Definition: lumiFormat.py:85
CaloNoise_fillDB.dt
dt
Definition: CaloNoise_fillDB.py:58
LArWaveHelper::getSumSquareRegion
double getSumSquareRegion(const LArWave &theWave, unsigned iFirst, unsigned iLast) const
Definition: LArWaveHelper.cxx:392
baseline
@ baseline
Definition: SUSYToolsTester.cxx:99
python.xAODType.dummy
dummy
Definition: xAODType.py:4
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:467
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:483
LArWaveHelper::getZeroCross
unsigned getZeroCross(const LArWave &theWave) const
Definition: LArWaveHelper.cxx:421
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
ReadOfcFromCool.nsamples
nsamples
Definition: ReadOfcFromCool.py:115
Threshold
Threshold
Definition: TRIGGERidentity.h:18
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:252
python.CaloScaleNoiseConfig.default
default
Definition: CaloScaleNoiseConfig.py:79
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
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:475
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
Generated on Sun Dec 22 2024 21:38:58 for ATLAS Offline Software by   doxygen 1.8.18