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

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