SiStripCommissioningAnalysis/src/Utility.cc

0001 #include "DQM/SiStripCommissioningAnalysis/src/Utility.h"
0002 #include <algorithm>
0003
0004 // ----------------------------------------------------------------------------
0005 //
0006 sistrip::LinearFit::LinearFit() : x_(), y_(), e_(), ss_(0.), sx_(0.), sy_(0.) { ; }
0007
0008 // ----------------------------------------------------------------------------
0009 //
0010 void sistrip::LinearFit::add(const float& x, const float& y) {
0011   float e = 1.;  // default
0012   x_.push_back(x);
0013   y_.push_back(y);
0014   e_.push_back(e);
0015   float wt = 1. / sqrt(e);
0016   ss_ += wt;
0017   sx_ += x * wt;
0018   sy_ += y * wt;
0019 }
0020
0021 // ----------------------------------------------------------------------------
0022 //
0023 void sistrip::LinearFit::add(const float& x, const float& y, const float& e) {
0024   if (e > 0.) {
0025     x_.push_back(x);
0026     y_.push_back(y);
0027     e_.push_back(e);
0028     float wt = 1. / sqrt(e);
0029     ss_ += wt;
0030     sx_ += x * wt;
0031     sy_ += y * wt;
0032   }
0033 }
0034
0035 // ----------------------------------------------------------------------------
0036 //
0037 void sistrip::LinearFit::fit(Params& params) {
0038   float s2 = 0.;
0039   float b = 0;
0040   for (uint16_t i = 0; i < x_.size(); i++) {
0041     float t = (x_[i] - sx_ / ss_) / e_[i];
0042     s2 += t * t;
0043     b += t * y_[i] / e_[i];
0044   }
0045
0046   // Set parameters
0047   params.n_ = x_.size();
0048   params.b_ = b / s2;
0049   params.a_ = (sy_ - sx_ * params.b_) / ss_;
0050   params.erra_ = sqrt((1. + (sx_ * sx_) / (ss_ * s2)) / ss_);
0051   params.errb_ = sqrt(1. / s2);
0052
0053   /*
0054     params.chi2_ = 0.;
0055     *q=1.0;
0056     if (mwt == 0) {
0057     for (i=1;i<=ndata;i++)
0058     *chi2 += SQR(y[i]-(*a)-(*b)*x[i]);
0059     sigdat=sqrt((*chi2)/(ndata-2));
0060     *sigb *= sigdat;
0061     */
0062 }
0063
0064 // ----------------------------------------------------------------------------
0065 //
0066 sistrip::MeanAndStdDev::MeanAndStdDev() : s_(0.), x_(0.), xx_(0.), vec_() { ; }
0067
0068 // ----------------------------------------------------------------------------
0069 //
0070 void sistrip::MeanAndStdDev::add(const float& x, const float& e) {
0071   if (e > 0.) {
0072     float wt = 1. / sqrt(e);
0073     s_ += wt;
0074     x_ += x * wt;
0075     xx_ += x * x * wt;
0076   } else {
0077     s_++;
0078     x_ += x;
0079     xx_ += x * x;
0080   }
0081   vec_.push_back(x);
0082 }
0083
0084 // ----------------------------------------------------------------------------
0085 //
0086 void sistrip::MeanAndStdDev::fit(Params& params) {
0087   if (s_ > 0.) {
0088     float m = x_ / s_;
0089     float t = xx_ / s_ - m * m;
0090     if (t > 0.) {
0091       t = sqrt(t);
0092     } else {
0093       t = 0.;
0094     }
0095     params.mean_ = m;
0096     params.rms_ = t;
0097   }
0098   if (!vec_.empty()) {
0099     sort(vec_.begin(), vec_.end());
0100     uint16_t index = vec_.size() % 2 ? vec_.size() / 2 : vec_.size() / 2 - 1;
0101     params.median_ = vec_[index];
0102   }
0103 }