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

Last change on this file since 4789 was 4771, checked in by gaug, 21 years ago
*** empty log message ***
File size: 18.3 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 << "First Hi Gain slice " << (int)fHiGainFirst << " out of range [0,";
270 *fLog << lastavailable << "]... start at slice " << diff << " of the Lo Gain " << endl;
271
272 fHiLoFirst = diff;
273 }
274
275 const Int_t lastdesired = (Int_t)fHiGainLast;
276
277 if (lastdesired > lastavailable)
278 {
279 Int_t diff = lastdesired - lastavailable;
280 lastavailable += (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
281
282 if (lastdesired > lastavailable)
283 {
284 *fLog << endl;
285 *fLog << "Last Hi Gain slice " << (int)fHiGainLast << " out of range [0,";
286 *fLog << lastavailable << "]... reduce upper limit by " << diff << endl;
287 diff = (Int_t)fRunHeader->GetNumSamplesLoGain();
288 }
289
290 fHiGainLast = (Int_t)fRunHeader->GetNumSamplesHiGain() - 1;
291 fHiLoLast = diff;
292 }
293
294 const Int_t range = fHiLoFirst ? fHiLoLast - fHiLoFirst + 1 : fHiGainLast - fHiGainFirst + fHiLoLast + 1;
295
296 fHiGainSignal = new Float_t[range];
297 memset(fHiGainSignal,0,range*sizeof(Float_t));
298 fHiGainFirstDeriv = new Float_t[range];
299 memset(fHiGainFirstDeriv,0,range*sizeof(Float_t));
300 fHiGainSecondDeriv = new Float_t[range];
301 memset(fHiGainSecondDeriv,0,range*sizeof(Float_t));
302
303 *fLog << endl;
304 *fLog << inf << ": Extracting "
305 << Form("%s",IsExtractionType(kAmplitude) ? "Amplitude " : " Integral ")
306 << " using " << range << " FADC samples from "
307 << Form("%s%2i",fHiLoFirst ? "Low Gain slice " : " High Gain slice ",
308 fHiLoFirst ? (Int_t)fHiLoFirst : (Int_t)fHiGainFirst)
309 << " to (including) "
310 << Form("%s%2i",fHiLoLast ? "Low Gain slice " : " High Gain slice ",
311 fHiLoLast ? (Int_t)fHiLoLast-1 : (Int_t)fHiGainLast )
312 << endl;
313
314 if (IsExtractionType(kFilter))
315 *fLog << inf << GetDescriptor() << ": Will use Filter using "
316 << (Int_t)(fLoGainLast-fLoGainFirst+1) << " FADC slices "
317 << "from High Gain slice " << (Int_t)fLoGainFirst
318 << " to High Gain slice " << (Int_t)fLoGainLast << endl;
319
320 fBlindPixel->SetUsedFADCSlices(fHiGainFirst, range);
321
322 return kTRUE;
323
324}
325
326// --------------------------------------------------------------------------
327//
328// FindSignalHiGain:
329//
330// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
331// - Sum up contents of *ptr
332// - If *ptr is greater than fSaturationLimit, raise sat by 1
333// - If fHiLoLast is set, loop from logain to (logain+fHiLoLast)
334// - Add contents of *logain to sum
335//
336void MExtractBlindPixel::FindIntegral(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat)
337{
338
339 Int_t summ = 0;
340 Byte_t *p = ptr;
341 Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
342
343 if (fHiLoFirst == 0)
344 {
345
346 while (p<end)
347 {
348 summ += *ptr;
349
350 if (*p++ >= fSaturationLimit)
351 sat++;
352 }
353
354 }
355
356 p = logain + fHiLoFirst;
357 end = logain + fHiLoLast;
358 while (p<end)
359 {
360 summ += *p;
361
362 if (*p++ >= fSaturationLimit)
363 sat++;
364 }
365
366 sum = (Float_t)summ;
367}
368
369// --------------------------------------------------------------------------
370//
371// FindSignalPhe:
372//
373// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
374// - Sum up contents of *ptr
375// - If *ptr is greater than fSaturationLimit, raise sat by 1
376// - If fHiLoLast is set, loop from logain to (logain+fHiLoLast)
377// - Add contents of *logain to sum
378//
379void MExtractBlindPixel::FindAmplitude(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat)
380{
381
382 Int_t range = 0;
383 Int_t count = 0;
384 Float_t abmaxpos = 0.;
385 Byte_t *p = ptr;
386 Byte_t *end;
387 Byte_t max = 0;
388 Byte_t maxpos = 0;
389 Int_t summ = 0;
390
391 if (fHiLoFirst == 0)
392 {
393
394 range = fHiGainLast - fHiGainFirst + 1;
395
396 end = ptr + range;
397 //
398 // Check for saturation in all other slices
399 //
400 while (p++<end)
401 {
402
403 fHiGainSignal[count] = (Float_t)*p;
404 summ += *p;
405
406 if (*p > max)
407 {
408 max = *p;
409 maxpos = count;
410 }
411
412 count++;
413
414 if (*p >= fSaturationLimit)
415 sat++;
416 }
417 }
418
419 if (fHiLoLast != 0)
420 {
421
422 p = logain + fHiLoFirst;
423 end = logain + fHiLoLast;
424
425 while (p<end)
426 {
427
428 fHiGainSignal[count] = (Float_t)*p;
429 summ += *p;
430
431 if (*p > max)
432 {
433 max = *p;
434 maxpos = count;
435 }
436
437 range++;
438 count++;
439
440 if (*p++ >= fSaturationLimit)
441 sat++;
442 }
443 }
444
445 //
446 // allow one saturated slice
447 //
448 if (sat > 1)
449 {
450 sum = gkOverflow;
451 return;
452 }
453
454 //
455 // Don't start if the maxpos is too close to the left limit.
456 //
457 if (maxpos < 2)
458 {
459 sum = (Float_t)max;
460 return;
461 }
462
463 Float_t pp;
464
465 for (Int_t i=1;i<range-1;i++)
466 {
467 pp = fHiGainSecondDeriv[i-1] + 4.;
468 fHiGainSecondDeriv[i] = -1.0/pp;
469 fHiGainFirstDeriv [i] = fHiGainSignal[i+1] - fHiGainSignal[i] - fHiGainSignal[i] + fHiGainSignal[i-1];
470 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
471 p++;
472 }
473
474 fHiGainSecondDeriv[range-1] = 0.;
475 for (Int_t k=range-2;k>=0;k--)
476 fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
477
478 //
479 // Now find the maximum
480 //
481 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
482 Float_t lower = (Float_t)maxpos-1.;
483 Float_t upper = (Float_t)maxpos;
484 Float_t x = lower;
485 Float_t y = 0.;
486 Float_t a = 1.;
487 Float_t b = 0.;
488 Int_t klo = maxpos-1;
489 Int_t khi = maxpos;
490 Float_t klocont = fHiGainSignal[klo];
491 Float_t khicont = fHiGainSignal[khi];
492 sum = (Float_t)khicont;
493 abmaxpos = lower;
494
495 //
496 // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
497 // interval maxpos+1.
498 //
499 while (x<upper-0.3)
500 {
501
502 x += step;
503 a -= step;
504 b += step;
505
506 y = a*klocont
507 + b*khicont
508 + (a*a*a-a)*fHiGainSecondDeriv[klo]
509 + (b*b*b-b)*fHiGainSecondDeriv[khi];
510
511 if (y > sum)
512 {
513 sum = y;
514 abmaxpos = x;
515 }
516 }
517
518 if (abmaxpos > upper-0.1)
519 {
520
521 upper = (Float_t)maxpos+1;
522 lower = (Float_t)maxpos;
523 x = lower;
524 a = 1.;
525 b = 0.;
526 khi = maxpos+1;
527 klo = maxpos;
528 klocont = fHiGainSignal[klo];
529 khicont = fHiGainSignal[khi];
530
531 while (x<upper-0.3)
532 {
533
534 x += step;
535 a -= step;
536 b += step;
537
538 y = a* klocont
539 + b* khicont
540 + (a*a*a-a)*fHiGainSecondDeriv[klo]
541 + (b*b*b-b)*fHiGainSecondDeriv[khi];
542
543 if (y > sum)
544 {
545 sum = y;
546 abmaxpos = x;
547 }
548 }
549 }
550
551 const Float_t up = abmaxpos+step-0.055;
552 const Float_t lo = abmaxpos-step+0.055;
553 const Float_t maxpossave = abmaxpos;
554
555 x = abmaxpos;
556 a = upper - x;
557 b = x - lower;
558
559 step = 0.04; // step size of 83 ps
560
561 while (x<up)
562 {
563
564 x += step;
565 a -= step;
566 b += step;
567
568 y = a* klocont
569 + b* khicont
570 + (a*a*a-a)*fHiGainSecondDeriv[klo]
571 + (b*b*b-b)*fHiGainSecondDeriv[khi];
572
573 if (y > sum)
574 {
575 sum = y;
576 abmaxpos = x;
577 }
578 }
579
580 if (abmaxpos < klo + 0.02)
581 {
582 klo--;
583 khi--;
584 klocont = fHiGainSignal[klo];
585 khicont = fHiGainSignal[khi];
586 upper--;
587 lower--;
588 }
589
590 x = maxpossave;
591 a = upper - x;
592 b = x - lower;
593
594 while (x>lo)
595 {
596
597 x -= step;
598 a += step;
599 b -= step;
600
601 y = a* klocont
602 + b* khicont
603 + (a*a*a-a)*fHiGainSecondDeriv[klo]
604 + (b*b*b-b)*fHiGainSecondDeriv[khi];
605
606 if (y > sum)
607 {
608 sum = y;
609 abmaxpos = x;
610 }
611 }
612
613}
614
615// --------------------------------------------------------------------------
616//
617// FindSignalFilter:
618//
619// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
620// - Sum up contents of *ptr
621// - If *ptr is greater than fSaturationLimit, raise sat by 1
622//
623void MExtractBlindPixel::FindSignalFilter(Byte_t *ptr, Int_t &sum, Byte_t &sat) const
624{
625
626 Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1;
627
628 while (ptr<end)
629 {
630 sum += *ptr;
631
632 if (*ptr++ >= fSaturationLimit)
633 sat++;
634 }
635}
636
637// --------------------------------------------------------------------------
638//
639// Calculate the integral of the FADC time slices and store them as a new
640// pixel in the MExtractedBlindPixel container.
641//
642Int_t MExtractBlindPixel::Process()
643{
644
645 MRawEvtPixelIter pixel(fRawEvt);
646
647 fBlindPixel->Clear();
648
649 for (Int_t id=0;id<fBlindPixelIdx.GetSize();id++)
650 {
651
652 pixel.Jump(fBlindPixelIdx[id]);
653
654 Int_t sum = 0;
655 Byte_t sat = 0;
656
657 if (IsExtractionType(kFilter))
658 {
659
660 FindSignalFilter(pixel.GetHiGainSamples()+fLoGainFirst, sum, sat);
661
662 if (sum > fNSBFilterLimit)
663 {
664 fBlindPixel->SetExtractedSignal(-1.,id);
665 fBlindPixel->SetNumSaturated(sat,id);
666 fBlindPixel->SetReadyToSave();
667 continue;
668 }
669
670 sum = 0;
671 FindSignalFilter(pixel.GetLoGainSamples(), sum, sat);
672
673 /*
674 if (fModified)
675 {
676 if (sum > fNSBFilterLimit)
677 {
678 fBlindPixel->SetExtractedSignal(-1.,id);
679 fBlindPixel->SetNumSaturated(sat,id);
680 fBlindPixel->SetReadyToSave();
681 continue;
682 }
683 }
684 */
685 }
686
687 Float_t newsum = 0.;
688 sat = 0;
689
690 if (IsExtractionType(kAmplitude))
691 FindAmplitude (pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), newsum, sat);
692 else
693 FindIntegral (pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), newsum, sat);
694
695
696 fBlindPixel->SetExtractedSignal(newsum,id);
697 fBlindPixel->SetNumSaturated(sat,id);
698 }
699
700 fBlindPixel->SetReadyToSave();
701 return kTRUE;
702}
703
704// ------------------------------------------------------------------------------------
705//
706// Returns true if the extraction type. Available are: kAmplitude, kIntegral and kFilter
707// The flags kIntegral and kFilter may be set both.
708//
709Bool_t MExtractBlindPixel::IsExtractionType( const ExtractionType_t typ )
710{
711
712 return TESTBIT( fExtractionType, typ );
713
714}
715
716// --------------------------------------------------------------------------
717//
718// Sets the extraction type. Available are: kAmplitude and kIntegral
719//
720void MExtractBlindPixel::SetExtractionType( const ExtractionType_t typ )
721{
722 SETBIT( fExtractionType, typ );
723}
Note: See TracBrowser for help on using the repository browser.