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

Last change on this file since 8462 was 8339, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 20.0 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MBadPixelsTreat.cc,v 1.38 2007-03-01 17:42:27 tbretz Exp $
3! --------------------------------------------------------------------------
4!
5! *
6! * This file is part of MARS, the MAGIC Analysis and Reconstruction
7! * Software. It is distributed to you in the hope that it can be a useful
8! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
9! * It is distributed WITHOUT ANY WARRANTY.
10! *
11! * Permission to use, copy, modify and distribute this software and its
12! * documentation for any purpose is hereby granted without fee,
13! * provided that the above copyright notice appear in all copies and
14! * that both that copyright notice and this permission notice appear
15! * in supporting documentation. It is provided "as is" without express
16! * or implied warranty.
17! *
18!
19!
20! Author(s): Oscar Blanch 12/2001 <mailto:blanch@ifae.es>
21! Author(s): Thomas Bretz 08/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
22! Author(s): Stefan Ruegamer <mailto:snruegam@astro.uni-wuerzburg.de>
23!
24! Copyright: MAGIC Software Development, 2000-2006
25!
26!
27\* ======================================================================== */
28
29/////////////////////////////////////////////////////////////////////////////
30//
31// MBadPixelsTreat
32//
33// You can use MBadPixelsTreat::SetUseInterpolation to replaced the
34// bad pixels by the average of its neighbors instead of unmapping
35// them. If you want to include the central pixel use
36// MBadPixelsTreat::SetUseCentralPixel. The bad pixels are taken from
37// an existing MBadPixelsCam.
38//
39// It check if there are enough neighbors to calculate the mean
40// If not, unmap the pixel. The minimum number of good neighbors
41// should be set using SetNumMinNeighbors
42//
43// If you want to interpolate unreliable pixels and unsuitable
44// (broken) pixels use SetHardTreatment().
45//
46//
47// Options:
48// --------
49// SetHardTreatment: Also interpolate unreliable pixels not only unsuitable
50// SetUseInterpolation: Interpolate pixels instead of unmapping them
51// SetUseCentralPixel: also use the pixel itself for interpolation
52// SetProcessPedestalRun: process the pedestals once per run/file
53// SetProcessPedestalEvt: process the pedestal once per event
54// SetProcessTimes: do interpolation of the arrival time
55//
56// If the arrival time treatment is switched on and "MPedPhotFromExtractor"
57// and "MPedPhotFromExtractorRndm" are found the pixel is filled with
58// a random gaus calculated from these two MPedPhotCams in the case
59// the pixels is detected as background.
60//
61//
62// Input Containers:
63// MSignalCam
64// MPedPhotCam
65// MBadPixelsCam
66// [MGeomCam]
67//
68// Output Containers:
69// MSignalCam
70//
71/////////////////////////////////////////////////////////////////////////////
72#include "MBadPixelsTreat.h"
73
74#include <fstream>
75
76#include <TEnv.h>
77#include <TRandom.h>
78#include <TObjString.h>
79
80//#include "MArrayD.h" // Used instead of TArrayD because operator[] has no range check
81
82#include "MLog.h"
83#include "MLogManip.h"
84
85#include "MParList.h"
86
87#include "MGeomPix.h"
88#include "MGeomCam.h"
89
90#include "MSignalPix.h"
91#include "MSignalCam.h"
92
93#include "MPedPhotPix.h"
94#include "MPedPhotCam.h"
95
96#include "MBadPixelsPix.h"
97#include "MBadPixelsCam.h"
98
99ClassImp(MBadPixelsTreat);
100
101using namespace std;
102
103static const TString gsDefName = "MBadPixelsTreat";
104static const TString gsDefTitle = "Task to treat bad pixels (interpolation, unmapping)";
105
106// --------------------------------------------------------------------------
107//
108// Default constructor.
109//
110MBadPixelsTreat::MBadPixelsTreat(const char *name, const char *title)
111 : fGeomCam(0), fEvt(0), fBadPixels(0), fPedPhot1(0), fPedPhot2(0),
112 fFlags(0), fNumMinNeighbors(3), fMaxArrivalTimeDiff(0.9)
113{
114 fName = name ? name : gsDefName.Data();
115 fTitle = title ? title : gsDefTitle.Data();
116
117 SetUseInterpolation();
118 SetProcessPedestal();
119 SetProcessTimes();
120}
121
122// --------------------------------------------------------------------------
123//
124// Returns the status of a pixel. If kHardTreatment is set a pixel must
125// be unsuitable or uriliable to be treated. If not it is treated only if
126// it is unsuitable
127// (IsBad() checks for any flag)
128//
129Bool_t MBadPixelsTreat::IsPixelBad(Int_t idx) const
130{
131 return TESTBIT(fFlags, kHardTreatment) ? (*fBadPixels)[idx].IsBad():(*fBadPixels)[idx].IsUnsuitable();
132}
133
134void MBadPixelsTreat::AddNamePedPhotCam(const char *name)
135{
136 fNamePedPhotCams.Add(new TObjString(name));
137}
138
139// --------------------------------------------------------------------------
140//
141// - Try to find or create MBlindPixels in parameter list.
142// - get the MSignalCam from the parlist (abort if missing)
143// - if no pixels are given by the user try to determin the starfield
144// from the monte carlo run header.
145//
146Int_t MBadPixelsTreat::PreProcess (MParList *pList)
147{
148 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
149 if (!fBadPixels)
150 {
151 *fLog << err << AddSerialNumber("MBadPixelsCam") << " not found... aborting." << endl;
152 return kFALSE;
153 }
154
155 fEvt = (MSignalCam*)pList->FindObject(AddSerialNumber("MSignalCam"));
156 if (!fEvt)
157 {
158 *fLog << err << AddSerialNumber("MSignalCam") << " not found... aborting." << endl;
159 return kFALSE;
160 }
161
162 fGeomCam = 0;
163 if (!IsUseInterpolation())
164 return kTRUE;
165
166 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
167 if (!fGeomCam)
168 {
169 *fLog << err << AddSerialNumber("MGeomCam") << " not found - can't use interpolation... abort." << endl;
170 *fLog << " Use MBadPixelsTreat::SetUseInterpolation(kFALSE) to switch interpolation" << endl;
171 *fLog << " off and use unmapping instead." << endl;
172 return kFALSE;
173 }
174
175 const Bool_t proc = IsProcessPedestalEvt() || IsProcessPedestalRun();
176
177 if (fNamePedPhotCams.GetSize()>0 && !proc)
178 {
179 *fLog << err << "Pedestal list contains entries, but pedestal treatment is switched off... abort." << endl;
180 return kFALSE;
181 }
182
183 if (proc)
184 {
185 if (fNamePedPhotCams.GetSize()==0)
186 {
187 *fLog << inf << "No container names specified... using default: MPedPhotCam." << endl;
188 AddNamePedPhotCam();
189 }
190
191 fPedPhotCams.Clear();
192
193 TIter Next(&fNamePedPhotCams);
194 TObject *o=0;
195 while ((o=Next()))
196 {
197 TObject *p = pList->FindObject(AddSerialNumber(o->GetName()), "MPedPhotCam");
198 if (!p)
199 {
200 *fLog << err << AddSerialNumber(o->GetName()) << " [MPedPhotCam] not found... aborting." << endl;
201 return kFALSE;
202 }
203
204 fPedPhotCams.Add(p);
205 }
206 }
207
208 if (IsProcessTimes())
209 {
210 fPedPhot1 = (MPedPhotCam*)pList->FindObject("MPedPhotFromExtractor", "MPedPhotCam");
211 fPedPhot2 = (MPedPhotCam*)pList->FindObject("MPedPhotFromExtractorRndm", "MPedPhotCam");
212
213 *fLog << inf << "Additional no-signal-interpolation switched ";
214 *fLog << (fPedPhot1 && fPedPhot2 ? "on" : "off");
215 *fLog << "." << endl;
216
217 if (fPedPhot1 && fPedPhot2)
218 *fLog << "Maximum arrival time difference: " << fMaxArrivalTimeDiff << endl;
219
220 }
221
222 *fLog << inf << "Minimum number of neighbor pixels: " << (int)fNumMinNeighbors << endl;
223
224 if (IsProcessPedestalEvt())
225 *fLog << "Processing Pedestals once per event..." << endl;
226 if (IsProcessPedestalRun())
227 *fLog << "Processing Pedestals once per run..." << endl;
228 if (IsProcessTimes())
229 *fLog << "Processing Arrival Times once per event..." << endl;
230
231 return kTRUE;
232}
233
234// --------------------------------------------------------------------------
235//
236// Replaces each pixel (signal, signal error, pedestal, pedestal rms)
237// by the average of its surrounding pixels.
238// If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
239// included.
240//
241void MBadPixelsTreat::InterpolateSignal() const
242{
243 const UShort_t entries = fGeomCam->GetNumPixels();
244
245 //
246 // Loop over all pixels
247 //
248 for (UShort_t i=0; i<entries; i++)
249 {
250 //
251 // Check whether pixel with idx i is blind
252 //
253 if (!IsPixelBad(i))
254 continue;
255
256 //
257 // Get the corresponding geometry and pedestal
258 //
259 MSignalPix &pix = (*fEvt)[i];
260 const MGeomPix &gpix = (*fGeomCam)[i];
261
262 // Do Not-Use-Central-Pixel
263 const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
264
265 Int_t num = nucp ? 0 : 1;
266
267 Double_t nphot = nucp ? 0 : pix.GetNumPhotons();
268 Double_t perr = nucp ? 0 : Pow2(pix.GetErrorPhot());
269
270 //
271 // The values are rescaled to the small pixels area for the right comparison
272 //
273 const Double_t ratio = fGeomCam->GetPixRatio(i);
274
275 nphot *= ratio;
276 perr *= ratio;
277
278 //
279 // Loop over all its neighbors
280 //
281 const Int_t n = gpix.GetNumNeighbors();
282 for (int j=0; j<n; j++)
283 {
284 const UShort_t nidx = gpix.GetNeighbor(j);
285
286 //
287 // Do not use blind neighbors
288 //
289 if (IsPixelBad(nidx))
290 continue;
291
292 //
293 // Get the geometry for the neighbor
294 //
295 const Double_t nratio = fGeomCam->GetPixRatio(nidx);
296
297 //
298 //The error is calculated as the quadratic sum of the errors
299 //
300 const MSignalPix &evtpix = (*fEvt)[nidx];
301
302 nphot += nratio*evtpix.GetNumPhotons();
303 perr += nratio*Pow2(evtpix.GetErrorPhot());
304
305 num++;
306 }
307
308 // Check if there are enough neighbors to calculate the mean
309 // If not, unmap the pixel. The maximum number of blind neighbors
310 // should be 2
311 if (num<fNumMinNeighbors)
312 {
313 pix.SetPixelUnmapped();
314 continue;
315 }
316
317 //
318 // Now the mean is calculated and the values rescaled back
319 // to the pixel area
320 //
321 nphot /= num*ratio;
322 perr = TMath::Sqrt(perr/(num*ratio));
323
324 pix.Set(nphot, perr);
325 }
326}
327
328// --------------------------------------------------------------------------
329//
330void MBadPixelsTreat::InterpolatePedestals(MPedPhotCam &pedphot) const
331{
332 const Int_t entries = pedphot.GetSize();
333
334 //
335 // Loop over all pixels
336 //
337 for (UShort_t i=0; i<entries; i++)
338 {
339 //
340 // Check whether pixel with idx i is blind
341 //
342 if (!IsPixelBad(i))
343 continue;
344
345 //
346 // Get the corresponding geometry and pedestal
347 //
348 const MGeomPix &gpix = (*fGeomCam)[i];
349 const MPedPhotPix &ppix = pedphot[i];
350
351 // Do Not-Use-Central-Pixel
352 const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
353
354 Int_t num = nucp ? 0 : 1;
355
356 Double_t ped = nucp ? 0 : ppix.GetMean();
357 Double_t rms = nucp ? 0 : Pow2(ppix.GetRms());
358
359 //
360 // The values are rescaled to the small pixels area for the right comparison
361 //
362 const Double_t ratio = fGeomCam->GetPixRatio(i);
363
364 ped *= ratio;
365 rms *= ratio;
366
367 //
368 // Loop over all its neighbors
369 //
370 const Int_t n = gpix.GetNumNeighbors();
371 for (int j=0; j<n; j++)
372 {
373 const UShort_t nidx = gpix.GetNeighbor(j);
374
375 //
376 // Do not use blind neighbors
377 //
378 if (IsPixelBad(nidx))
379 continue;
380
381 //
382 // Get the geometry for the neighbor
383 //
384 const Double_t nratio = fGeomCam->GetPixRatio(nidx);
385 const MPedPhotPix &nppix = pedphot[nidx];
386
387 //
388 //The error is calculated as the quadratic sum of the errors
389 //
390 ped += nratio*nppix.GetMean();
391 rms += nratio*Pow2(nppix.GetRms());
392
393 num++;
394 }
395
396 // Check if there are enough neighbors to calculate the mean
397 // If not, unmap the pixel. The minimum number of good neighbors
398 // should be fNumMinNeighbors
399 if (num<fNumMinNeighbors)
400 {
401 (*fEvt)[i].SetPixelUnmapped();
402 continue;
403 }
404
405 //
406 // Now the mean is calculated and the values rescaled back
407 // to the pixel area
408 //
409 ped /= num*ratio;
410 rms = TMath::Sqrt(rms/(num*ratio));
411
412 pedphot[i].Set(ped, rms);
413 }
414 pedphot.SetReadyToSave();
415}
416
417// --------------------------------------------------------------------------
418//
419// loop over all MPedPhotCam and interpolate them
420//
421void MBadPixelsTreat::InterpolatePedestals() const
422{
423 TIter Next(&fPedPhotCams);
424 MPedPhotCam *cam=0;
425 while ((cam=(MPedPhotCam*)Next()))
426 {
427 InterpolatePedestals(*cam);
428 cam->ReCalc(*fGeomCam, fBadPixels);
429 }
430}
431
432// --------------------------------------------------------------------------
433//
434void MBadPixelsTreat::InterpolateTimes() const
435{
436 const Int_t n = fEvt->GetNumPixels();
437 for (int i=0; i<n; i++)
438 {
439 // Check whether pixel with idx i is bad
440 if (!IsPixelBad(i))
441 continue;
442
443 // Geometry of bad pixel
444 const MGeomPix &gpix = (*fGeomCam)[i];
445
446 // Number of neighbor pixels
447 const Int_t n2 = gpix.GetNumNeighbors();
448
449 // Copy the arrival time of all neighboring bad pixels
450 // to a new array for simplicity
451 Double_t time[6];
452 Int_t cnt = 0;
453 for (Int_t j=0; j<n2; j++)
454 {
455 const Int_t idx = gpix.GetNeighbor(j);
456 if (!IsPixelBad(idx))
457 time[cnt++] = (*fEvt)[idx].GetArrivalTime();
458 }
459
460 // if there are too few neighbours, don't interpolate the pixel
461 //if ((cnt < 3 && n2 > 3) || (cnt < 2 && n2 == 3))
462 if (cnt<fNumMinNeighbors)
463 {
464 (*fEvt)[i].SetPixelUnmapped();
465 continue;
466 }
467
468 Double_t min = FLT_MAX; // Find minimum arrival time
469 Double_t max = -FLT_MAX; // Find maximum arrival time
470
471 Double_t sum2 = 0; // Sum of arrival times of the pixels
472 Int_t cnt2 = 0; // Number of pixels summed in sum2
473
474 for (Int_t j=0; j<cnt; j++)
475 {
476 const Double_t tm1 = time[j]; // time of one neighbor pixel
477 const Double_t tm2 = time[(j+1)%cnt]; // time of its neighbor pixel
478
479 // Calculate mean arrival time of pixel probably inside the shower
480 if (TMath::Abs(tm1 - tm2)<fMaxArrivalTimeDiff)
481 {
482 sum2 += tm1+tm2;
483 cnt2++;
484 }
485
486 // Find minimum arrival time
487 if (tm1<min)
488 min = tm1;
489
490 // Find maximum arrival time
491 if (tm1>max)
492 max = tm1;
493 }
494
495 // If less than two nbeighbors belong to a shower the pixel doesn't
496 // belong to the shower, too. Set the arrival time to a uniform
497 // random value, otherwise use the mean value of the pixels belonging
498 // to the shower.
499 if (cnt2<=2)
500 {
501 sum2 = gRandom->Uniform(max-min)+min; // FIXME? Set Seed value?
502
503 // Proceed with a treatment of the signal of empty pixels
504 // better than the interpolation. (FIXME: Maybe a function
505 // different from a gaussian could be a better choice...)
506 if (fPedPhot1 && fPedPhot2)
507 {
508 const Int_t aidx = gpix.GetAidx();
509 // This is to which bias level the signal fluctuates
510 const Double_t mean = fPedPhot1->GetArea(aidx).GetMean();
511 // This is how the signal fluctuates
512 const Double_t rms = fPedPhot2->GetArea(aidx).GetRms();
513 const Double_t phe = gRandom->Gaus(mean, rms);
514
515 (*fEvt)[i].SetNumPhotons(phe);
516 }
517 }
518 else
519 sum2 /= cnt2*2;
520
521 (*fEvt)[i].SetArrivalTime(sum2);
522 }
523}
524
525// --------------------------------------------------------------------------
526//
527// Removes all blind pixels from the analysis by setting their state
528// to unused.
529//
530void MBadPixelsTreat::Unmap() const
531{
532 const UShort_t entries = fEvt->GetNumPixels();
533
534 //
535 // remove the pixels in fPixelsIdx if they are set to be used,
536 // (set them to 'unused' state)
537 //
538 for (UShort_t i=0; i<entries; i++)
539 {
540 if (IsPixelBad(i))
541 (*fEvt)[i].SetPixelUnmapped();
542 }
543}
544
545// --------------------------------------------------------------------------
546//
547// Interpolate Pedestals if kProcessPedestal not set
548//
549Bool_t MBadPixelsTreat::ReInit(MParList *pList)
550{
551 if (IsUseInterpolation() && IsProcessPedestalRun())
552 InterpolatePedestals();
553 return kTRUE;
554}
555
556// --------------------------------------------------------------------------
557//
558// Treat the blind pixels
559//
560Int_t MBadPixelsTreat::Process()
561{
562 if (IsUseInterpolation())
563 {
564 InterpolateSignal();
565 if (IsProcessPedestalEvt())
566 InterpolatePedestals();
567 if (IsProcessTimes())
568 InterpolateTimes();
569 }
570 else
571 Unmap();
572
573 return kTRUE;
574}
575
576void MBadPixelsTreat::StreamPrimitive(ostream &out) const
577{
578 out << " MBadPixelsTreat " << GetUniqueName();
579 if (fName!=gsDefName || fTitle!=gsDefTitle)
580 {
581 out << "(\"" << fName << "\"";
582 if (fTitle!=gsDefTitle)
583 out << ", \"" << fTitle << "\"";
584 out <<")";
585 }
586 out << ";" << endl;
587
588 if (IsUseInterpolation())
589 out << " " << GetUniqueName() << ".SetUseInterpolation();" << endl;
590 if (IsUseCentralPixel())
591 out << " " << GetUniqueName() << ".SetUseCentralPixel();" << endl;
592 if (IsProcessPedestalRun())
593 out << " " << GetUniqueName() << ".SetProcessPedestalRun();" << endl;
594 if (IsProcessPedestalEvt())
595 out << " " << GetUniqueName() << ".SetProcessPedestalEvt();" << endl;
596 if (IsProcessTimes())
597 out << " " << GetUniqueName() << ".SetProcessTimes();" << endl;
598 if (IsHardTreatment())
599 out << " " << GetUniqueName() << ".SetHardTreatment();" << endl;
600 if (fNumMinNeighbors!=3)
601 out << " " << GetUniqueName() << ".SetNumMinNeighbors(" << (int)fNumMinNeighbors << ");" << endl;
602}
603
604// --------------------------------------------------------------------------
605//
606// Read the setup from a TEnv, eg:
607// MBadPixelsTreat.UseInterpolation: no
608// MBadPixelsTreat.UseCentralPixel: no
609// MBadPixelsTreat.HardTreatment: no
610// MBadPixelsTreat.ProcessPedestalRun: no
611// MBadPixelsTreat.ProcessPedestalEvt: no
612// MBadPixelsTreat.NumMinNeighbors: 3
613// MBadPixelsTreat.MaxArrivalTimeDiff: 0.9
614//
615Int_t MBadPixelsTreat::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
616{
617 Bool_t rc = kFALSE;
618 if (IsEnvDefined(env, prefix, "UseInterpolation", print))
619 {
620 rc = kTRUE;
621 SetUseInterpolation(GetEnvValue(env, prefix, "UseInterpolation", IsUseInterpolation()));
622 }
623 if (IsEnvDefined(env, prefix, "UseCentralPixel", print))
624 {
625 rc = kTRUE;
626 SetUseCentralPixel(GetEnvValue(env, prefix, "UseCentralPixel", IsUseCentralPixel()));
627 }
628 if (IsEnvDefined(env, prefix, "HardTreatment", print))
629 {
630 rc = kTRUE;
631 SetHardTreatment(GetEnvValue(env, prefix, "HardTreatment", IsHardTreatment()));
632 }
633 if (IsEnvDefined(env, prefix, "ProcessPedestalRun", print))
634 {
635 rc = kTRUE;
636 SetProcessPedestalRun(GetEnvValue(env, prefix, "ProcessPedestalRun", IsProcessPedestalRun()));
637 }
638 if (IsEnvDefined(env, prefix, "ProcessPedestalEvt", print))
639 {
640 rc = kTRUE;
641 SetProcessPedestalEvt(GetEnvValue(env, prefix, "ProcessPedestalEvt", IsProcessPedestalEvt()));
642 }
643 if (IsEnvDefined(env, prefix, "ProcessTimes", print))
644 {
645 rc = kTRUE;
646 SetProcessTimes(GetEnvValue(env, prefix, "ProcessTimes", IsProcessTimes()));
647 }
648 if (IsEnvDefined(env, prefix, "NumMinNeighbors", print))
649 {
650 rc = kTRUE;
651 SetNumMinNeighbors(GetEnvValue(env, prefix, "NumMinNeighbors", fNumMinNeighbors));
652 }
653 if (IsEnvDefined(env, prefix, "MaxArrivalTimeDiff", print))
654 {
655 rc = kTRUE;
656 SetMaxArrivalTimeDiff(GetEnvValue(env, prefix, "MaxArrivalTimeDiff", fMaxArrivalTimeDiff));
657 }
658 return rc;
659}
Note: See TracBrowser for help on using the repository browser.