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

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