source: branches/Mars_McMismatchStudy/mbadpixels/MBadPixelsCalc.cc@ 20115

Last change on this file since 20115 was 9374, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 10.5 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, 02/2004 <mailto:tbretz@astro.uni.wuerzburg.de>
19! Author(s): Stefan Ruegamer, 08/2005 <mailto:snruegam@astro.uni.wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2005
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MBadPixelsCalc
29//
30//
31// The job of the task is to determin bad pixels event-wise. This must be
32// redone for each event. This particular task is for what is explained
33// below.
34// New algorithms may enter new tasks.
35//
36//
37// Check the pedestal RMS of every pixel with respect to the mean
38// pedestal RMS of the camera.
39//
40// The pixels can be set as blind if the pedestalRMS is too big or 0.
41//
42// If you don't want to use this option set the PedestalLevel<=0;
43//
44// MBadPixelsCalc calc;
45// calc.SetPedestalLevel(-1);
46//
47//
48// Input Containers:
49// [MPedPhotCam]
50// [MGeomCam]
51//
52// Output Containers:
53// MBadPixels
54//
55/////////////////////////////////////////////////////////////////////////////
56#include "MBadPixelsCalc.h"
57
58#include <TEnv.h>
59
60#include "MArrayI.h"
61#include "MArrayD.h"
62
63#include "MLog.h"
64#include "MLogManip.h"
65
66#include "MParList.h"
67
68#include "MGeomCam.h"
69#include "MGeom.h"
70
71#include "MPedPhotCam.h"
72#include "MPedPhotPix.h"
73
74#include "MBadPixelsCam.h"
75#include "MBadPixelsPix.h"
76
77ClassImp(MBadPixelsCalc);
78
79using namespace std;
80
81static const TString gsDefName = "MBadPixelsCalc";
82static const TString gsDefTitle = "Find hot spots (star, broken pixels, etc)";
83
84// --------------------------------------------------------------------------
85//
86// Default constructor.
87//
88MBadPixelsCalc::MBadPixelsCalc(const char *name, const char *title)
89 : fPedestalLevel(3), fPedestalLevelVarianceLo(5),
90 fPedestalLevelVarianceHi(5), fNamePedPhotCam("MPedPhotCam"),
91 fCheckInProcess(kTRUE), fCheckInPostProcess(kFALSE)
92{
93 fName = name ? name : gsDefName.Data();
94 fTitle = title ? title : gsDefTitle.Data();
95}
96
97// --------------------------------------------------------------------------
98//
99//
100Int_t MBadPixelsCalc::PreProcess (MParList *pList)
101{
102 fBadPixels = (MBadPixelsCam*)pList->FindCreateObj(AddSerialNumber("MBadPixelsCam"));
103 if (!fBadPixels)
104 return kFALSE;
105
106 if (fPedestalLevel>0)
107 {
108 fPedPhotCam = (MPedPhotCam*)pList->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam");
109 if (!fPedPhotCam)
110 {
111 *fLog << err << fNamePedPhotCam << "[MPedPhotCam] not found... aborting." << endl;
112 return kFALSE;
113 }
114
115 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
116 if (!fGeomCam)
117 {
118 *fLog << err << "MGeomCam not found... aborting." << endl;
119 return kFALSE;
120 }
121 }
122
123 *fLog << inf << "Name of MPedPhotCam to get pedestal rms from: " << fNamePedPhotCam << endl;
124 if (fPedestalLevel>0)
125 *fLog << "Checking mean 'pedestal rms' against absolute value with level " << fPedestalLevel << endl;
126 if (fPedestalLevelVarianceLo>0)
127 *fLog << "Checking mean 'pedestal rms' against its lower variance with level " << fPedestalLevelVarianceLo << endl;
128 if (fPedestalLevelVarianceHi>0)
129 *fLog << "Checking mean 'pedestal rms' against its upper variance with level " << fPedestalLevelVarianceHi << endl;
130
131 return kTRUE;
132}
133
134// --------------------------------------------------------------------------
135//
136// Check the pedestal Rms of the pixels: compute with 2 iterations the mean
137// for inner and outer pixels. Set as blind the pixels with too small or
138// too high pedestal Rms with respect to the mean.
139//
140Bool_t MBadPixelsCalc::CheckPedestalRms(MBadPixelsPix::UnsuitableType_t type) const
141{
142 const Bool_t checklo = fPedestalLevelVarianceLo>0;
143 const Bool_t checkhi = fPedestalLevelVarianceHi>0;
144
145 if (fPedestalLevel<=0 && !checklo && !checkhi)
146 return kTRUE;
147
148 if (!fGeomCam || !fPedPhotCam || !fBadPixels)
149 {
150 *fLog << err << "MBadPixelsCalc::CheckPedestalRms: ERROR - One of the necessary container are not initialized..." << endl;
151 return kFALSE;
152 }
153
154 const Int_t entries = fPedPhotCam->GetSize();
155
156 const Int_t na = fGeomCam->GetNumAreas();
157
158 MArrayD meanrms(na);
159 MArrayI npix(na);
160
161 for (Int_t i=0; i<entries; i++)
162 {
163 if ((*fBadPixels)[i].IsUnsuitable())
164 continue;
165
166 const Double_t rms = (*fPedPhotCam)[i].GetRms();
167
168 if (rms<=0 || rms>=200*fGeomCam->GetPixRatioSqrt(i))
169 continue;
170
171 const Byte_t aidx = (*fGeomCam)[i].GetAidx();
172
173 meanrms[aidx] += rms;
174 npix[aidx]++;
175 }
176
177 //if no pixel has a minimum signal, return
178 Int_t counter=0;
179 for (int i=0; i<na; i++)
180 {
181 if (npix[i]==0 || meanrms[i]==0) //no valid Pedestals Rms
182 {
183 counter++;
184 continue;
185 }
186
187 meanrms[i] /= npix[i];
188 npix[i]=0;
189 }
190
191 if (counter==na)
192 {
193 *fLog << err << "MBadPixelsCalc::CheckPedestalRms: ERROR - No pixel seems to contain a valid pedestal RMS..." << endl;
194 return kFALSE;
195 }
196
197 MArrayD meanrms2(na);
198 MArrayD varrms2(na);
199 for (Int_t i=0; i<entries; i++)
200 {
201 if ((*fBadPixels)[i].IsUnsuitable())
202 continue;
203
204 const Double_t rms = (*fPedPhotCam)[i].GetRms();
205 const Byte_t aidx = (*fGeomCam)[i].GetAidx();
206
207 //Calculate the corrected means:
208 if (rms<=meanrms[aidx]/1.5 || rms>=meanrms[aidx]*1.5)
209 continue;
210
211 meanrms2[aidx] += rms;
212 varrms2 [aidx] += rms*rms;
213 npix[aidx]++;
214 }
215
216 //if no pixel has a minimum signal, return
217 MArrayD lolim1(na), lolim2(na); // Precalcualtion of limits
218 MArrayD uplim1(na), uplim2(na); // for speeed reasons
219 counter = 0;
220 for (int i=0; i<na; i++)
221 {
222 if (npix[i]==0 || meanrms2[i]==0)
223 {
224 counter++;
225 continue;
226 }
227
228 meanrms2[i] /= npix[i];
229
230 if (fPedestalLevel>0)
231 {
232 lolim1[i] = meanrms2[i]/fPedestalLevel;
233 uplim1[i] = meanrms2[i]*fPedestalLevel;
234 }
235
236 if (checklo || checkhi)
237 {
238 varrms2[i] /= npix[i];
239 varrms2[i] = TMath::Sqrt(varrms2[i]-meanrms2[i]*meanrms2[i]);
240
241 lolim2[i] = meanrms2[i]-fPedestalLevelVarianceLo*varrms2[i];
242 uplim2[i] = meanrms2[i]+fPedestalLevelVarianceHi*varrms2[i];
243 }
244 }
245
246 if (counter==na)
247 {
248 *fLog << err << "MBadPixelsCalc::CheckPedestalRms: ERROR - No pixel seems to contain a valid pedestal RMS anymore..." << endl;
249 return kFALSE;
250 }
251
252 Int_t bads = 0;
253
254 //Blind the Bad Pixels
255 for (Int_t i=0; i<entries; i++)
256 {
257 if ((*fBadPixels)[i].IsUnsuitable())
258 continue;
259
260 const Double_t rms = (*fPedPhotCam)[i].GetRms();
261 const Byte_t aidx = (*fGeomCam)[i].GetAidx();
262
263 if ((fPedestalLevel<=0 || (rms> lolim1[aidx] && rms<=uplim1[aidx])) &&
264 (!checklo || rms> lolim2[aidx]) &&
265 (!checkhi || rms<=uplim2[aidx])
266 )
267 continue;
268
269 (*fBadPixels)[i].SetUnsuitable(type);
270 if (type==MBadPixelsPix::kUnsuitableRun)
271 (*fBadPixels)[i].SetUncalibrated(MBadPixelsPix::kDeadPedestalRms);
272
273 bads++;
274 }
275
276 // Check if the number of pixels to blind is > 60% of total number of pixels
277 //
278 if (bads>0.5*entries)
279 {
280 *fLog << err << "ERROR - More than 50% unsuitable pixels (" << bads << "/" << entries << ")... something must be wrong!" << endl;
281 return kFALSE;
282 }
283
284 return kTRUE;
285}
286
287Bool_t MBadPixelsCalc::CheckPedestalRms(MBadPixelsCam &cam, const MPedPhotCam &ped, MBadPixelsPix::UnsuitableType_t t)
288{
289 MBadPixelsCam *store1 = fBadPixels;
290 const MPedPhotCam *store2 = fPedPhotCam;
291
292 fBadPixels = &cam;
293 fPedPhotCam = &ped;
294
295 const Bool_t rc = CheckPedestalRms(t);
296
297 fBadPixels = store1;
298 fPedPhotCam = store2;
299
300 return rc;
301}
302
303Bool_t MBadPixelsCalc::CheckPedestalRms(MBadPixelsCam &cam, const MPedestalCam &ped, MBadPixelsPix::UnsuitableType_t t)
304{
305 return CheckPedestalRms(cam, MPedPhotCam(ped), t);
306}
307
308// --------------------------------------------------------------------------
309//
310//
311Int_t MBadPixelsCalc::Process()
312{
313 if (!fCheckInProcess)
314 return kTRUE;
315
316 return CheckPedestalRms(MBadPixelsPix::kUnsuitableEvt) ? kTRUE : kERROR;
317}
318
319// --------------------------------------------------------------------------
320//
321//
322Int_t MBadPixelsCalc::PostProcess()
323{
324 return fCheckInPostProcess ? CheckPedestalRms(MBadPixelsPix::kUnsuitableRun) : kTRUE;
325}
326
327// --------------------------------------------------------------------------
328//
329// Read the setup from a TEnv, eg:
330// MBadPixelsCalc.PedestalLevel: 3.0
331//
332// MBadPixelsCalc.PedestalLevelVariance: 5.0
333// overwrites
334// MBadPixelsCalc.PedestalLevelVarianceLo: 5.0
335// and
336// MBadPixelsCalc.PedestalLevelVarianceHi: 5.0
337//
338Int_t MBadPixelsCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
339{
340 Bool_t rc = kFALSE;
341 if (IsEnvDefined(env, prefix, "PedestalLevel", print))
342 {
343 rc = kTRUE;
344 SetPedestalLevel(GetEnvValue(env, prefix, "PedestalLevel", fPedestalLevel));
345 }
346 if (IsEnvDefined(env, prefix, "PedestalLevelVarianceLo", print))
347 {
348 rc = kTRUE;
349 SetPedestalLevelVarianceLo(GetEnvValue(env, prefix, "PedestalLevelVarianceLo", fPedestalLevelVarianceLo));
350 }
351 if (IsEnvDefined(env, prefix, "PedestalLevelVarianceHi", print))
352 {
353 rc = kTRUE;
354 SetPedestalLevelVarianceHi(GetEnvValue(env, prefix, "PedestalLevelVarianceHi", fPedestalLevelVarianceHi));
355 }
356
357 if (IsEnvDefined(env, prefix, "PedestalLevelVariance", print))
358 {
359 rc = kTRUE;
360 SetPedestalLevelVariance(GetEnvValue(env, prefix, "PedestalLevelVariance", -1));
361 }
362 return rc;
363}
Note: See TracBrowser for help on using the repository browser.