source: branches/removing_cpp11_features/mdrs/MHDrsCalibration.cc

Last change on this file was 16891, checked in by tbretz, 11 years ago
Added a check for underflows, for the case the gains are read as offsets... only for checking purpose.
File size: 10.7 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 2013 <mailto:thomas.bretz@epfl.ch>
19!
20! Copyright: MAGIC Software Development, 2000-2013
21!
22!
23\* ======================================================================== */
24
25///////////////////////////////////////////////////////////////////////
26//
27// MHDrsCalibration
28//
29// This class contains a list of MHFadcPix.
30//
31///////////////////////////////////////////////////////////////////////
32#include "MHDrsCalibration.h"
33
34#include <fstream>
35#include <sstream>
36#include <limits.h>
37
38#include "MLog.h"
39#include "MLogManip.h"
40
41#include "MParList.h"
42#include "MGeomCam.h"
43#include "MBinning.h"
44
45#include "MRawEvtData.h"
46#include "MRawRunHeader.h"
47#include "MHCamera.h"
48#include "MDrsCalibration.h"
49
50#include <TH2.h>
51#include <TMath.h>
52
53ClassImp(MHDrsCalibration);
54
55using namespace std;
56
57// --------------------------------------------------------------------------
58//
59// default constructor
60// creates an a list of histograms for all pixels and both gain channels
61//
62MHDrsCalibration::MHDrsCalibration(const char *name, const char *title) : fRun(0),
63fEvt(0), fResult(0), fBuffer(1440*1024*6+160*1024*2+4),
64fStep(-1), fScale(0), fNumPixels(0), fNumSamples(0)
65{
66 //
67 // set the name and title of this object
68 //
69 fName = name ? name : "MHDrsCalibration";
70 fTitle = title ? title : "Container for ADC spectra histograms";
71
72 for (int i=1024*1440*2+4; i<1024*1440*3+4; i++)
73 fBuffer[i] = 2000./4096; // Set mean to 0.5
74}
75
76// --------------------------------------------------------------------------
77//
78// To setup the object we get the number of pixels from a MGeomCam object
79// in the Parameter list.
80//
81Bool_t MHDrsCalibration::SetupFill(const MParList *pList)
82{
83 fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
84 if (!fRun)
85 {
86 *fLog << err << "MRawRunHeader not found... aborting." << endl;
87 return kFALSE;
88 }
89
90 switch (++fStep)
91 {
92 case 0:
93 SetTitle("DRS average baseline (physical pipeline);;ADC signal [mV];");
94 break;
95 case 1:
96 SetTitle("DRS average gain (baseline subtracted, physical pipeline);;ADC signal [mV];");
97 break;
98 case 2:
99 SetTitle("DRS average trigger offset (baseline subtracted, logical pipeline);;ADC signal [mV];");
100 break;
101 }
102
103 if (!MHCamEvent::SetupFill(pList))
104 return kFALSE;
105
106 fSum->ResetBit(MHCamera::kProfile);
107
108 fData.Reset();
109
110 fScale = 1;
111
112 *fLog << warn << "WARNING --- oulier detection needed!" << endl;
113
114 return kTRUE;
115}
116
117Bool_t MHDrsCalibration::ReInit(MParList *pList)
118{
119/* fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
120 if (!fRun)
121 {
122 *fLog << err << "MRawRunHeader not found... aborting." << endl;
123 return kFALSE;
124 }*/
125
126 fEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
127 if (!fEvt)
128 {
129 *fLog << err << "MRawEvtData not found... aborting." << endl;
130 return kFALSE;
131 }
132
133 if (!fEvt->HasStartCells())
134 {
135 *fLog << err << "MRawEvtData has no start cells stored... aborting." << endl;
136 return kFALSE;
137 }
138
139 // FIXME: check consistency...?
140
141 if (fEvt->GetNumPixels()==0 || fEvt->GetNumSamples()==0)
142 {
143 *fLog << err << "MRawEvtData contains either no pixels or no samples... aborting." << endl;
144 return kFALSE;
145 }
146
147 if (fStep>=3)
148 {
149 *fLog << err << "DRS calibration already finished." << endl;
150 return kFALSE;
151 }
152
153
154 fData.InitSize(fEvt->GetNumPixels(), fEvt->GetNumSamples());
155
156 fResult = (MDrsCalibration*)pList->FindCreateObj("MDrsCalibration");
157 if (!fResult)
158 return kFALSE;
159
160 if (fStep==1 || fStep==2)
161 fScale = fResult->fNumOffset;
162
163 fNumPixels = fEvt->GetNumPixels();
164 fNumSamples = fEvt->GetNumSamples();
165
166 if (fStep<2 && fNumSamples!=1024)
167 {
168 *fLog << err << "Only 1024 samples are supported for step " << fStep << "... file has " << fNumSamples << endl;
169 return kFALSE;
170 }
171
172 return kTRUE;
173}
174
175// --------------------------------------------------------------------------
176Int_t MHDrsCalibration::Fill(const MParContainer *par, const Stat_t w)
177{
178
179 if (fStep==0)
180 {
181 fData.AddRel(reinterpret_cast<int16_t*>(fEvt->GetSamples()),
182 reinterpret_cast<int16_t*>(fEvt->GetStartCells()));
183 }
184 if (fStep==1)
185 {
186 fData.AddRel(reinterpret_cast<int16_t*>(fEvt->GetSamples()),
187 reinterpret_cast<int16_t*>(fEvt->GetStartCells()),
188 fResult->fOffset.data(),
189 fResult->fNumOffset);
190 }
191 if (fStep==2)
192 {
193 fData.AddAbs(reinterpret_cast<int16_t*>(fEvt->GetSamples()),
194 reinterpret_cast<int16_t*>(fEvt->GetStartCells()),
195 fResult->fOffset.data(),
196 fResult->fNumOffset);
197 }
198
199 /*
200 if (fStep==3)
201 {
202 fData.AddAbs(reinterpret_cast<int16_t*>(fEvt->GetSamples()),
203 reinterpret_cast<int16_t*>(fEvt->GetStartCells()),
204 fResult->fOffset.data(),
205 fResult->fNumOffset);
206 }*/
207
208 return kTRUE;
209}
210
211void MHDrsCalibration::InitHistogram()
212{
213 if (!fEvt)
214 return;
215
216 // Be carefull: the contents of fData are not cloned!
217 pair<vector<double>,vector<double> > p = fData.GetSampleStats();
218
219 vector<double> &mean = p.first;
220 vector<double> &rms = p.second;
221
222 if (mean.size()==0)
223 return;
224
225 fSum->Reset();
226
227 for (size_t ch=0; ch<fNumPixels; ch++)
228 {
229 double m = 0;
230 double r = 0;
231 for (size_t i=0; i<fNumSamples; i++)
232 {
233 m += mean[ch*fNumSamples+i];
234 r += rms [ch*fNumSamples+i];
235 }
236
237 fSum->AddBinContent(fEvt->GetPixelIds()[ch]+1, m);
238 fSum->SetBinError(fEvt->GetPixelIds()[ch]+1, r);
239 }
240
241 fSum->Scale(1./fNumSamples/fScale*2000/4096); // 0.5mV/ADC count
242 fSum->SetEntries(fData.GetNumEntries());
243 fSum->SetAllUsed();
244}
245
246template<typename T>
247Bool_t MHDrsCalibration::CopyData(vector<T> &dest) const
248{
249 const vector<int64_t> &sum = fData.GetSum();
250 dest.resize(sum.size());
251
252 for (size_t ch=0; ch<fNumPixels; ch++)
253 for (size_t i=0; i<fNumSamples; i++)
254 {
255 if (sum[ch*fNumSamples+i]>UINT_MAX)
256 {
257 *fLog << err << "ERROR - Sum over slices exceeded maximum allowed range." << endl;
258 return kFALSE;
259 }
260 dest[ch*fNumSamples+i] = sum[ch*fNumSamples+i];
261 }
262
263 return kTRUE;
264}
265
266Bool_t MHDrsCalibration::Finalize()
267{
268 if (!fResult)
269 return kTRUE;
270
271 fResult->fRoi = fNumSamples;
272 fResult->fStep = fStep;
273 fResult->fNumTm = 0;
274
275 // ==========================================================================
276 if (fStep==0)
277 {
278 fResult->fNumOffset = fData.GetNumEntries();
279 if (!CopyData(fResult->fOffset))
280 return kFALSE;
281
282 //--------------------------------------------
283 for (int i=0; i<1024*1440; i++)
284 fResult->fGain[i] = 4096*fData.GetNumEntries();
285
286 // Scale ADC data from 12bit to 2000mV
287 fData.GetSampleStats(fBuffer.data()+4, 2000./4096);
288 reinterpret_cast<uint32_t*>(fBuffer.data())[1] = fRun->GetFileNumber();
289 }
290
291 // ==========================================================================
292 if (fStep==1)
293 {
294 fResult->fNumGain = fData.GetNumEntries();//*fResult->fNumCellOffset;
295 if (!CopyData(fResult->fGain))
296 return kFALSE;
297
298 //---------------------------------------------
299 // DAC: 0..2.5V == 0..65535 2500*50000 625*50000 625*3125
300 // V-mV: 1000 ---------- --------- --------
301 //fNumGain *= 2500*50000; 65536 16384 1024
302 //for (int i=0; i<1024*1440; i++)
303 // fGain[i] *= 65536;
304 fResult->fNumGain *= 1953125;
305 for (unsigned int i=0; i<fNumPixels*fNumSamples; i++)
306 fResult->fGain[i] *= 1024;
307
308 // Scale ADC data from 12bit to 2000mV
309 fData.GetSampleStats(fBuffer.data()+1024*1440*2+4, 2000./4096/fResult->fNumOffset);//0.5);
310 reinterpret_cast<uint32_t*>(fBuffer.data())[2] = fRun->GetFileNumber();
311 }
312
313 // ==========================================================================
314 if (fStep==2)
315 {
316 fResult->fNumTrgOff = fData.GetNumEntries();//*fResult->fNumCellOffset;
317 if (!CopyData(fResult->fTrgOff))
318 return kFALSE;
319
320 //--------------------------------------------
321 // Scale ADC data from 12bit to 2000mV
322 fData.GetSampleStats(fBuffer.data()+1024*1440*4+4, 2000./4096/fResult->fNumOffset);//0.5);
323 reinterpret_cast<uint32_t*>(fBuffer.data())[0] = fNumSamples;
324 reinterpret_cast<uint32_t*>(fBuffer.data())[3] = fRun->GetFileNumber();
325 }
326
327 // --------------------------------------
328
329 InitHistogram();
330
331 if (fOutputPath.IsNull())
332 return kTRUE;
333
334 const std::string fname(Form("%s/%d_%03d.drs.fits", fOutputPath.Data(),
335 fRun->GetRunNumber(), fRun->GetFileNumber()));
336
337 const Bool_t exists = !gSystem->AccessPathName(fname.c_str(), kFileExists);
338 if (exists)
339 {
340 *fLog << err << "File '" << fname << "' already exists." << endl;
341 return kFALSE;
342 }
343
344 const std::string msg = fResult->WriteFitsImp(fname, fBuffer);
345
346 if (msg.empty())
347 {
348 *fLog << inf << "Wrote DRS calibration file " << fname << endl;
349 return kTRUE;
350 }
351
352 *fLog << err << "Error writing to " << fname << ": " << msg << endl;
353 return kFALSE;
354}
355
356void MHDrsCalibration::Paint(Option_t *o)
357{
358 InitHistogram();
359 MHCamEvent::Paint(o);
360}
361
362Short_t MHDrsCalibration::GetNumUnderflows(float lvl) const
363{
364 if (!fResult)
365 return -1;
366
367 short cnt = 0;
368
369 for (int drs=0; drs<160; drs++)
370 {
371 const int p = drs*9*1024;
372
373 double avg = 0;
374 for (int i=0; i<1024*9; i++)
375 avg += fResult->fOffset[p+i];
376
377 avg /= 1024*9*fResult->fNumOffset;
378
379 if (avg/1000<lvl)
380 cnt++;
381 }
382
383 return cnt;
384
385}
Note: See TracBrowser for help on using the repository browser.