source: trunk/Mars/mpedestal/MPedestalSubtract.cc@ 13408

Last change on this file since 13408 was 11449, checked in by tbretz, 13 years ago
Improved some output.
File size: 8.0 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, 10/2006 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2008
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MPedestalSubtract
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 MPedestalSubtractedEvt.
32//
33// Input Containers:
34// MRawEvtData
35// MRawRunHeader
36// MPedestalCam
37//
38// Output Containers:
39// MPedestalSubtractedEvt
40//
41//////////////////////////////////////////////////////////////////////////////
42#include "MPedestalSubtract.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 "MPedestalSubtractedEvt.h"
59
60#include "MExtractedSignalCam.h"
61#include "MExtractedSignalPix.h"
62
63ClassImp(MPedestalSubtract);
64
65using namespace std;
66
67const char *MPedestalSubtract::fgNamePedestalCam = "MPedestalCam";
68const char *MPedestalSubtract::fgNamePedestalSubtractedEvt = "MPedestalSubtractedEvt";
69
70// --------------------------------------------------------------------------
71//
72// Default constructor.
73//
74MPedestalSubtract::MPedestalSubtract(const char *name, const char *title)
75 : fRawEvt(NULL), fPedestals(NULL), fSignal(NULL)
76{
77 fName = name ? name : "MPedestalSubtract";
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// - MPedestalSubtractedEvt
92//
93Int_t MPedestalSubtract::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 fSignal = (MPedestalSubtractedEvt*)pList->FindCreateObj("MPedestalSubtractedEvt");//, AddSerialNumber(fNamePedestalSubtractedEvt));
103 if (!fSignal)
104 return kFALSE;
105
106 if (fPedestals)
107 {
108 *fLog << inf << "Pedestals given by pointer will be subtracted." << endl;
109 return kTRUE;
110 }
111
112 if (fNamePedestalCam.IsNull())
113 {
114 *fLog << inf << "No name for MPedestalCam given, pedestal subtraction will be skipped." << endl;
115 return kTRUE;
116 }
117
118 fPedestals = (MPedestalCam*)pList->FindObject(AddSerialNumber(fNamePedestalCam), "MPedestalCam");
119 if (!fPedestals)
120 {
121 *fLog << err << AddSerialNumber(fNamePedestalCam) << " [MPedestalCam] not found... aborting" << endl;
122 return kFALSE;
123 }
124
125 *fLog << inf << "Pedestals " << fNamePedestalCam << " will be subtracted." << endl;
126
127 return kTRUE;
128}
129
130void MPedestalSubtract::Memcpy(void *dest, void *src, Int_t cnt) const
131{
132 if (fRawEvt->GetNumBytesPerSample()==2)
133 memcpy(dest, src, cnt*2);
134 else
135 {
136 const Byte_t *b = (Byte_t*)src;
137 for (USample_t *ptr=(USample_t*)dest; ptr<(USample_t*)dest+cnt; ptr++)
138 *ptr = *b++;
139 }
140}
141
142Bool_t MPedestalSubtract::ReInit(MParList *pList)
143{
144 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
145 if (!fRunHeader)
146 {
147 *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
148 return kFALSE;
149 }
150 return kTRUE;
151}
152
153// --------------------------------------------------------------------------
154//
155//
156Int_t MPedestalSubtract::Process()
157{
158 // Please note:
159 // - for data with only hi-gain samples numl is 0
160 // - for data with hi- and lo-gain samples
161 // numl is 0 if read from a raw-data file or a new MC root-file(?)
162 // numl is not 0 if read from an old MC root-file
163
164 // Total number of samples
165 const Int_t numh = fRawEvt->GetNumHiGainSamples();
166 const Int_t numl = fRawEvt->GetNumLoGainSamples();
167
168 // Check if event is empty (presumably MC event -- sanity check)
169 if (numh+numl==0)
170 return kCONTINUE;
171
172 // Check for consistency (our simulation can do weird things!)
173 if (numh+numl!=fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain())
174 {
175 *fLog << err << "MPedestalSubtract::Process: ERROR - Number of samples in event ";
176 *fLog << "(hi+lo=" << numh+numl << ")" << endl << " doesn't match number in run-header ";
177 *fLog << "(" << fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain() << ")" << endl;
178 *fLog << " In case you are processing real data you have just found a bug." << endl;
179 *fLog << " If you process Monte Carlo data, it means probably that the length" << endl;
180 *fLog << " of the arrays in MRawEvtData are inconsistent with the run-header." << endl;
181 return kERROR;
182 }
183
184 // Get scale between FADC units and 256 ;-)
185 const UInt_t scale = fRawEvt->GetScale();
186
187 // initialize fSignal
188 fSignal->InitSamples(numh+numl);//, fRawEvt->GetNumPixels(), numh+numl);
189
190 // iterate over all pixels
191 MRawEvtPixelIter pixel(fRawEvt);
192 while (pixel.Next())
193 {
194 // Get index ofthis pixel
195 const Int_t pixidx = pixel.GetPixelId();
196
197 if (pixidx>=fSignal->GetNumPixels())
198 {
199 *fLog << err << "ERROR - Pixel index " << pixidx << " out of bounds... abort." << endl;
200 *fLog << " Index from MRawEvtPixelIter is " << pixidx << ", max in MPedestalSubtractedEvt is " << fSignal->GetNumPixels() << endl;
201 return kERROR;
202 }
203 // Get pointer were to store merged raw data
204 USample_t *sample = fSignal->GetSamplesRaw(pixidx);
205
206 // copy hi- and lo-gains samples together
207 Memcpy(sample, pixel.GetHiGainSamples(), numh);
208 Memcpy(sample+numh, pixel.GetLoGainSamples(), numl);
209
210 // start of destination array, end of hi-gain destination array
211 // and start of hi-gain samples
212 Float_t *beg = fSignal->GetSamples(pixidx);
213 Float_t *end = beg + fSignal->GetNumSamples();
214
215 const USample_t *src = sample;
216
217 // if no pedestals are given just convert the data into
218 // floats and we are finished
219 if (!fPedestals)
220 {
221 while (beg<end)
222 *beg++ = *src++;//Float_t(*src++)/scale;
223 continue;
224 }
225
226 // get pedestal information for this pixel
227 const MPedestalPix &pedpix = (*fPedestals)[pixidx];
228
229 // pedestal information
230 const Int_t ab = pixel.HasABFlag() ? 1 : 0;
231 const Float_t ped = pedpix.GetPedestal();
232
233 // determine with which pedestal (+/- AB offset) to start
234 const Bool_t swap = (ab&1)==1;
235 const Float_t offh = swap ? -pedpix.GetPedestalABoffset() : pedpix.GetPedestalABoffset();
236 const Float_t mean[2] = { ped + offh, ped - offh };
237
238 // Copy hi-gains into array and substract pedestal
239 // FIXME: Shell we really subtract the pedestal from saturating slices???
240 for (Float_t *ptr=beg; ptr<end; ptr++)
241 *ptr = Float_t(*src++)/scale - mean[(ptr-beg)&1];
242 }
243
244 return kTRUE;
245}
Note: See TracBrowser for help on using the repository browser.