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

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