source: trunk/MagicSoft/Mars/msignal/MExtractBlindPixel.cc@ 4964

Last change on this file since 4964 was 4899, checked in by gaug, 21 years ago
*** empty log message ***
File size: 18.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): Markus Gaug, 02/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MExtractBlindPixel
28//
29// Extracts the signal from a fixed window in a given range.
30//
31// Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast)
32// to modify the ranges. The "low-gain" ranges are used for the NSB rejection
33// whereas the high-gain ranges for blind pixel signal extraction. "High-gain"
34// ranges can extend to the slices stored as "low-gain" in MRawEvtPixelIter
35//
36// Defaults are:
37//
38// fHiGainFirst = fgHiGainFirst = 10
39// fHiGainLast = fgHiGainLast = 29
40// fLoGainFirst = fgLoGainFirst = 0
41// fLoGainLast = fgLoGainLast = 7
42//
43// The switches:
44// - SetExtractionType ( kAmplitude ) and SetExtractionType ( kIntegral )
45// can be used to choose between amplitude extraction (using a spline) and
46// summed integral.
47// - SetExtractionType ( kFilter )
48// can be used to apply a filter discarding events passing over a threshold
49// defined in fNSBFilterLimit
50//
51//////////////////////////////////////////////////////////////////////////////
52#include "MExtractBlindPixel.h"
53
54#include "MLog.h"
55#include "MLogManip.h"
56
57#include "MParList.h"
58
59#include "MRawEvtData.h"
60#include "MRawRunHeader.h"
61#include "MRawEvtPixelIter.h"
62
63#include "MExtractedSignalBlindPixel.h"
64
65#include "MPedestalCam.h"
66#include "MPedestalPix.h"
67
68ClassImp(MExtractBlindPixel);
69
70using namespace std;
71
72const UInt_t MExtractBlindPixel::fgBlindPixelIdx = 559;
73const Byte_t MExtractBlindPixel::fgHiGainFirst = 10;
74const Byte_t MExtractBlindPixel::fgHiGainLast = 19;
75const Byte_t MExtractBlindPixel::fgLoGainFirst = 0;
76const Byte_t MExtractBlindPixel::fgLoGainLast = 7;
77const Int_t MExtractBlindPixel::fgNSBFilterLimit = 70;
78const Float_t MExtractBlindPixel::fgResolution = 0.003;
79const Float_t MExtractBlindPixel::gkOverflow = 300.;
80// --------------------------------------------------------------------------
81//
82// Default constructor.
83//
84// Initializes:
85// - fBlindPixelIdx to fgBlindPixelIdx
86// - fNSBFilterLimit to fgNSBFilterLimit
87// - fResolution to fgResolution
88// - fExtractionType to 0.
89//
90// Calls:
91// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
92//
93MExtractBlindPixel::MExtractBlindPixel(const char *name, const char *title)
94 : fHiGainSignal(NULL),
95 fHiGainFirstDeriv(NULL),
96 fHiGainSecondDeriv(NULL)
97{
98
99 fName = name ? name : "MExtractBlindPixel";
100 fTitle = title ? title : "Task to extract the signal from the FADC slices";
101
102 AddToBranchList("MRawEvtData.*");
103
104 SetResolution();
105 SetNSBFilterLimit();
106 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
107
108 SetNumBlindPixels();
109
110 Clear();
111
112}
113
114// --------------------------------------------------------------------------
115//
116// Destructor:
117//
118// - Deletes, (if not NULL):
119// fHiGainSignal;
120// fHiGainFirstDeriv;
121// fHiGainSecondDeriv;
122//
123MExtractBlindPixel::~MExtractBlindPixel()
124{
125
126 if (fHiGainSignal)
127 delete [] fHiGainSignal;
128 if (fHiGainFirstDeriv)
129 delete [] fHiGainFirstDeriv;
130 if (fHiGainSecondDeriv)
131 delete [] fHiGainSecondDeriv;
132
133}
134
135// --------------------------------------------------------------------------
136//
137// Clear
138//
139// Initializes:
140// - fBlindPixelIdx to 0
141// - fExtractionType to 0
142//
143// Calls:
144// - SetBlindPixelIdx()
145//
146// Deletes and sets to NULL (if exists):
147// - fHiGainSignal
148// - fHiGainFirstDeriv
149// - fHiGainSecondDeriv
150//
151void MExtractBlindPixel::Clear( const Option_t *o)
152{
153
154 fExtractionType = 0;
155
156 fBlindPixelIdx.Set(0);
157 SetBlindPixelIdx();
158
159 if (fHiGainSignal)
160 {
161 delete [] fHiGainSignal;
162 fHiGainSignal = NULL;
163 }
164 if (fHiGainFirstDeriv)
165 {
166 delete [] fHiGainFirstDeriv;
167 fHiGainFirstDeriv = NULL;
168 }
169 if (fHiGainSecondDeriv)
170 {
171 delete [] fHiGainSecondDeriv;
172 fHiGainSecondDeriv = NULL;
173 }
174
175}
176
177void MExtractBlindPixel::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
178{
179
180 MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
181
182 fNumHiGainSamples = (Float_t)(fHiGainLast-fHiGainFirst+1);
183 if (lolast)
184 fNumLoGainSamples = (Float_t)(fLoGainLast-fLoGainFirst+1);
185 else
186 fNumLoGainSamples = 0.;
187
188 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
189 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
190
191 fHiLoFirst = 0;
192 fHiLoLast = 0;
193}
194
195// --------------------------------------------------------------------------
196//
197// Calls:
198// - MExtractor::PreProcess(pList)
199//
200// The following output containers are also searched and created if
201// they were not found:
202//
203// - MExtractedBlindPixel
204//
205Int_t MExtractBlindPixel::PreProcess(MParList *pList)
206{
207
208 if (!MExtractor::PreProcess(pList))
209 return kFALSE;
210
211 fBlindPixel = (MExtractedSignalBlindPixel*)pList->FindCreateObj(AddSerialNumber("MExtractedSignalBlindPixel"));
212 if (!fBlindPixel)
213 return kFALSE;
214
215
216 return kTRUE;
217}
218
219// -------------------------------------------------------------------------- //
220//
221// The ReInit searches for:
222// - MRawRunHeader::GetNumSamplesHiGain()
223// - MRawRunHeader::GetNumSamplesLoGain()
224//
225// In case that the variables fHiGainLast and fLoGainLast are smaller than
226// the even part of the number of samples obtained from the run header, a
227// warning is given an the range is set back accordingly. A call to:
228// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or
229// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff)
230// is performed in that case. The variable diff means here the difference
231// between the requested range (fHiGainLast) and the available one. Note that
232// the functions SetRange() are mostly overloaded and perform more checks,
233// modifying the ranges again, if necessary.
234//
235Bool_t MExtractBlindPixel::ReInit(MParList *pList)
236{
237
238 if (fHiGainSignal)
239 delete [] fHiGainSignal;
240 if (fHiGainFirstDeriv)
241 delete [] fHiGainFirstDeriv;
242 if (fHiGainSecondDeriv)
243 delete [] fHiGainSecondDeriv;
244
245 for (Int_t i=0;i<fNumBlindPixels;i++)
246 fBlindPixel->SetBlindPixelIdx(fBlindPixelIdx.At(i),i);
247
248 fBlindPixel->SetExtractionType(fExtractionType);
249
250 for (Int_t i=0;i<fBlindPixelIdx.GetSize();i++)
251 {
252
253 MPedestalPix &pedpix = (*fPedestals)[fBlindPixelIdx.At(i)];
254
255 if (&pedpix)
256 {
257 fBlindPixel->SetPed ( pedpix.GetPedestal() * fNumLoGainSamples, i );
258 fBlindPixel->SetPedErr ( pedpix.GetPedestalRms()* fNumLoGainSamples
259 / TMath::Sqrt((Float_t)fPedestals->GetTotalEntries()), i );
260 fBlindPixel->SetPedRms ( pedpix.GetPedestalRms()* TMath::Sqrt((Float_t)fNumLoGainSamples), i );
261 fBlindPixel->SetPedRmsErr( fBlindPixel->GetPedErr()/2., i );
262 }
263 }
264
265 const Int_t higainsamples = fRunHeader->GetNumSamplesHiGain();
266 const Int_t logainsamples = fRunHeader->GetNumSamplesLoGain();
267 Int_t lastavailable = higainsamples-1;
268
269 if (logainsamples)
270 {
271 //
272 // If the signal is searched entirely in the low-gain range, have
273 // to skip the higain completely. This is steered by the variable fHiLoFirst
274 //
275 const Int_t firstdesired = (Int_t)fHiGainFirst;
276
277 if (firstdesired > lastavailable)
278 {
279 const Int_t diff = firstdesired - lastavailable;
280 *fLog << endl;
281 *fLog << warn << "First Hi Gain slice " << (int)fHiGainFirst << " out of range [0,";
282 *fLog << lastavailable << "]... start at slice " << diff << " of the Lo Gain " << endl;
283
284 fHiLoFirst = diff;
285 }
286 }
287
288 const Int_t lastdesired = (Int_t)fHiGainLast;
289
290 if (lastdesired > lastavailable)
291 {
292 Int_t diff = lastdesired - lastavailable;
293 lastavailable += logainsamples ? logainsamples-1 : 0;
294
295 if (lastdesired > lastavailable)
296 {
297 *fLog << endl;
298 *fLog << "Last Hi Gain slice " << (int)fHiGainLast << " out of range [0,";
299 *fLog << lastavailable << "]... reduce upper limit by " << diff << endl;
300 diff = logainsamples;
301 }
302
303 fHiGainLast = higainsamples - 1;
304 fHiLoLast = logainsamples ? diff : 0;
305 }
306
307 const Int_t range = fHiLoFirst ? fHiLoLast - fHiLoFirst + 1 : fHiGainLast - fHiGainFirst + fHiLoLast + 1;
308
309 fHiGainSignal = new Float_t[range];
310 memset(fHiGainSignal,0,range*sizeof(Float_t));
311 fHiGainFirstDeriv = new Float_t[range];
312 memset(fHiGainFirstDeriv,0,range*sizeof(Float_t));
313 fHiGainSecondDeriv = new Float_t[range];
314 memset(fHiGainSecondDeriv,0,range*sizeof(Float_t));
315
316 *fLog << endl;
317 *fLog << inf << ": Extracting "
318 << Form("%s",IsExtractionType(kAmplitude) ? "Amplitude " : " Integral ")
319 << " using " << range << " FADC samples from "
320 << Form("%s%2i",fHiLoFirst ? "Low Gain slice " : " High Gain slice ",
321 fHiLoFirst ? (Int_t)fHiLoFirst : (Int_t)fHiGainFirst)
322 << " to (including) "
323 << Form("%s%2i",fHiLoLast ? "Low Gain slice " : " High Gain slice ",
324 fHiLoLast ? (Int_t)fHiLoLast-1 : (Int_t)fHiGainLast )
325 << endl;
326
327 if (IsExtractionType(kFilter))
328 *fLog << inf << GetDescriptor() << ": Will use Filter using "
329 << (Int_t)(fLoGainLast-fLoGainFirst+1) << " FADC slices "
330 << "from High Gain slice " << (Int_t)fLoGainFirst
331 << " to High Gain slice " << (Int_t)fLoGainLast << endl;
332
333 fBlindPixel->SetUsedFADCSlices(fHiGainFirst, range);
334
335 return kTRUE;
336
337}
338
339// --------------------------------------------------------------------------
340//
341// FindSignalHiGain:
342//
343// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
344// - Sum up contents of *ptr
345// - If *ptr is greater than fSaturationLimit, raise sat by 1
346// - If fHiLoLast is set, loop from logain to (logain+fHiLoLast)
347// - Add contents of *logain to sum
348//
349void MExtractBlindPixel::FindIntegral(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat)
350{
351
352 Int_t summ = 0;
353 Byte_t *p = ptr;
354 Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
355
356 if (fHiLoFirst == 0)
357 {
358
359 while (p<end)
360 {
361 summ += *ptr;
362
363 if (*p++ >= fSaturationLimit)
364 sat++;
365 }
366
367 }
368
369 p = logain + fHiLoFirst;
370 end = logain + fHiLoLast;
371 while (p<end)
372 {
373 summ += *p;
374
375 if (*p++ >= fSaturationLimit)
376 sat++;
377 }
378
379 sum = (Float_t)summ;
380}
381
382// --------------------------------------------------------------------------
383//
384// FindSignalPhe:
385//
386// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
387// - Sum up contents of *ptr
388// - If *ptr is greater than fSaturationLimit, raise sat by 1
389// - If fHiLoLast is set, loop from logain to (logain+fHiLoLast)
390// - Add contents of *logain to sum
391//
392void MExtractBlindPixel::FindAmplitude(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat)
393{
394
395 Int_t range = 0;
396 Int_t count = 0;
397 Float_t abmaxpos = 0.;
398 Byte_t *p = ptr;
399 Byte_t *end;
400 Byte_t max = 0;
401 Byte_t maxpos = 0;
402 Int_t summ = 0;
403
404 if (fHiLoFirst == 0)
405 {
406
407 range = fHiGainLast - fHiGainFirst + 1;
408
409 end = ptr + range;
410 //
411 // Check for saturation in all other slices
412 //
413 while (p++<end)
414 {
415
416 fHiGainSignal[count] = (Float_t)*p;
417 summ += *p;
418
419 if (*p > max)
420 {
421 max = *p;
422 maxpos = count;
423 }
424
425 count++;
426
427 if (*p >= fSaturationLimit)
428 sat++;
429 }
430 }
431
432 if (fHiLoLast != 0)
433 {
434
435 p = logain + fHiLoFirst;
436 end = logain + fHiLoLast;
437
438 while (p<end)
439 {
440
441 fHiGainSignal[count] = (Float_t)*p;
442 summ += *p;
443
444 if (*p > max)
445 {
446 max = *p;
447 maxpos = count;
448 }
449
450 range++;
451 count++;
452
453 if (*p++ >= fSaturationLimit)
454 sat++;
455 }
456 }
457
458 //
459 // allow one saturated slice
460 //
461 if (sat > 1)
462 {
463 sum = gkOverflow;
464 return;
465 }
466
467 //
468 // Don't start if the maxpos is too close to the left limit.
469 //
470 if (maxpos < 2)
471 {
472 sum = (Float_t)max;
473 return;
474 }
475
476 Float_t pp;
477
478 for (Int_t i=1;i<range-1;i++)
479 {
480 pp = fHiGainSecondDeriv[i-1] + 4.;
481 fHiGainSecondDeriv[i] = -1.0/pp;
482 fHiGainFirstDeriv [i] = fHiGainSignal[i+1] - fHiGainSignal[i] - fHiGainSignal[i] + fHiGainSignal[i-1];
483 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
484 p++;
485 }
486
487 fHiGainSecondDeriv[range-1] = 0.;
488 for (Int_t k=range-2;k>=0;k--)
489 fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
490
491 //
492 // Now find the maximum
493 //
494 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
495 Float_t lower = (Float_t)maxpos-1.;
496 Float_t upper = (Float_t)maxpos;
497 Float_t x = lower;
498 Float_t y = 0.;
499 Float_t a = 1.;
500 Float_t b = 0.;
501 Int_t klo = maxpos-1;
502 Int_t khi = maxpos;
503 Float_t klocont = fHiGainSignal[klo];
504 Float_t khicont = fHiGainSignal[khi];
505 sum = (Float_t)khicont;
506 abmaxpos = lower;
507
508 //
509 // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
510 // interval maxpos+1.
511 //
512 while (x<upper-0.3)
513 {
514
515 x += step;
516 a -= step;
517 b += step;
518
519 y = a*klocont
520 + b*khicont
521 + (a*a*a-a)*fHiGainSecondDeriv[klo]
522 + (b*b*b-b)*fHiGainSecondDeriv[khi];
523
524 if (y > sum)
525 {
526 sum = y;
527 abmaxpos = x;
528 }
529 }
530
531 if (abmaxpos > upper-0.1)
532 {
533
534 upper = (Float_t)maxpos+1;
535 lower = (Float_t)maxpos;
536 x = lower;
537 a = 1.;
538 b = 0.;
539 khi = maxpos+1;
540 klo = maxpos;
541 klocont = fHiGainSignal[klo];
542 khicont = fHiGainSignal[khi];
543
544 while (x<upper-0.3)
545 {
546
547 x += step;
548 a -= step;
549 b += step;
550
551 y = a* klocont
552 + b* khicont
553 + (a*a*a-a)*fHiGainSecondDeriv[klo]
554 + (b*b*b-b)*fHiGainSecondDeriv[khi];
555
556 if (y > sum)
557 {
558 sum = y;
559 abmaxpos = x;
560 }
561 }
562 }
563
564 const Float_t up = abmaxpos+step-0.055;
565 const Float_t lo = abmaxpos-step+0.055;
566 const Float_t maxpossave = abmaxpos;
567
568 x = abmaxpos;
569 a = upper - x;
570 b = x - lower;
571
572 step = 0.04; // step size of 83 ps
573
574 while (x<up)
575 {
576
577 x += step;
578 a -= step;
579 b += step;
580
581 y = a* klocont
582 + b* khicont
583 + (a*a*a-a)*fHiGainSecondDeriv[klo]
584 + (b*b*b-b)*fHiGainSecondDeriv[khi];
585
586 if (y > sum)
587 {
588 sum = y;
589 abmaxpos = x;
590 }
591 }
592
593 if (abmaxpos < klo + 0.02)
594 {
595 klo--;
596 khi--;
597 klocont = fHiGainSignal[klo];
598 khicont = fHiGainSignal[khi];
599 upper--;
600 lower--;
601 }
602
603 x = maxpossave;
604 a = upper - x;
605 b = x - lower;
606
607 while (x>lo)
608 {
609
610 x -= step;
611 a += step;
612 b -= step;
613
614 y = a* klocont
615 + b* khicont
616 + (a*a*a-a)*fHiGainSecondDeriv[klo]
617 + (b*b*b-b)*fHiGainSecondDeriv[khi];
618
619 if (y > sum)
620 {
621 sum = y;
622 abmaxpos = x;
623 }
624 }
625
626}
627
628// --------------------------------------------------------------------------
629//
630// FindSignalFilter:
631//
632// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
633// - Sum up contents of *ptr
634// - If *ptr is greater than fSaturationLimit, raise sat by 1
635//
636void MExtractBlindPixel::FindSignalFilter(Byte_t *ptr, Int_t &sum, Byte_t &sat) const
637{
638
639 Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1;
640
641 while (ptr<end)
642 {
643 sum += *ptr;
644
645 if (*ptr++ >= fSaturationLimit)
646 sat++;
647 }
648}
649
650// --------------------------------------------------------------------------
651//
652// Calculate the integral of the FADC time slices and store them as a new
653// pixel in the MExtractedBlindPixel container.
654//
655Int_t MExtractBlindPixel::Process()
656{
657
658 MRawEvtPixelIter pixel(fRawEvt);
659
660 fBlindPixel->Clear();
661
662 for (Int_t id=0;id<fBlindPixelIdx.GetSize();id++)
663 {
664
665 pixel.Jump(fBlindPixelIdx[id]);
666
667 Int_t sum = 0;
668 Byte_t sat = 0;
669
670 if (IsExtractionType(kFilter))
671 {
672
673 FindSignalFilter(pixel.GetHiGainSamples()+fLoGainFirst, sum, sat);
674
675 if (sum > fNSBFilterLimit)
676 {
677 fBlindPixel->SetExtractedSignal(-1.,id);
678 fBlindPixel->SetNumSaturated(sat,id);
679 fBlindPixel->SetReadyToSave();
680 continue;
681 }
682
683 sum = 0;
684 if (pixel.HasLoGain())
685 FindSignalFilter(pixel.GetLoGainSamples(), sum, sat);
686
687 /*
688 if (fModified)
689 {
690 if (sum > fNSBFilterLimit)
691 {
692 fBlindPixel->SetExtractedSignal(-1.,id);
693 fBlindPixel->SetNumSaturated(sat,id);
694 fBlindPixel->SetReadyToSave();
695 continue;
696 }
697 }
698 */
699 }
700
701 Float_t newsum = 0.;
702 sat = 0;
703
704 if (IsExtractionType(kAmplitude))
705 FindAmplitude (pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), newsum, sat);
706 else
707 FindIntegral (pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), newsum, sat);
708
709
710 fBlindPixel->SetExtractedSignal(newsum,id);
711 fBlindPixel->SetNumSaturated(sat,id);
712 }
713
714 fBlindPixel->SetReadyToSave();
715 return kTRUE;
716}
717
718// ------------------------------------------------------------------------------------
719//
720// Returns true if the extraction type. Available are: kAmplitude, kIntegral and kFilter
721// The flags kIntegral and kFilter may be set both.
722//
723Bool_t MExtractBlindPixel::IsExtractionType( const ExtractionType_t typ )
724{
725
726 return TESTBIT( fExtractionType, typ );
727
728}
729
730// --------------------------------------------------------------------------
731//
732// Sets the extraction type. Available are: kAmplitude and kIntegral
733//
734void MExtractBlindPixel::SetExtractionType( const ExtractionType_t typ )
735{
736 SETBIT( fExtractionType, typ );
737}
Note: See TracBrowser for help on using the repository browser.