source: trunk/MagicSoft/Mars/mimage/MNewImagePar.cc@ 8462

Last change on this file since 8462 was 8210, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 11.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): 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 = FLT_MAX; //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 const MSignalPix &pix = evt[i];
169 if (pix.IsPixelUnmapped())
170 continue;
171
172 // Get geometry of pixel
173 const MGeomPix &gpix = geom[i];
174
175 // Find the three pixels which are next to the COG
176 const Double_t dx = gpix.GetX() - hillas.GetMeanX(); // [mm]
177 const Double_t dy = gpix.GetY() - hillas.GetMeanY(); // [mm]
178
179 const Double_t dist0 = dx*dx+dy*dy;
180
181 if (dist0<dist[0])
182 {
183 dist[2] = dist[1];
184 dist[1] = dist[0];
185 dist[0] = dist0;
186
187 idx[2] = idx[1];
188 idx[1] = idx[0];
189 idx[0] = i;
190 }
191 else
192 if (dist0<dist[1])
193 {
194 dist[2] = dist[1];
195 dist[1] = dist0;
196
197 idx[2] = idx[1];
198 idx[1] = i;
199 }
200 else
201 if (dist0<dist[2])
202 {
203 dist[2] = dist0;
204 idx[2] = i;
205 }
206
207 if (!pix.IsPixelUsed())
208 continue;
209
210 // Check for requested islands
211 if (island>=0 && pix.GetIdxIsland()!=island)
212 continue;
213
214 // count used and core pixels
215 if (pix.IsPixelCore())
216 {
217 fNumCorePixels++;
218 fCoreArea += gpix.GetA();
219 }
220
221 // count used pixels
222 fNumUsedPixels++;
223 fUsedArea += gpix.GetA();
224
225 // signal in pixel
226 Double_t nphot = pix.GetNumPhotons();
227
228 //
229 // Calculate signal contained inside ellipse
230 //
231 const Double_t dzx = hillas.GetCosDelta()*dx + hillas.GetSinDelta()*dy; // [mm]
232 const Double_t dzy = -hillas.GetSinDelta()*dx + hillas.GetCosDelta()*dy; // [mm]
233 const Double_t dz = gpix.GetD()*gpix.GetD()/4;
234 const Double_t tana = dzy*dzy/(dzx*dzx);
235 const Double_t distr = (1+tana)/(rl + tana*rw);
236 if (distr>dist0-dz || dzx==0)
237 fConcCore += nphot;
238
239 //
240 // count photons in outer rings of camera
241 //
242 if (gpix.IsInOutermostRing())
243 edgepix1 += nphot;
244 if (gpix.IsInOuterRing())
245 edgepix2 += nphot;
246
247 const Double_t ratio = geom.GetPixRatio(i);
248 if (TMath::Nint(ratio)==1) // Means this is a small (= inner pixel)
249 {
250 fInnerSize += nphot;
251
252 // Do calculation of "inner leakage" only for MAGIC-like geometry,
253 // i.e., 577 pixels, pixels of 2 different areas, inner part
254 // from pixel 0 to 396:
255 if (ismagiclike)
256 {
257 if(i > 270) // last two "rings" of inner pixels
258 {
259 edgepixin2 += nphot;
260 if(i > 330) // last "ring" of inner pixels
261 edgepixin1 += nphot;
262 }
263 }
264 }
265
266 //
267 // Now convert nphot from absolute number of photons or phe to signal
268 // density (divide by pixel area), to find the pixel with highest signal
269 // density:
270 //
271 nphot *= ratio;
272
273 // Look for signal density in two highest pixels:
274 if (nphot>maxpix1)
275 {
276 maxpix2 = maxpix1;
277 maxpix1 = nphot; // [1]
278 }
279 else
280 if (nphot>maxpix2)
281 maxpix2 = nphot; // [1]
282 }
283
284 fInnerLeakage1 = edgepixin1 / fInnerSize;
285 fInnerLeakage2 = edgepixin2 / fInnerSize;
286
287 fLeakage1 = edgepix1 / hillas.GetSize();
288 fLeakage2 = edgepix2 / hillas.GetSize();
289
290 // FIXME?: in case the pixel with highest signal density is an outer
291 // pixel, the value of fConc (ratio of signal in two highest pixels
292 // to SIZE) should rather be 2*fConc1, under the simplest assumption
293 // that the light density inside the outer (large) pixel is uniform.
294 fConc = (maxpix1+maxpix2)/hillas.GetSize(); // [ratio]
295 fConc1 = maxpix1/hillas.GetSize(); // [ratio]
296
297 //
298 // Concentration around COG (it is calculated here, because the
299 // distance of the pixel to COG is calculated anyhow)
300 //
301 fConcCOG = 0;
302 for (UInt_t i=0; i<TMath::Min(3U, npix); i++)
303 fConcCOG += idx[i]<0 ? 0 : evt[idx[i]].GetNumPhotons()*geom.GetPixRatio(idx[i]);
304 fConcCOG /= hillas.GetSize(); // [ratio]
305
306 // Concentration of signal contained in ellipse
307 fConcCore /= hillas.GetSize();
308
309 SetReadyToSave();
310}
311
312// --------------------------------------------------------------------------
313//
314void MNewImagePar::Print(Option_t *) const
315{
316 *fLog << all;
317 *fLog << GetDescriptor() << endl;
318 *fLog << " - Leakage1 [1] = " << fLeakage1 << endl;
319 *fLog << " - Leakage2 [1] = " << fLeakage2 << endl;
320 *fLog << " - InnerLeakage1 [1] = " << fInnerLeakage1 << endl;
321 *fLog << " - InnerLeakage2 [1] = " << fInnerLeakage2 << endl;
322 *fLog << " - InnerSize [phe] = " << fInnerSize << endl;
323 *fLog << " - Conc [1] = " << fConc << endl;
324 *fLog << " - Conc1 [1] = " << fConc1 << endl;
325 *fLog << " - ConcCOG [1] = " << fConcCOG << endl;
326 *fLog << " - ConcCore [1] = " << fConcCore << endl;
327 *fLog << " - Num Used Pixels [#] = " << fNumUsedPixels << endl;
328 *fLog << " - Num Core Pixels [#] = " << fNumCorePixels << endl;
329 *fLog << " - Used Area [mm^2] = " << fUsedArea << endl;
330 *fLog << " - Core Area [mm^2] = " << fCoreArea << endl;
331}
332
333// -------------------------------------------------------------------------
334//
335// Print contents of MNewImagePar to *fLog, depending on the geometry in
336// units of deg.
337//
338void MNewImagePar::Print(const MGeomCam &geom) const
339{
340 *fLog << all;
341 *fLog << GetDescriptor() << endl;
342 *fLog << " - Leakage1 [1] = " << fLeakage1 << endl;
343 *fLog << " - Leakage2 [1] = " << fLeakage2 << endl;
344 *fLog << " - InnerLeakage1 [1] = " << fInnerLeakage1 << endl;
345 *fLog << " - InnerLeakage2 [1] = " << fInnerLeakage2 << endl;
346 *fLog << " - InnerSize [phe] = " << fInnerSize << endl;
347 *fLog << " - Conc [1] = " << fConc << endl;
348 *fLog << " - Conc1 [1] = " << fConc1 << endl;
349 *fLog << " - ConcCOG [1] = " << fConcCOG << endl;
350 *fLog << " - ConcCore [1] = " << fConcCore << endl;
351 *fLog << " - Num Used Pixels [#] = " << fNumUsedPixels << endl;
352 *fLog << " - Num Core Pixels [#] = " << fNumCorePixels << endl;
353 *fLog << " - Used Area [deg^2] = " << fUsedArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
354 *fLog << " - Core Area [deg^2] = " << fCoreArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
355}
Note: See TracBrowser for help on using the repository browser.