HcalAlgos/src/HcalSiPMnonlinearity.cc

0001 #include <iostream>
0002 #include <gsl/gsl_poly.h>
0003 #include <gsl/gsl_complex.h>
0004 #include "CalibCalorimetry/HcalAlgos/interface/HcalSiPMnonlinearity.h"
0005
0006 // Assume parameters come to us from the reco side; i.e.,
0007 // true pes = corfun(pixelsfired). But we want to invert that.
0008 //
0009 int HcalSiPMnonlinearity::getPixelsFired(int inpes) const {
0010   gsl_complex z[3];
0011   double w = -inpes;
0012   // normalize params
0013   double a = a2 / w;
0014   double b = b1 / w;
0015   double c = c0 / w;
0016   int nroots = gsl_poly_complex_solve_cubic(a, b, c, &z[1], &z[2], &z[3]);
0017   assert(nroots);
0018
0019   // all use cases tested over the full range of anticipated values;
0020   // the first root is always the right one.
0021   double realpix = 0;
0022   // find real roots
0023   for (int i = 0; i < 3; ++i) {
0024     if (z[i].dat[1] == 0) {
0025       realpix = z[i].dat[0];
0026       break;
0027     }
0028   }
0029
0030   return realpix > 0 ? (int)(realpix + 0.5) : 0;
0031 }