source: trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc@ 8359

Last change on this file since 8359 was 8358, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 15.0 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.63 2007-03-03 22:46:10 tbretz Exp $
3! --------------------------------------------------------------------------
4!
5! *
6! * This file is part of MARS, the MAGIC Analysis and Reconstruction
7! * Software. It is distributed to you in the hope that it can be a useful
8! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
9! * It is distributed WITHOUT ANY WARRANTY.
10! *
11! * Permission to use, copy, modify and distribute this software and its
12! * documentation for any purpose is hereby granted without fee,
13! * provided that the above copyright notice appear in all copies and
14! * that both that copyright notice and this permission notice appear
15! * in supporting documentation. It is provided "as is" without express
16! * or implied warranty.
17! *
18!
19!
20! Author(s): Markus Gaug, 05/2004 <mailto:markus@ifae.es>
21! Author(s): Thomas Bretz, 05/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
22!
23! Copyright: MAGIC Software Development, 2000-2006
24!
25!
26\* ======================================================================== */
27
28//////////////////////////////////////////////////////////////////////////////
29//
30// MExtractTimeAndCharge
31//
32// Base class for the signal extractors which extract the arrival time
33// and the signal at the same time. Uses the functions
34// FindTimeAndChargeHiGain() and FindTimeAndChargeLoGain() to extract
35// the signal.
36//
37// The following figure gives and example of possible inheritance trees.
38// An extractor class can inherit from each of the following base classes:
39// - MExtractor
40// - MExtractTime
41// - MExtractTimeAndCharge
42//
43//Begin_Html
44/*
45<img src="images/ExtractorClasses.gif">
46*/
47//End_Html
48//
49// The following variables have to be set by the derived class and
50// do not have defaults:
51// - fNumHiGainSamples
52// - fNumLoGainSamples
53// - fSqrtHiGainSamples
54// - fSqrtLoGainSamples
55//
56//
57// Class Version 3:
58// ----------------
59// - Byte_t fMaxBinContent;
60//
61//
62// Input Containers:
63// MRawEvtData
64// MRawRunHeader
65// MPedestalCam
66//
67// Output Containers:
68// MArrivalTimeCam
69// MExtractedSignalCam
70//
71// For weired events check: Run 94127 Event 672, 1028
72//
73//////////////////////////////////////////////////////////////////////////////
74#include "MExtractTimeAndCharge.h"
75
76#include <TRandom.h>
77#include <TVector3.h>
78
79#include "MMath.h"
80
81#include "MLog.h"
82#include "MLogManip.h"
83
84#include "MParList.h"
85
86#include "MRawRunHeader.h"
87#include "MRawEvtPixelIter.h"
88
89#include "MArrivalTimeCam.h"
90#include "MArrivalTimePix.h"
91
92#include "MExtractedSignalCam.h"
93#include "MExtractedSignalPix.h"
94
95#include "MPedestalSubtractedEvt.h"
96
97ClassImp(MExtractTimeAndCharge);
98
99using namespace std;
100
101const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -2.5;
102const Byte_t MExtractTimeAndCharge::fgLoGainSwitch = 120;
103
104// --------------------------------------------------------------------------
105//
106// Default constructor.
107//
108// Sets:
109// - fWindowSizeHiGain and fWindowSizeLoGain to 0
110// - fLoGainStartShift to fgLoGainStartShift
111// - fLoGainSwitch to fgLoGainSwitch
112//
113MExtractTimeAndCharge::MExtractTimeAndCharge(const char *name, const char *title)
114 : fWindowSizeHiGain(0), fWindowSizeLoGain(0)
115{
116 fName = name ? name : "MExtractTimeAndCharge";
117 fTitle = title ? title : "Base class for signal and time extractors";
118
119 SetLoGainStartShift();
120 SetLoGainSwitch();
121}
122
123// --------------------------------------------------------------------------
124//
125// The PreProcess searches for the following input containers:
126// - MRawEvtData
127// - MRawRunHeader
128// - MPedestalCam
129//
130// The following output containers are also searched and created if
131// they were not found:
132//
133// - MExtractedSignalCam
134// - MArrivalTimeCam
135//
136Int_t MExtractTimeAndCharge::PreProcess(MParList *pList)
137{
138 if (!MExtractTime::PreProcess(pList))
139 return kFALSE;
140
141 fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam));
142 if (!fSignals)
143 return kFALSE;
144
145 *fLog << flush << inf;
146 return kTRUE;
147}
148
149// --------------------------------------------------------------------------
150//
151// The ReInit calls:
152// - MExtractor::ReInit()
153//
154// Call:
155// - MArrivalTimeCam::SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
156// fLoGainFirst, fLoGainLast, fNumLoGainSamples);
157// - InitArrays();
158//
159Bool_t MExtractTimeAndCharge::ReInit(MParList *pList)
160{
161 if (!MExtractTime::ReInit(pList))
162 return kFALSE;
163
164 if (!InitArrays(fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain()))
165 return kFALSE;
166
167 if (fSignals)
168 fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
169 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
170
171 if (!HasLoGain())
172 {
173 *fLog << inf << "No lo-gains... resetting lo-gain switch." << endl;
174 fLoGainSwitch=0xff;
175 }
176
177 return kTRUE;
178}
179
180// --------------------------------------------------------------------------
181//
182// Return the x-value lower than sat0 at which the signal has been
183// fallen bwlow maxcont/2. This time is determined using a simple second
184// order polynomial interpolation.
185//
186Double_t MExtractTimeAndCharge::GetSaturationTime(Int_t sat0, const Float_t *sig, Int_t maxcont) const
187{
188 const Int_t p = sat0>1 ? sat0-2 : sat0-1;
189 if (sat0<=0)
190 return 0;
191 if (sat0==1)
192 return sig[0]>maxcont/2 ? 0 : 0.5;
193
194 if (sig[p]>sig[p+1] || sig[p+1]>sig[p+2])
195 return sig[p+1]>maxcont/2 ? sat0-1 : sat0-0.5;
196
197 // Find the place at which the signal is maxcont/2
198 const TVector3 vx(sig[p], sig[p+1], sig[p+2]);
199 const TVector3 vy(p, p+1, p+2);
200
201 return MMath::InterpolParabLin(vx, vy, maxcont/2);
202}
203
204// --------------------------------------------------------------------------
205//
206// Calculate the integral of the FADC time slices and store them as a new
207// pixel in the MArrivalTimeCam container.
208// Calculate the integral of the FADC time slices and store them as a new
209// pixel in the MExtractedSignalCam container.
210// The functions FindTimeAndChargeHiGain() and FindTimeAndChargeLoGain() are
211// supposed to extract the signal themselves.
212//
213Int_t MExtractTimeAndCharge::Process()
214{
215 const Int_t numh = fRunHeader->GetNumSamplesHiGain();
216 const Int_t numl = fRunHeader->GetNumSamplesLoGain();
217
218 MRawEvtPixelIter pixel(fRawEvt);
219 while (pixel.Next())
220 {
221 const Int_t pixidx = pixel.GetPixelId();
222
223 const Float_t *sig = fSignal->GetSamples(pixidx);
224
225 Int_t sathi0 = fHiGainFirst; // First slice to extract and first saturating slice
226 Int_t sathi1 = fHiGainLast; // Last slice to extract and last saturating slice
227
228 Int_t maxcont;
229 Int_t maxposhi = fSignal->GetMax(pixidx, sathi0, sathi1, maxcont);
230 // Would it be better to take lastsat-firstsat?
231 Int_t numsathi = fSignal->GetSaturation(pixidx, fSaturationLimit, sathi0, sathi1);
232
233 Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid
234 Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid
235
236 // Do not even try to extract the hi-gain if we have
237 // more than one saturating slice
238 const Int_t rangehi = fHiGainLast - fHiGainFirst + 1;
239
240 if (numsathi<2)
241 FindTimeAndChargeHiGain2(sig+fHiGainFirst, rangehi,
242 sumhi, deltasumhi, timehi, deltatimehi,
243 numsathi, maxposhi);
244
245 // If we have saturating slices try to get a better estimate
246 // of the arrival time than timehi or sathi0. This is
247 // usefull to know where to start lo-gain extraction.
248 if (numsathi>1)
249 {
250 timehi = GetSaturationTime(sathi0, sig, maxcont)-fHiGainFirst;
251 deltatimehi = 0;
252 }
253
254 // Make sure that in cases the time couldn't be correctly determined
255 // more meaningfull default values are assigned.
256 // For extractors like the digital filter and the spline
257 // we allow extracpolation by one slice.
258 // FIXME: Defined Out-Of-Range better so that the extractors
259 // know what to return!
260 if (deltatimehi>-0.5 && (timehi<-1 || timehi>=rangehi))
261 {
262 // Flag this as unreliable!
263 timehi = gRandom->Uniform(rangehi+1)-1;
264 // deltatimehi=-1; // Set PIXEL to UNRELIABLE?
265 }
266
267 timehi += fHiGainFirst;
268
269 Float_t sumlo =0, deltasumlo =-1; // invalidate logain of MExtractedSignalPix
270 Float_t timelo=0, deltatimelo=-1; // invalidate logain of MArrivalTimePix
271 Int_t numsatlo=0;
272
273 //
274 // Adapt the low-gain extraction range from the obtained high-gain time
275 //
276
277 // IN THIS CASE THE PIXEL SHOULD BE MARKED BAD!!!!
278 // MEANS: Hi gain has saturated, but the signal is to dim
279 // to extract the lo-gain properly
280 // This could happen because the maxcont was not extracted from
281 // all slices!
282 // THIS produces pulse positions ~= -1
283 // The signal might be handled in MCalibrateData, but hwat's about
284 // the arrival times in MCalibrateRelTime
285 if (numsathi>0 && maxcont<=fLoGainSwitch)
286 deltasumlo=deltasumhi=deltatimelo=deltatimehi=-1;
287
288 // If more than 8 hi-gain slices have saturated this is
289 // no physical event. We just assume that something with
290 // the extraction is wrong
291 if (numsathi>8) // FIXME: Should be something like 2?
292 deltasumhi=deltatimehi=-1;
293
294 // FIXME: What to do with the pixel if it saturates too early???
295 if (numl>0 && maxcont>fLoGainSwitch)
296 {
297 Int_t first = numh+fLoGainFirst;
298 Int_t last = numh+fLoGainLast;
299
300 // To determin the window in which the lo-gain is extracted
301 // clearly more information about the relation between the
302 // extraction window and the reslting time is necessary.
303 //
304 // numh + fLoGainStartShift == 14 / fLoGainStartShift=-1
305 //
306 // The lo-gain is expected to have its raising edge
307 // at timehi+numh+fOffsetLoGain. We start to search for the
308 // lo-gain fLoGainStartShift slices earlier.
309 //
310 // Instead of fLoGainStartShift the extractor should now how many
311 // slices before the expected raising edge the start of the
312 // search must be placed and from there we can step 1.5 slices
313 // back to be on the safe side.
314 //
315 // The jitter in the hi-/lo-gain offset ssems to be around +/-0.5
316 const Float_t tm = deltatimehi<0 ? -1.+fHiGainFirst : timehi;
317 first = TMath::FloorNint(tm+numh+fOffsetLoGain+fLoGainStartShift);
318
319 if (first<0)
320 first = 0;
321 if (first>last)
322 first=last;
323
324 /*
325 // Currently we have to stick to this check because at least
326 // the spline has arrays of this size...
327 if (first>last)
328 first = last;
329 if (first<numh+fLoGainFirst)
330 first = numh+fLoGainFirst;
331 */
332 Int_t satlo0 = first; // First slice to extract and first saturating slice
333 Int_t satlo1 = last; // Last slice to extract and last saturating slice
334
335 // Would it be better to take lastsat-firstsat?
336 Int_t maxlo;
337 Int_t maxposlo = fSignal->GetMax(pixidx, satlo0, satlo1, maxlo);
338 numsatlo = fSignal->GetSaturation(pixidx, fSaturationLimit, satlo0, satlo1);
339
340 //if (satlo0>first && satlo1<last && numsatlo>2)
341 //{
342 // fSignal->InterpolateSaturation(pixidx, fSaturationLimit, satlo0, satlo1);
343 // numsatlo = 0;
344 //}
345
346 const Int_t rangelo = last-first+1;
347 FindTimeAndChargeLoGain2(sig+first, rangelo,
348 sumlo, deltasumlo, timelo, deltatimelo,
349 numsatlo, maxposlo);
350
351 // If we have saturating slices try to get a better estimate
352 // of the arrival time than timehi or sathi0. This is
353 // usefull to know where to start lo-gain extraction.
354 if (numsatlo>1)
355 {
356 timelo = GetSaturationTime(satlo0, sig, maxlo)-numh-first;
357 deltatimelo = 0;
358 }
359
360 // Make sure that in cases the time couldn't be correctly determined
361 // more meaningfull default values are assigned
362 // For extractors like the digital filter and the spline
363 // we allow extracpolation by one slice.
364 if (deltatimelo>-0.5 && (timelo<-1 || timelo>=rangelo))
365 {
366 // Flag this as unreliable!
367 timelo = gRandom->Uniform(rangelo+1)-1;
368 //deltatimelo=-1; // Set PIXEL to UNRELIABLE?
369 }
370
371 timelo += first-numh;
372
373 // If more than 6 lo-gain slices have saturated this is
374 // no physical event. We just assume that something with
375 // the extraction is wrong
376 if (numsatlo>6)
377 deltasumlo=deltatimelo=-1;
378
379 // The extracted lo-gain signal cannot be zero or
380 // negative at all, so it must be wrong
381 if (sumlo<=0)
382 deltasumlo=-1;
383
384 //if (TMath::Abs(timelo-fOffsetLoGain - timehi)>1.0)
385 // deltatimelo = -1;
386 }
387
388 // Now store the result in the corresponding containers
389 MExtractedSignalPix &pix = (*fSignals)[pixidx];
390 MArrivalTimePix &tix = (*fArrTime)[pixidx];
391 pix.SetExtractedSignal(sumhi, deltasumhi, sumlo, deltasumlo);
392 pix.SetGainSaturation(numsathi, numsatlo);
393
394 tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
395 tix.SetGainSaturation(numsathi, numsatlo);
396 } /* while (pixel.Next()) */
397
398 fArrTime->SetReadyToSave();
399 fSignals->SetReadyToSave();
400
401 return kTRUE;
402}
403
404// --------------------------------------------------------------------------
405//
406// In addition to the resources of the base-class MExtractor:
407// MJPedestal.MExtractor.LoGainStartShift: -2.8
408//
409Int_t MExtractTimeAndCharge::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
410{
411 Bool_t rc = MExtractTime::ReadEnv(env, prefix, print);
412
413 if (IsEnvDefined(env, prefix, "LoGainStartShift", print))
414 {
415 fLoGainStartShift = GetEnvValue(env, prefix, "LoGainStartShift", fgLoGainStartShift);
416 SetLoGainStartShift(fLoGainStartShift);
417 rc = kTRUE;
418 }
419
420 if (IsEnvDefined(env, prefix, "LoGainSwitch", print))
421 {
422 fLoGainSwitch = GetEnvValue(env, prefix, "LoGainSwitch", fLoGainSwitch);
423 rc = kTRUE;
424 }
425
426 return rc;
427}
428
429void MExtractTimeAndCharge::Print(Option_t *o) const
430{
431 MExtractTime::Print(o);
432
433 if (HasLoGain())
434 {
435 *fLog << dec;
436 *fLog << " LoGainStartShift: " << fLoGainStartShift << endl;
437 *fLog << " LoGainSwitch: " << (Int_t)fLoGainSwitch << endl;
438 }
439}
Note: See TracBrowser for help on using the repository browser.