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

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