CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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 ( )

Definition at line 16 of file PulseFitWithShape.cc.

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

16  {
17  fNsamples = 0;
18  fNsamplesShape = 0;
19  fNum_samp_bef_max = 0;
20  fNum_samp_after_max = 0;
21  fNoise = 0.0;
22 }
PulseFitWithShape::~PulseFitWithShape ( )
override

Definition at line 25 of file PulseFitWithShape.cc.

25 {}

Member Function Documentation

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

Definition at line 58 of file PulseFitWithShape.cc.

References isotrackTrainRegressor::a1, isotrackTrainRegressor::a2, c, HLT_FULL_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, gpuVertexFinder::noise, and pshape.

Referenced by EcalLaserAnalyzer2::endJob().

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

Definition at line 28 of file PulseFitWithShape.h.

References fAmp_fitted_max.

Referenced by EcalLaserAnalyzer2::endJob().

28 { return fAmp_fitted_max; }
double PulseFitWithShape::getTime ( )
inline

Definition at line 29 of file PulseFitWithShape.h.

References fTim_fitted_max.

Referenced by EcalLaserAnalyzer2::endJob().

29 { return fTim_fitted_max; }
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 29 of file PulseFitWithShape.cc.

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().

35  {
36  fNsamples = n_samples;
37  fNsamplesShape = n_samplesShape;
38  fNb_iter = niter;
39  fNum_samp_bef_max = samplb;
40  fNum_samp_after_max = sampla;
41 
42  if (fNsamplesShape < fNum_samp_bef_max + fNum_samp_after_max + 1) {
43  std::cout << "PulseFitWithShape::init: Error! Configuration samples in fit greater than total number of samples!"
44  << std::endl;
45  }
46 
47  for (int i = 0; i < fNsamplesShape; i++) {
48  pshape.push_back(shape[i]);
49  dshape.push_back(0.0);
50  }
51 
52  fNoise = nois;
53  return;
54 }
mps_fire.i
i
Definition: mps_fire.py:428
PulseFitWithShape::dshape
std::vector< double > dshape
Definition: PulseFitWithShape.h:37
PulseFitWithShape::pshape
std::vector< double > pshape
Definition: PulseFitWithShape.h:36
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144

Member Data Documentation

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

Definition at line 37 of file PulseFitWithShape.h.

Referenced by doFit(), and init().

double PulseFitWithShape::fAmp_fitted_max

Definition at line 25 of file PulseFitWithShape.h.

Referenced by doFit(), and getAmpl().

int PulseFitWithShape::fNb_iter
private

Definition at line 39 of file PulseFitWithShape.h.

Referenced by doFit(), and init().

double PulseFitWithShape::fNoise
private

Definition at line 34 of file PulseFitWithShape.h.

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

int PulseFitWithShape::fNsamples
private

Definition at line 32 of file PulseFitWithShape.h.

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

int PulseFitWithShape::fNsamplesShape
private

Definition at line 33 of file PulseFitWithShape.h.

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

int PulseFitWithShape::fNum_samp_after_max
private

Definition at line 41 of file PulseFitWithShape.h.

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

int PulseFitWithShape::fNum_samp_bef_max
private

Definition at line 40 of file PulseFitWithShape.h.

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

double PulseFitWithShape::fTim_fitted_max

Definition at line 26 of file PulseFitWithShape.h.

Referenced by doFit(), and getTime().

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

Definition at line 36 of file PulseFitWithShape.h.

Referenced by doFit(), and init().

Generated for CMSSW Reference Manual by doxygen 1.8.5