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

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