00001
00007 #include <math.h>
00008 #include <iostream>
00009
00010 #include "CondFormats/CastorObjects/interface/CastorPedestalWidth.h"
00011
00012 namespace {
00013 int offset (int fCapId1, int fCapId2) {
00014 static int offsets [4] = {0, 1, 3, 6};
00015 if (fCapId1 < fCapId2) {
00016 int tmp = fCapId1; fCapId1 = fCapId2; fCapId2 = tmp;
00017 }
00018 return offsets [fCapId1] + fCapId2;
00019 }
00020 }
00021
00022 CastorPedestalWidth::CastorPedestalWidth (int fId) : mId (fId) {
00023 for (int i = 10; --i >= 0; *(&mSigma00 + i) = 0) {}
00024 }
00025
00026 float CastorPedestalWidth::getWidth (int fCapId) const {
00027 return sqrt (*(getValues () + offset (fCapId, fCapId)));
00028 }
00029
00030 float CastorPedestalWidth::getSigma (int fCapId1, int fCapId2) const {
00031 return *(getValues () + offset (fCapId1, fCapId2));
00032 }
00033
00034 void CastorPedestalWidth::setSigma (int fCapId1, int fCapId2, float fSigma) {
00035 *(&mSigma00 + offset (fCapId1, fCapId2)) = fSigma;
00036 }
00037
00038
00039 void CastorPedestalWidth::makeNoise (unsigned fFrames, const double* fGauss, double* fNoise) const {
00040 double s_xx_mean = (getSigma (0,0) + getSigma (1,1) + getSigma (2,2) + getSigma (3,3)) / 4;
00041 double s_xy_mean = (getSigma (1,0) + getSigma (2,1) + getSigma (3,2) + getSigma (3,0)) / 4;
00042 double sigma = sqrt (0.5 * (s_xx_mean + sqrt (s_xx_mean*s_xx_mean - 2*s_xy_mean*s_xy_mean)));
00043 double corr = sigma == 0 ? 0 : 0.5*s_xy_mean / sigma;
00044 for (unsigned i = 0; i < fFrames; i++) {
00045 fNoise [i] = fGauss[i]*sigma;
00046 if (i > 0) fNoise [i] += fGauss[i-1]*corr;
00047 if (i < fFrames-1) fNoise [i] += fGauss[i+1]*corr;
00048 }
00049 }