source: trunk/Mars/msignal/MExtractFACT.cc@ 19852

Last change on this file since 19852 was 18272, checked in by Daniela Dorner, 10 years ago
fixed bug in setting values for pixel where signal cannot be extracted
File size: 7.7 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): Thomas Bretz, 2014 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2014
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MExtractFACT
28//
29//////////////////////////////////////////////////////////////////////////////
30#include "MExtractFACT.h"
31
32#include <algorithm>
33
34#include <TRandom.h>
35
36#include "MLog.h"
37#include "MLogManip.h"
38
39#include "MParList.h"
40
41#include "MRawRunHeader.h"
42
43#include "MArrivalTimeCam.h"
44#include "MArrivalTimePix.h"
45
46#include "MExtractedSignalCam.h"
47#include "MExtractedSignalPix.h"
48
49#include "MPedestalSubtractedEvt.h"
50
51ClassImp(MExtractFACT);
52
53using namespace std;
54
55// --------------------------------------------------------------------------
56//
57// Default constructor.
58//
59// Sets:
60// - fWindowSizeHiGain and fWindowSizeLoGain to 0
61// - fLoGainStartShift to fgLoGainStartShift
62// - fLoGainSwitch to fgLoGainSwitch
63//
64MExtractFACT::MExtractFACT(const char *name, const char *title)
65{
66 fName = name ? name : "MExtractFACT";
67 fTitle = title ? title : "Base class for signal and time extractors";
68}
69
70// --------------------------------------------------------------------------
71//
72// The PreProcess searches for the following input containers:
73// - MRawEvtData
74// - MRawRunHeader
75// - MPedestalCam
76//
77// The following output containers are also searched and created if
78// they were not found:
79//
80// - MExtractedSignalCam
81// - MArrivalTimeCam
82//
83Int_t MExtractFACT::PreProcess(MParList *pList)
84{
85 if (!MExtractTime::PreProcess(pList))
86 return kFALSE;
87
88 fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam));
89 if (!fSignals)
90 return kFALSE;
91
92 *fLog << flush << inf;
93 return kTRUE;
94}
95
96// --------------------------------------------------------------------------
97//
98// The ReInit calls:
99// - MExtractor::ReInit()
100//
101// Call:
102// - MArrivalTimeCam::SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
103// fLoGainFirst, fLoGainLast, fNumLoGainSamples);
104// - InitArrays();
105//
106Bool_t MExtractFACT::ReInit(MParList *pList)
107{
108 if (!MExtractTime::ReInit(pList))
109 return kFALSE;
110
111 //if (fSignals)
112 // fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
113 // fLoGainFirst, fLoGainLast, fNumLoGainSamples);
114
115 return kTRUE;
116}
117
118// --------------------------------------------------------------------------
119//
120// Calculate the integral of the FADC time slices and store them as a new
121// pixel in the MArrivalTimeCam container.
122// Calculate the integral of the FADC time slices and store them as a new
123// pixel in the MExtractedSignalCam container.
124// The functions FindTimeAndChargeHiGain() and FindTimeAndChargeLoGain() are
125// supposed to extract the signal themselves.
126//
127Int_t MExtractFACT::Process()
128{
129 // FIXME: Check ranges
130 //const UInt_t nroi = fSignal->GetNumSamples();
131 const UInt_t npix = fSignal->GetNumPixels();
132
133 const UInt_t rangehi = fHiGainLast - fHiGainFirst - 1;
134
135 for (UInt_t ch=0; ch<npix; ch++)
136 {
137 const Float_t *ptr = fSignal->GetSamples(ch);
138
139 const Float_t *pbeg = ptr+fHiGainFirst+1;
140 const Float_t *pend = ptr+fHiGainLast-1;
141
142 // Find maximum
143 const Float_t *pmax = pbeg;
144
145 if (!IsNoiseCalculation())
146 {
147 for (const Float_t *p=pbeg; p<pend; p++)
148 if (*p>*pmax)
149 pmax = p;
150 }
151
152 //Float_t max = *pmax;
153
154 // -------------------------------------------------------------------------
155 // Calculate a parabola through this and the surrounding points
156 // on logarithmic values (that's a gaussian)
157 //
158 // Quadratic interpolation
159 //
160 // calculate the parameters of a parabula such that
161 // y(i) = a + b*x(i) + c*x(i)^2
162 // y'(i) = b + 2*c*x(i)
163 //
164 // y = exp( - (x-k)^2 / s^2 / 2 )
165 //
166 // -2*s^2 * log(y) = x^2 - 2*k*x + k^2
167 //
168 // a = (k/s0)^2/2
169 // b = k/s^2
170 // c = -1/(2s^2) --> s = sqrt(-1/2c)
171 //
172 Float_t max = -1;
173 Float_t tmax = -1;
174 if (pmax>pbeg && pmax<pend-1)
175 {
176 const Float_t &v1 = pmax[-1];
177 const Float_t &v2 = pmax[ 0];
178 const Float_t &v3 = pmax[+1];
179 if (v1>0 && v2>0 && v3>0)
180 {
181 const double y1 = log(v1);
182 const double y2 = log(v2);
183 const double y3 = log(v3);
184
185 const double a = y2;
186 const double b = (y3-y1)/2;
187 const double c = (y1+y3)/2 - y2;
188 if (c<0)
189 {
190 const double w = -1./(2*c);
191 const double dx = b*w;
192
193 if (dx>=-1 && dx<=1)
194 {
195 max = exp(a + b*dx + c*dx*dx);
196
197 // tmax is position of maximum
198 tmax = dx;
199 tmax += Int_t(pmax-ptr);
200 }
201 }
202 }
203 }
204
205 Float_t time = -1;
206 Float_t slope = -1;
207 if (max>=0)
208 {
209 // -10: maximum search window
210
211 // Time is position of half hight leading edge
212 pend = std::max(pbeg, pmax-15);
213 for (;pmax>=pend; pmax--)
214 if (*pmax<max/2)
215 break;
216
217 if (pmax>pend && pmax[0]!=pmax[1])
218 {
219 time = (max/2-pmax[0])/(pmax[1]-pmax[0]);
220 time += Int_t(pmax-ptr);
221 slope = tmax - time;
222 }
223 }
224
225 if (time<0)
226 {
227 time = gRandom->Uniform(rangehi)+fHiGainFirst+1;
228 max = ptr[uint16_t(time)];
229 }
230
231 // Now store the result in the corresponding containers
232 MExtractedSignalPix &pix = (*fSignals)[ch];
233 MArrivalTimePix &tix = (*fArrTime)[ch];
234 pix.SetExtractedSignal(max);
235 pix.SetGainSaturation(0);
236
237 tix.SetArrivalTime(time, slope);
238 tix.SetGainSaturation(0);
239 }
240
241 fArrTime->SetReadyToSave();
242 fSignals->SetReadyToSave();
243
244 return kTRUE;
245}
246
247// --------------------------------------------------------------------------
248//
249// In addition to the resources of the base-class MExtractor:
250// MJPedestal.MExtractor.LoGainStartShift: -2.8
251//
252Int_t MExtractFACT::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
253{
254 Bool_t rc = MExtractTime::ReadEnv(env, prefix, print);
255 /*
256 if (IsEnvDefined(env, prefix, "LoGainStartShift", print))
257 {
258 fLoGainStartShift = GetEnvValue(env, prefix, "LoGainStartShift", fgLoGainStartShift);
259 SetLoGainStartShift(fLoGainStartShift);
260 rc = kTRUE;
261 }
262
263 if (IsEnvDefined(env, prefix, "LoGainSwitch", print))
264 {
265 fLoGainSwitch = GetEnvValue(env, prefix, "LoGainSwitch", (Int_t)fLoGainSwitch);
266 rc = kTRUE;
267 }
268 */
269 return rc;
270}
271
Note: See TracBrowser for help on using the repository browser.