source: trunk/MagicSoft/Mars/manalysis/MBlindPixelsCalc2.cc@ 4094

Last change on this file since 4094 was 3415, checked in by tonello, 21 years ago
*** empty log message ***
File size: 18.1 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): Nadia Tonello 02/2004 <mailto:tonello@mppmu.mpg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2003
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28//
29// MBlindPixelsCalc2
30//
31// NT:: WARNING: This class will change soon to fit into mbadpixels!!
32// It does all MBlindPixelCalc did, plus check the pedestal Rms
33// to define bad pixels.
34//
35// This is the specific image cleaning for a list of pixels. This task
36// sets the pixels listed in fPixelsIdx to unused so they should not be
37// used for analysis (eg calculation of hillas parameters).
38//
39// NT::You can use MBlindPixelsCalc2::SetCheckPedestalRms to see which
40// pixels to blind. If the PedestalRms in one pixel is more than
41// 3 times the mean or less than 1/3 the mean, they are set as blind.
42//
43// You can use MBlindPixelsCalc2::SetUseInterpolation to replace the
44// blind pixels by the average of its neighbors instead of unmapping
45// them. If you want to include the central pixel use
46// MBlindPixelsCalc2::SetUseCentralPixel.
47// NT:: The interpolation will be done only if there are at least
48// 3 non-blind pixels.
49//
50// Example:
51// MBadPixelCalc2 blind;
52// blind.SetCheckPedestalRms();
53// blind.SetUseInterpolation();
54//
55// If you do not use SetUseInterpolation(), the pixels selected with
56// SetCheckPedestalRms will be Unmapped and not used for the image
57// analysis.
58//
59// Input Containers:
60// MCerPhotEvt
61// [MBlindPixels]
62//
63// Output Containers:
64// MCerPhotEvt
65// MBlindPixels
66//
67/////////////////////////////////////////////////////////////////////////////
68#include "MBlindPixelsCalc2.h"
69
70#include <fstream>
71
72#include "MLog.h"
73#include "MLogManip.h"
74
75#include "MParList.h"
76
77#include "MGeomPix.h"
78#include "MGeomCam.h"
79#include "MCerPhotPix.h"
80#include "MCerPhotEvt.h"
81#include "MPedPhotPix.h"
82#include "MPedPhotCam.h"
83#include "MPedestalPix.h"
84#include "MPedestalCam.h"
85#include "MBlindPixels.h"
86
87#include "MMcRunHeader.hxx"
88
89ClassImp(MBlindPixelsCalc2);
90
91using namespace std;
92
93static const TString gsDefName = "MBlindPixelsCalc2";
94static const TString gsDefTitle = "Find hot spots (star, broken pixels, etc)";
95
96// --------------------------------------------------------------------------
97//
98// Default constructor.
99//
100MBlindPixelsCalc2::MBlindPixelsCalc2(const char *name, const char *title)
101 : fFlags(0), fErrors(3)
102{
103 fName = name ? name : gsDefName.Data();
104 fTitle = title ? title : gsDefTitle.Data();
105}
106
107// --------------------------------------------------------------------------
108//
109// - Try to find or create MBlindPixels in parameter list.
110// - get the MCerPhotEvt from the parlist (abort if missing)
111//
112Int_t MBlindPixelsCalc2::PreProcess (MParList *pList)
113{
114 if (TESTBIT(fFlags, kUseBlindPixels))
115 fPixels = (MBlindPixels*)pList->FindObject(AddSerialNumber("MBlindPixels"));
116 else
117 fPixels = (MBlindPixels*)pList->FindCreateObj(AddSerialNumber("MBlindPixels"));
118
119 if (!fPixels)
120 {*fLog << err << dbginf << "Try to remove -SetUseBlindPixels- from the task list to go on" << endl;
121 return kFALSE;}
122
123 fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
124 if (!fEvt)
125 {
126 *fLog << err << dbginf << "MCerPhotEvt not found... aborting." << endl;
127 return kFALSE;
128 }
129
130 fPed = (MPedPhotCam*)pList->FindObject(AddSerialNumber("MPedPhotCam"));
131 if (!fPed)
132 {
133 *fLog << err << dbginf << "MPedPhotCam not found... aborting." << endl;
134 return kFALSE;
135 }
136
137 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
138 if (!fGeomCam)
139 *fLog << warn << dbginf << "No camera geometry available... can't use interpolation." << endl;
140
141 const UShort_t size = fPixelsIdx.GetSize();
142 if (size == 0)
143 {
144 if (!pList->FindObject("MMcRunHeader") && !fPixels)
145 {
146 *fLog << warn << "Warning - Neither blind pixels are given nor a MMcRunHeader was found... removing MBlindPixelsCalc2 from list." << endl;
147 return kSKIP;
148 }
149 return kTRUE;
150 }
151
152 // Set as blind pixels the global blind pixels, which are given
153 // through the macros
154
155 UShort_t numids = fPixelsIdx.GetSize();
156
157 for(Int_t i = 0; i<numids; i++)
158 fPixels->SetPixelBlind(fPixelsIdx[i]);
159
160
161 memset(fErrors.GetArray(), 0, sizeof(Char_t)*fErrors.GetSize());
162
163
164 fErrors[0]= 0;
165 fErrors[1]= 0;
166 fErrors[2]= 0;
167
168 return kTRUE;
169}
170
171// --------------------------------------------------------------------------
172//
173// Replaces each pixel (signal, signal error, pedestal, pedestal rms)
174// by the average of its surrounding pixels.
175// If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
176// included.
177//
178void MBlindPixelsCalc2::Interpolate() const
179{
180 const UShort_t entries = fGeomCam->GetNumPixels();
181
182 //
183 // Create arrays
184 //
185 Double_t *nphot = new Double_t[entries];
186 Double_t *perr = new Double_t[entries];
187
188 //
189 // Loop over all pixels
190 //
191 for (UShort_t i=0; i<entries; i++)
192 {
193 //
194 // Check whether pixel with idx i is blind
195 //
196 if (!fPixels->IsBlind(i))
197 continue;
198
199 //
200 // Get a pointer to this pixel. If it is not yet existing
201 // create a new entry for this pixel in MCerPhotEvt
202 //
203
204 MCerPhotPix *pix = fEvt->GetPixById(i);
205 if (!pix)
206 pix = fEvt->AddPixel(i, 0, 0);
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 (fPixels->IsBlind(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 //
265 // Now the mean is calculated and the values rescaled back to the pixel area
266 //
267 nphot[i] /= (num*ratio);
268 perr[i] = TMath::Sqrt(perr[i]/(num*ratio));
269
270
271 // Check if there are enough neighbors to calculate the mean
272 // If not, unmap the pixel. The maximum number of blind neighbors should be 2
273
274 if (num < 3)
275 { pix->SetPixelUnmapped();
276 // pix->Set(nphot[i], perr[i]);
277 continue;
278 }
279
280 pix->Set(nphot[i], perr[i]);
281 }
282 //
283 // Delete the intermediat arrays
284 //
285 delete nphot;
286 delete perr;
287}
288
289// --------------------------------------------------------------------------
290//
291
292void MBlindPixelsCalc2::InterpolatePedestals() const
293{
294 const Int_t entries = fPed->GetSize();
295 Double_t *ped = new Double_t[entries];
296 Double_t *pedrms = new Double_t[entries];
297
298 //
299 // Loop over all pixels
300 //
301 for (UShort_t i=0; i<entries; i++)
302 {
303 //
304 // Check whether pixel with idx i is blind
305 //
306 if (!fPixels->IsBlind(i))
307 continue;
308
309 //
310 // Get the corresponding geometry and pedestal
311 //
312 const MGeomPix &gpix = (*fGeomCam)[i];
313 const MPedPhotPix &ppix = (*fPed)[i];
314
315 // Do Not-Use-Central-Pixel
316 const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
317
318 Int_t num = nucp ? 0 : 1;
319
320 ped[i] = nucp ? 0 : ppix.GetMean();
321 pedrms[i] = nucp ? 0 : Pow2(ppix.GetRms());
322
323 //
324 // The values are rescaled to the small pixels area for the right comparison
325 //
326 const Double_t ratio = fGeomCam->GetPixRatio(i);
327
328 ped[i] *= ratio;
329 pedrms[i] *= ratio;
330 //
331 // Loop over all its neighbors
332 //
333 const Int_t n = gpix.GetNumNeighbors();
334 for (int j=0; j<n; j++)
335 {
336 const UShort_t nidx = gpix.GetNeighbor(j);
337
338 //
339 // Do not use blind neighbors
340 //
341 if (fPixels->IsBlind(nidx))
342 continue;
343
344 //
345 // Get the geometry for the neighbor
346 //
347 const Double_t nratio = fGeomCam->GetPixRatio(nidx);
348 MPedPhotPix &nppix = (*fPed)[nidx];
349
350 //
351 //The error is calculated as the quadratic sum of the errors
352 //
353 ped[i] += (nratio*nppix.GetMean());
354 pedrms[i] += nratio*Pow2(nppix.GetRms());
355
356 num++;
357 }
358
359 // Check if there are enough neighbors to calculate the mean
360 // If not, unmap the pixel. The minimum number of good neighbors
361 // should be 3
362
363 if (num < 3)
364 {
365 MCerPhotPix *pix = fEvt->GetPixById(i);
366 if (!pix)
367 pix = fEvt->AddPixel(i, 0, 0);
368 pix->SetPixelUnmapped();
369 continue;
370 }
371
372 //
373 // Now the mean is calculated and the values rescaled back to the pixel area
374 //
375 ped[i] /= (num*ratio);
376 pedrms[i] = TMath::Sqrt(pedrms[i]/(num*ratio));
377
378 (*fPed)[i].Set(ped[i], pedrms[i]);
379 }
380
381 //
382 // Delete the arrays
383 //
384 delete ped;
385 delete pedrms;
386
387 return;
388}
389
390// --------------------------------------------------------------------------
391//
392// Removes all blind pixels from the analysis by setting their state
393// to unmapped.
394//
395void MBlindPixelsCalc2::Unmap() const
396{
397 const UShort_t entries = fEvt->GetNumPixels();
398
399 //
400 // remove the pixels in fPixelsIdx if they are set to be used,
401 // (set them to 'unused' state)
402 //
403 for (UShort_t i=0; i<entries; i++)
404 {
405 MCerPhotPix &pix = (*fEvt)[i];
406
407 if (fPixels->IsBlind(pix.GetPixId()))
408 pix.SetPixelUnmapped();
409 }
410}
411
412// --------------------------------------------------------------------------
413// Check the pedestal Rms of the pixels: compute with 2 iterations the mean
414// for inner and outer pixels. Set as blind the pixels with too small or
415// too high pedestal Rms with respect to the mean.
416//
417Bool_t MBlindPixelsCalc2::CheckPedestalRms()
418
419{
420 const Int_t entries = fPed->GetSize();
421 Int_t pixtoblind = 0;
422 Int_t pixin = 0;
423 Int_t pixout = 0;
424
425 Float_t meanPedRms =0.;
426 Float_t meanPedRmsout = 0.;
427 Float_t meanPedRms2 = 0.;
428 Float_t meanPedRmsout2 = 0.;
429
430 for (Int_t i=0; i<entries; i++)
431 {
432 //
433 // get pixel as entry from list
434 //
435 const MPedPhotPix &ppix = (*fPed)[i];
436 const Double_t nratio = fGeomCam->GetPixRatioSqrt(i);
437 const Double_t pixPedRms = ppix.GetRms();
438
439 //Calculate the corrected means:
440
441 if (nratio == 1 )
442 {
443 if( pixPedRms >0. && pixPedRms <200. )
444 {meanPedRms += pixPedRms;
445 pixin++;}
446 }
447 else
448 {
449 if( pixPedRms >0. && pixPedRms <500. )
450 {meanPedRmsout += pixPedRms;
451 pixout++;}
452 }
453 }
454
455 //if no pixel has a minimum signal, return
456 if( pixin==0 || pixout==0 )
457 {fErrors[1]++; //no valid Pedestals Rms
458 return kFALSE;
459 }
460
461 meanPedRms = meanPedRms / pixin ;
462 meanPedRmsout = meanPedRmsout / pixout ;
463
464 if ( meanPedRms == 0. || meanPedRmsout ==0.)
465 {fErrors[1]++; //no valid Pedestals Rms
466 return kFALSE;
467 }
468
469 pixin = 0;
470 pixout = 0;
471
472 for (Int_t i=0; i<entries; i++)
473 {
474 //
475 // get pixel as entry from list
476 //
477 const MPedPhotPix &ppix = (*fPed)[i];
478 const Double_t nratio = fGeomCam->GetPixRatioSqrt(i);
479 const Double_t pixPedRms = ppix.GetRms();
480
481 //Calculate the corrected means:
482
483 if (nratio == 1)
484 {if (pixPedRms > 0.5 * meanPedRms && pixPedRms < 1.5 * meanPedRms)
485 {meanPedRms2 += pixPedRms;
486 pixin++;}
487 }
488 else
489 {if (pixPedRms > 0.5 * meanPedRmsout && pixPedRms < 1.5 * meanPedRmsout)
490 {meanPedRmsout2 += pixPedRms;
491 pixout++;}
492 }
493 }
494
495 //if no pixel has a minimum signal, return
496 if( pixin==0 || pixout==0 )
497 {fErrors[1]++; //no valid Pedestals Rms
498 return kFALSE;}
499
500 meanPedRms = meanPedRms2 / pixin ;
501 meanPedRmsout = meanPedRmsout2 / pixout ;
502
503 //Blind the Bad Pixels
504 for (Int_t i=0; i<entries; i++)
505 {
506 const MPedPhotPix &ppix = (*fPed)[i];
507 const Float_t nratio = fGeomCam->GetPixRatioSqrt(i);
508 const Double_t pixPedRms = ppix.GetRms();
509
510 if (nratio == 1)
511 { if (pixPedRms > 3.* meanPedRms || pixPedRms <= meanPedRms/3.)
512 { fPixels->SetPixelBlind(i);
513 pixtoblind++;}
514 }
515 else
516 { if (pixPedRms > 3.* meanPedRmsout || pixPedRms <= meanPedRmsout/3.)
517 { fPixels->SetPixelBlind(i);
518 pixtoblind++; }
519 }
520 }
521
522 // Check if the number of pixels to blind is > 60% of total number of pixels
523 //
524 if (pixtoblind>0.6*entries)
525 { fErrors[2]++;
526 return kFALSE; }
527
528
529 return kTRUE;
530
531}
532
533// --------------------------------------------------------------------------
534//
535// Treat the blind pixels
536//
537Int_t MBlindPixelsCalc2::Process()
538{
539 if (TESTBIT(fFlags, kCheckPedestalRms))
540 {
541 // if the number of blind pixels is too high, do not interpolate
542 if (CheckPedestalRms()== kFALSE)
543 return kTRUE;
544 else
545 {if (TESTBIT(fFlags, kUseInterpolation))
546 InterpolatePedestals();
547 }
548 }
549
550 if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam)
551 { Interpolate();
552 }
553 else
554 { Unmap();
555 }
556 fErrors[0]++;
557 return kTRUE;
558}
559
560
561// --------------------------------------------------------------------------
562//
563// - Check whether pixels to disable are available. If pixels are
564// given by the user nothing more is done.
565//
566Bool_t MBlindPixelsCalc2::ReInit(MParList *pList)
567{
568 if (TESTBIT(fFlags, kUseBlindPixels))
569 return kTRUE;
570
571 //
572 // If pixels are given by the user, we are already done
573 //
574 if (fPixelsIdx.GetSize() > 0)
575 return kTRUE;
576
577 //
578 // Delete the old array holding the blind pixels for the last file
579 //
580 fPixels->Clear();
581
582 if (!fGeomCam->InheritsFrom("MGeomCamMagic"))
583 {
584 *fLog << warn << "MBlindPixelsCalc2::ReInit: Warning - Starfield only implemented for Magic standard Camera... no action." << endl;
585 return kTRUE;
586 }
587 return kTRUE;
588}
589
590void MBlindPixelsCalc2::StreamPrimitive(ofstream &out) const
591{
592 out << " MBlindPixelsCalc2 " << GetUniqueName();
593 if (fName!=gsDefName || fTitle!=gsDefTitle)
594 {
595 out << "(\"" << fName << "\"";
596 if (fTitle!=gsDefTitle)
597 out << ", \"" << fTitle << "\"";
598 out <<")";
599 }
600 out << ";" << endl;
601
602 if (TESTBIT(fFlags, kUseInterpolation))
603 out << " " << GetUniqueName() << ".SetUseInterpolation();" << endl;
604 if (TESTBIT(fFlags, kUseCentralPixel))
605 out << " " << GetUniqueName() << ".SetUseCentralPixel();" << endl;
606 if (TESTBIT(fFlags, kCheckPedestalRms))
607 out << " " << GetUniqueName() << ".SetCheckPedestalRms();" << endl;
608
609 if (fPixelsIdx.GetSize()==0)
610 return;
611
612 out << " {" << endl;
613 out << " TArrayS dummy;" << endl;
614 for (int i=0; i<fPixelsIdx.GetSize(); i++)
615 out << " dummy[" << i << "]=" << ((TArrayS)fPixelsIdx)[i] << ";" << endl;
616 out << " " << GetUniqueName() << ".SetPixels(dummy);" << endl;
617 out << " }" << endl;
618}
619
620
621// --------------------------------------------------------------------------
622//
623// This is used to print the output in the PostProcess. Later (having
624// millions of events) we may want to improve the output.
625//
626void MBlindPixelsCalc2::PrintSkipped(int i, const char *str) const
627{
628 *fLog << " " << setw(7) << fErrors[i] << " (";
629 *fLog << setw(3) << (int)(100.*fErrors[i]/GetNumExecutions());
630 *fLog << "%) Evts skipped due to: " << str << endl;
631}
632
633// --------------------------------------------------------------------------
634//
635// Prints some statistics about the calculation. The percentage
636// is calculated with respect to the number of executions of this task.
637//
638Int_t MBlindPixelsCalc2::PostProcess()
639{
640 if (GetNumExecutions()==0)
641 return kTRUE;
642
643 *fLog << inf << endl;
644 *fLog << GetDescriptor() << " execution statistics:" << endl;
645 *fLog << dec << setfill(' ');
646
647 // Number of events where the computation of the mean pedestal Rms failed
648 // On these events no blind pixels have been selected
649
650 if (TESTBIT(fFlags, kCheckPedestalRms))
651 PrintSkipped(1, "Big fluctuations in Pedestal Rms, no CHECK");
652
653 // Number of events with too many bad pixels have been found
654 // On these events, no interpolation has been applied
655
656 if (TESTBIT(fFlags, kUseInterpolation))
657 PrintSkipped(2, " >60% of pixels are bad, no INTERPOLATION");
658
659 //Number of events that have been succesfully checked and corrected
660
661 *fLog << " " << (int)fErrors[0] << " (" << (int)(100.*fErrors[0]/GetNumExecutions()) << "%) Evts well treated!" << endl;
662 *fLog << endl;
663
664 return kTRUE;
665}
666
667
668
669
670
671
Note: See TracBrowser for help on using the repository browser.