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

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