Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
matchedfilter.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file matchedfilter.C
1 
2 // ROOT macro for simulating the signal of different shapes, superimposed with the
3 // gaussian noise and to process resulting sequence using matching filter,
4 // extracting amplitude and time of the signal.
5 //
6 // Usage:
7 // To load the tools
8 // .L matchedfilter.C
9 // To generate signal with unipolar shape with 16 samples and calculate corrections
10 // init(16,0)
11 // The filter coefficients can be manually modified using, for example:
12 // g[2] = 0.5; // to set filter coeffcient 2 to 0.5
13 // To generate signal and noise,
14 // for example signal positioned at sample 50 with amplitude 100 and rms noise 10:
15 // series(50.,10.,100);
16 //
17 #define PRECISION_COEFF 1./64. // Precision of the filter coefficients. 1/64 is good enough
18 //#define PRECISION_SAMPLING 1.e-9 // High precision of the ADC sampling
19 //#define PRECISION_SAMPLING 1./10. // 14-bit precision of the ADC sampling
20 
21 #define SMAL_UNDEF 1.0e-99
22 #define ALevel 100.0 // Amplitude level
23 
24 #define NNX 400
25 Int_t NX = NNX;
26 Double_t x[NNX],y[NNX],yma[NNX],dy[NNX],dyma[NNX],s[NNX],v[NNX],ni[NNX];
27 //Double_t y_s[NNX];
28 Double_t gmax, gmaxx;
29 Double_t tau =0.;
30 Double_t ymax =0.;
31 Int_t quiet = 0;
32 #define NN 128;
33 Int_t N ; // filter length, set during initialization
34 Double_t fN; // widely used (Double_t)N;
35 Double_t g[NN], gma[NN];
36 Double_t norm, norm_ma;
37 Double_t gpos, gnoise;
38 Int_t npars = 2;
39 Double_t pars[] = {1.,0.,0,0};
40 #define grf_NPTS 15
41 #define grf_NY 8
42 Double_t serY[grf_NY];
43 Double_t grf_x[grf_NPTS] = {2000.,1000.,500.,200.,80.,60.,40.,30.,20.,10.,7.,5.,4.,2.5,2.0};
44 Double_t grf_y[grf_NY][grf_NPTS];
45 Double_t grf_ey[grf_NY][grf_NPTS];
46 Double_t grf_ex[grf_NPTS];
47 Char_t * grf_name[grf_NY] =
48  {
49  "T/SigmaT, T/ErrorT",
50  "A/SigmaA, A/ErrorA",
51  "Mean(Time)",
52  "Mean(A)",
53  "1/(Systematic Error) of Time",
54  "1/(Systematic Error) of Amplitude",
55  "T/StDev(Time)",
56  "A/StDev(A)"
57  };
58 TGraph * gr = 0;
59 TF1 * gf;
60 
61 #define NCorr 256 //number of correction points
62 Double_t serCorr[2][NCorr];
63 Double_t aCorr[NCorr],tCorr[NCorr];
64 TCanvas * gCnv=NULL;
65 
66 //shapes
67 Int_t SigShape;
68 #define UNIPOLAR 0
69 #define BIPOLAR 1
70 #define RECTANGULAR 2
71 #define GAUSS 3
72 #define STRIANGLE 4 //symmetric triangle
73 #define RTRIANGLE 5 //rigth triangle
74 #define UTRIANGLE 6
75 #define UNKNOWN_SHAPE 7
76 
77 #define NEED_FILTER_NOISE // define when you are interested in the noise output of the filter
78 
79 //#define TRY_WEIGHTED // to try weighted average finding of the max of the filtered signal
80 #ifdef TRY_WEIGHTED
81 Double_t tauw=0.;
82 #endif
83 
84 // **********************************************************************************************
85 
86 TCanvas *gCanv[5];
87 void init(UChar_t nsamples=16, UChar_t sigshape=0)
88 {
89 // This macro loads necessary functions and initialize the system.
90 // To simulate the series of measurements:
91 // root[] doSeries(pos, ntry);
92 // where the 'pos' is the position of the pulse in sample numbers
93 // the ntry is number of tries in each series
94  Int_t ii;
95  TCanvas *canv;
96 
97  gCanv[0] = new TCanvas();
98  gCanv[1] = new TCanvas();
99  gCnv=gCanv[1]; //standard canvas for show
100  gCanv[2] = new TCanvas();
101  gCanv[3] = new TCanvas();
102  gCanv[4] = new TCanvas();
103  gCanv[5] = new TCanvas();
104 
105  // consistency check
106  if(sigshape >= UNKNOWN_SHAPE)
107  {cout<<"Unknown shape, >= "<<UNKNOWN_SHAPE<<"\n"; exit(1);}
108  SigShape = (Int_t) sigshape;
109  if(nsamples > NN)
110  {cout<<"The filter length too long, >= "<<NN<<"\n"; exit(1);}
111  N = (Int_t)nsamples;
112  fN = (Double_t)N;
113 
114  cout<<"Preparing simulation of ";
115  switch (SigShape)
116  {
117  case UNIPOLAR: cout<<"UNIPOLAR"; break;
118  case BIPOLAR: cout<<"BIPOLAR"; break;
119  case RECTANGULAR: cout<<"RECTANGULAR"; break;
120  case GAUSS: cout<<"GAUSS"; break;
121  }
122  cout<<" shape["<<N<<"]\n";
123 
124  for(ii=0;ii<NX;ii++) ni[ii] = (Double_t) ii;
125 
126  getCoeffs();
127 
128  cout<<"Round coefficients to precision "<<PRECISION_COEFF<<endl;
129  roundCoeffs(PRECISION_COEFF);
130 
131  getCorrections();
132 
133  showShape();
134 
135  showAutocorrelation();
136 
137  #ifdef NEED_FILTER_NOISE
138  cout<<"Attention, NEED_FILTER_NOISE is active, it slows down the calculations\n";
139  #endif
140  cout<<"The filter coefficients g[i] can be modified now."<<endl;
141  cout<<"If they have been modified, execute: showShape(); getCorrections();"<<endl;
142  cout<<"To simulate, execute: doSeries(pos,ntry)\n";
143 }
144 
145 // **********************************************************************************************
146 
147 Double_t shape(Double_t xx)
148 {
149  return fshape(&xx,pars);
150 }
151 
152 // ----------------------------------------------------------------------------------------------
153 Double_t fshape(Double_t *xx, Double_t *par)
154 {
155  //Double_t r;
156  Double_t x1 = xx[0] - par[1],v1,v;
157  if(x1 < 0.) x1 = 0.;
158  switch (SigShape)
159  {
160  case UNIPOLAR:
161  // Landau-like
162  return par[0]*pow(x1,4.)*(exp(-x1*16./fN)); // fN is set to (float)NSAMPLES
163 0 case BIPOLAR:
164  // bipolar shape, simply derivative of the UNIPOLAR
165  return par[0]*pow(x1,3.)*exp(-x1*16./fN) * (4.-16./fN*x1);
166  case RECTANGULAR:// symmetrical trapezoidal
167  v1 = 1.;
168  if(x1<v1) v = par[0]*x1/v1;
169  else if(x1<=fN-v1) v = par[0];
170  else v = par[0]*(fN-x1)/v1;
171  return v;
172  case GAUSS:
173  x1 -= fN/2.;
174  return par[0]*exp(-(x1*x1)/1.);
175  case STRIANGLE: // symmetrical triangle
176  return (x1<fN/2.) ? par[0]*x1 : par[0]*(fN-x1);
177  case RTRIANGLE:
178  v1 = 1.;
179  return (x1<v1) ? par[0]*(fN-v1)*x1/v1 : par[0]*(fN-x1);
180  case UTRIANGLE: // asimmetrical triangle
181  v1 = 1.1;
182  return (x1<v1) ? par[0]*(fN-v1)*x1/v1 : par[0]*(fN-x1);
183  default:
184  cout<<"Unknown shape\n";
185  exit(1);
186  }
187 }
188 
189 // ----------------------------------------------------------------------------------------------
190 
191 Double_t fshapeACorr(Double_t *xx, Double_t *par)
192 {// AutoCorrelation at shift par[0]
193  Double_t par0[10],vv;
194  for(Int_t ii=0;ii<10;ii++) par0[ii] = par[ii];
195  par0[2] = par[0]; // the shift parameter for fshape
196  vv = fshape(xx,par0+1)*fshape(xx,par+1);
197  //cout<<"fshapeACorr("<<*xx<<","<<par[0]<<")="<<vv<<"\n";
198  return vv;
199 }
200 
201 // **********************************************************************************************
202 
203 // ----------------------------------------------------------------------------------------------
204 
205 // first derivative of the unipolar shape
206 Double_t difup(Double_t *xx, Double_t *par)
207 {
208  Double_t x1 = xx[0] - par[1];
209  if(x1 < 0.) x1 = 0.;
210  return par[0]*pow(x1,3.)*exp(-x1*16./fN) * (4.-16./fN*x1);
211 }
212 
213 // ----------------------------------------------------------------------------------------------
214 
215 // **********************************************************************************************
216 
217 void show(Double_t *xx, Int_t nn=200, Char_t *xname=NULL, Char_t *yname=NULL, Int_t lcolor=kBlue, Int_t mstyle=kStar, Double_t msize=0.75)
218 {
219  gr = new TGraph(nn,ni,xx);
220  gr->SetLineColor(lcolor);
221  gr->SetLineWidth(2);
222  gr->SetMarkerColor(kBlack);
223  gr->SetMarkerStyle(mstyle);
224  gr->SetMarkerSize(msize);
225  gr->SetTitle("");
226  if (xname) {gr->GetXaxis()->SetTitle(xname);}
227  else gr->GetXaxis()->SetTitle("Sample");
228  //if (yname) {gr->GetYaxis()->SetTitle(yname);}
229  if (yname) {gr->SetTitle(yname);}
230  //gr->Draw("ACP");
231  gr->Draw("APL");
232  if(gCnv == NULL) gCnv = c1;
233  gCnv->SetGridx(); gCnv->SetGridy();
234  gCnv->SetCrosshair();gCnv->SetTickx(); gCnv->SetTicky();
235  gCnv = NULL;
236 }
237 
238 // ----------------------------------------------------------------------------------------------
239 
240 void showsig(Double_t ampl, Double_t position, Double_t length, Char_t *fname=NULL)
241 {
242  if(fname==NULL) fname="f1_";
243  cout<<"drawing "<<fname<<"\n";
244  TF1 *tf = new TF1(fname,fshape,position,position+length,2);
245  tf->SetParameter(0,ampl);
246  tf->SetParameter(1,position);
247  tf->SetLineStyle(1);
248  tf->SetLineWidth(2);
249  tf->Draw("same");
250 }
251 
252 // ----------------------------------------------------------------------------------------------
253 
254 Double_t rms(Double_t *xx, Int_t nn=NX)
255 {
256  Double_t sum=0., sum2=0., xmax=0.,fN=(Double_t)nn;
257  for(Int_t ii=0;ii<nn;ii++)
258  {
259  if(xx[ii]>xmax) xmax = xx[ii];
260  sum += xx[ii];
261  sum2 += xx[ii]*xx[ii];
262  }
263  sum2 = sqrt(sum2/fN);
264  cout<<"srsq="<<sum2*sqrt(fN)<<", rms="<<sum2<<", mean="<<sum/fN<<", max="<<xmax<<"\n";
265  return sum2;
266 }
267 
268 
269 // ----------------------------------------------------------------------------------------------
270 
271 Double_t meanw(Double_t *xx, Int_t pos, Int_t nn)
272 {
273  Double_t sumw=0.,sumwin=0.;
274  for(Int_t ii=pos-nn+1;ii<=pos;ii++)
275  {
276  sumwin += xx[ii];
277  sumw += xx[ii]*((Double_t)ii);
278  }
279  return sumw/sumwin;
280 }
281 
282 // ----------------------------------------------------------------------------------------------
283 
284 TRandom grand;
285 void setSample(Double_t shift, Double_t noise_rms=10.)
286 {
287  Double_t vv;
288  Int_t ishift = (Int_t) shift;
289  Int_t ii;
290  gpos = shift; gnoise = noise_rms;
291  for(ii=0;ii<NX;ii++)
292  {
293  v[ii] = grand.Gaus(0.,noise_rms);
294  x[ii] = v[ii];
295  s[ii] = 0.;
296  if(ii>=shift && ii<shift+N)
297  {
298  s[ii] = shape((Double_t) (ii-shift))*ALevel/gmax;
299  x[ii] += s[ii];
300  }
301  }
302  if(!quiet)
303  {
304  cout<<"Signal: ";
305  vv = rms((s+ishift),N);
306  cout<<"Noise: ";
307  vv = rms(v+ishift,N);
308  cout<<"S+N: ";
309  vv = rms(x+ishift,N);
310  }
311 }
312 
313 // ----------------------------------------------------------------------------------------------
314 
315 //Double_t round(Double_t val, Double_t precision = PRECISION_SAMPLING)
316 //{
317 // return = floor(val/precision)*precision;
318 //}
319 Double_t round(Double_t x) {return x;}
320 void go()
321 {
322  // Calculate filter output and t/T
323  Int_t ii,i1,imax;
324  Int_t timewin = 10;
325  ymax = 0.;
326  for (ii=N;ii<NX;ii++) {y[ii] = 0.;yma[ii] = 0.;}//y_s[ii] = 0.;}
327  #ifdef NEED_FILTER_NOISE
328  for (ii=N;ii<NX;ii++)
329  #else
330  for (ii=N;ii<=((Int_t)gpos+1+3*N);ii++)
331  #endif
332  {
333  for(i1=0;i1<N;i1++)
334  {
335  y[ii] += round(x[ii-i1])*g[N-i1]; // calculate filter output
336  yma[ii] += x[ii-i1]*gma[N-i1]; // the same for moving average, for comparison
337  }
338  if(ymax<y[ii]) {ymax=y[ii]; imax=ii;}
339  dy[ii] = round(y[ii]) - round(y[ii-1]);
340  dyma[ii] = round(yma[ii]) - round(yma[ii-1]);
341  }
342  if(dy[imax]-dy[imax+1] <= 0.) {cout<<"Computational error!\n";exit();}
343  //tau = imax + dy[imax]/(dy[imax]-dy[imax+1]);
344  // parabolic approximation
345  tau = round(dy[imax])/(round(dy[imax])-round(dy[imax+1])) - 0.5;
346  //cout<<tau<<" - tau\n";
347  //correction line for parabolic approximation of amplitude is broken
348  //at point tau=1/2. The corrected max(Y) estimator is more stable and
349  //it is more easy implement in hardware.
350  //ymax += (tau/4.)*(dy[imax]+dy[imax+1]);;
351  tau += (Double_t)imax - fN;
352  #ifdef TRY_WEIGHTED
353  tauw = meanw(y,imax+timewin/2.,timewin);
354  #endif
355  if(!quiet)
356  {
357  cout<<"ymax["<<tau<<"]= "<<ymax<<"\n";
358  #ifdef TRY_WEIGHTED
359  cout<<"weigthed mean = "<<tauw<<"\n";
360  #endif
361  }
362 }
363 
364 // ----------------------------------------------------------------------------------------------
365 
366 
367 // Generate set of 'ntry' sequences with signal positioned at 'pos',
368 // superimposed with 'noise'
369 // If nscan > 1 then collect the correction coeffcients
370 // If acorr > 0 - do the amplitude correction
371 // Time correseries(ction is always performed.
372 // Evidently, if pos is integral number then corrections would be zero
373 Double_t series(Double_t pos, Double_t noise, Int_t ntry, Int_t nscans=1, Int_t acorr=1)
374 {
375  cout<<"series("<<pos<<","<<noise<<","<<ntry<<","<<nscans<<","<<acorr<<")"<<endl;
376  Int_t ii,iscan,i1;//,previ;
377  Double_t sumt,sumt2,ty[4];
378  #ifdef TRY_WEIGHTED
379  Double_t sumw=0.,sumw2=0.;
380  #endif
381  Double_t suma, suma2;
382  Double_t sumdevt2,sumdeva2;
383  //Double_t meana=ALevel/gmax*(1.+aCorr[(Int_t)(((pos-floor(pos))*NCorr))]);
384  Double_t meana=ALevel/gmax;
385  Double_t mean,sigma,ns,rc,poss;
386  Double_t v1,v2;
387  TCanvas *CnvAcc[2];
388  quiet = 1;
389  //cout<<"Attention. quiet set to 1\n";
390  //if(ntry>1000) return;
391  for(iscan=0;iscan<nscans;iscan++)
392  {
393  poss = pos + (Double_t)iscan / (Double_t)nscans;
394  suma = 0.; suma2=0.; sumt=0.; sumt2=0.; sumdevt2=0.; sumdeva2=0.;
395  for(ii = 0; ii<ntry; ii++)
396  {
397  // generate signal and noise
398  setSample(poss,noise);
399  // process the signal
400  go();
401  // the rest is to get the deviations of the corrected paremeters from the actual ones
402  sumt += tau; sumt2 += tau*tau;
403  #ifdef TRY_WEIGHTED
404  sumw += tauw; sumw2 += tauw*tauw;
405  #endif
406  suma += ymax; suma2 += ymax*ymax;
407  // first, get the index from the correction table
408  i1 = TMath::Nint((tau-floor(tau))*NCorr);
409  if(i1>=NCorr) {cout<<"Index error "<<i1<<endl; continue;}
410  //cout<<"poss="<<poss<<", tau="<<tau<<", Corr["<<i1<<"]="<<tCorr[i1]<<","<<aCorr[i1]<<"\n";
411  v2 = tau - tCorr[i1];
412  v1 = v2 - pos;
413  sumdevt2 += v1*v1;
414  // for amplitude use corrected index
415  v1 = ymax/(1.+aCorr[i1]) - meana;
416  //v1 = ymax - meana;
417  //cout<<ii<<": nc="<<v1<<", corrected="<<ymax*(1.-aCorr[i1]) - meana<<"\n";
418  sumdeva2 += v1*v1;
419  // periodical print
420  if(ii%1000==999)
421  {
422  //cout<<"["<<ii<<"]"<<tau<<", "<<tauw<<"\n";
423  ns=(Double_t)(ii+1);
424  mean = sumt/ns; sigma = sqrt(fabs(sumt2/ns - mean*mean));
425  cout<<"["<<ii<<"] mean="<<mean<<", s="<<sigma;
426  mean = suma/ns; sigma = sqrt(fabs(suma2/ns - mean*mean));
427  cout<<", F="<<mean<<" sF="<<sigma
428  <<", A="<<mean*gmax<<" sA="<<sigma*gmax<<"\n";
429  //cout<<"srt="<<sqrt(sumdevt2/ns)<<","<<1./sqrt(sumdevt2/ns)
430  // <<". sra="<<sqrt(sumdeva2/ns)<<","<<1./sqrt(sumdeva2/ns)<<"\n";
431  }
432  }
433  ns=(Double_t)ntry;
434  //cout<<"lindif: ";
435  ty[2] = sumt/ns;
436  ty[3] = suma/ns;
437  ty[0] = sqrt(fabs(sumt2/ns-ty[2]*ty[2]));
438  //ty[1] = gmax*sqrt(fabs(suma2/ns - (ty[3]*ty[3])));
439  ty[1] = sqrt(fabs(suma2/ns - (ty[3]*ty[3])));
440  if(ty[0]!=0.) ty[0] = 1./ty[0]; else ty[0] = 1.e100;
441  if(ty[1]!=0.) ty[1] = 1./ty[1]; else ty[1] = 1.e100;
442  ty[1] *= ty[3];
443  ty[3] *= gmax;
444  if(nscans>1)
445  {// collect corrections
446  if(iscan==0) v1= 0.; // often having floor problem here
447  else v1 = (ty[2] - floor(ty[2]))*(Double_t)NCorr;
448  i1 = TMath::Nint(v1); if(i1>=NCorr) i1=NCorr-1; if(i1<0) i1=0;
449  //if(i1==previ) cout<<"Replaced element "<<i1<<", try to increase corr. table\n"; previ = i1;
450  //i1 = (Int_t)v1;
451  //cout<<"p["<<iscan<<"]="<<poss<<" t="<<ty[2]<<","<<v1<<" i="<<i1<<" tc="<<(ty[2] - poss)<<" ac="<<(ty[3]/ALevel - 1.)<<"\n";
452  serCorr[0][i1] = ty[2] - poss;
453  serCorr[1][i1] = ty[3]/ALevel - 1.;
454  }
455  //record only the first scan
456  if(iscan==0)
457  {
458  serY[0]=ty[0]; serY[1]=ty[1]; serY[2]=ty[2]; serY[3]=ty[3];
459  serY[6] = 1./sqrt(sumdevt2/ns);
460  serY[7] = meana/sqrt(sumdeva2/ns);
461  // make correction if requested
462  if(acorr!=0)
463  {
464  // index of correction element calculated from the time/period
465  //i1 = (Int_t)((ty[2]-floor(ty[2]))*NCorr);
466  i1 = TMath::Nint((ty[2]-floor(ty[2]))*NCorr);
467  cout<<"Corr["<<i1<<"]="<<tCorr[i1]<<","<<aCorr[i1]<<"\n";
468  serY[2] -= tCorr[i1];
469  cout<<"time correction["<<i1<<"]="<<tCorr[i1]<<" 1/stdev="<<1./(gpos-serY[2])<<" from "<<ty[2]<<" to "<<serY[2]<<"\n";
470  serY[3] /= (1.+aCorr[i1]);// correct amplitudes
471  cout<<"ampl. correction="<<aCorr[i1]<<" from "<<ty[3]<<" to "<<serY[3]<<"\n";
472  }
473  }
474  /*
475  * cout<<"time="<<ty[2]<<", 1/s="<<ty[0]<<", Acc="<<serCorr[0][iscan]<<"\n";
476  * cout<<"ampl="<<ty[3]<<", M/s="<<ty[1]<<", Acc="<<serCorr[1][iscan]<<"\n";
477  */
478  }
479  #ifdef TRY_WEIGHTED
480  cout<<"weighted:";
481  mean = sumw/ns; sigma = sqrt(fabs(sumw2/ns-mean*mean));
482  cout<<"mean="<<mean<<", "<<sigma<<"\n";
483  #endif
484  quiet = 0;
485  return rc;
486 }
487 
488 // ----------------------------------------------------------------------------------------------
489 
490 void nfits(Double_t pos, Double_t noise, Int_t ntry)
491 {
492  Double_t sumt=0.,sumt2=0.,sum=0.,sum2=0.;
493  Double_t v,vi,ca=gmax/pars[0];
494  Double_t mean,sigma,ns=(Double_t)ntry;
495  quiet = 1;
496  for(Int_t ii=0;ii<ntry;ii++)
497  {
498  setSample(pos,noise);
499  gf->SetParameter(0,ALevel/gmax);
500  gf->SetParameter(1,pos);
501  if(gr) delete gr;
502  gr = new TGraph(200,ni,x);
503  gr->Fit(gf,"Q");
504  v = gf->GetParameter(0);
505  sum += v; sum2 += v*v;
506  v = gf->GetParameter(1);
507  sumt += v; sumt2 += v*v;
508  if(ii%100==99)
509  {
510  vi = (Double_t)(ii+1.);
511  mean = sumt/vi; sigma = sqrt(fabs(sumt2/vi - mean*mean));
512  cout<<"["<<ii<<"] Pos="<<mean<<", sigma="<<sigma
513  <<". ParError="<<gf->GetParError(1)<<"\n";
514  }
515  }
516  mean = sumt/ns; sigma = sqrt(fabs(sumt2/ns - mean*mean));
517  cout<<"Pos="<<mean<<", sPos="<<sigma<<", ePos="<<gf->GetParError(1);
518  mean = sum/ns; sigma = sqrt(fabs(sum2/ns - mean*mean));
519  cout<<". A="<<mean*ca<<", sA="<<sigma*ca<<", eA="<<ca*gf->GetParError(0)<<"\n";
520  quiet = 0;
521 }
522 
523 // ----------------------------------------------------------------------------------------------
524 
525 void saveCorrections()
526 {
527  Double_t v1,v2;
528  Int_t ii,kk,i1;
529  for(kk=0;kk<NCorr;kk++)
530  {
531  //if((serCorr[0][kk])==SMAL_UNDEF)
532  if((serCorr[0][kk])==SMAL_UNDEF)
533  {// some correction could be undefined, then use average from neighbors
534  //linear interpolation till next defined index
535  for(ii=1;ii<10;ii++) if((serCorr[0][kk+ii])!=SMAL_UNDEF) break;
536  //cout<<"kk="<<kk<<" ii="<<ii<<"\n";
537  for(i1=0;i1<ii;i1++)
538  {
539  tCorr[kk+i1] = serCorr[0][kk-1] + (serCorr[0][kk+ii] - serCorr[0][kk-1])*(Double_t)(i1+1.)/(Double_t)ii;
540  aCorr[kk+i1] = serCorr[1][kk-1] + (serCorr[1][kk+ii] - serCorr[1][kk-1])*(Double_t)(i1+1.)/(Double_t)ii;
541  }
542  kk += ii-1;
543  }
544  else
545  {
546  aCorr[kk] = serCorr[1][kk];
547  tCorr[kk] = serCorr[0][kk];
548  }
549  }
550  v1 = fabs(TMath::MinElement(NCorr,serCorr[0]));
551  v2 = fabs(TMath::MaxElement(NCorr,serCorr[0]));
552  serY[4] = 1./TMath::Max(v1,v2);
553  v1 = fabs(TMath::MinElement(NCorr,serCorr[1]));
554  v2 = fabs(TMath::MaxElement(NCorr,serCorr[1]));
555  serY[5] = 1./TMath::Max(v1,v2);
556  cout<<"1/syserr for t,A = "<<serY[4]<<", "<<serY[5]<<"\n";
557 }
558 
559 // ----------------------------------------------------------------------------------------------
560 
561 TCanvas *fc[grf_NY];
562 TGraphErrors *fg[grf_NY];
563 void doSeries(Double_t pos=50., Int_t ntry=100)
564 {
565  Int_t ii,kk;
566 
567  //for(ii=0;ii<4;ii++) if(fc[ii]) {cout<<"Deleting canvas "<<ii<<"\n";delete fc[ii];}
568  for(ii=0;ii<grf_NPTS;ii++)
569  {
570  grf_ex[ii] = 0.;
571  cout<<"Series "<<ii<<" out of "<<grf_NPTS
572  <<" for "<<ntry<<" points with S/N="<<grf_x[ii]<<"\n";
573  series(pos, ALevel/grf_x[ii],ntry,1,1);
574  for(kk=0;kk<grf_NY;kk++)
575  {
576  grf_y[kk][ii] = serY[kk];
577  if(kk==2) grf_ey[kk][ii] = 1./serY[kk-2];
578  else if(kk==3) grf_ey[kk][ii] = serY[kk]/serY[kk-2];
579  else grf_ey[kk][ii] = 0.;
580  }
581  }
582  // draw graphs
583  for(kk=0;kk<grf_NY;kk++)
584  {
585  Int_t lcolor,lwidth,lstyle,lmcolor,lmstyle;
586  Double_t lmsize;
587  //if(kk==5)continue; // skip "1/(Systematic Error) of Amplitude"
588  if(kk<4)
589  {
590  //cout<<"Creating canvas "<<kk<<"\n";
591  if(fc[kk]==NULL) fc[kk] = new TCanvas();
592  else fc[kk]->Clear();
593  //lstyle = (kk<2) ? kBlack : kBlack;
594  lstyle = 1;
595  lcolor = (kk<2) ? kBlack : kBlue;
596  lmstyle = kStar;
597  lwidth=2; lmcolor=lcolor; lmsize=1.;
598  fc[kk]->SetLogx(1);
599  if(kk<2)
600  {
601  fc[kk]->SetLogy(1);
602  }
603  }
604  else if(kk<6)
605  {//graph 4 and 5
606  // basic graphs
607  fc[kk-4]->cd();
608  lcolor=kRed; lwidth=4; lstyle=7; lmcolor=2; lmstyle=kDot; lmsize=.75;
609  }
610  else
611  {//graph 6 and 7
612  fc[kk-6]->cd();
613  lcolor=kBlue; lwidth=2; lstyle=2; lmcolor=lcolor; lmstyle=kCircle; lmsize=1.;
614  }
615  fg[kk] = new TGraphErrors(grf_NPTS,grf_x,grf_y[kk],grf_ex,grf_ey[kk]);
616  //gr = new TGraph(nn,ni,xx);
617  fg[kk]->SetLineColor(lcolor);
618  fg[kk]->SetLineWidth(lwidth);
619  fg[kk]->SetLineStyle(lstyle);
620  fg[kk]->SetMarkerColor(lmcolor);
621  fg[kk]->SetMarkerStyle(lmstyle);
622  fg[kk]->SetMarkerSize(lmsize);
623  //fg[kk]->SetTitle("");
624  (fg[kk]->GetXaxis())->SetTitle("A/RMS(noise)");
625  //fg[kk]->Draw("ACP");
626  if(kk<4)
627  {
628  fg[kk]->Draw("APL");
629  fg[kk]->SetTitle(grf_name[kk]);
630 
631  fc[kk]->SetGridx();
632  fc[kk]->SetGridy();
633  fc[kk]->SetCrosshair();
634  fc[kk]->SetTickx();
635  fc[kk]->SetTicky();
636 
637  }
638  else fg[kk]->Draw("PL");
639  }
640 }
641 
642 // ----------------------------------------------------------------------------------------------
643 
644 void getCoeffs()
645 {
646  Int_t ii;
647  pars[0] = 1.;
648 
649  // calculate maximum of the shape
650  gf = new TF1("gf",fshape,0,N,npars);
651  gf->SetParameter(0,pars[0]);
652  gmax = gf->GetMaximum();
653  gmaxx = gf->GetMaximumX();
654  cout<<"Filter max["<<gmaxx<<"]="<<gmax<<"\n";
655  // set up the matching filter
656  // Use L2 normalization
657  norm = 0.; norm_ma = 0.;
658  for(ii=0;ii<N;ii++)
659  {
660  g[ii] = shape(ii);
661  norm += g[ii]*g[ii];
662  gma[ii] = 1.;
663  norm_ma += gma[ii]*gma[ii];
664  }
665  //norm = sqrt(norm/N);
666  norm = sqrt(norm);
667  //norm_ma = sqrt(norm_ma/N);
668  norm_ma = sqrt(norm_ma);
669  cout<<"normalization = "<<norm<<", "<<norm_ma<<"\n";
670  for(ii=0;ii<N;ii++)
671  {
672  g[ii] = g[ii]/norm;
673  gma[ii] /= norm_ma;
674  y[ii]=0.; dy[ii]=0.; dyma[ii]=0.;
675  }
676  // Lets have the function normalized the same way
677  pars[0] = 1./norm;
678  gf->SetParameter(0,pars[0]);
679  gmax = gf->GetMaximum();
680  gmaxx = gf->GetMaximumX();
681  cout<<"Filter max["<<gmaxx<<"]="<<gmax<<"\n";
682 }
683 
684 // ----------------------------------------------------------------------------------------------
685 
686 void roundCoeffs(Double_t precision=1./64.)
687 {
688  for(Int_t ii=0;ii<N;ii++)
689  {
690  g[ii] = floor(g[ii]/precision)*precision;
691  }
692 }
693 
694 // ----------------------------------------------------------------------------------------------
695 
696 void showShape()
697 {
698  // show the shape
699  gCanv[0]->cd();
700  show(g,N,0,0,0,21,2.); gr->GetYaxis()->SetTitle("Filter Coefficients");
701  gr->Draw(); showsig(pars[0],0.,N,"f1_shape");
702 
703  // show the input signal superposed with noise
704  gCanv[1]->cd();
705  Double_t rmsnoise = ALevel/10.;
706  series(50.,rmsnoise,100);
707  TString txt = "Signal (A=";
708  txt+=ALevel; txt+=" plus noise (RMS="; txt+=rmsnoise; txt+=")";
709  show(x);
710  //gr->GetYaxis()->SetTitle("Signal (A=100) plus noise (RMS=10)");
711  gr->GetYaxis()->SetTitle(txt);
712  gr->Draw();
713  showsig(pars[0]*ALevel/gmax,gpos,N,"f1_signal");
714 
715  gCanv[2]->cd();
716  show(y); gr->GetYaxis()->SetTitle("Filter output");
717  gr->Draw();
718 }
719 
720 // ----------------------------------------------------------------------------------------------
721 
722 void clearCorrections()
723 {
724  Int_t ii;
725  for(ii=0;ii<NCorr;ii++){ aCorr[ii]=0.; tCorr[ii]=0.;}
726  for(ii=0;ii<NCorr;ii++) {serCorr[0][ii] = SMAL_UNDEF; serCorr[1][ii] = SMAL_UNDEF;}
727 }
728 
729 // ----------------------------------------------------------------------------------------------
730 
731 void getCorrections()
732 {
733  Int_t ii;
734  // Calculate the systematic errors and amplitude corrections
735  // by running series with zero noise and 100-point scan inside one clock
736  // period
737  clearCorrections();
738  cout<<"Collecting corrections...\n";
739  gCanv[3]->cd();
740  series(NN,0.,1,NCorr,0);
741  saveCorrections();
742  // show corrections
743  Double_t t1[NCorr];
744  TGraph *corrGraph0,*corrGraph1;
745  for(ii=0;ii<NCorr;ii++) t1[ii] = (Double_t)ii/(Double_t)NCorr;
746  corrGraph0 = new TGraph(NCorr,t1,tCorr);
747  corrGraph0->SetLineColor(kBlue);corrGraph0->SetLineWidth(2);
748  corrGraph0->SetMarkerColor(kBlack);corrGraph0->SetMarkerStyle(kDot);
749  corrGraph0->GetXaxis()->SetTitle("t/T");
750  corrGraph0->SetTitle("Time Corrections");
751  gCanv[3]->SetGridx(); gCanv[3]->SetGridy();
752  gCanv[3]->SetTickx(); gCanv[3]->SetTicky();
753  corrGraph0->Draw("APL");
754  gCanv[4]->cd();
755  canv = gCanv[4];
756  corrGraph1 = new TGraph(NCorr,t1,aCorr);
757  corrGraph1->SetLineColor(kRed);corrGraph1->SetLineWidth(2);
758  corrGraph1->SetMarkerColor(kBlack);corrGraph1->SetMarkerStyle(kDot);
759  corrGraph1->GetXaxis()->SetTitle("t/T");
760  corrGraph1->SetTitle("Amplitude Corrections");
761  canv->SetGridx(); canv->SetGridy();
762  canv->SetTickx(); canv->SetTicky();
763  corrGraph1->Draw("APL");
764  //corrGraph1->Draw("SAME");
765 }
766 
767 // ----------------------------------------------------------------------------------------------
768 
769 void showAutocorrelation()
770 {
771  Int_t ii;
772  cout<<"Drawing the autocorrelation function in the interval [0:1]\n";
773  gCanv[5]->cd();
774  Double_t xacrl[NCorr],yacrl[NCorr],vg;
775  TF1 *facrl = new TF1("facrl",fshapeACorr,0.,NCorr,3);
776  facrl->SetParameter(1,1.);
777  facrl->SetParameter(0,0.);
778  vg = facrl->Integral(0.,fN);
779  for(ii=0;ii<NCorr;ii++)
780  {
781  xacrl[ii] = (Double_t)ii/(Double_t)NCorr;
782  facrl->SetParameter(0,xacrl[ii]);
783  yacrl[ii] = facrl->Integral(0.,fN)/vg;
784  //cout<<"aCorr("<<xacrl[ii]<<"="<<yacrl[ii]<<"\n";
785  }
786  TGraph *gACorr = new TGraph(NCorr,xacrl,yacrl);
787  gACorr->SetLineStyle(1);
788  gACorr->SetLineWidth(2);
789  gACorr->SetTitle("AutoCorrelation");
790  gACorr->GetXaxis()->SetTitle("t/T");
791  gACorr->Draw("AL");
792 }
793