source: trunk/MagicSoft/Mars/msignal/MArrivalTimeCalc2.cc@ 9620

Last change on this file since 9620 was 7810, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 7.3 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, 02/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Hendrik Bartko, 02/2004 <mailto:hbartko@mppmu.mpg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MArrivalTimeCalc2
29//
30// Calculates the arrival time as the mean time of the fWindowSize time slices
31// which have the highest integral content.
32//
33//
34//////////////////////////////////////////////////////////////////////////////
35#include "MArrivalTimeCalc2.h"
36
37#include "MLog.h"
38#include "MLogManip.h"
39
40#include "MParList.h"
41
42#include "MRawEvtData.h"
43#include "MRawEvtPixelIter.h"
44
45#include "MPedestalCam.h"
46#include "MPedestalPix.h"
47
48#include "MArrivalTimeCam.h"
49#include "MArrivalTimePix.h"
50
51
52ClassImp(MArrivalTimeCalc2);
53
54using namespace std;
55
56const Byte_t MArrivalTimeCalc2::fgSaturationLimit = 254;
57const Byte_t MArrivalTimeCalc2::fgFirst = 0;
58const Byte_t MArrivalTimeCalc2::fgLast = 14;
59const Byte_t MArrivalTimeCalc2::fgWindowSize = 6;
60
61// --------------------------------------------------------------------------
62//
63// Default constructor.
64//
65MArrivalTimeCalc2::MArrivalTimeCalc2(const char *name, const char *title)
66 : fSaturationLimit(fgSaturationLimit)
67{
68
69 fName = name ? name : "MArrivalTimeCalc2";
70 fTitle = title ? title : "Task to extract the signal from the FADC slices";
71
72 AddToBranchList("MRawEvtData.*");
73
74 SetRange();
75}
76
77void MArrivalTimeCalc2::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast, Byte_t windowsize)
78{
79 fNumHiGainSamples = hilast-hifirst+1;
80 fNumLoGainSamples = lolast-lofirst+1;
81
82 fHiGainFirst = hifirst;
83 fLoGainFirst = lofirst;
84
85 fWindowSize = windowsize & ~1;
86
87 if (fWindowSize != windowsize)
88 *fLog << warn << "MArrivalTimeCalc2::SetRange - window size has to be even, set to: " << int(fWindowSize) << " samples " << endl;
89
90 if (fWindowSize<2)
91 {
92 fWindowSize = 2;
93 *fLog << warn << "MArrivalTimeCalc2::SetRange - window size set to two samples" << endl;
94 }
95
96 if (fWindowSize > fNumHiGainSamples)
97 {
98 fWindowSize = fNumLoGainSamples & ~1;
99 *fLog << warn << "MArrivalTimeCalc2::SetRange - window size set to " << int(fWindowSize) << " samples " << endl;
100 }
101
102 if (fWindowSize > fNumLoGainSamples)
103 {
104 fWindowSize = fNumLoGainSamples & ~1;
105 *fLog << warn << "MArrivalTimeCalc2::SetRange - window size set to " << int(fWindowSize) << " samples " << endl;
106 }
107
108 fWindowSizeSqrt = TMath::Sqrt((Float_t)fWindowSize);
109}
110
111// --------------------------------------------------------------------------
112//
113// The PreProcess searches for the following input containers:
114// - MRawEvtData
115// - MPedestalCam
116//
117// The following output containers are also searched and created if
118// they were not found:
119//
120// - MExtractedSignalCam
121//
122Int_t MArrivalTimeCalc2::PreProcess(MParList *pList)
123{
124 fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData"));
125 if (!fRawEvt)
126 {
127 *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl;
128 return kFALSE;
129 }
130
131 fPedestals = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
132 if (!fPedestals)
133 {
134 *fLog << err << AddSerialNumber("MPedestalCam") << " not found... aborting" << endl;
135 return kFALSE;
136 }
137
138 fArrivalTime = (MArrivalTimeCam*)pList->FindCreateObj(AddSerialNumber("MArrivalTimeCam"));
139 if (!fArrivalTime)
140 return kFALSE;
141
142 return kTRUE;
143}
144
145void MArrivalTimeCalc2::FindSignalTime(Byte_t *ptr, Byte_t size, Float_t &time, Float_t &deltatime, Int_t &sat, Float_t pedes, Float_t pedrms) const
146{
147 const Byte_t *end = ptr + size;
148
149 Int_t sum=0; // integral content of the actual window
150 Int_t max = 0; // highest integral content of all windows
151
152 //
153 // Calculate the sum of the first fWindowSize slices
154 //
155 sat = 0;
156 Byte_t *p = ptr;
157
158 while (p<ptr+fWindowSize)
159 {
160 sum += *p;
161 if (*p++ >= fSaturationLimit)
162 sat++;
163 }
164
165 //
166 // Check for saturation in all other slices
167 //
168 while (p<end)
169 if (*p++ >= fSaturationLimit)
170 sat++;
171
172 //
173 // Calculate the i-th sum as
174 // sum_i+1 = sum_i + slice[i+8] - slice[i]
175 // This is fast and accurate (because we are using int's)
176 //
177 max=sum;
178 Byte_t *ptrmax=ptr; // pointer to the first slice of the maximum window
179
180 for (p=ptr; p+fWindowSize<end; p++)
181 {
182 sum += *(p+fWindowSize) - *p;
183
184 if (sum>max)
185 {
186 max = sum;
187 ptrmax = p+1;
188 }
189 }
190
191 // now calculate the time for the maximum window
192 Int_t timesignalsum = 0;
193 Int_t timesquaredsum =0;
194 Int_t timesum =0;
195
196 for (p=ptrmax; p < ptrmax + fWindowSize; p++)
197 {
198 timesignalsum += *p*(p-ptr);
199 timesum += p-ptr;
200 timesquaredsum += (p-ptr)*(p-ptr);
201 }
202
203 const Float_t pedsubsum = max - fWindowSize*pedes;
204 const Float_t pedsubtimesignalsum = timesignalsum - timesum*pedes;
205
206 time = pedsubsum != 0 ? pedsubtimesignalsum / pedsubsum : 1;
207 deltatime = pedsubsum != 0 ? pedrms / pedsubsum * sqrt(timesquaredsum - fWindowSize*time) : 1;
208}
209
210// --------------------------------------------------------------------------
211//
212// Calculates the arrival time as the mean time of the fWindowSize time slices
213// which have the highest integral content.
214//
215Int_t MArrivalTimeCalc2::Process()
216{
217 MRawEvtPixelIter pixel(fRawEvt);
218
219 Int_t sat=0;
220 while (pixel.Next())
221 {
222 //
223 // Find signal in hi- and lo-gain
224 //
225 Float_t timehi, timelo, deltatimehi, deltatimelo;
226 Int_t sathi, satlo;
227
228 //
229 // Take correspodning pedestal
230 //
231 const Int_t pixid = pixel.GetPixelId();
232
233 const MPedestalPix &ped = (*fPedestals)[pixid];
234
235 MArrivalTimePix &pix = (*fArrivalTime)[pixid];
236
237 const Float_t pedes = ped.GetPedestal();
238 const Float_t pedrms = ped.GetPedestalRms();
239
240 FindSignalTime(pixel.GetHiGainSamples()+fHiGainFirst, fNumHiGainSamples, timehi, deltatimehi, sathi, pedes, pedrms);
241 FindSignalTime(pixel.GetLoGainSamples()+fLoGainFirst, fNumLoGainSamples, timelo, deltatimelo, satlo, pedes, pedrms);
242
243 if (satlo)
244 sat++;
245
246 pix.SetArrivalTime(timehi+ Float_t(fHiGainFirst), deltatimehi, timelo + Float_t(fLoGainFirst), deltatimelo);
247 pix.SetGainSaturation(sathi, satlo);
248 }
249
250 fArrivalTime->SetReadyToSave();
251
252 //
253 // Print a warning if event has saturationg lo-gains
254 //
255 // if (sat)
256 // *fLog << warn << "WARNING - Lo Gain saturated in " << sat << " pixels." << endl;
257
258 return kTRUE;
259}
Note: See TracBrowser for help on using the repository browser.