CMS 3D CMS Logo

List of all members | Public Member Functions | Public Attributes | Private Attributes
PulseFitWithShape Class Reference

#include <PulseFitWithShape.h>

Inheritance diagram for PulseFitWithShape:

Public Member Functions

virtual double doFit (double *, double *cova=nullptr)
 
double getAmpl ()
 
double getTime ()
 
virtual void init (int, int, int, int, int, const std::vector< double > &, double)
 
 PulseFitWithShape ()
 
 ~PulseFitWithShape () override
 

Public Attributes

double fAmp_fitted_max
 
double fTim_fitted_max
 

Private Attributes

std::vector< double > dshape
 
int fNb_iter
 
double fNoise
 
int fNsamples
 
int fNsamplesShape
 
int fNum_samp_after_max
 
int fNum_samp_bef_max
 
std::vector< double > pshape
 

Detailed Description

Definition at line 10 of file PulseFitWithShape.h.

Constructor & Destructor Documentation

◆ PulseFitWithShape()

PulseFitWithShape::PulseFitWithShape ( )

Definition at line 18 of file PulseFitWithShape.cc.

18  {
19  fNsamples = 0;
20  fNsamplesShape = 0;
23  fNoise = 0.0;
24 }

References fNoise, fNsamples, fNsamplesShape, fNum_samp_after_max, and fNum_samp_bef_max.

◆ ~PulseFitWithShape()

PulseFitWithShape::~PulseFitWithShape ( )
override

Definition at line 27 of file PulseFitWithShape.cc.

27 {}

Member Function Documentation

◆ doFit()

double PulseFitWithShape::doFit ( double *  adc,
double *  cova = nullptr 
)
virtual

Definition at line 60 of file PulseFitWithShape.cc.

60  {
61  // xpar = fit paramaters
62  // [0] = signal amplitude
63  // [1] = residual pedestal
64  // [2] = clock phase
65 
66  bool useCova = true;
67  if (cova == nullptr)
68  useCova = false;
69 
70  double xpar[3];
71  double chi2;
72 
73  fAmp_fitted_max = 0.;
74  fTim_fitted_max = 0.;
75 
76  // for now don't fit pedestal
77 
78  xpar[1] = 0.0;
79 
80  // Sample noise. If the cova matrix is defined, use it :
81 
82  double noise = fNoise;
83  //if(cova[0] > 0.) noise=1./sqrt(cova[0]);
84 
85  xpar[0] = 0.;
86  xpar[2] = 0.;
87 
88  // first locate max:
89 
90  int imax = 0;
91  double amax = 0.0;
92  for (int i = 0; i < fNsamples; i++) {
93  if (adc[i] > amax) {
94  amax = adc[i];
95  imax = i;
96  }
97  }
98 
99  // Shift pulse shape and calculate its derivative:
100 
101  double qms = 0.;
102  int ims = 0;
103 
104  // 1) search for maximum
105 
106  for (int is = 0; is < fNsamplesShape; is++) {
107  if (pshape[is] > qms) {
108  qms = pshape[is];
109  ims = is;
110  }
111 
112  // 2) compute shape derivative :
113 
114  if (is < fNsamplesShape - 2)
115  dshape[is] = (pshape[is + 2] - pshape[is]) * 12.5;
116  else
117  dshape[is] = dshape[is - 1];
118  }
119 
120  // 3) compute pol2 max
121 
122  double sq1 = pshape[ims - 1];
123  double sq2 = pshape[ims];
124  double sq3 = pshape[ims + 1];
125 
126  double a2 = (sq3 + sq1) / 2.0 - sq2;
127  double a1 = sq2 - sq1 + a2 * (1 - 2 * ims);
128 
129  double t_ims = 0;
130  if (a2 != 0)
131  t_ims = -a1 / (2.0 * a2);
132 
133  // Starting point of the fit : qmax and tmax given by a
134  // P2 fit on 3 max samples.
135 
136  double qm = 0.;
137  int im = 0;
138 
139  int nsamplef = fNum_samp_bef_max + fNum_samp_after_max + 1; // number of samples used in the fit
140  int nsampleo = imax - fNum_samp_bef_max; // first sample number = sample max-fNum_samp_bef_max
141 
142  for (int is = 0; is < nsamplef; is++) {
143  if (adc[is + nsampleo] > qm) {
144  qm = adc[is + nsampleo];
145  im = nsampleo + is;
146  }
147  }
148 
149  double tm;
150  if (qm > 5. * noise) {
151  if (im >= nsamplef + nsampleo)
152  im = nsampleo + nsamplef - 1;
153  double q1 = adc[im - 1];
154  double q2 = adc[im];
155  double q3 = adc[im + 1];
156  double y2 = (q1 + q3) / 2. - q2;
157  double y1 = q2 - q1 + y2 * (1 - 2 * im);
158  double y0 = q2 - y1 * (double)im - y2 * (double)(im * im);
159  tm = -y1 / 2. / y2;
160  xpar[0] = y0 + y1 * tm + y2 * tm * tm;
161  xpar[2] = (double)ims / 25. - tm;
162  }
163 
164  double chi2old = 999999.;
165  chi2 = 99999.;
166  int nsfit = nsamplef;
167  int iloop = 0;
168  int nloop = fNb_iter;
169  if (qm <= 5 * noise)
170  nloop = 1; // Just one iteration for very low signal
171 
172  std::vector<double> resi(fNsamples, 0.0);
173 
174  while (std::fabs(chi2old - chi2) > 0.1 && iloop < nloop) {
175  iloop++;
176  chi2old = chi2;
177 
178  double c = 0.;
179  double y1 = 0.;
180  double s1 = 0.;
181  double s2 = 0.;
182  double ys1 = 0.;
183  double sp1 = 0.;
184  double sp2 = 0.;
185  double ssp = 0.;
186  double ysp = 0.;
187 
188  for (int is = 0; is < nsfit; is++) {
189  int iis = is;
190  iis = is + nsampleo;
191 
192  double t1 = (double)iis + xpar[2];
193  double xbin = t1 * 25.;
194  int ibin1 = (int)xbin;
195 
196  if (ibin1 < 0)
197  ibin1 = 0;
198 
199  double amp1, amp11, amp12, der1, der11, der12;
200 
201  if (ibin1 <= fNsamplesShape - 2) { // Interpolate shape values to get the right number :
202 
203  int ibin2 = ibin1 + 1;
204  double xfrac = xbin - ibin1;
205  amp11 = pshape[ibin1];
206  amp12 = pshape[ibin2];
207  amp1 = (1. - xfrac) * amp11 + xfrac * amp12;
208  der11 = dshape[ibin1];
209  der12 = dshape[ibin2];
210  der1 = (1. - xfrac) * der11 + xfrac * der12;
211 
212  } else { // Or extraoplate if we are outside the array :
213 
214  amp1 = pshape[fNsamplesShape - 1] + dshape[fNsamplesShape - 1] * (xbin - double(fNsamplesShape - 1)) / 25.;
215  der1 = dshape[fNsamplesShape - 1];
216  }
217 
218  if (useCova) { // Use covariance matrix:
219  for (int js = 0; js < nsfit; js++) {
220  int jjs = js;
221  jjs += nsampleo;
222 
223  t1 = (double)jjs + xpar[2];
224  xbin = t1 * 25.;
225  ibin1 = (int)xbin;
226  if (ibin1 < 0)
227  ibin1 = 0;
228  if (ibin1 > fNsamplesShape - 2)
229  ibin1 = fNsamplesShape - 2;
230  int ibin2 = ibin1 + 1;
231  double xfrac = xbin - ibin1;
232  amp11 = pshape[ibin1];
233  amp12 = pshape[ibin2];
234  double amp2 = (1. - xfrac) * amp11 + xfrac * amp12;
235  der11 = dshape[ibin1];
236  der12 = dshape[ibin2];
237  double der2 = (1. - xfrac) * der11 + xfrac * der12;
238  c = c + cova[js * fNsamples + is];
239  y1 = y1 + adc[iis] * cova[js * fNsamples + is];
240  s1 = s1 + amp1 * cova[js * fNsamples + is];
241  s2 = s2 + amp1 * amp2 * cova[js * fNsamples + is];
242  ys1 = ys1 + adc[iis] * amp2 * cova[js * fNsamples + is];
243  sp1 = sp1 + der1 * cova[is * fNsamples + js];
244  sp2 = sp2 + der1 * der2 * cova[js * fNsamples + is];
245  ssp = ssp + amp1 * der2 * cova[js * fNsamples + is];
246  ysp = ysp + adc[iis] * der2 * cova[js * fNsamples + is];
247  }
248  } else { // Don't use covariance matrix:
249  c++;
250  y1 = y1 + adc[iis];
251  s1 = s1 + amp1;
252  s2 = s2 + amp1 * amp1;
253  ys1 = ys1 + amp1 * adc[iis];
254  sp1 = sp1 + der1;
255  sp2 = sp2 + der1 * der1;
256  ssp = ssp + der1 * amp1;
257  ysp = ysp + der1 * adc[iis];
258  }
259  }
260  xpar[0] = (ysp * ssp - ys1 * sp2) / (ssp * ssp - s2 * sp2);
261  xpar[2] += (ysp / xpar[0] / sp2 - ssp / sp2);
262 
263  for (int is = 0; is < nsfit; is++) {
264  int iis = is;
265  iis = is + nsampleo;
266 
267  double t1 = (double)iis + xpar[2];
268  double xbin = t1 * 25.;
269  int ibin1 = (int)xbin;
270  if (ibin1 < 0)
271  ibin1 = 0;
272 
273  if (ibin1 < 0)
274  ibin1 = 0;
275  if (ibin1 > fNsamplesShape - 2)
276  ibin1 = fNsamplesShape - 2;
277  int ibin2 = ibin1 + 1;
278  double xfrac = xbin - ibin1;
279  double amp11 = xpar[0] * pshape[ibin1];
280  double amp12 = xpar[0] * pshape[ibin2];
281  double amp1 = xpar[1] + (1. - xfrac) * amp11 + xfrac * amp12;
282  resi[iis] = adc[iis] - amp1;
283  }
284 
285  chi2 = 0.;
286  for (int is = 0; is < nsfit; is++) {
287  int iis = is;
288  iis += nsampleo;
289 
290  if (useCova) {
291  for (int js = 0; js < nsfit; js++) {
292  int jjs = js;
293  jjs += nsampleo;
294  chi2 += resi[iis] * resi[jjs] * cova[js * fNsamples + is];
295  }
296 
297  } else
298  chi2 += resi[iis] * resi[iis];
299  }
300  }
301 
302  fAmp_fitted_max = xpar[0];
303  fTim_fitted_max = (double)t_ims / 25. - xpar[2];
304 
305  return chi2;
306 }

References testProducerWithPsetDescEmpty_cfi::a2, ecalMGPA::adc(), HltBtagPostValidation_cff::c, hltPixelTracks_cff::chi2, dshape, fAmp_fitted_max, fNb_iter, fNoise, fNsamples, fNsamplesShape, fNum_samp_after_max, fNum_samp_bef_max, fTim_fitted_max, mps_fire::i, createfilelist::int, hgchebackDigitizer_cfi::noise, pshape, q1, q2, indexGen::s2, RandomServiceHelper::t1, testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

Referenced by EcalLaserAnalyzer2::endJob().

◆ getAmpl()

double PulseFitWithShape::getAmpl ( )
inline

Definition at line 28 of file PulseFitWithShape.h.

28 { return fAmp_fitted_max; }

References fAmp_fitted_max.

Referenced by EcalLaserAnalyzer2::endJob().

◆ getTime()

double PulseFitWithShape::getTime ( )
inline

Definition at line 29 of file PulseFitWithShape.h.

29 { return fTim_fitted_max; }

References fTim_fitted_max.

Referenced by EcalLaserAnalyzer2::endJob().

◆ init()

void PulseFitWithShape::init ( int  n_samples,
int  samplb,
int  sampla,
int  niter,
int  n_samplesShape,
const std::vector< double > &  shape,
double  nois 
)
virtual

Definition at line 31 of file PulseFitWithShape.cc.

37  {
38  fNsamples = n_samples;
39  fNsamplesShape = n_samplesShape;
40  fNb_iter = niter;
41  fNum_samp_bef_max = samplb;
42  fNum_samp_after_max = sampla;
43 
45  std::cout << "PulseFitWithShape::init: Error! Configuration samples in fit greater than total number of samples!"
46  << std::endl;
47  }
48 
49  for (int i = 0; i < fNsamplesShape; i++) {
50  pshape.push_back(shape[i]);
51  dshape.push_back(0.0);
52  }
53 
54  fNoise = nois;
55  return;
56 }

References gather_cfg::cout, dshape, fNb_iter, fNoise, fNsamples, fNsamplesShape, fNum_samp_after_max, fNum_samp_bef_max, mps_fire::i, and pshape.

Referenced by EcalLaserAnalyzer2::endJob().

Member Data Documentation

◆ dshape

std::vector<double> PulseFitWithShape::dshape
private

Definition at line 37 of file PulseFitWithShape.h.

Referenced by doFit(), and init().

◆ fAmp_fitted_max

double PulseFitWithShape::fAmp_fitted_max

Definition at line 25 of file PulseFitWithShape.h.

Referenced by doFit(), and getAmpl().

◆ fNb_iter

int PulseFitWithShape::fNb_iter
private

Definition at line 39 of file PulseFitWithShape.h.

Referenced by doFit(), and init().

◆ fNoise

double PulseFitWithShape::fNoise
private

Definition at line 34 of file PulseFitWithShape.h.

Referenced by doFit(), init(), and PulseFitWithShape().

◆ fNsamples

int PulseFitWithShape::fNsamples
private

Definition at line 32 of file PulseFitWithShape.h.

Referenced by doFit(), init(), and PulseFitWithShape().

◆ fNsamplesShape

int PulseFitWithShape::fNsamplesShape
private

Definition at line 33 of file PulseFitWithShape.h.

Referenced by doFit(), init(), and PulseFitWithShape().

◆ fNum_samp_after_max

int PulseFitWithShape::fNum_samp_after_max
private

Definition at line 41 of file PulseFitWithShape.h.

Referenced by doFit(), init(), and PulseFitWithShape().

◆ fNum_samp_bef_max

int PulseFitWithShape::fNum_samp_bef_max
private

Definition at line 40 of file PulseFitWithShape.h.

Referenced by doFit(), init(), and PulseFitWithShape().

◆ fTim_fitted_max

double PulseFitWithShape::fTim_fitted_max

Definition at line 26 of file PulseFitWithShape.h.

Referenced by doFit(), and getTime().

◆ pshape

std::vector<double> PulseFitWithShape::pshape
private

Definition at line 36 of file PulseFitWithShape.h.

Referenced by doFit(), and init().

ecalMGPA::adc
constexpr int adc(sample_type sample)
get the ADC sample (12 bits)
Definition: EcalMGPASample.h:11
mps_fire.i
i
Definition: mps_fire.py:355
PulseFitWithShape::fNum_samp_bef_max
int fNum_samp_bef_max
Definition: PulseFitWithShape.h:40
gather_cfg.cout
cout
Definition: gather_cfg.py:144
indexGen.s2
s2
Definition: indexGen.py:107
PulseFitWithShape::fTim_fitted_max
double fTim_fitted_max
Definition: PulseFitWithShape.h:26
hltPixelTracks_cff.chi2
chi2
Definition: hltPixelTracks_cff.py:25
PulseFitWithShape::fAmp_fitted_max
double fAmp_fitted_max
Definition: PulseFitWithShape.h:25
PulseFitWithShape::fNum_samp_after_max
int fNum_samp_after_max
Definition: PulseFitWithShape.h:41
testProducerWithPsetDescEmpty_cfi.a2
a2
Definition: testProducerWithPsetDescEmpty_cfi.py:35
PulseFitWithShape::fNsamples
int fNsamples
Definition: PulseFitWithShape.h:32
testProducerWithPsetDescEmpty_cfi.y1
y1
Definition: testProducerWithPsetDescEmpty_cfi.py:29
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
PulseFitWithShape::dshape
std::vector< double > dshape
Definition: PulseFitWithShape.h:37
q2
double q2[4]
Definition: TauolaWrapper.h:88
testProducerWithPsetDescEmpty_cfi.y2
y2
Definition: testProducerWithPsetDescEmpty_cfi.py:30
q1
double q1[4]
Definition: TauolaWrapper.h:87
PulseFitWithShape::fNsamplesShape
int fNsamplesShape
Definition: PulseFitWithShape.h:33
PulseFitWithShape::fNb_iter
int fNb_iter
Definition: PulseFitWithShape.h:39
createfilelist.int
int
Definition: createfilelist.py:10
PulseFitWithShape::fNoise
double fNoise
Definition: PulseFitWithShape.h:34
PulseFitWithShape::pshape
std::vector< double > pshape
Definition: PulseFitWithShape.h:36
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
hgchebackDigitizer_cfi.noise
noise
Definition: hgchebackDigitizer_cfi.py:42