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
PulseFitWithShape Class Reference

#include <PulseFitWithShape.h>

Inheritance diagram for PulseFitWithShape:

Public Member Functions

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

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

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

21 {
22 
23  fNsamples = 0;
24  fNsamplesShape = 0;
25  fNum_samp_bef_max = 0;
26  fNum_samp_after_max = 0;
27  fNoise = 0.0;
28 }
PulseFitWithShape::~PulseFitWithShape ( )
virtual

Definition at line 31 of file PulseFitWithShape.cc.

32 {
33 }

Member Function Documentation

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

Definition at line 62 of file PulseFitWithShape.cc.

References EnergyCorrector::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, i, j, HLT_FULL_cff::noise, pshape, q1, q2, and indexGen::s2.

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

Definition at line 29 of file PulseFitWithShape.h.

References fAmp_fitted_max.

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

Definition at line 30 of file PulseFitWithShape.h.

References fTim_fitted_max.

30 { 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 37 of file PulseFitWithShape.cc.

References gather_cfg::cout, dshape, fNb_iter, fNoise, fNsamples, fNsamplesShape, fNum_samp_after_max, fNum_samp_bef_max, i, pshape, and reco::return().

38 {
39 
40  fNsamples = n_samples ;
41  fNsamplesShape = n_samplesShape ;
42  fNb_iter = niter ;
43  fNum_samp_bef_max = samplb ;
44  fNum_samp_after_max = sampla ;
45 
46 
47  if( fNsamplesShape < fNum_samp_bef_max+fNum_samp_after_max+1){
48  std::cout<<"PulseFitWithShape::init: Error! Configuration samples in fit greater than total number of samples!" << std::endl;
49  }
50 
51  for(int i=0;i<fNsamplesShape;i++){
52  pshape.push_back(shape[i]);
53  dshape.push_back(0.0);
54  }
55 
56  fNoise=nois;
57  return ;
58  }
PulseFitWithShape::pshape
std::vector< double > pshape
Definition: PulseFitWithShape.h:39
i
int i
Definition: DBlmapReader.cc:9
PulseFitWithShape::dshape
std::vector< double > dshape
Definition: PulseFitWithShape.h:40
reco::return
return(e1-e2)*(e1-e2)+dp *dp
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:145

Member Data Documentation

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

Definition at line 40 of file PulseFitWithShape.h.

Referenced by doFit(), and init().

double PulseFitWithShape::fAmp_fitted_max

Definition at line 26 of file PulseFitWithShape.h.

Referenced by doFit(), and getAmpl().

int PulseFitWithShape::fNb_iter
private

Definition at line 42 of file PulseFitWithShape.h.

Referenced by doFit(), and init().

double PulseFitWithShape::fNoise
private

Definition at line 37 of file PulseFitWithShape.h.

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

int PulseFitWithShape::fNsamples
private

Definition at line 35 of file PulseFitWithShape.h.

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

int PulseFitWithShape::fNsamplesShape
private

Definition at line 36 of file PulseFitWithShape.h.

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

int PulseFitWithShape::fNum_samp_after_max
private

Definition at line 44 of file PulseFitWithShape.h.

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

int PulseFitWithShape::fNum_samp_bef_max
private

Definition at line 43 of file PulseFitWithShape.h.

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

double PulseFitWithShape::fTim_fitted_max

Definition at line 27 of file PulseFitWithShape.h.

Referenced by doFit(), and getTime().

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

Definition at line 39 of file PulseFitWithShape.h.

Referenced by doFit(), and init().

Generated for CMSSW Reference Manual by doxygen 1.8.5