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

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