source: trunk/Mars/msignal/MExtractTimeAndCharge.cc@ 10113

Last change on this file since 10113 was 9069, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 16.1 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.70 2008-08-02 11:11:21 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// Class Version 4:
62// ----------------
63// - Byte_t fLoGainSwitch
64// + UInt_t fLoGainSwitch
65//
66//
67// Input Containers:
68// MRawEvtData
69// MRawRunHeader
70// MPedestalCam
71//
72// Output Containers:
73// MArrivalTimeCam
74// MExtractedSignalCam
75//
76// For weired events check: Run 94127 Event 672, 1028
77//
78//////////////////////////////////////////////////////////////////////////////
79#include "MExtractTimeAndCharge.h"
80
81#include <TRandom.h>
82#include <TVector3.h>
83
84#include "MMath.h"
85
86#include "MLog.h"
87#include "MLogManip.h"
88
89#include "MParList.h"
90
91#include "MRawRunHeader.h"
92#include "MRawEvtPixelIter.h"
93
94#include "MArrivalTimeCam.h"
95#include "MArrivalTimePix.h"
96
97#include "MExtractedSignalCam.h"
98#include "MExtractedSignalPix.h"
99
100#include "MPedestalSubtractedEvt.h"
101
102ClassImp(MExtractTimeAndCharge);
103
104using namespace std;
105
106const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -1.0; // was -2.5
107const UInt_t MExtractTimeAndCharge::fgLoGainSwitch = 120;
108
109// --------------------------------------------------------------------------
110//
111// Default constructor.
112//
113// Sets:
114// - fWindowSizeHiGain and fWindowSizeLoGain to 0
115// - fLoGainStartShift to fgLoGainStartShift
116// - fLoGainSwitch to fgLoGainSwitch
117//
118MExtractTimeAndCharge::MExtractTimeAndCharge(const char *name, const char *title)
119 : fWindowSizeHiGain(0), fWindowSizeLoGain(0)
120{
121 fName = name ? name : "MExtractTimeAndCharge";
122 fTitle = title ? title : "Base class for signal and time extractors";
123
124 SetLoGainStartShift();
125 SetLoGainSwitch();
126}
127
128// --------------------------------------------------------------------------
129//
130// The PreProcess searches for the following input containers:
131// - MRawEvtData
132// - MRawRunHeader
133// - MPedestalCam
134//
135// The following output containers are also searched and created if
136// they were not found:
137//
138// - MExtractedSignalCam
139// - MArrivalTimeCam
140//
141Int_t MExtractTimeAndCharge::PreProcess(MParList *pList)
142{
143 if (!MExtractTime::PreProcess(pList))
144 return kFALSE;
145
146 fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam));
147 if (!fSignals)
148 return kFALSE;
149
150 *fLog << flush << inf;
151 return kTRUE;
152}
153
154// --------------------------------------------------------------------------
155//
156// The ReInit calls:
157// - MExtractor::ReInit()
158//
159// Call:
160// - MArrivalTimeCam::SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
161// fLoGainFirst, fLoGainLast, fNumLoGainSamples);
162// - InitArrays();
163//
164Bool_t MExtractTimeAndCharge::ReInit(MParList *pList)
165{
166 if (!MExtractTime::ReInit(pList))
167 return kFALSE;
168
169 if (!InitArrays(fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain()))
170 return kFALSE;
171
172 if (fSignals)
173 fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
174 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
175
176 if (!HasLoGain())
177 {
178 *fLog << inf << "No lo-gains... resetting lo-gain switch." << endl;
179 fLoGainSwitch=0xff;
180 }
181
182 return kTRUE;
183}
184
185// --------------------------------------------------------------------------
186//
187// Return the x-value lower than sat0 at which the signal has been
188// fallen bwlow maxcont/2. This time is determined using a simple second
189// order polynomial interpolation.
190//
191Double_t MExtractTimeAndCharge::GetSaturationTime(Int_t sat0, const Float_t *sig, Int_t maxconthalf) const
192{
193 const Int_t p = sat0>1 ? sat0-2 : sat0-1;
194 if (sat0<=0)
195 return 0;
196
197 if (sat0==1)
198 return sig[0]>maxconthalf ? 0 : 0.5;
199
200 if (sig[p]>sig[p+1] || sig[p+1]>sig[p+2])
201 return sig[p+1]>maxconthalf ? sat0-1 : sat0-0.5;
202
203 // Find the place at which the signal is maxcont/2
204 const TVector3 vx(sig[p], sig[p+1], sig[p+2]);
205 const TVector3 vy(p, p+1, p+2);
206
207 return MMath::InterpolParabLin(vx, vy, maxconthalf);
208}
209
210// --------------------------------------------------------------------------
211//
212// Calculate the integral of the FADC time slices and store them as a new
213// pixel in the MArrivalTimeCam container.
214// Calculate the integral of the FADC time slices and store them as a new
215// pixel in the MExtractedSignalCam container.
216// The functions FindTimeAndChargeHiGain() and FindTimeAndChargeLoGain() are
217// supposed to extract the signal themselves.
218//
219Int_t MExtractTimeAndCharge::Process()
220{
221 const Int_t nums = fSignal->GetNumSamples();
222
223 const Int_t numh = fRunHeader->GetNumSamplesHiGain();
224 const Int_t numl = fRunHeader->GetNumSamplesLoGain();
225
226 // Some sanity checks to get rid of some weird behaviour of the simulation
227 if (nums!=numh+numl)
228 {
229 *fLog << err << "MExtractTimeAndCharge::Process: ERROR - Number of samples in event ";
230 *fLog << "(" << nums << ")" << endl << " doesn't match number in run-header ";
231 *fLog << "(" << numh+numl << ")" << endl;
232 *fLog << " In case you are processing real data you have just found a bug." << endl;
233 *fLog << " If you process Monte Carlo data, it means probably that the length" << endl;
234 *fLog << " of the arrays in MRawEvtData are inconsistent with the run-header." << endl;
235 return kERROR;
236 }
237
238 // Start extraction
239 const UInt_t satlim = fSaturationLimit*fRunHeader->GetScale();
240
241 MRawEvtPixelIter pixel(fRawEvt);
242 while (pixel.Next())
243 {
244 const Int_t pixidx = pixel.GetPixelId();
245
246 const Float_t *sig = fSignal->GetSamples(pixidx);
247
248 // Would it be better to take lastsat-firstsat?
249 Int_t sathi0 = fHiGainFirst; // First slice to extract and first saturating slice
250 Int_t sathi1 = fHiGainLast; // Last slice to extract and last saturating slice
251 Int_t numsathi = fSignal->GetSaturation(pixidx, satlim, sathi0, sathi1);
252
253 Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid
254 Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid
255
256 // Do not even try to extract the hi-gain if we have
257 // more than one saturating slice
258 const Int_t rangehi = fHiGainLast - fHiGainFirst + 1;
259
260 if (numsathi<2)
261 {
262 const Int_t maxposhi = fSignal->GetMaxPos(pixidx, fHiGainFirst, fHiGainLast);
263 FindTimeAndChargeHiGain2(sig+fHiGainFirst, rangehi,
264 sumhi, deltasumhi, timehi, deltatimehi,
265 numsathi, maxposhi);
266 }
267
268 // If we have saturating slices try to get a better estimate
269 // of the arrival time than timehi or sathi0. This is
270 // usefull to know where to start lo-gain extraction.
271 const UInt_t maxcont = fSignal->GetRawMaxVal(pixidx, fHiGainFirst, fHiGainLast);
272 if (numsathi>1)
273 {
274 timehi = GetSaturationTime(sathi0, sig, maxcont/2)-fHiGainFirst;
275 deltatimehi = 0;
276 }
277
278 // Make sure that in cases the time couldn't be correctly determined
279 // more meaningfull default values are assigned.
280 // For extractors like the digital filter and the spline
281 // we allow extracpolation by one slice.
282 // FIXME: Defined Out-Of-Range better so that the extractors
283 // know what to return!
284 if (deltatimehi>-0.5 && (timehi<-1 || timehi>=rangehi))
285 {
286 // Flag this as unreliable!
287 timehi = gRandom->Uniform(rangehi+1)-1;
288 // deltatimehi=-1; // Set PIXEL to UNRELIABLE?
289 }
290
291 timehi += fHiGainFirst;
292
293 Float_t sumlo =0, deltasumlo =-1; // invalidate logain of MExtractedSignalPix
294 Float_t timelo=0, deltatimelo=-1; // invalidate logain of MArrivalTimePix
295 Int_t numsatlo=0;
296
297 //
298 // Adapt the low-gain extraction range from the obtained high-gain time
299 //
300
301 // IN THIS CASE THE PIXEL SHOULD BE MARKED BAD!!!!
302 // MEANS: Hi gain has saturated, but the signal is to dim
303 // to extract the lo-gain properly
304 // This could happen because the maxcont was not extracted from
305 // all slices!
306 // THIS produces pulse positions ~= -1
307 // The signal might be handled in MCalibrateData, but hwat's about
308 // the arrival times in MCalibrateRelTime
309 if (numsathi>0 && maxcont<=fLoGainSwitch)
310 deltasumlo=deltasumhi=deltatimelo=deltatimehi=-1;
311
312 // If more than 8 hi-gain slices have saturated this is
313 // no physical event. We just assume that something with
314 // the extraction is wrong
315 if (numsathi>8) // FIXME: Should be something like 2?
316 deltasumhi=deltatimehi=-1;
317
318 // FIXME: What to do with the pixel if it saturates too early???
319 if (numl>0 && maxcont>fLoGainSwitch)
320 {
321 Int_t first = numh+fLoGainFirst;
322 Int_t last = numh+fLoGainLast;
323
324 // To determin the window in which the lo-gain is extracted
325 // clearly more information about the relation between the
326 // extraction window and the reslting time is necessary.
327 //
328 // numh + fLoGainStartShift == 14 / fLoGainStartShift=-1
329 //
330 // The lo-gain is expected to have its raising edge
331 // at timehi+numh+fOffsetLoGain. We start to search for the
332 // lo-gain fLoGainStartShift slices earlier.
333 //
334 // Instead of fLoGainStartShift the extractor should now how many
335 // slices before the expected raising edge the start of the
336 // search must be placed and from there we can step 1.5 slices
337 // back to be on the safe side.
338 //
339 // The jitter in the hi-/lo-gain offset ssems to be around +/-0.5
340 const Float_t tm = deltatimehi<0 ? -1.+fHiGainFirst : timehi;
341 first = TMath::FloorNint(tm+numh+fOffsetLoGain+fLoGainStartShift);
342
343 if (first<0)
344 first = 0;
345 if (first>last)
346 first=last;
347
348 /*
349 // Currently we have to stick to this check because at least
350 // the spline has arrays of this size...
351 if (first>last)
352 first = last;
353 if (first<numh+fLoGainFirst)
354 first = numh+fLoGainFirst;
355 */
356 // Would it be better to take lastsat-firstsat?
357 Int_t satlo0 = first; // First slice to extract and first saturating slice
358 Int_t satlo1 = last; // Last slice to extract and last saturating slice
359 numsatlo = fSignal->GetSaturation(pixidx, satlim, satlo0, satlo1);
360
361 //if (satlo0>first && satlo1<last && numsatlo>2)
362 //{
363 // fSignal->InterpolateSaturation(pixidx, fSaturationLimit, satlo0, satlo1);
364 // numsatlo = 0;
365 //}
366
367 const Int_t rangelo = last-first+1;
368 const Int_t maxposlo = fSignal->GetMaxPos(pixidx, first, last);
369 FindTimeAndChargeLoGain2(sig+first, rangelo,
370 sumlo, deltasumlo, timelo, deltatimelo,
371 numsatlo, maxposlo);
372
373 // If we have saturating slices try to get a better estimate
374 // of the arrival time than timehi or sathi0. This is
375 // usefull to know where to start lo-gain extraction.
376 if (numsatlo>1)
377 {
378 const UInt_t maxcontlo = fSignal->GetRawMaxVal(pixidx, fHiGainFirst, fHiGainLast);
379 timelo = GetSaturationTime(satlo0, sig, maxcontlo/2)-numh-first;
380 deltatimelo = 0;
381 }
382
383 // Make sure that in cases the time couldn't be correctly determined
384 // more meaningfull default values are assigned
385 // For extractors like the digital filter and the spline
386 // we allow extracpolation by one slice.
387 if (deltatimelo>-0.5 && (timelo<-1 || timelo>=rangelo))
388 {
389 // Flag this as unreliable!
390 timelo = gRandom->Uniform(rangelo+1)-1;
391 //deltatimelo=-1; // Set PIXEL to UNRELIABLE?
392 }
393
394 timelo += first-numh;
395
396 // If more than 6 lo-gain slices have saturated this is
397 // no physical event. We just assume that something with
398 // the extraction is wrong
399 if (numsatlo>6)
400 deltasumlo=deltatimelo=-1;
401
402 // The extracted lo-gain signal cannot be zero or
403 // negative at all, so it must be wrong
404 if (sumlo<=0)
405 deltasumlo=-1;
406
407 //if (TMath::Abs(timelo-fOffsetLoGain - timehi)>1.0)
408 // deltatimelo = -1;
409 }
410
411 // Now store the result in the corresponding containers
412 MExtractedSignalPix &pix = (*fSignals)[pixidx];
413 MArrivalTimePix &tix = (*fArrTime)[pixidx];
414 pix.SetExtractedSignal(sumhi, deltasumhi, sumlo, deltasumlo);
415 pix.SetGainSaturation(numsathi, numsatlo);
416
417 tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
418 tix.SetGainSaturation(numsathi, numsatlo);
419 } /* while (pixel.Next()) */
420
421 fArrTime->SetReadyToSave();
422 fSignals->SetReadyToSave();
423
424 return kTRUE;
425}
426
427// --------------------------------------------------------------------------
428//
429// In addition to the resources of the base-class MExtractor:
430// MJPedestal.MExtractor.LoGainStartShift: -2.8
431//
432Int_t MExtractTimeAndCharge::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
433{
434 Bool_t rc = MExtractTime::ReadEnv(env, prefix, print);
435
436 if (IsEnvDefined(env, prefix, "LoGainStartShift", print))
437 {
438 fLoGainStartShift = GetEnvValue(env, prefix, "LoGainStartShift", fgLoGainStartShift);
439 SetLoGainStartShift(fLoGainStartShift);
440 rc = kTRUE;
441 }
442
443 if (IsEnvDefined(env, prefix, "LoGainSwitch", print))
444 {
445 fLoGainSwitch = GetEnvValue(env, prefix, "LoGainSwitch", (Int_t)fLoGainSwitch);
446 rc = kTRUE;
447 }
448
449 return rc;
450}
451
452void MExtractTimeAndCharge::Print(Option_t *o) const
453{
454 MExtractTime::Print(o);
455
456 if (HasLoGain())
457 {
458 *fLog << dec;
459 *fLog << " LoGainStartShift: " << fLoGainStartShift << endl;
460 *fLog << " LoGainSwitch: " << fLoGainSwitch << endl;
461 }
462}
Note: See TracBrowser for help on using the repository browser.