source: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc@ 6461

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