source: trunk/Mars/mimage/MNewImagePar.cc@ 14623

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