SiStripCommissioningAnalysis/src/SiStripPulseShape.cc

0001 #include "DQM/SiStripCommissioningAnalysis/interface/SiStripPulseShape.h"
0002 #include <TF1.h>
0003 #include <TMath.h>
0004
0005 /* New Pulse Shape Fit by G. Auzinger June 2017
0006 following methods are used to describe a pulse shape for a 3 stage CR-CR-RC pre-amplifier (CR)
0007 + shaper (CR + RC) with time constants x, y, z respectively
0008 some special cases (x=0. y=z, x=0 && y=z) are considered (divergence of equation) and
0009 the shaper CR y is approximated via a pol(6) of x which gives a resonable respresentation and
0010 reduces the number of free parameters -- this shape is derived from first principle and can
0011 be used to fit peak mode signals and deconvolution mode signals
0012 following parameters are added in addition to the time constants x, z:
0013 a_0 ... baseline offset amplitude
0014 s   ... scale parameter
0015 t_0 ... turn_on time of the pulse
0016 par[5] ... scale parameter for the time slices in deco mode
0017 */
0018
0019 double pulse_raw(double x, double y, double z, double t) {
0020   double result1, result2, result3;
0021
0022   result1 = z * y * exp(-t / y);
0023   result1 /= pow(y, 2) - (x + z) * y + z * x;
0024
0025   result2 = z * x * exp(-t / x);
0026   result2 /= pow(x, 2) - (x - z) * y - z * x;
0027
0028   result3 = pow(z, 2) * exp(-t / z);
0029   result3 /= pow(z, 2) + (x - z) * y - z * x;
0030
0031   return result1 + result2 + result3;
0032 }
0033
0034 double pulse_x0(double y, double z, double t) { return z / (y - z) * (exp(-t / y) - exp(-t / z)); }
0035
0036 double pulse_yz(double x, double z, double t) {
0037   double result1, result2;
0038
0039   result1 = exp(-t / x) - exp(-t / z);
0040   result1 *= z * x / (z - x);
0041
0042   result2 = t * exp(-t / z);
0043
0044   return (result1 + result2) / (z - x);
0045 }
0046
0047 double pulse_x0_yz(double z, double t) { return t / z * exp(-t / z); }
0048
0049 double pulse(double x, double y, double z, double t) {
0050   if (x > y) {
0051     double pivot = x;
0052     x = y;
0053     y = pivot;
0054   }
0055
0056   if ((x == 0) && (y == z))
0057     return pulse_x0_yz(z, t);
0058   else if (x == 0)
0059     return pulse_x0(y, z, t);
0060   else if (y == z)
0061     return pulse_yz(x, z, t);
0062   else
0063     return pulse_raw(x, y, z, t);
0064 }
0065
0066 double fpeak(double *x, double *par) {
0067   double xx = par[0];
0068   double y = par[1];
0069   double z = par[2];
0070   double a_0 = par[3];
0071   double s = par[4];
0072   double t_0 = par[5];
0073   double t = x[0] - t_0;
0074
0075   // below turn-on time return just a constant
0076   if (x[0] < t_0)
0077     return a_0;
0078   // elswhere return the pulse
0079   return a_0 + s * pulse(xx, y, z, t);
0080 }
0081
0082 double fturnOn(double *x, double *par) {
0083   double a_0 = par[0];
0084   double s = par[1];
0085   double t_0 = par[2];
0086   double width = par[3];
0087
0088   return a_0 + s * TMath::Erf((x[0] - t_0) / width);
0089 }
0090
0091 double fdecay(double *x, double *par) {
0092   double s = par[0];
0093   double c_exp = par[1];
0094   double c_pow = par[2];
0095
0096   return s * TMath::Exp(x[0] * c_exp) * (1 + x[0] * c_pow);
0097 }
0098
0099 double fdecay(double *x, double *par);
0100
0101 double fdeconv(double *x, double *par) {
0102   double xm = par[6] * (x[0] - 25);
0103   double xp = par[6] * (x[0] + 25);
0104   double xz = par[6] * x[0];
0105   return 1.2131 * fpeak(&xp, par) - 1.4715 * fpeak(&xz, par) + 0.4463 * fpeak(&xm, par);
0106 }
0107
0108 double fpeak_convoluted(double *x, double *par) {
0109   TF1 f("peak_convoluted", fpeak, 0, 250, 4);
0110   return f.IntegralError(x[0] - par[4] / 2., x[0] + par[4] / 2., par, nullptr, 1.) / (par[4]);
0111 }
0112
0113 double fdeconv_convoluted(double *x, double *par) {
0114   double xm = (x[0] - 25);
0115   double xp = (x[0] + 25);
0116   double xz = x[0];
0117   return 1.2131 * fpeak_convoluted(&xp, par) - 1.4715 * fpeak_convoluted(&xz, par) +
0118          0.4463 * fpeak_convoluted(&xm, par);
0119 }