source: tags/Mars-V0.9.5/mimage/MNewImagePar.cc

Last change on this file was 7543, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 11.6 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): Wolfgang Wittek 03/2003 <mailto:wittek@mppmu.mpg.de>
19! Author(s): Thomas Bretz <mailto:tbretz@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MNewImagePar
29//
30// Storage Container for new image parameters
31//
32// Float_t fLeakage1; // (photons in most outer ring of pixels) over fSize
33// Float_t fLeakage2; // (photons in the 2 outer rings of pixels) over fSize
34// Float_t fInnerLeakage1; // (photons in most outer rings of inner pixels) over fInnerSize
35// Float_t fInnerLeakage2; // (photons in the 2 outer rings of inner pixels) over fInnerSize
36// Float_t fInnerSize; //
37//
38// Float_t fConc; // [ratio] concentration ratio: sum of the two highest pixels / fSize
39// Float_t fConc1; // [ratio] concentration ratio: sum of the highest pixel / fSize
40// Float_t fConcCOG; // [ratio] concentration of the three pixels next to COG
41// Float_t fConcCore; // [ratio] concentration of signals inside or touching the ellipse
42//
43// Float_t fUsedArea; // Area of pixels which survived the image cleaning
44// Float_t fCoreArea; // Area of core pixels
45// Short_t fNumUsedPixels; // Number of pixels which survived the image cleaning
46// Short_t fNumCorePixels; // number of core pixels
47// Short_t fNumHGSaturatedPixels; // number of pixels with saturating hi-gains
48// Short_t fNumSaturatedPixels; // number of pixels with saturating lo-gains
49//
50// Version 2:
51// ----------
52// - added fNumSaturatedPixels
53//
54// Version 3:
55// ----------
56// - added fNumHGSaturatedPixels
57// - added fInnerLeakage1
58// - added fInnerLeakage2
59// - added fInnerSize
60// - added fUsedArea
61// - added fCoreArea
62//
63// Version 4:
64// ----------
65// - moved cleaning/island independant parameters to MImagePar:
66// + removed fNumHGSaturatedPixels
67// + removed fNumSaturatedPixels
68//
69// Version 5:
70// ----------
71// - added fConcCOG
72// - added fConcCore
73//
74//
75/////////////////////////////////////////////////////////////////////////////
76#include "MNewImagePar.h"
77
78#include "MLog.h"
79#include "MLogManip.h"
80
81#include "MHillas.h"
82
83#include "MGeomCam.h"
84#include "MGeomPix.h"
85
86#include "MSignalCam.h"
87#include "MSignalPix.h"
88
89ClassImp(MNewImagePar);
90
91using namespace std;
92
93// --------------------------------------------------------------------------
94//
95// Default constructor.
96//
97MNewImagePar::MNewImagePar(const char *name, const char *title)
98{
99 fName = name ? name : "MNewImagePar";
100 fTitle = title ? title : "New image parameters";
101
102 Reset();
103}
104
105// --------------------------------------------------------------------------
106//
107void MNewImagePar::Reset()
108{
109 fLeakage1 = -1;
110 fLeakage2 = -1;
111
112 fInnerLeakage1 = -1;
113 fInnerLeakage2 = -1;
114 fInnerSize = -1;
115
116 fConc = -1;
117 fConc1 = -1;
118 fConcCOG = -1;
119 fConcCore = -1;
120
121 fNumUsedPixels = -1;
122 fNumCorePixels = -1;
123
124 fUsedArea = -1;
125 fCoreArea = -1;
126}
127
128// --------------------------------------------------------------------------
129//
130// Calculation of new image parameters
131//
132void MNewImagePar::Calc(const MGeomCam &geom, const MSignalCam &evt,
133 const MHillas &hillas, Int_t island)
134{
135 fNumUsedPixels = 0;
136 fNumCorePixels = 0;
137
138 fUsedArea = 0;
139 fCoreArea = 0;
140
141 fInnerSize = 0;
142 fConcCore = 0;
143
144 Double_t edgepix1 = 0;
145 Double_t edgepix2 = 0;
146
147 Double_t edgepixin1 = 0;
148 Double_t edgepixin2 = 0;
149
150 Float_t maxpix1 = 0; // [#phot]
151 Float_t maxpix2 = 0; // [#phot]
152
153 const Double_t d = geom.GetMaxRadius()*geom.GetMaxRadius();
154 Double_t dist[3] = { d, d, d };
155 Int_t idx[3] = { -1, -1, -1};
156
157 const Double_t rl = 1./(hillas.GetLength()*hillas.GetLength());
158 const Double_t rw = 1./(hillas.GetWidth() *hillas.GetWidth());
159
160 const Bool_t ismagiclike =
161 geom.GetNumPixels() == 577 &&
162 geom.GetNumAreas() == 2 &&
163 geom.GetPixRatio(396) > geom.GetPixRatio(397);
164
165 UInt_t npix = evt.GetNumPixels();
166 for (UInt_t i=0; i<npix; i++)
167 {
168 // Get geometry of pixel
169 const MGeomPix &gpix = geom[i];
170
171 // Find the three pixels which are next to the COG
172 const Double_t dx = gpix.GetX() - hillas.GetMeanX(); // [mm]
173 const Double_t dy = gpix.GetY() - hillas.GetMeanY(); // [mm]
174
175 const Double_t dist0 = dx*dx+dy*dy;
176
177 if (dist0<dist[0])
178 {
179 dist[2] = dist[1];
180 dist[1] = dist[0];
181 dist[0] = dist0;
182
183 idx[2] = idx[1];
184 idx[1] = idx[0];
185 idx[0] = i;
186 }
187 else
188 if (dist0<dist[1])
189 {
190 dist[2] = dist[1];
191 dist[1] = dist0;
192
193 idx[2] = idx[1];
194 idx[1] = i;
195 }
196 else
197 if (dist0<dist[2])
198 {
199 dist[2] = dist0;
200 idx[2] = i;
201 }
202
203 const MSignalPix &pix = evt[i];
204 if (!pix.IsPixelUsed())
205 continue;
206
207 // Check for requested islands
208 if (island>=0 && pix.GetIdxIsland()!=island)
209 continue;
210
211 // count used and core pixels
212 if (pix.IsPixelCore())
213 {
214 fNumCorePixels++;
215 fCoreArea += gpix.GetA();
216 }
217
218 // count used pixels
219 fNumUsedPixels++;
220 fUsedArea += gpix.GetA();
221
222 // signal in pixel
223 Double_t nphot = pix.GetNumPhotons();
224
225 //
226 // Calculate signal contained inside ellipse
227 //
228 const Double_t dzx = hillas.GetCosDelta()*dx + hillas.GetSinDelta()*dy; // [mm]
229 const Double_t dzy = -hillas.GetSinDelta()*dx + hillas.GetCosDelta()*dy; // [mm]
230 const Double_t dz = gpix.GetD()*gpix.GetD()/4;
231 const Double_t tana = dzy*dzy/(dzx*dzx);
232 const Double_t distr = (1+tana)/(rl + tana*rw);
233 if (distr>dist0-dz || dzx==0)
234 fConcCore += nphot;
235
236 //
237 // count photons in outer rings of camera
238 //
239 if (gpix.IsInOutermostRing())
240 edgepix1 += nphot;
241 if (gpix.IsInOuterRing())
242 edgepix2 += nphot;
243
244 const Double_t ratio = geom.GetPixRatio(i);
245 if (TMath::Nint(ratio)==1) // Means this is a small (= inner pixel)
246 {
247 fInnerSize += nphot;
248
249 // Do calculation of "inner leakage" only for MAGIC-like geometry,
250 // i.e., 577 pixels, pixels of 2 different areas, inner part
251 // from pixel 0 to 396:
252 if (ismagiclike)
253 {
254 if(i > 270) // last two "rings" of inner pixels
255 {
256 edgepixin2 += nphot;
257 if(i > 330) // last "ring" of inner pixels
258 edgepixin1 += nphot;
259 }
260 }
261 }
262
263 //
264 // Now convert nphot from absolute number of photons or phe to signal
265 // density (divide by pixel area), to find the pixel with highest signal
266 // density:
267 //
268 nphot *= ratio;
269
270 // Look for signal density in two highest pixels:
271 if (nphot>maxpix1)
272 {
273 maxpix2 = maxpix1;
274 maxpix1 = nphot; // [1]
275 }
276 else
277 if (nphot>maxpix2)
278 maxpix2 = nphot; // [1]
279 }
280
281 fInnerLeakage1 = edgepixin1 / fInnerSize;
282 fInnerLeakage2 = edgepixin2 / fInnerSize;
283
284 fLeakage1 = edgepix1 / hillas.GetSize();
285 fLeakage2 = edgepix2 / hillas.GetSize();
286
287 // FIXME?: in case the pixel with highest signal density is an outer
288 // pixel, the value of fConc (ratio of signal in two highest pixels
289 // to SIZE) should rather be 2*fConc1, under the simplest assumption
290 // that the light density inside the outer (large) pixel is uniform.
291 fConc = (maxpix1+maxpix2)/hillas.GetSize(); // [ratio]
292 fConc1 = maxpix1/hillas.GetSize(); // [ratio]
293
294 //
295 // Concentration around COG (it is calculated here, because the
296 // distance of the pixel to COG is calculated anyhow)
297 //
298 fConcCOG = 0;
299 for (UInt_t i=0; i<TMath::Min(3U, npix); i++)
300 fConcCOG += idx[i]<0 ? 0 : evt[idx[i]].GetNumPhotons()*geom.GetPixRatio(idx[i]);
301 fConcCOG /= hillas.GetSize(); // [ratio]
302
303 // Concentration of signal contained in ellipse
304 fConcCore /= hillas.GetSize();
305
306 SetReadyToSave();
307}
308
309// --------------------------------------------------------------------------
310//
311void MNewImagePar::Print(Option_t *) const
312{
313 *fLog << all;
314 *fLog << GetDescriptor() << endl;
315 *fLog << " - Leakage1 [1] = " << fLeakage1 << endl;
316 *fLog << " - Leakage2 [1] = " << fLeakage2 << endl;
317 *fLog << " - InnerLeakage1 [1] = " << fInnerLeakage1 << endl;
318 *fLog << " - InnerLeakage2 [1] = " << fInnerLeakage2 << endl;
319 *fLog << " - InnerSize [phe] = " << fInnerSize << endl;
320 *fLog << " - Conc [1] = " << fConc << endl;
321 *fLog << " - Conc1 [1] = " << fConc1 << endl;
322 *fLog << " - ConcCOG [1] = " << fConcCOG << endl;
323 *fLog << " - ConcCore [1] = " << fConcCore << endl;
324 *fLog << " - Num Used Pixels [#] = " << fNumUsedPixels << endl;
325 *fLog << " - Num Core Pixels [#] = " << fNumCorePixels << endl;
326 *fLog << " - Used Area [mm^2] = " << fUsedArea << endl;
327 *fLog << " - Core Area [mm^2] = " << fCoreArea << endl;
328}
329
330// -------------------------------------------------------------------------
331//
332// Print contents of MNewImagePar to *fLog, depending on the geometry in
333// units of deg.
334//
335void MNewImagePar::Print(const MGeomCam &geom) const
336{
337 *fLog << all;
338 *fLog << GetDescriptor() << endl;
339 *fLog << " - Leakage1 [1] = " << fLeakage1 << endl;
340 *fLog << " - Leakage2 [1] = " << fLeakage2 << endl;
341 *fLog << " - InnerLeakage1 [1] = " << fInnerLeakage1 << endl;
342 *fLog << " - InnerLeakage2 [1] = " << fInnerLeakage2 << endl;
343 *fLog << " - InnerSize [phe] = " << fInnerSize << endl;
344 *fLog << " - Conc [1] = " << fConc << endl;
345 *fLog << " - Conc1 [1] = " << fConc1 << endl;
346 *fLog << " - ConcCOG [1] = " << fConcCOG << endl;
347 *fLog << " - ConcCore [1] = " << fConcCore << endl;
348 *fLog << " - Num Used Pixels [#] = " << fNumUsedPixels << endl;
349 *fLog << " - Num Core Pixels [#] = " << fNumCorePixels << endl;
350 *fLog << " - Used Area [deg^2] = " << fUsedArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
351 *fLog << " - Core Area [deg^2] = " << fCoreArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
352}
Note: See TracBrowser for help on using the repository browser.