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

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