Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:11:29

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 }