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

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