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

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