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

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