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

Last change on this file since 4693 was 4669, checked in by gaug, 20 years ago
*** empty log message ***
File size: 18.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): 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 fNumLoGainSamples = (Float_t)(fLoGainLast-fLoGainFirst+1);
184
185 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
186 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
187
188 fHiLoFirst = 0;
189 fHiLoLast = 0;
190}
191
192// --------------------------------------------------------------------------
193//
194// Calls:
195// - MExtractor::PreProcess(pList)
196//
197// The following output containers are also searched and created if
198// they were not found:
199//
200// - MExtractedBlindPixel
201//
202Int_t MExtractBlindPixel::PreProcess(MParList *pList)
203{
204
205 if (!MExtractor::PreProcess(pList))
206 return kFALSE;
207
208 fBlindPixel = (MExtractedSignalBlindPixel*)pList->FindCreateObj(AddSerialNumber("MExtractedSignalBlindPixel"));
209 if (!fBlindPixel)
210 return kFALSE;
211
212
213 return kTRUE;
214}
215
216// -------------------------------------------------------------------------- //
217//
218// The ReInit searches for:
219// - MRawRunHeader::GetNumSamplesHiGain()
220// - MRawRunHeader::GetNumSamplesLoGain()
221//
222// In case that the variables fHiGainLast and fLoGainLast are smaller than
223// the even part of the number of samples obtained from the run header, a
224// warning is given an the range is set back accordingly. A call to:
225// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or
226// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff)
227// is performed in that case. The variable diff means here the difference
228// between the requested range (fHiGainLast) and the available one. Note that
229// the functions SetRange() are mostly overloaded and perform more checks,
230// modifying the ranges again, if necessary.
231//
232Bool_t MExtractBlindPixel::ReInit(MParList *pList)
233{
234
235 if (fHiGainSignal)
236 delete [] fHiGainSignal;
237 if (fHiGainFirstDeriv)
238 delete [] fHiGainFirstDeriv;
239 if (fHiGainSecondDeriv)
240 delete [] fHiGainSecondDeriv;
241
242 for (Int_t i=0;i<fNumBlindPixels;i++)
243 fBlindPixel->SetBlindPixelIdx(fBlindPixelIdx.At(i),i);
244
245 fBlindPixel->SetExtractionType(fExtractionType);
246
247 for (Int_t i=0;i<fBlindPixelIdx.GetSize();i++)
248 {
249
250 MPedestalPix &pedpix = (*fPedestals)[fBlindPixelIdx.At(i)];
251
252 if (&pedpix)
253 {
254 fBlindPixel->SetPed ( pedpix.GetPedestal() * fNumLoGainSamples, i );
255 fBlindPixel->SetPedErr ( pedpix.GetPedestalRms()* fNumLoGainSamples
256 / TMath::Sqrt((Float_t)fPedestals->GetTotalEntries()), i );
257 fBlindPixel->SetPedRms ( pedpix.GetPedestalRms()* TMath::Sqrt((Float_t)fNumLoGainSamples), i );
258 fBlindPixel->SetPedRmsErr( fBlindPixel->GetPedErr()/2., i );
259 }
260 }
261
262 const Int_t firstdesired = (Int_t)fHiGainFirst;
263 Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
264
265 if (firstdesired > lastavailable)
266 {
267 const Int_t diff = firstdesired - lastavailable;
268 *fLog << endl;
269 *fLog << warn << GetDescriptor()
270 << Form("%s%2i%s%2i%s",": Selected First Hi Gain FADC slice ",
271 (int)fHiGainFirst,
272 " ranges out of the available limits: [0,",lastavailable,"].") << endl;
273 *fLog << warn << GetDescriptor()
274 << Form("%s%2i%s",": Will start with slice ",diff," of the Low-Gain for the High-Gain extraction")
275 << endl;
276
277 fHiLoFirst = diff;
278 }
279
280 const Int_t lastdesired = (Int_t)fHiGainLast;
281
282 if (lastdesired > lastavailable)
283 {
284 Int_t diff = lastdesired - lastavailable;
285 lastavailable += (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
286
287 if (lastdesired > lastavailable)
288 {
289 *fLog << endl;
290 *fLog << warn << GetDescriptor()
291 << Form("%s%2i%s%2i%s",": Selected Last Hi Gain FADC slice ",
292 (int)fHiGainLast,
293 " ranges out of the available limits: [0,",lastavailable,"].") << endl;
294 *fLog << warn << GetDescriptor()
295 << Form("%s%2i",": Will reduce upper limit by ",diff)
296 << endl;
297 diff = (Int_t)fRunHeader->GetNumSamplesLoGain();
298 }
299
300 fHiGainLast = (Int_t)fRunHeader->GetNumSamplesHiGain() - 1;
301 fHiLoLast = diff;
302 }
303
304 const Int_t range = fHiLoFirst ? fHiLoLast - fHiLoFirst + 1 : fHiGainLast - fHiGainFirst + fHiLoLast + 1;
305
306 fHiGainSignal = new Float_t[range];
307 memset(fHiGainSignal,0,range*sizeof(Float_t));
308 fHiGainFirstDeriv = new Float_t[range];
309 memset(fHiGainFirstDeriv,0,range*sizeof(Float_t));
310 fHiGainSecondDeriv = new Float_t[range];
311 memset(fHiGainSecondDeriv,0,range*sizeof(Float_t));
312
313 *fLog << endl;
314 *fLog << inf << GetDescriptor() << ": Extracting "
315 << Form("%s",IsExtractionType(kAmplitude) ? "Amplitude " : " Integral ")
316 << " using " << range << " FADC samples from "
317 << Form("%s%2i",fHiLoFirst ? "Low Gain slice " : " High Gain slice ",
318 fHiLoFirst ? (Int_t)fHiLoFirst : (Int_t)fHiGainFirst)
319 << " to (including) "
320 << Form("%s%2i",fHiLoLast ? "Low Gain slice " : " High Gain slice ",
321 fHiLoLast ? (Int_t)fHiLoLast-1 : (Int_t)fHiGainLast )
322 << endl;
323
324 if (IsExtractionType(kFilter))
325 *fLog << inf << GetDescriptor() << ": Will use Filter using "
326 << (Int_t)(fLoGainLast-fLoGainFirst+1) << " FADC slices "
327 << "from High Gain slice " << (Int_t)fLoGainFirst
328 << " to High Gain slice " << (Int_t)fLoGainLast << endl;
329
330 fBlindPixel->SetUsedFADCSlices(fHiGainFirst, range);
331
332 return kTRUE;
333
334}
335
336// --------------------------------------------------------------------------
337//
338// FindSignalHiGain:
339//
340// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
341// - Sum up contents of *ptr
342// - If *ptr is greater than fSaturationLimit, raise sat by 1
343// - If fHiLoLast is set, loop from logain to (logain+fHiLoLast)
344// - Add contents of *logain to sum
345//
346void MExtractBlindPixel::FindIntegral(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat)
347{
348
349 Int_t summ = 0;
350 Byte_t *p = ptr;
351 Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
352
353 if (fHiLoFirst == 0)
354 {
355
356 while (p<end)
357 {
358 summ += *ptr;
359
360 if (*p++ >= fSaturationLimit)
361 sat++;
362 }
363
364 }
365
366 p = logain + fHiLoFirst;
367 end = logain + fHiLoLast;
368 while (p<end)
369 {
370 summ += *p;
371
372 if (*p++ >= fSaturationLimit)
373 sat++;
374 }
375
376 sum = (Float_t)summ;
377}
378
379// --------------------------------------------------------------------------
380//
381// FindSignalPhe:
382//
383// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
384// - Sum up contents of *ptr
385// - If *ptr is greater than fSaturationLimit, raise sat by 1
386// - If fHiLoLast is set, loop from logain to (logain+fHiLoLast)
387// - Add contents of *logain to sum
388//
389void MExtractBlindPixel::FindAmplitude(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat)
390{
391
392 Int_t range = 0;
393 Int_t count = 0;
394 Float_t abmaxpos = 0.;
395 Byte_t *p = ptr;
396 Byte_t *end;
397 Byte_t max = 0;
398 Byte_t maxpos = 0;
399 Int_t summ = 0;
400
401 if (fHiLoFirst == 0)
402 {
403
404 range = fHiGainLast - fHiGainFirst + 1;
405
406 end = ptr + range;
407 //
408 // Check for saturation in all other slices
409 //
410 while (p++<end)
411 {
412
413 fHiGainSignal[count] = (Float_t)*p;
414 summ += *p;
415
416 if (*p > max)
417 {
418 max = *p;
419 maxpos = count;
420 }
421
422 count++;
423
424 if (*p >= fSaturationLimit)
425 sat++;
426 }
427 }
428
429 if (fHiLoLast != 0)
430 {
431
432 p = logain + fHiLoFirst;
433 end = logain + fHiLoLast;
434
435 while (p<end)
436 {
437
438 fHiGainSignal[count] = (Float_t)*p;
439 summ += *p;
440
441 if (*p > max)
442 {
443 max = *p;
444 maxpos = count;
445 }
446
447 range++;
448 count++;
449
450 if (*p++ >= fSaturationLimit)
451 sat++;
452 }
453 }
454
455 //
456 // allow one saturated slice
457 //
458 if (sat > 1)
459 {
460 sum = gkOverflow;
461 return;
462 }
463
464 //
465 // Don't start if the maxpos is too close to the left limit.
466 //
467 if (maxpos < 2)
468 {
469 sum = (Float_t)max;
470 return;
471 }
472
473 Float_t pp;
474
475 for (Int_t i=1;i<range-1;i++)
476 {
477 pp = fHiGainSecondDeriv[i-1] + 4.;
478 fHiGainSecondDeriv[i] = -1.0/pp;
479 fHiGainFirstDeriv [i] = fHiGainSignal[i+1] - fHiGainSignal[i] - fHiGainSignal[i] + fHiGainSignal[i-1];
480 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
481 p++;
482 }
483
484 fHiGainSecondDeriv[range-1] = 0.;
485 for (Int_t k=range-2;k>=0;k--)
486 fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
487
488 //
489 // Now find the maximum
490 //
491 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
492 Float_t lower = (Float_t)maxpos-1.;
493 Float_t upper = (Float_t)maxpos;
494 Float_t x = lower;
495 Float_t y = 0.;
496 Float_t a = 1.;
497 Float_t b = 0.;
498 Int_t klo = maxpos-1;
499 Int_t khi = maxpos;
500 Float_t klocont = fHiGainSignal[klo];
501 Float_t khicont = fHiGainSignal[khi];
502 sum = (Float_t)khicont;
503 abmaxpos = lower;
504
505 //
506 // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
507 // interval maxpos+1.
508 //
509 while (x<upper-0.3)
510 {
511
512 x += step;
513 a -= step;
514 b += step;
515
516 y = a*klocont
517 + b*khicont
518 + (a*a*a-a)*fHiGainSecondDeriv[klo]
519 + (b*b*b-b)*fHiGainSecondDeriv[khi];
520
521 if (y > sum)
522 {
523 sum = y;
524 abmaxpos = x;
525 }
526 }
527
528 if (abmaxpos > upper-0.1)
529 {
530
531 upper = (Float_t)maxpos+1;
532 lower = (Float_t)maxpos;
533 x = lower;
534 a = 1.;
535 b = 0.;
536 khi = maxpos+1;
537 klo = maxpos;
538 klocont = fHiGainSignal[klo];
539 khicont = fHiGainSignal[khi];
540
541 while (x<upper-0.3)
542 {
543
544 x += step;
545 a -= step;
546 b += step;
547
548 y = a* klocont
549 + b* khicont
550 + (a*a*a-a)*fHiGainSecondDeriv[klo]
551 + (b*b*b-b)*fHiGainSecondDeriv[khi];
552
553 if (y > sum)
554 {
555 sum = y;
556 abmaxpos = x;
557 }
558 }
559 }
560
561 const Float_t up = abmaxpos+step-0.055;
562 const Float_t lo = abmaxpos-step+0.055;
563 const Float_t maxpossave = abmaxpos;
564
565 x = abmaxpos;
566 a = upper - x;
567 b = x - lower;
568
569 step = 0.04; // step size of 83 ps
570
571 while (x<up)
572 {
573
574 x += step;
575 a -= step;
576 b += step;
577
578 y = a* klocont
579 + b* khicont
580 + (a*a*a-a)*fHiGainSecondDeriv[klo]
581 + (b*b*b-b)*fHiGainSecondDeriv[khi];
582
583 if (y > sum)
584 {
585 sum = y;
586 abmaxpos = x;
587 }
588 }
589
590 if (abmaxpos < klo + 0.02)
591 {
592 klo--;
593 khi--;
594 klocont = fHiGainSignal[klo];
595 khicont = fHiGainSignal[khi];
596 upper--;
597 lower--;
598 }
599
600 x = maxpossave;
601 a = upper - x;
602 b = x - lower;
603
604 while (x>lo)
605 {
606
607 x -= step;
608 a += step;
609 b -= step;
610
611 y = a* klocont
612 + b* khicont
613 + (a*a*a-a)*fHiGainSecondDeriv[klo]
614 + (b*b*b-b)*fHiGainSecondDeriv[khi];
615
616 if (y > sum)
617 {
618 sum = y;
619 abmaxpos = x;
620 }
621 }
622
623}
624
625// --------------------------------------------------------------------------
626//
627// FindSignalFilter:
628//
629// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
630// - Sum up contents of *ptr
631// - If *ptr is greater than fSaturationLimit, raise sat by 1
632//
633void MExtractBlindPixel::FindSignalFilter(Byte_t *ptr, Int_t &sum, Byte_t &sat) const
634{
635
636 Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1;
637
638 while (ptr<end)
639 {
640 sum += *ptr;
641
642 if (*ptr++ >= fSaturationLimit)
643 sat++;
644 }
645}
646
647// --------------------------------------------------------------------------
648//
649// Calculate the integral of the FADC time slices and store them as a new
650// pixel in the MExtractedBlindPixel container.
651//
652Int_t MExtractBlindPixel::Process()
653{
654
655 MRawEvtPixelIter pixel(fRawEvt);
656
657 fBlindPixel->Clear();
658
659 for (Int_t id=0;id<fBlindPixelIdx.GetSize();id++)
660 {
661
662 pixel.Jump(fBlindPixelIdx[id]);
663
664 Int_t sum = 0;
665 Byte_t sat = 0;
666
667 if (IsExtractionType(kFilter))
668 {
669
670 FindSignalFilter(pixel.GetHiGainSamples()+fLoGainFirst, sum, sat);
671
672 if (sum > fNSBFilterLimit)
673 {
674 fBlindPixel->SetExtractedSignal(-1.,id);
675 fBlindPixel->SetNumSaturated(sat,id);
676 fBlindPixel->SetReadyToSave();
677 continue;
678 }
679
680 sum = 0;
681 FindSignalFilter(pixel.GetLoGainSamples(), sum, sat);
682
683 /*
684 if (fModified)
685 {
686 if (sum > fNSBFilterLimit)
687 {
688 fBlindPixel->SetExtractedSignal(-1.,id);
689 fBlindPixel->SetNumSaturated(sat,id);
690 fBlindPixel->SetReadyToSave();
691 continue;
692 }
693 }
694 */
695 }
696
697 Float_t newsum = 0.;
698 sat = 0;
699
700 if (IsExtractionType(kAmplitude))
701 FindAmplitude (pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), newsum, sat);
702 else
703 FindIntegral (pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), newsum, sat);
704
705
706 fBlindPixel->SetExtractedSignal(newsum,id);
707 fBlindPixel->SetNumSaturated(sat,id);
708 }
709
710 fBlindPixel->SetReadyToSave();
711 return kTRUE;
712}
713
714// ------------------------------------------------------------------------------------
715//
716// Returns true if the extraction type. Available are: kAmplitude, kIntegral and kFilter
717// The flags kIntegral and kFilter may be set both.
718//
719Bool_t MExtractBlindPixel::IsExtractionType( const ExtractionType_t typ )
720{
721
722 return TESTBIT( fExtractionType, typ );
723
724}
725
726// --------------------------------------------------------------------------
727//
728// Sets the extraction type. Available are: kAmplitude and kIntegral
729//
730void MExtractBlindPixel::SetExtractionType( const ExtractionType_t typ )
731{
732 SETBIT( fExtractionType, typ );
733}
Note: See TracBrowser for help on using the repository browser.