CastorObjects/src/CastorPedestalWidth.cc

0001 /**
0002 \class CastorPedestalWidth
0003 \author Fedor Ratnikov (UMd)
0004 correlation matrix for pedestals
0005 $Author: ratnikov
0006 $Date: 2009/03/24 16:05:35 $
0007 $Revision: 1.9 $
0008 Adapted for CASTOR by L. Mundim
0009 */
0010
0011 #include <cmath>
0012 #include <iostream>
0013
0014 #include "CondFormats/CastorObjects/interface/CastorPedestalWidth.h"
0015
0016 namespace {
0017   int offset(int fCapId1, int fCapId2) {
0018     //    static int offsets [4] = {0, 1, 3, 6};
0019     //    if (fCapId1 < fCapId2) { // swap
0020     //      int tmp = fCapId1; fCapId1 = fCapId2; fCapId2 = tmp;
0021     //    }
0022     //    return offsets [fCapId1] + fCapId2;
0023     return fCapId1 * 4 + fCapId2;
0024   }
0025 }  // namespace
0026
0027 CastorPedestalWidth::CastorPedestalWidth(int fId) : mId(fId) {
0028   for (int i = 16; --i >= 0; *(&mSigma00 + i) = 0) {
0029   }
0030 }
0031
0032 float CastorPedestalWidth::getWidth(int fCapId) const { return sqrt(*(getValues() + offset(fCapId, fCapId))); }
0033
0034 float CastorPedestalWidth::getSigma(int fCapId1, int fCapId2) const {
0035   return *(getValues() + offset(fCapId1, fCapId2));
0036 }
0037
0038 void CastorPedestalWidth::setSigma(int fCapId1, int fCapId2, float fSigma) {
0039   *(&mSigma00 + offset(fCapId1, fCapId2)) = fSigma;
0040 }
0041
0042 // produces pedestal noise in assumption of near correlations and small variations
0043 void CastorPedestalWidth::makeNoise(unsigned fFrames, const double* fGauss, double* fNoise) const {
0044   double s_xx_mean = (getSigma(0, 0) + getSigma(1, 1) + getSigma(2, 2) + getSigma(3, 3)) / 4;
0045   double s_xy_mean = (getSigma(1, 0) + getSigma(2, 1) + getSigma(3, 2) + getSigma(3, 0)) / 4;
0046   double sigma = sqrt(0.5 * (s_xx_mean + sqrt(s_xx_mean * s_xx_mean - 2 * s_xy_mean * s_xy_mean)));
0047   double corr = sigma == 0 ? 0 : 0.5 * s_xy_mean / sigma;
0048   for (unsigned i = 0; i < fFrames; i++) {
0049     fNoise[i] = fGauss[i] * sigma;
0050     if (i > 0)
0051       fNoise[i] += fGauss[i - 1] * corr;
0052     if (i < fFrames - 1)
0053       fNoise[i] += fGauss[i + 1] * corr;
0054   }
0055 }