source: trunk/MagicSoft/Mars/mimage/MHillas.cc@ 3464

Last change on this file since 3464 was 3183, checked in by tbretz, 21 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): Thomas Bretz 12/2000 <mailto:tbretz@uni-sw.gwdg.de>
19! Author(s): Harald Kornmayer 1/2001
20! Author(s): Rudolf Bock 10/2001 <mailto:Rudolf.Bock@cern.ch>
21! Author(s): Wolfgang Wittek 6/2002 <mailto:wittek@mppmu.mpg.de>
22!
23! Copyright: MAGIC Software Development, 2000-2003
24!
25!
26\* ======================================================================== */
27
28/////////////////////////////////////////////////////////////////////////////
29//
30// MHillas
31//
32// Storage Container for image parameters
33//
34// basic image parameters
35//
36// Version 1:
37// ----------
38// fLength [mm] major axis of ellipse
39// fWidth [mm] minor axis
40// fDelta [rad] angle of major axis with x-axis
41// by definition the major axis is pointing into
42// the hemisphere x>0, thus -pi/2 < delta < pi/2
43// fSize [#CerPhot] total sum of pixels
44// fMeanX [mm] x of center of ellipse
45// fMeanY [mm] y of center of ellipse
46//
47// Version 2:
48// ----------
49// fNumCorePixels number of pixels called core
50// fNumUsedPixels number of pixels which survived the cleaning
51//
52// Version 3:
53// ----------
54// fNumCorePixels moved to MNewImagePar
55// fNumUsedPixels moved to MNewImagePar
56// fCosDelta added
57// fSinDelte added
58//
59/////////////////////////////////////////////////////////////////////////////
60#include "MHillas.h"
61
62#include <fstream>
63
64#include <TArrayF.h>
65#include <TEllipse.h>
66
67#include "MGeomPix.h"
68#include "MGeomCam.h"
69
70#include "MCerPhotPix.h"
71#include "MCerPhotEvt.h"
72
73#include "MLog.h"
74#include "MLogManip.h"
75
76ClassImp(MHillas);
77
78using namespace std;
79
80// --------------------------------------------------------------------------
81//
82// Default constructor.
83//
84MHillas::MHillas(const char *name, const char *title)
85{
86 fName = name ? name : "MHillas";
87 fTitle = title ? title : "Storage container for image parameters of one event";
88
89 Reset();
90}
91
92// --------------------------------------------------------------------------
93//
94// Destructor. Deletes the TEllipse if one exists.
95//
96MHillas::~MHillas()
97{
98}
99
100// --------------------------------------------------------------------------
101//
102// Initializes the values with defaults. For the default values see the
103// source code.
104//
105void MHillas::Reset()
106{
107 fLength = -1;
108 fWidth = -1;
109 fDelta = 0;
110
111 fSize = -1;
112 fMeanX = 0;
113 fMeanY = 0;
114}
115
116// --------------------------------------------------------------------------
117//
118// Print the hillas Parameters to *fLog
119//
120void MHillas::Print(Option_t *) const
121{
122 Double_t atg = atan2(fMeanY, fMeanX)*kRad2Deg;
123
124 if (atg<0)
125 atg += 180;
126
127 *fLog << all;
128 *fLog << "Basic Image Parameters (" << GetName() << ")" << endl;
129 *fLog << " - Length [mm] = " << fLength << endl;
130 *fLog << " - Width [mm] = " << fWidth << endl;
131 *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl;
132 *fLog << " - Size [1] = " << fSize << " #CherPhot" << endl;
133 *fLog << " - Meanx [mm] = " << fMeanX << endl;
134 *fLog << " - Meany [mm] = " << fMeanY << endl;
135 *fLog << " - atg(y/x) [deg] = " << atg << endl;
136}
137
138// --------------------------------------------------------------------------
139//
140// Print the hillas Parameters to *fLog, depending on the geometry in
141// units of deg
142//
143void MHillas::Print(const MGeomCam &geom) const
144{
145 Double_t atg = atan2(fMeanY, fMeanX)*kRad2Deg;
146
147 if (atg<0)
148 atg += 180;
149
150 *fLog << all;
151 *fLog << "Basic Image Parameters (" << GetName() << ")" << endl;
152 *fLog << " - Length [deg] = " << fLength*geom.GetConvMm2Deg() << endl;
153 *fLog << " - Width [deg] = " << fWidth *geom.GetConvMm2Deg() << endl;
154 *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl;
155 *fLog << " - Size [1] = " << fSize << " #CherPhot" << endl;
156 *fLog << " - Meanx [deg] = " << fMeanX *geom.GetConvMm2Deg() << endl;
157 *fLog << " - Meany [deg] = " << fMeanY *geom.GetConvMm2Deg() << endl;
158 *fLog << " - atg(y/x) [deg] = " << atg << endl;
159}
160
161// --------------------------------------------------------------------------
162//
163// Paint the ellipse corresponding to the parameters
164//
165void MHillas::Paint(Option_t *opt)
166{
167 if (fLength<0 || fWidth<0)
168 return;
169
170 TEllipse e(fMeanX, fMeanY, fLength, fWidth, 0, 360, fDelta*kRad2Deg+180);
171 e.SetLineWidth(2);
172 e.Paint();
173}
174
175// --------------------------------------------------------------------------
176//
177// Calculate the image parameters from a Cherenkov photon event
178// assuming Cher.photons/pixel and pixel coordinates are given
179// In case you don't call Calc from within an eventloop make sure, that
180// you call the Reset member function before.
181// Returns:
182// 0 no error
183// 1 number of pixels < 3
184// 2 size==0
185// 3 number of used pixel < 3
186// 4 CorrXY == 0
187//
188Int_t MHillas::Calc(const MGeomCam &geom, const MCerPhotEvt &evt)
189{
190 //
191 // sanity check 1
192 //
193 if (evt.GetNumPixels()<3)
194 return 1;
195
196 //
197 // calculate mean value of pixel coordinates and fSize
198 // -----------------------------------------------------
199 //
200 // Because this are only simple sums of roughly 600 values
201 // with an accuracy less than three digits there should not
202 // be any need to use double precision (Double_t) for the
203 // calculation of fMeanX, fMeanY and fSize.
204 //
205 fMeanX = 0;
206 fMeanY = 0;
207 fSize = 0;
208
209 MCerPhotPix *pix = 0;
210
211 TIter Next(evt);
212 UInt_t numused = 0;
213 while ((pix=(MCerPhotPix*)Next()))
214 {
215 const MGeomPix &gpix = geom[pix->GetPixId()];
216
217 const Float_t nphot = pix->GetNumPhotons();
218
219 fSize += nphot; // [counter]
220 fMeanX += nphot * gpix.GetX(); // [mm]
221 fMeanY += nphot * gpix.GetY(); // [mm]
222
223 numused++;
224 }
225
226 //
227 // sanity check 2
228 //
229 if (fSize==0)
230 return 2;
231
232 fMeanX /= fSize; // [mm]
233 fMeanY /= fSize; // [mm]
234
235 //
236 // sanity check 3
237 //
238 if (numused<3)
239 return 3;
240
241 //
242 // calculate 2nd moments
243 // ---------------------
244 //
245 // To make sure that the more complicated sum is evaluated as
246 // accurate as possible (because it is needed for more
247 // complicated calculations (see below) we calculate it using
248 // double prcision.
249 //
250 Double_t corrxx=0; // [m^2]
251 Double_t corrxy=0; // [m^2]
252 Double_t corryy=0; // [m^2]
253
254 Next.Reset();
255 while ((pix=(MCerPhotPix*)Next()))
256 {
257 const MGeomPix &gpix = geom[pix->GetPixId()];
258
259 const Float_t dx = gpix.GetX() - fMeanX; // [mm]
260 const Float_t dy = gpix.GetY() - fMeanY; // [mm]
261
262 const Float_t nphot = pix->GetNumPhotons(); // [#phot]
263
264 corrxx += nphot * dx*dx; // [mm^2]
265 corrxy += nphot * dx*dy; // [mm^2]
266 corryy += nphot * dy*dy; // [mm^2]
267 }
268
269 //
270 // If corrxy=0 (which should happen not very often, because fSize>0)
271 // we cannot calculate Length and Width.
272 // In reallity it is almost impossible to have a distribution of
273 // cerenkov photons in the used pixels which is exactly symmetric
274 // along one of the axis.
275 // It seems to be less than 0.1% of all events.
276 //
277 if (corrxy==0)
278 return 4;
279
280 //
281 // calculate the basic Hillas parameters: orientation and size of axes
282 // -------------------------------------------------------------------
283 //
284 // fDelta is the angle between the major axis of the ellipse and
285 // the x axis. It is independent of the position of the ellipse
286 // in the camera it has values between -pi/2 and pi/2 degrees
287 //
288 // Rem: I tested replacing sqrt() by hypot() but they exactly
289 // consume the same amount of time
290 //
291 const Double_t d0 = corryy - corrxx;
292 const Double_t d1 = corrxy*2;
293 const Double_t d2 = d0 + TMath::Sqrt(d0*d0 + d1*d1);
294 const Double_t tand = d2 / d1;
295 const Double_t tand2 = tand*tand;
296
297 fDelta = TMath::ATan(tand);
298
299 const Double_t s2 = tand2+1;
300 const Double_t s = TMath::Sqrt(s2);
301
302 fCosDelta = 1.0/s; // need these in derived classes
303 fSinDelta = tand/s; // like MHillasExt
304
305 const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/fSize;
306 const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/fSize;
307
308 //
309 // fLength^2 is the second moment along the major axis of the ellipse
310 // fWidth^2 is the second moment along the minor axis of the ellipse
311 //
312 // From the algorithm we get: fWidth <= fLength is always true
313 //
314 // very small numbers can get negative by rounding
315 //
316 fLength = axis1<0 ? 0 : TMath::Sqrt(axis1); // [mm]
317 fWidth = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm]
318
319 SetReadyToSave();
320
321 return 0;
322}
323
324// --------------------------------------------------------------------------
325//
326// This function is ment for special usage, please never try to set
327// values via this function
328//
329void MHillas::Set(const TArrayF &arr)
330{
331 if (arr.GetSize() != 6)
332 return;
333
334 fLength = arr.At(0); // [mm] major axis of ellipse
335 fWidth = arr.At(1); // [mm] minor axis of ellipse
336 fDelta = arr.At(2); // [rad] angle of major axis with x-axis
337 fSize = arr.At(3); // [#CerPhot] sum of content of all pixels (number of Cherenkov photons)
338 fMeanX = arr.At(4); // [mm] x-coordinate of center of ellipse
339 fMeanY = arr.At(5); // [mm] y-coordinate of center of ellipse
340}
Note: See TracBrowser for help on using the repository browser.