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

Last change on this file since 19894 was 19158, checked in by Daniela Dorner, 6 years ago
ignore negative values
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 // signal in pixel
162 Double_t nphot = pix.GetNumPhotons();
163 if (nphot<0)
164 continue;
165
166 // Get geometry of pixel
167 const MGeom &gpix = geom[i];
168
169 // Find the three pixels which are next to the COG
170 const Double_t dx = gpix.GetX() - hillas.GetMeanX(); // [mm]
171 const Double_t dy = gpix.GetY() - hillas.GetMeanY(); // [mm]
172
173 const Double_t dist0 = dx*dx+dy*dy;
174
175 if (dist0<dist[0])
176 {
177 dist[2] = dist[1];
178 dist[1] = dist[0];
179 dist[0] = dist0;
180
181 idx[2] = idx[1];
182 idx[1] = idx[0];
183 idx[0] = i;
184 }
185 else
186 if (dist0<dist[1])
187 {
188 dist[2] = dist[1];
189 dist[1] = dist0;
190
191 idx[2] = idx[1];
192 idx[1] = i;
193 }
194 else
195 if (dist0<dist[2])
196 {
197 dist[2] = dist0;
198 idx[2] = i;
199 }
200
201 if (!pix.IsPixelUsed())
202 continue;
203
204 // Check for requested islands
205 if (island>=0 && pix.GetIdxIsland()!=island)
206 continue;
207
208 // count used and core pixels
209 if (pix.IsPixelCore())
210 {
211 fNumCorePixels++;
212 fCoreArea += gpix.GetA();
213 }
214
215 // count used pixels
216 fNumUsedPixels++;
217 fUsedArea += gpix.GetA();
218
219 //
220 // Calculate signal contained inside ellipse
221 //
222 const Double_t dzx = hillas.GetCosDelta()*dx + hillas.GetSinDelta()*dy; // [mm]
223 const Double_t dzy = -hillas.GetSinDelta()*dx + hillas.GetCosDelta()*dy; // [mm]
224 const Double_t dz = gpix.GetT()/2;
225 const Double_t distr = (dzy*dzy+dzx*dzx)/(dzx*dzx*rl + dzy*dzy*rw);
226 if ((dzx==0 && dzy==0) || sqrt(distr)>sqrt(dist0)-dz)
227 fConcCore += nphot;
228
229 //
230 // count photons in outer rings of camera
231 //
232 if (gpix.IsInOutermostRing())
233 edgepix1 += nphot;
234 if (gpix.IsInOuterRing())
235 edgepix2 += nphot;
236
237 //
238 // Now convert nphot from absolute number of photons or phe to signal
239 // density (divide by pixel area), to find the pixel with highest signal
240 // density:
241 //
242 nphot *= geom.GetPixRatio(i);
243
244 // Look for signal density in two highest pixels:
245 if (nphot>maxpix1)
246 {
247 maxpix2 = maxpix1;
248 maxpix1 = nphot; // [1]
249 }
250 else
251 if (nphot>maxpix2)
252 maxpix2 = nphot; // [1]
253 }
254
255 fLeakage1 = edgepix1 / hillas.GetSize();
256 fLeakage2 = edgepix2 / hillas.GetSize();
257
258 // FIXME?: in case the pixel with highest signal density is an outer
259 // pixel, the value of fConc (ratio of signal in two highest pixels
260 // to SIZE) should rather be 2*fConc1, under the simplest assumption
261 // that the light density inside the outer (large) pixel is uniform.
262 fConc = (maxpix1+maxpix2)/hillas.GetSize(); // [ratio]
263 fConc1 = maxpix1/hillas.GetSize(); // [ratio]
264
265 //
266 // Concentration around COG (it is calculated here, because the
267 // distance of the pixel to COG is calculated anyhow)
268 //
269 fConcCOG = 0;
270 for (UInt_t i=0; i<TMath::Min(3U, npix); i++)
271 fConcCOG += idx[i]<0 ? 0 : evt[idx[i]].GetNumPhotons()*geom.GetPixRatio(idx[i]);
272 fConcCOG /= hillas.GetSize(); // [ratio]
273
274 // This can for example happen in case of Muon Rings
275 if (fConcCOG<0)
276 fConcCOG=0;
277
278 // Concentration of signal contained in ellipse
279 fConcCore /= hillas.GetSize();
280
281 SetReadyToSave();
282}
283
284// --------------------------------------------------------------------------
285//
286void MNewImagePar::Print(Option_t *) const
287{
288 *fLog << all;
289 *fLog << GetDescriptor() << endl;
290 *fLog << " - Leakage1 [1] = " << fLeakage1 << endl;
291 *fLog << " - Leakage2 [1] = " << fLeakage2 << endl;
292 *fLog << " - Conc [1] = " << fConc << endl;
293 *fLog << " - Conc1 [1] = " << fConc1 << endl;
294 *fLog << " - ConcCOG [1] = " << fConcCOG << endl;
295 *fLog << " - ConcCore [1] = " << fConcCore << endl;
296 *fLog << " - Num Used Pixels [#] = " << fNumUsedPixels << endl;
297 *fLog << " - Num Core Pixels [#] = " << fNumCorePixels << endl;
298 *fLog << " - Used Area [mm" << UTF8::kSquare << "] = " << fUsedArea << endl;
299 *fLog << " - Core Area [mm" << UTF8::kSquare << "] = " << fCoreArea << endl;
300}
301
302// -------------------------------------------------------------------------
303//
304// Print contents of MNewImagePar to *fLog, depending on the geometry in
305// units of deg.
306//
307void MNewImagePar::Print(const MGeomCam &geom) const
308{
309 *fLog << all;
310 *fLog << GetDescriptor() << endl;
311 *fLog << " - Leakage1 [1] = " << fLeakage1 << endl;
312 *fLog << " - Leakage2 [1] = " << fLeakage2 << endl;
313 *fLog << " - Conc [1] = " << fConc << endl;
314 *fLog << " - Conc1 [1] = " << fConc1 << endl;
315 *fLog << " - ConcCOG [1] = " << fConcCOG << endl;
316 *fLog << " - ConcCore [1] = " << fConcCore << endl;
317 *fLog << " - Num Used Pixels [#] = " << fNumUsedPixels << endl;
318 *fLog << " - Num Core Pixels [#] = " << fNumCorePixels << endl;
319 *fLog << " - Used Area [deg" << UTF8::kSquare << "] = " << fUsedArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
320 *fLog << " - Core Area [deg" << UTF8::kSquare << "] = " << fCoreArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
321}
Note: See TracBrowser for help on using the repository browser.