00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "pdscalibdata/CsPadCommonModeSubV1.h"
00017
00018
00019
00020
00021 #include <algorithm>
00022 #include <stdexcept>
00023 #include <fstream>
00024 #include <cmath>
00025
00026
00027
00028
00029 #include "MsgLogger/MsgLogger.h"
00030 #include "pdscalibdata/CsPadPixelStatusV1.h"
00031
00032
00033
00034
00035
00036 namespace {
00037
00038 const char logger[] = "CsPadCommonModeSubV1";
00039
00040 }
00041
00042
00043
00044
00045
00046
00047 namespace pdscalibdata {
00048
00049
00050
00051
00052 CsPadCommonModeSubV1::CsPadCommonModeSubV1 ()
00053 : m_mode(uint32_t(None))
00054 {
00055 std::fill_n(m_data, int(DataSize), 0.0);
00056 }
00057
00058 CsPadCommonModeSubV1::CsPadCommonModeSubV1 (const std::string& fname)
00059 : m_mode(uint32_t(None))
00060 {
00061 std::fill_n(m_data, int(DataSize), 0.0);
00062
00063 m_mode =1;
00064 m_data[0]=50;
00065 m_data[1]=10;
00066 m_data[2]=100;
00067
00068
00069 std::ifstream in(fname.c_str());
00070 if (in.good()) {
00071
00072
00073 if (not (in >> m_mode)) {
00074 const std::string msg = "Common mode file does not have enough data: "+fname;
00075 MsgLogRoot(error, msg);
00076 throw std::runtime_error(msg);
00077 }
00078
00079
00080
00081 double* it = m_data;
00082 size_t count = 0;
00083 while(in and count != DataSize) {
00084 in >> *it++;
00085 ++ count;
00086 }
00087 }
00088
00089 MsgLog(logger, trace, "CsPadCommonModeSubV1: mode=" << m_mode << " data=" << m_data[0] << "," << m_data[1] << "," << m_data[2]);
00090 }
00091
00092 CsPadCommonModeSubV1::CsPadCommonModeSubV1 (CommonMode mode, const double data[DataSize])
00093 : m_mode(uint32_t(mode))
00094 {
00095 std::copy(data, data+DataSize, m_data);
00096 MsgLog(logger, trace, "CsPadCommonModeSubV1: mode=" << m_mode << " data=" << m_data[0] << "," << m_data[1] << "," << m_data[2]);
00097 }
00098
00099
00100
00101
00102
00103 CsPadCommonModeSubV1::~CsPadCommonModeSubV1 ()
00104 {
00105 }
00106
00107 float
00108 CsPadCommonModeSubV1::findCommonMode(const int16_t* sdata,
00109 const float* peddata,
00110 const uint16_t *pixStatus,
00111 unsigned ssize,
00112 int stride) const
00113 {
00114
00115 if (m_mode == None) return float(UnknownCM);
00116
00117
00118
00119 if (not peddata) return float(UnknownCM);
00120
00121
00122 const int low = -1000;
00123 const int high = 2000;
00124 const unsigned hsize = high-low;
00125 int hist[hsize];
00126 std::fill_n(hist, hsize, 0);
00127 unsigned long count = 0;
00128
00129
00130 for (unsigned c = 0, p = 0; c != ssize; ++ c, p += stride) {
00131
00132
00133
00134 if (pixStatus and pixStatus[p]) continue;
00135
00136
00137 double dval = sdata[p] - peddata[p];
00138 int val = dval < 0 ? int(dval - 0.5) : int(dval + 0.5);
00139
00140
00141 unsigned bin = unsigned(val - low);
00142
00143
00144 if (bin < hsize) {
00145 ++hist[bin] ;
00146 ++ count;
00147 }
00148
00149 }
00150
00151 MsgLog(logger, debug, "histo filled count = " << count);
00152
00153
00154
00155
00156
00157 int peakPos = -1;
00158 int peakCount = -1;
00159 int hmRight = hsize;
00160 int thresh = int(m_data[2]);
00161 if(thresh<=0) thresh=100;
00162 for (unsigned i = 0; i < hsize; ++ i ) {
00163 if (hist[i] > peakCount and hist[i] > thresh) {
00164 peakPos = i;
00165 peakCount = hist[i];
00166 } else if (peakCount > 0 and hist[i] <= peakCount/2) {
00167 hmRight = i;
00168 break;
00169 }
00170 }
00171
00172
00173 if (peakPos < 0) {
00174 MsgLog(logger, debug, "peakPos = " << peakPos);
00175 return float(UnknownCM);
00176 }
00177
00178
00179 int hmLeft = -1;
00180 for (int i = peakPos; hmLeft < 0 and i >= 0; -- i) {
00181 if(hist[i] <= peakCount/2) hmLeft = i;
00182 }
00183 MsgLog(logger, debug, "peakPos = " << peakPos << " peakCount = " << peakCount
00184 << " hmLeft = " << hmLeft << " hmRight = " << hmRight);
00185
00186
00187 int fwhm = hmRight - hmLeft;
00188 double sigma = fwhm / 2.36;
00189
00190
00191 double mean = peakPos;
00192 for (int j = 0; j < 2; ++j) {
00193 int s0 = 0;
00194 double s1 = 0;
00195 double s2 = 0;
00196 int d = int(sigma*2+0.5);
00197 for (int i = std::max(0,peakPos-d); i < int(hsize) and i <= peakPos+d; ++ i) {
00198 s0 += hist[i];
00199 s1 += (i-mean)*hist[i];
00200 s2 += (i-mean)*(i-mean)*hist[i];
00201 }
00202 mean = mean + s1/s0;
00203 sigma = std::sqrt(s2/s0 - (s1/s0)*(s1/s0));
00204 }
00205 mean += low;
00206
00207 MsgLog(logger, debug, "mean = " << mean << " sigma = " << sigma);
00208
00209
00210 if (mean > m_data[0] or sigma > m_data[1]) return float(UnknownCM);
00211
00212 return mean;
00213 }
00214
00215 }