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

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