CSCCondDB/src/CSCFitAFEBThr.cc

0001 #include "Minuit2/FunctionMinimum.h"
0002 #include "Minuit2/MnStrategy.h"
0003
0004 #include <OnlineDB/CSCCondDB/interface/CSCFitAFEBThr.h>
0005 #include <OnlineDB/CSCCondDB/interface/CSCThrTurnOnFcn.h>
0006
0007 #include <cmath>
0008 #include <vector>
0009 #include <iostream>
0010
0011 using namespace ROOT::Minuit2;
0012 using namespace std;
0013
0014 CSCFitAFEBThr::CSCFitAFEBThr() {
0015   theOBJfun = new CSCThrTurnOnFcn();
0016   theFitter = new VariableMetricMinimizer();
0017 }
0018
0019 CSCFitAFEBThr::~CSCFitAFEBThr() {
0020   delete theFitter;
0021   delete theOBJfun;
0022 }
0023
0024 bool CSCFitAFEBThr::ThresholdNoise(const std::vector<float>& inputx,
0025                                    const std::vector<float>& inputy,
0026                                    const int& npulses,
0027                                    std::vector<int>& dacoccup,
0028                                    std::vector<float>& mypar,
0029                                    std::vector<float>& ermypar,
0030                                    float& ercorr,
0031                                    float& chisq,
0032                                    int& ndf,
0033                                    int& niter,
0034                                    float& edm) const {
0035   bool status = false;
0036
0037   std::vector<double> parinit(2, 0);
0038   std::vector<double> erparinit(2, 0);
0039
0040   /// initial parameters, parinit[0]-threshold,  parinit[1]-noise
0041   parinit[0] = 30.0;
0042   parinit[1] = 2.0;
0043
0044   erparinit[0] = 20;
0045   erparinit[1] = 0.5;
0046
0047   /// do not fit input[y]==max and input[y]==0.0; calculate binom. error;
0048   std::vector<float> x;
0049   std::vector<float> y;
0050   std::vector<float> ynorm;
0051   std::vector<float> ery;
0052   x.clear();
0053   y.clear();
0054   ynorm.clear();
0055   ery.clear();
0056
0057   /// ndf > 0 if there is input data,number of points to fit > 2 and
0058   ///         fit did not fail.
0059   /// ndf = 0 if number of points to fit = 2
0060   /// ndf =-1 .......................... = 1
0061   /// ndf =-2 .......................... = 0
0062   /// ndf =-3 fit failed (number of points to fit was > 2)
0063   /// ndf =-4 no input data
0064
0065   int sum = 0;
0066   float r;
0067   for (size_t i = 0; i < inputx.size(); i++) {
0068     if (inputy[i] > 0.0)
0069       sum++;
0070     r = inputy[i] / (float)dacoccup[i];
0071     ynorm.push_back(r);
0072     //     std::cout<<" "<<ynorm[i];
0073   }
0074   //  std::cout<<std::endl;
0075   if (sum == 0) {
0076     ndf = -4;
0077     return status;
0078   }
0079
0080   int nbeg = inputx.size();
0081   // for(size_t i=inputx.size()-1;i>=0;i--) // Wrong.
0082   // Because i is unsigned, i>=0 is always true,
0083   // and the loop termination condition  is never reached.
0084   // We offset by 1.
0085   for (size_t i = inputx.size(); i > 0; i--) {
0086     if (ynorm[i - 1] < 1.0)
0087       nbeg--;
0088     if (ynorm[i - 1] == 1.0)
0089       break;
0090   }
0091
0092   for (size_t i = nbeg; i < inputx.size(); i++) {
0093     if (ynorm[i] > 0.0) {
0094       x.push_back(inputx[i]);
0095       y.push_back(ynorm[i]);
0096
0097       float p = inputy[i] / (float)dacoccup[i];
0098       float s = (float)dacoccup[i] * p * (1.0 - p);
0099       s = sqrt(s) / (float)dacoccup[i];
0100       ery.push_back(s);
0101     }
0102   }
0103
0104   /// do not fit data with less than 3 points
0105   ndf = x.size() - 2;
0106   if (ndf <= 0)
0107     return status;
0108
0109   /// Calculate approximate initial threshold par[0]
0110   float half = 0.5;
0111   float dmin = 999999.0;
0112   float diff;
0113   for (size_t i = 0; i < x.size(); i++) {
0114     diff = y[i] - half;
0115     if (diff < 0.0)
0116       diff = -diff;
0117     if (diff < dmin) {
0118       dmin = diff;
0119       parinit[0] = x[i];
0120     }  // par[0] from data
0121     //std::cout<<i+1<<" "<<x[i]<<" "<<y[i]<<" "<<ery[i]<<std::endl;
0122   }
0123
0124   /// store data, errors and npulses for fit
0125   theOBJfun->setData(x, y);
0126   theOBJfun->setErrors(ery);
0127   theOBJfun->setNorm(1.0);
0128
0129   // for(size_t int i=0;i<x.size();i++) std::cout<<" "<<x[i]<<" "<<y[i]
0130   //                                               <<" "<<ery[i]<<std::endl;
0131
0132   /// Fit  as 1D, <=500 iterations, edm=10**-5 (->0.1)
0133   FunctionMinimum fmin =
0134       theFitter->Minimize(*theOBJfun, MnUserParameterState{parinit, erparinit}, MnStrategy{1}, 500, 0.1);
0135
0136   status = fmin.IsValid();
0137
0138   if (status) {
0139     mypar[0] = (float)fmin.Parameters().Vec()(0);
0140     mypar[1] = (float)fmin.Parameters().Vec()(1);
0141     ermypar[0] = (float)sqrt(fmin.Error().Matrix()(0, 0));
0142     ermypar[1] = (float)sqrt(fmin.Error().Matrix()(1, 1));
0143     ercorr = 0;
0144     if (ermypar[0] != 0.0 && ermypar[1] != 0.0)
0145       ercorr = (float)fmin.Error().Matrix()(0, 1) / (ermypar[0] * ermypar[1]);
0146
0147     chisq = fmin.Fval();
0148     ndf = y.size() - mypar.size();
0149     niter = fmin.NFcn();
0150     edm = fmin.Edm();
0151   } else
0152     ndf = -3;
0153   return status;
0154 }