source: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.cc@ 3490

Last change on this file since 3490 was 3490, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 15.2 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): Oscar Blanch 12/2001 <mailto:blanch@ifae.es>
19! Author(s): Thomas Bretz 08/2002 <mailto:tbretz@astro.uni.wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MBadPixelsTreat
29//
30// You can use MBadPixelsTreat::SetUseInterpolation to replaced the
31// bad pixels by the average of its neighbors instead of unmapping
32// them. If you want to include the central pixel use
33// MBadPixelsTreat::SetUseCentralPixel. The bad pixels are taken from
34// an existing MBadPixelsCam.
35//
36// It check if there are enough neighbors to calculate the mean
37// If not, unmap the pixel. The minimum number of good neighbors
38// should be set using SetNumMinNeighbors
39//
40// Input Containers:
41// MCerPhotEvt
42// MPedPhotCam
43// MBadPixelsCam
44// [MGeomCam]
45//
46// Output Containers:
47// MCerPhotEvt
48//
49/////////////////////////////////////////////////////////////////////////////
50#include "MBadPixelsTreat.h"
51
52#include <fstream>
53
54#include <TArrayD.h>
55
56#include "MLog.h"
57#include "MLogManip.h"
58
59#include "MParList.h"
60
61#include "MGeomPix.h"
62#include "MGeomCam.h"
63
64#include "MCerPhotPix.h"
65#include "MCerPhotEvt.h"
66
67#include "MPedPhotPix.h"
68#include "MPedPhotCam.h"
69
70#include "MBadPixelsPix.h"
71#include "MBadPixelsCam.h"
72
73ClassImp(MBadPixelsTreat);
74
75using namespace std;
76
77static const TString gsDefName = "MBadPixelsTreat";
78static const TString gsDefTitle = "Task to treat bad pixels (interpolation, unmapping)";
79
80// --------------------------------------------------------------------------
81//
82// Default constructor.
83//
84MBadPixelsTreat::MBadPixelsTreat(const char *name, const char *title)
85 : fFlags(0), fNumMinNeighbors(3)
86{
87 fName = name ? name : gsDefName.Data();
88 fTitle = title ? title : gsDefTitle.Data();
89}
90
91// --------------------------------------------------------------------------
92//
93// - Try to find or create MBlindPixels in parameter list.
94// - get the MCerPhotEvt from the parlist (abort if missing)
95// - if no pixels are given by the user try to determin the starfield
96// from the monte carlo run header.
97//
98Int_t MBadPixelsTreat::PreProcess (MParList *pList)
99{
100 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
101 if (!fBadPixels)
102 {
103 *fLog << err << "MBadPixelsCam not found... aborting." << endl;
104 return kFALSE;
105 }
106
107 fPedPhot = (MPedPhotCam*)pList->FindObject(AddSerialNumber("MPedPhotCam"));
108 if (!fPedPhot)
109 {
110 *fLog << err << "MPedPhotCam not found... aborting." << endl;
111 return kFALSE;
112 }
113
114 fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
115 if (!fEvt)
116 {
117 *fLog << err << "MCerPhotEvt not found... aborting." << endl;
118 return kFALSE;
119 }
120
121 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
122 if (!fGeomCam && TESTBIT(fFlags, kUseInterpolation))
123 {
124 *fLog << err << "MGeomCam not found... can't use interpolation." << endl;
125 return kFALSE;
126 }
127
128 return kTRUE;
129}
130
131// --------------------------------------------------------------------------
132//
133// Replaces each pixel (signal, signal error, pedestal, pedestal rms)
134// by the average of its surrounding pixels.
135// If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
136// included.
137//
138void MBadPixelsTreat::InterpolateSignal() const
139{
140 const UShort_t entries = fGeomCam->GetNumPixels();
141
142 //
143 // Create arrays
144 //
145 TArrayD nphot(entries);
146 TArrayD perr(entries);
147
148 //
149 // Loop over all pixels
150 //
151 for (UShort_t i=0; i<entries; i++)
152 {
153 MCerPhotPix *pix = fEvt->GetPixById(i);
154
155 //
156 // Check whether pixel with idx i is blind
157 //
158 if (pix && (*fBadPixels)[i].IsOK())
159 continue;
160
161 //
162 // Get a pointer to this pixel. If it is not yet existing
163 // create a new entry for this pixel in MCerPhotEvt
164 //
165 if (!pix)
166 {
167 pix = fEvt->AddPixel(i, 0, 0);
168 (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
169 }
170
171 //
172 // Get the corresponding geometry and pedestal
173 //
174 const MGeomPix &gpix = (*fGeomCam)[i];
175
176 // Do Not-Use-Central-Pixel
177 const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
178
179 Int_t num = nucp ? 0 : 1;
180
181 nphot[i] = nucp ? 0 : pix->GetNumPhotons();
182 perr[i] = nucp ? 0 : Pow2(pix->GetErrorPhot());
183
184 //
185 // The values are rescaled to the small pixels area for the right comparison
186 //
187 const Double_t ratio = fGeomCam->GetPixRatio(i);
188
189 nphot[i] *= ratio;
190 perr[i] *= ratio;
191
192 //
193 // Loop over all its neighbors
194 //
195 const Int_t n = gpix.GetNumNeighbors();
196 for (int j=0; j<n; j++)
197 {
198 const UShort_t nidx = gpix.GetNeighbor(j);
199
200 //
201 // Do not use blind neighbors
202 //
203 if ((*fBadPixels)[i].IsBad())
204 continue;
205
206 //
207 // Check whether the neighbor has a signal stored
208 //
209 const MCerPhotPix *evtpix = fEvt->GetPixById(nidx);
210 if (!evtpix)
211 continue;
212
213 //
214 // Get the geometry for the neighbor
215 //
216 const Double_t nratio = fGeomCam->GetPixRatio(nidx);
217
218 //
219 //The error is calculated as the quadratic sum of the errors
220 //
221 nphot[i] += nratio*evtpix->GetNumPhotons();
222 perr[i] += nratio* Pow2(evtpix->GetErrorPhot());
223
224 num++;
225 }
226
227 //
228 // Now the mean is calculated and the values rescaled back
229 // to the pixel area
230 //
231 nphot[i] /= (num*ratio);
232 perr[i] = TMath::Sqrt(perr[i]/(num*ratio));
233
234 // Check if there are enough neighbors to calculate the mean
235 // If not, unmap the pixel. The maximum number of blind neighbors
236 // should be 2
237 if (num<3)
238 {
239 pix->SetPixelUnmapped();
240 continue;
241 }
242
243 pix->Set(nphot[i], perr[i]);
244 }
245}
246
247// --------------------------------------------------------------------------
248//
249void MBadPixelsTreat::InterpolatePedestals() const
250{
251 const Int_t entries = fPedPhot->GetSize();
252
253 TArrayD ped(entries);
254 TArrayD rms(entries);
255
256 //
257 // Loop over all pixels
258 //
259 for (UShort_t i=0; i<entries; i++)
260 {
261 //
262 // Check whether pixel with idx i is blind
263 //
264 if ((*fBadPixels)[i].IsOK())
265 continue;
266
267 //
268 // Get the corresponding geometry and pedestal
269 //
270 const MGeomPix &gpix = (*fGeomCam)[i];
271 const MPedPhotPix &ppix = (*fPedPhot)[i];
272
273 // Do Not-Use-Central-Pixel
274 const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
275
276 Int_t num = nucp ? 0 : 1;
277
278 ped[i] = nucp ? 0 : ppix.GetMean();
279 rms[i] = nucp ? 0 : Pow2(ppix.GetRms());
280
281 //
282 // The values are rescaled to the small pixels area for the right comparison
283 //
284 const Double_t ratio = fGeomCam->GetPixRatio(i);
285
286 ped[i] *= ratio;
287 rms[i] *= ratio;
288
289 //
290 // Loop over all its neighbors
291 //
292 const Int_t n = gpix.GetNumNeighbors();
293 for (int j=0; j<n; j++)
294 {
295 const UShort_t nidx = gpix.GetNeighbor(j);
296
297 //
298 // Do not use blind neighbors
299 //
300 if ((*fBadPixels)[i].IsBad())
301 continue;
302
303 //
304 // Get the geometry for the neighbor
305 //
306 const Double_t nratio = fGeomCam->GetPixRatio(nidx);
307 const MPedPhotPix &nppix = (*fPedPhot)[nidx];
308
309 //
310 //The error is calculated as the quadratic sum of the errors
311 //
312 ped[i] += nratio*nppix.GetMean();
313 rms[i] += nratio*Pow2(nppix.GetRms());
314
315 num++;
316 }
317
318 // Check if there are enough neighbors to calculate the mean
319 // If not, unmap the pixel. The minimum number of good neighbors
320 // should be fNumMinNeighbors
321 if (num < fNumMinNeighbors)
322 {
323 MCerPhotPix *pix =fEvt->GetPixById(i);
324 if (!pix)
325 pix = fEvt->AddPixel(i, 0, 0);
326 pix->SetPixelUnmapped();
327 continue;
328 }
329
330 //
331 // Now the mean is calculated and the values rescaled back
332 // to the pixel area
333 //
334 ped[i] /= (num*ratio);
335 rms[i] = TMath::Sqrt(rms[i]/(num*ratio));
336
337 (*fPedPhot)[i].Set(ped[i], rms[i]);
338 }
339}
340
341// --------------------------------------------------------------------------
342//
343// Replaces each pixel (signal, signal error, pedestal, pedestal rms)
344// by the average of its surrounding pixels.
345// If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
346// included.
347//
348// NT: Informations about the interpolated pedestals are added.
349// When the option Interpolated is called, the blind pixel with the new
350// values of signal and fluttuation is included in the calculation of
351// the Image Parameters.
352//
353void MBadPixelsTreat::Interpolate() const
354{
355 const UShort_t entries = fGeomCam->GetNumPixels();
356
357 //
358 // Create arrays
359 //
360 TArrayD nphot(entries);
361 TArrayD perr(entries);
362 TArrayD ped(entries);
363 TArrayD pedrms(entries);
364
365 //
366 // Loop over all pixels
367 //
368 for (UShort_t i=0; i<entries; i++)
369 {
370 MCerPhotPix *pix = fEvt->GetPixById(i);
371
372 //
373 // Check whether pixel with idx i is blind
374 //
375 if (pix && (*fBadPixels)[i].IsOK())
376 continue;
377
378 //
379 // Get a pointer to this pixel. If it is not yet existing
380 // create a new entry for this pixel in MCerPhotEvt
381 //
382 if (!pix)
383 {
384 pix = fEvt->AddPixel(i, 0, 0);
385 (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
386 }
387
388 //
389 // Get the corresponding geometry and pedestal
390 //
391 const MGeomPix &gpix = (*fGeomCam)[i];
392 const MPedPhotPix &ppix = (*fPedPhot)[i];
393
394 // Do Not-Use-Central-Pixel
395 const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
396
397 Int_t num = nucp ? 0 : 1;
398
399 nphot[i] = nucp ? 0 : pix->GetNumPhotons();
400 ped[i] = nucp ? 0 : ppix.GetMean();
401 perr[i] = nucp ? 0 : Pow2(pix->GetErrorPhot());
402 pedrms[i] = nucp ? 0 : Pow2(ppix.GetRms());
403
404 //
405 // The values are rescaled to the small pixels area for the right comparison
406 //
407 const Double_t ratio = fGeomCam->GetPixRatio(i);
408
409 if (nucp)
410 {
411 nphot[i] *= ratio;
412 perr[i] *= ratio;
413 ped[i] *= ratio;
414 pedrms[i] *= ratio;
415 }
416
417 //
418 // Loop over all its neighbors
419 //
420 const Int_t n = gpix.GetNumNeighbors();
421 for (int j=0; j<n; j++)
422 {
423 const UShort_t nidx = gpix.GetNeighbor(j);
424
425 //
426 // Do not use blind neighbors
427 //
428 if ((*fBadPixels)[nidx].IsBad())
429 continue;
430
431 //
432 // Check whether the neighbor has a signal stored
433 //
434 const MCerPhotPix *evtpix = fEvt->GetPixById(nidx);
435 if (!evtpix)
436 continue;
437
438 //
439 // Get the geometry for the neighbor
440 //
441 const Double_t nratio = fGeomCam->GetPixRatio(nidx);
442 MPedPhotPix &nppix = (*fPedPhot)[nidx];
443
444 //
445 // The error is calculated as the quadratic sum of the errors
446 //
447 nphot[i] += nratio*evtpix->GetNumPhotons();
448 ped[i] += nratio*nppix.GetMean();
449 perr[i] += nratio*Pow2(evtpix->GetErrorPhot());
450 pedrms[i] += nratio*Pow2(nppix.GetRms());
451
452 num++;
453 }
454
455 if (num<2)
456 {
457 pix->SetPixelUnmapped();
458 nphot[i] = 0;
459 ped[i] = 0;
460 perr[i] = 0;
461 pedrms[i] = 0;
462 continue;
463 }
464
465 //
466 // Now the mean is calculated and the values rescaled back to the pixel area
467 //
468 nphot[i] /= num*ratio;
469 ped[i] /= num*ratio;
470 perr[i] = TMath::Sqrt(perr[i]/(num*ratio));
471 pedrms[i] = TMath::Sqrt(pedrms[i]/(num*ratio));
472
473 }
474
475 //
476 // Now the new pixel values are calculated and can be replaced in
477 // the corresponding containers
478 //
479 for (UShort_t i=0; i<entries; i++)
480 {
481 //
482 // Do not use blind neighbors
483 //
484 if ((*fBadPixels)[i].IsOK())
485 continue;
486
487 //
488 // It must exist, we have created it in the loop before.
489 //
490 fEvt->GetPixById(i)->Set(nphot[i], perr[i]);
491 (*fPedPhot)[i].Set(ped[i], pedrms[i]);
492 }
493}
494
495// --------------------------------------------------------------------------
496//
497// Removes all blind pixels from the analysis by setting their state
498// to unused.
499//
500void MBadPixelsTreat::Unmap() const
501{
502 const UShort_t entries = fEvt->GetNumPixels();
503
504 //
505 // remove the pixels in fPixelsIdx if they are set to be used,
506 // (set them to 'unused' state)
507 //
508 for (UShort_t i=0; i<entries; i++)
509 {
510 MCerPhotPix &pix = (*fEvt)[i];
511
512 if ((*fBadPixels)[pix.GetPixId()].IsBad())
513 pix.SetPixelUnused();
514 }
515}
516
517// --------------------------------------------------------------------------
518//
519// Treat the blind pixels
520//
521Int_t MBadPixelsTreat::Process()
522{
523 /*
524 if (TESTBIT(fFlags, kCheckPedestalRms))
525 {
526 // if the number of blind pixels is too high, do not interpolate
527 if (CheckPedestalRms()==kFALSE)
528 return kTRUE;
529
530 if (TESTBIT(fFlags, kUseInterpolation))
531 InterpolatePedestals();
532 }
533 */
534
535 if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam)
536 Interpolate();
537 else
538 Unmap();
539
540
541 //fErrors[0]++;
542
543 return kTRUE;
544}
545
546void MBadPixelsTreat::StreamPrimitive(ofstream &out) const
547{
548 out << " MBadPixelsTreat " << GetUniqueName();
549 if (fName!=gsDefName || fTitle!=gsDefTitle)
550 {
551 out << "(\"" << fName << "\"";
552 if (fTitle!=gsDefTitle)
553 out << ", \"" << fTitle << "\"";
554 out <<")";
555 }
556 out << ";" << endl;
557
558 if (TESTBIT(fFlags, kUseInterpolation))
559 out << " " << GetUniqueName() << ".SetUseInterpolation();" << endl;
560 if (TESTBIT(fFlags, kUseCentralPixel))
561 out << " " << GetUniqueName() << ".SetUseCentralPixel();" << endl;
562}
Note: See TracBrowser for help on using the repository browser.