00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "cspad_mod/CsPadPedestals.h"
00017
00018
00019
00020
00021 #include <fstream>
00022 #include <cmath>
00023
00024
00025
00026
00027 #include "MsgLogger/MsgLogger.h"
00028
00029
00030
00031
00032
00033 using namespace Psana;
00034
00035
00036 using namespace cspad_mod;
00037 PSANA_MODULE_FACTORY(CsPadPedestals)
00038
00039 namespace {
00040
00041
00042 double stddev2(unsigned count, int64_t sum, int64_t sum2)
00043 {
00044
00045
00046
00047 int64_t offset = sum / count;
00048 int64_t sum2o = sum2 - 2*offset*sum + count*offset*offset;
00049 int64_t sumo = sum - offset*count;
00050 double avg = double(sumo) / count;
00051 return double(sum2o) / count - avg*avg;
00052 }
00053
00054 }
00055
00056
00057
00058
00059
00060 namespace cspad_mod {
00061
00062
00063
00064
00065 CsPadPedestals::CsPadPedestals (const std::string& name)
00066 : Module(name)
00067 , m_pedFile()
00068 , m_noiseFile()
00069 , m_src()
00070 , m_segMask()
00071 , m_sum()
00072 , m_sum2()
00073 {
00074 m_pedFile = configStr("output", "cspad-pedestals.dat");
00075 m_noiseFile = configStr("noise", "cspad-noise.dat");
00076
00077
00078 std::fill_n(&m_count[0], int(MaxQuads), 0UL);
00079 std::fill_n(&m_segMask[0], int(MaxQuads), 0U);
00080 std::fill_n(&m_sum[0][0][0][0], MaxQuads*MaxSectors*NumColumns*NumRows, int64_t(0));
00081 std::fill_n(&m_sum2[0][0][0][0], MaxQuads*MaxSectors*NumColumns*NumRows, int64_t(0));
00082 }
00083
00084
00085
00086
00087 CsPadPedestals::~CsPadPedestals ()
00088 {
00089 }
00090
00091
00092 void
00093 CsPadPedestals::beginRun(Event& evt, Env& env)
00094 {
00095
00096
00097
00098
00099 Source src(configStr("source", "DetInfo(:Cspad)"));
00100 int count = 0;
00101
00102
00103 shared_ptr<Psana::CsPad::ConfigV1> config1 = env.configStore().get(src, &m_src);
00104 if (config1.get()) {
00105 for (int i = 0; i < MaxQuads; ++i) {
00106 m_segMask[i] = config1->asicMask()==1 ? 0x3 : 0xff;
00107 }
00108 ++ count;
00109 }
00110
00111 shared_ptr<Psana::CsPad::ConfigV2> config2 = env.configStore().get(src, &m_src);
00112 if (config2.get()) {
00113 for (int i = 0; i < MaxQuads; ++i) {
00114 m_segMask[i] = config2->roiMask(i);
00115 }
00116 ++ count;
00117 }
00118
00119 shared_ptr<Psana::CsPad::ConfigV3> config3 = env.configStore().get(src, &m_src);
00120 if (config3.get()) {
00121 for (int i = 0; i < MaxQuads; ++i) {
00122 m_segMask[i] = config3->roiMask(i);
00123 }
00124 ++ count;
00125 }
00126
00127 shared_ptr<Psana::CsPad::ConfigV4> config4 = env.configStore().get(src, &m_src);
00128 if (config4.get()) {
00129 for (int i = 0; i < MaxQuads; ++i) {
00130 m_segMask[i] = config4->roiMask(i);
00131 }
00132 ++ count;
00133 }
00134
00135 shared_ptr<Psana::CsPad::ConfigV5> config5 = env.configStore().get(src, &m_src);
00136 if (config5.get()) {
00137 for (int i = 0; i < MaxQuads; ++i) {
00138 m_segMask[i] = config5->roiMask(i);
00139 }
00140 ++ count;
00141 }
00142
00143 if (not count) {
00144 MsgLog(name(), error, "No CsPad configuration objects found, terminating.");
00145 terminate();
00146 return;
00147 }
00148
00149 if (count > 1) {
00150 MsgLog(name(), error, "Multiple CsPad configuration objects found, use more specific source address. Terminating.");
00151 terminate();
00152 return;
00153 }
00154
00155 MsgLog(name(), info, "Found CsPad object with address " << m_src);
00156 if (m_src.level() != Pds::Level::Source) {
00157 MsgLog(name(), error, "Found Cspad configuration object with address not at Source level. Terminating.");
00158 terminate();
00159 return;
00160 }
00161
00162 const Pds::DetInfo& dinfo = static_cast<const Pds::DetInfo&>(m_src);
00163
00164
00165 if (dinfo.device() != Pds::DetInfo::Cspad) {
00166 MsgLog(name(), error, "Found Cspad configuration object with invalid address. Terminating.");
00167 terminate();
00168 return;
00169 }
00170 }
00171
00172
00173
00174 void
00175 CsPadPedestals::event(Event& evt, Env& env)
00176 {
00177
00178
00179
00180 shared_ptr<Psana::CsPad::DataV1> data1 = evt.get(m_src);
00181 if (data1.get()) {
00182
00183 int nQuads = data1->quads_shape()[0];
00184 for (int iq = 0; iq != nQuads; ++ iq) {
00185
00186
00187 const CsPad::ElementV1& quad = data1->quads(iq);
00188
00189
00190 const ndarray<const int16_t, 3>& data = quad.data();
00191 collectStat(quad.quad(), data.data());
00192
00193 ++ m_count[quad.quad()];
00194 }
00195
00196 }
00197
00198 shared_ptr<Psana::CsPad::DataV2> data2 = evt.get(m_src);
00199 if (data2.get()) {
00200
00201 int nQuads = data2->quads_shape()[0];
00202 for (int iq = 0; iq != nQuads; ++ iq) {
00203
00204
00205 const CsPad::ElementV2& quad = data2->quads(iq);
00206
00207
00208 const ndarray<const int16_t, 3>& data = quad.data();
00209 collectStat(quad.quad(), data.data());
00210
00211 ++ m_count[quad.quad()];
00212 }
00213
00214 }
00215
00216 }
00217
00218
00219
00220
00221 void
00222 CsPadPedestals::endJob(Event& evt, Env& env)
00223 {
00224
00225 MsgLog(name(), info, "collected total " << m_count[0] << " events");
00226
00227 if (not m_pedFile.empty()) {
00228
00229
00230 std::ofstream out(m_pedFile.c_str());
00231 for (int iq = 0; iq != MaxQuads; ++ iq) {
00232 for (int is = 0; is != MaxSectors; ++ is) {
00233 for (int ic = 0; ic != NumColumns; ++ ic) {
00234 for (int ir = 0; ir != NumRows; ++ ir) {
00235 double avg = m_count[iq] ? double(m_sum[iq][is][ic][ir]) / m_count[iq] : 0;
00236 out << avg << ' ';
00237 }
00238 out << '\n';
00239 }
00240 }
00241 }
00242
00243 out.close();
00244
00245 }
00246
00247 if (not m_noiseFile.empty()) {
00248
00249
00250 std::ofstream out(m_noiseFile.c_str());
00251 for (int iq = 0; iq != MaxQuads; ++ iq) {
00252 for (int is = 0; is != MaxSectors; ++ is) {
00253 for (int ic = 0; ic != NumColumns; ++ ic) {
00254 for (int ir = 0; ir != NumRows; ++ ir) {
00255 double stdev = 0;
00256 if (m_count[iq] > 1) {
00257 stdev = std::sqrt(::stddev2(m_count[iq], m_sum[iq][is][ic][ir], m_sum2[iq][is][ic][ir]));
00258 }
00259 out << stdev << ' ';
00260 }
00261 out << '\n';
00262 }
00263 }
00264 }
00265
00266 out.close();
00267
00268 }
00269
00270 }
00271
00272
00273
00274 void
00275 CsPadPedestals::collectStat(unsigned qNum, const int16_t* data)
00276 {
00277
00278
00279 int seg = 0;
00280 for (int is = 0; is < MaxSectors; ++ is) {
00281 if (m_segMask[qNum] & (1 << is)) {
00282
00283
00284 int64_t* sum = &m_sum[qNum][is][0][0];
00285 int64_t* sum2 = &m_sum2[qNum][is][0][0];
00286
00287 const int16_t* segData = data + seg*NumColumns*NumRows;
00288
00289
00290 for (int i = 0; i < NumColumns*NumRows ; ++ i) {
00291 sum[i] += segData[i];
00292 sum2[i] += segData[i]*segData[i];
00293 }
00294
00295 ++seg;
00296 }
00297 }
00298
00299 }
00300
00301 }