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, 10/2006 <mailto:tbretz@astro.uni-wuerzburg.de>
|
---|
19 | !
|
---|
20 | ! Copyright: MAGIC Software Development, 2000-2008
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 | //////////////////////////////////////////////////////////////////////////////
|
---|
26 | //
|
---|
27 | // MDrsCalibApply
|
---|
28 | //
|
---|
29 | // This class merges hi- and lo-gain samples into one array and
|
---|
30 | // subtracts the pedestal (including the AB-offset) from the
|
---|
31 | // data and stores the result in MDrsCalibApplyedEvt.
|
---|
32 | //
|
---|
33 | // Input Containers:
|
---|
34 | // MRawEvtData
|
---|
35 | // MRawRunHeader
|
---|
36 | // MPedestalCam
|
---|
37 | //
|
---|
38 | // Output Containers:
|
---|
39 | // MDrsCalibApplyedEvt
|
---|
40 | //
|
---|
41 | //////////////////////////////////////////////////////////////////////////////
|
---|
42 | #include "MDrsCalibApply.h"
|
---|
43 |
|
---|
44 | #include "MLog.h"
|
---|
45 | #include "MLogManip.h"
|
---|
46 |
|
---|
47 | #include "MParList.h"
|
---|
48 |
|
---|
49 | #include "MArrayB.h"
|
---|
50 |
|
---|
51 | #include "MRawRunHeader.h"
|
---|
52 | #include "MRawEvtData.h"
|
---|
53 | #include "MRawEvtPixelIter.h"
|
---|
54 |
|
---|
55 | #include "MPedestalCam.h"
|
---|
56 | #include "MPedestalPix.h"
|
---|
57 |
|
---|
58 | #include "../mhist/MHDrsCalib.h"
|
---|
59 | #include "MPedestalSubtractedEvt.h"
|
---|
60 |
|
---|
61 | #include "MExtractedSignalCam.h"
|
---|
62 | #include "MExtractedSignalPix.h"
|
---|
63 |
|
---|
64 | ClassImp(MDrsCalibResult);
|
---|
65 | ClassImp(MDrsCalibApply);
|
---|
66 |
|
---|
67 | using namespace std;
|
---|
68 |
|
---|
69 | // --------------------------------------------------------------------------
|
---|
70 | //
|
---|
71 | // Default constructor.
|
---|
72 | //
|
---|
73 | MDrsCalibApply::MDrsCalibApply(const char *name, const char *title)
|
---|
74 | : fRawEvt(NULL), fDrsCalib(NULL), fSignal(NULL), fResult(0),
|
---|
75 | fMaxNumPrevEvents(5), fRemoveSpikes(3), fSlidingAverage(0)
|
---|
76 | {
|
---|
77 | fName = name ? name : "MDrsCalibApply";
|
---|
78 | fTitle = title ? title : "Class to subtract pedestal";
|
---|
79 | }
|
---|
80 |
|
---|
81 | // --------------------------------------------------------------------------
|
---|
82 | //
|
---|
83 | // The PreProcess searches for the following input containers:
|
---|
84 | // - MRawEvtData
|
---|
85 | // - MRawRunHeader
|
---|
86 | // - MPedestalCam
|
---|
87 | //
|
---|
88 | // The following output containers are also searched and created if
|
---|
89 | // they were not found:
|
---|
90 | //
|
---|
91 | // - MDrsCalibApplyedEvt
|
---|
92 | //
|
---|
93 | Int_t MDrsCalibApply::PreProcess(MParList *pList)
|
---|
94 | {
|
---|
95 | fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData"));
|
---|
96 | if (!fRawEvt)
|
---|
97 | {
|
---|
98 | *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl;
|
---|
99 | return kFALSE;
|
---|
100 | }
|
---|
101 |
|
---|
102 | fDrsCalib = (MDrsCalibration*)pList->FindObject("MDrsCalibration");
|
---|
103 | if (!fDrsCalib)
|
---|
104 | *fLog << warn << "No DRS calibration container (MDrsCalibration) found... skipping." << endl;
|
---|
105 |
|
---|
106 | else
|
---|
107 | *fLog << inf << "DRS calibration will be applied." << endl;
|
---|
108 |
|
---|
109 | fResult = 0;
|
---|
110 | if (fMaxNumPrevEvents>0)
|
---|
111 | {
|
---|
112 | fPrevStart.clear();
|
---|
113 | fResult = (MDrsCalibResult*)pList->FindCreateObj("MDrsCalibResult");
|
---|
114 | if (!fResult)
|
---|
115 | return kFALSE;
|
---|
116 | }
|
---|
117 |
|
---|
118 | fSignal = (MPedestalSubtractedEvt*)pList->FindCreateObj("MPedestalSubtractedEvt");//, AddSerialNumber(fNamePedestalSubtractedEvt));
|
---|
119 | if (!fSignal)
|
---|
120 | return kFALSE;
|
---|
121 |
|
---|
122 | return kTRUE;
|
---|
123 | }
|
---|
124 |
|
---|
125 | Bool_t MDrsCalibApply::ReInit(MParList *pList)
|
---|
126 | {
|
---|
127 | fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
|
---|
128 | if (!fRunHeader)
|
---|
129 | {
|
---|
130 | *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
|
---|
131 | return kFALSE;
|
---|
132 | }
|
---|
133 |
|
---|
134 | if (!fDrsCalib)
|
---|
135 | return kTRUE;
|
---|
136 |
|
---|
137 | if (fDrsCalib->fRoi>0 && fDrsCalib->fRoi != fRunHeader->GetNumSamplesHiGain())
|
---|
138 | {
|
---|
139 | *fLog << warn << "Slice number of DRS calibration " << fDrsCalib->fRoi << " and run-header " << fRunHeader->GetNumSamplesHiGain() << " don't match...\n";
|
---|
140 | *fLog << "secondary baseline correction will not be applied." << endl;
|
---|
141 | return kFALSE;
|
---|
142 | }
|
---|
143 |
|
---|
144 | return kTRUE;
|
---|
145 | }
|
---|
146 |
|
---|
147 | // --------------------------------------------------------------------------
|
---|
148 | //
|
---|
149 | //
|
---|
150 | Int_t MDrsCalibApply::Process()
|
---|
151 | {
|
---|
152 | // Please note:
|
---|
153 | // - for data with only hi-gain samples numl is 0
|
---|
154 | // - for data with hi- and lo-gain samples
|
---|
155 | // numl is 0 if read from a raw-data file or a new MC root-file(?)
|
---|
156 | // numl is not 0 if read from an old MC root-file
|
---|
157 |
|
---|
158 | // Total number of samples
|
---|
159 | const Int_t roi = fRawEvt->GetNumHiGainSamples();
|
---|
160 |
|
---|
161 | // Check if event is empty (presumably MC event -- sanity check)
|
---|
162 | if (roi==0)
|
---|
163 | return kCONTINUE;
|
---|
164 |
|
---|
165 | // Check for consistency (our simulation can do weird things!)
|
---|
166 | if (roi!=fRunHeader->GetNumSamplesHiGain())
|
---|
167 | {
|
---|
168 | *fLog << err << "MDrsCalibApply::Process: ERROR - Number of samples in event ";
|
---|
169 | *fLog << "(hi=" << roi << ")" << endl << " doesn't match number in run-header ";
|
---|
170 | *fLog << "(" << fRunHeader->GetNumSamplesHiGain() << ")" << endl;
|
---|
171 | *fLog << " In case you are processing real data you have just found a bug." << endl;
|
---|
172 | *fLog << " If you process Monte Carlo data, it means probably that the length" << endl;
|
---|
173 | *fLog << " of the arrays in MRawEvtData are inconsistent with the run-header." << endl;
|
---|
174 | return kERROR;
|
---|
175 | }
|
---|
176 |
|
---|
177 | if (fRawEvt->GetNumLoGainSamples()+fRunHeader->GetNumSamplesLoGain()>0)
|
---|
178 | {
|
---|
179 | *fLog << err << "MDrsCalibApply::Process: ERROR - Data must not contain lo gain slices." << endl;
|
---|
180 | return kERROR;
|
---|
181 | }
|
---|
182 |
|
---|
183 | // initialize fSignal
|
---|
184 | fSignal->InitSamples(roi);
|
---|
185 |
|
---|
186 | if (!fDrsCalib) // FIXME: Do a fake (none) calibration
|
---|
187 | return kTRUE;
|
---|
188 |
|
---|
189 | const int16_t *val = reinterpret_cast<int16_t*>(fRawEvt->GetSamples());
|
---|
190 |
|
---|
191 | const int32_t *offset = fDrsCalib->fOffset.data();
|
---|
192 | const int64_t *gain = fDrsCalib->fGain.data();
|
---|
193 | const int64_t *trgoff = fDrsCalib->fTrgOff.data();
|
---|
194 |
|
---|
195 | const uint64_t scaleabs = fDrsCalib->fNumOffset;
|
---|
196 | const uint64_t scalegain = fDrsCalib->fNumGain;
|
---|
197 | const uint64_t scalerel = fDrsCalib->fNumTrgOff;
|
---|
198 |
|
---|
199 | const int16_t *start = reinterpret_cast<int16_t*>(fRawEvt->GetStartCells());
|
---|
200 |
|
---|
201 | Float_t *vec = fSignal->GetSamples(0);
|
---|
202 |
|
---|
203 | const UShort_t *idx = fRawEvt->GetPixelIds();
|
---|
204 | for (size_t ch=0; ch<fRawEvt->GetNumPixels(); ch++)
|
---|
205 | {
|
---|
206 | const size_t drs = ch*1024;
|
---|
207 | const size_t hw = ch*roi;
|
---|
208 | const size_t sw = idx[ch]*roi;
|
---|
209 |
|
---|
210 | DrsCalibrate::ApplyCh(vec+sw, val+hw, start[ch], roi,
|
---|
211 | offset+drs, scaleabs,
|
---|
212 | gain +drs, scalegain,
|
---|
213 | trgoff+hw, scalerel);
|
---|
214 | }
|
---|
215 |
|
---|
216 | if (fResult)
|
---|
217 | {
|
---|
218 | fResult->fData.resize(fPrevStart.size()*2);
|
---|
219 |
|
---|
220 | deque<vector<Short_t>>::const_iterator it = fPrevStart.begin();
|
---|
221 |
|
---|
222 | for (size_t i=0; i<fPrevStart.size()*2; i++)
|
---|
223 | {
|
---|
224 | fResult->fData[i] = DrsCalibrate::Step();
|
---|
225 |
|
---|
226 | const int32_t off = i%2 ? 3 : roi+10;
|
---|
227 |
|
---|
228 | fResult->fData[i] =
|
---|
229 | DrsCalibrate::CorrectStep(vec, fRawEvt->GetNumPixels(), roi,
|
---|
230 | it->data(), start, off, idx);
|
---|
231 |
|
---|
232 | if (i%2==1)
|
---|
233 | it++;
|
---|
234 | }
|
---|
235 |
|
---|
236 | fPrevStart.push_front(vector<Short_t>(start, start+fRawEvt->GetNumPixels()));
|
---|
237 | if (fPrevStart.size()>fMaxNumPrevEvents)
|
---|
238 | fPrevStart.pop_back();
|
---|
239 | }
|
---|
240 |
|
---|
241 | switch (fRemoveSpikes)
|
---|
242 | {
|
---|
243 | case 1:
|
---|
244 | DrsCalibrate::RemoveSpikes(vec, roi);
|
---|
245 | break;
|
---|
246 | case 2:
|
---|
247 | DrsCalibrate::RemoveSpikes2(vec, roi);
|
---|
248 | break;
|
---|
249 | case 3:
|
---|
250 | DrsCalibrate::RemoveSpikes3(vec, roi);
|
---|
251 | break;
|
---|
252 | }
|
---|
253 |
|
---|
254 |
|
---|
255 | if (fSlidingAverage)
|
---|
256 | DrsCalibrate::SlidingAverage(vec, roi, fSlidingAverage);
|
---|
257 |
|
---|
258 |
|
---|
259 | return kTRUE;
|
---|
260 | }
|
---|