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

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