source: branches/MarsMoreSimulationTruth/mbadpixels/MBadPixelsTreat.cc@ 20115

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