source: branches/removing_cpp11_features/mimage/MHillas.cc@ 18558

Last change on this file since 18558 was 9374, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 12.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): 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-2008
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 <TArrayF.h>
63
64#include <TLine.h>
65#include <TEllipse.h>
66#include <TVector2.h>
67
68#include "MGeom.h"
69#include "MGeomCam.h"
70
71#include "MSignalPix.h"
72#include "MSignalCam.h"
73
74#include "MLog.h"
75#include "MLogManip.h"
76
77ClassImp(MHillas);
78
79using namespace std;
80
81// --------------------------------------------------------------------------
82//
83// Default constructor.
84//
85MHillas::MHillas(const char *name, const char *title)
86{
87 fName = name ? name : "MHillas";
88 fTitle = title ? title : "Storage container for image parameters of one event";
89
90 Reset();
91}
92
93// --------------------------------------------------------------------------
94//
95// Initializes the values with defaults. For the default values see the
96// source code.
97//
98void MHillas::Reset()
99{
100 fLength = -1;
101 fWidth = -1;
102 fDelta = 0;
103
104 fSize = -1;
105 fMeanX = 0;
106 fMeanY = 0;
107}
108
109// --------------------------------------------------------------------------
110//
111// return distance to center (coordinate origin)
112//
113Double_t MHillas::GetDist0() const
114{
115 return TMath::Hypot(fMeanX, fMeanY);
116}
117
118// --------------------------------------------------------------------------
119//
120// return the mean position as TVector2(meanx, meany)
121//
122TVector2 MHillas::GetMean() const
123{
124 return TVector2(fMeanX, fMeanY);
125}
126
127// --------------------------------------------------------------------------
128//
129// return the direction of the main axis (normalized). The main axis
130// always has the direction --> x>0
131//
132TVector2 MHillas::GetNormAxis() const
133{
134 return TVector2(fCosDelta, fSinDelta);
135}
136
137// --------------------------------------------------------------------------
138//
139// Analytical estimation of length of border line:
140// U = pi*(a+b - sqrt(a*b))
141//
142// GetBorderLine() [mm]
143//
144Double_t MHillas::GetBorderLine() const
145{
146 const Double_t a = fWidth+fLength;
147 const Double_t s = fWidth*fLength;
148
149 return TMath::Pi()*(1.5*a - TMath::Sqrt(s));
150}
151
152// --------------------------------------------------------------------------
153//
154// Analytical estimation of length of border line:
155// A = pi*a*b
156//
157// GetArea() [mm^2]
158//
159Double_t MHillas::GetArea() const
160{
161 return TMath::Pi()*fWidth*fLength;
162}
163
164// --------------------------------------------------------------------------
165//
166// Print the hillas Parameters to *fLog
167//
168void MHillas::Print(Option_t *) const
169{
170 const Double_t atg = atan2(fMeanX, fMeanY)*kRad2Deg;
171
172 *fLog << all;
173 *fLog << GetDescriptor() << endl;
174 *fLog << " - Length [mm] = " << fLength << endl;
175 *fLog << " - Width [mm] = " << fWidth << endl;
176 *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl;
177 *fLog << " - Size [phe] = " << fSize << endl;
178 *fLog << " - Meanx [mm] = " << fMeanX << endl;
179 *fLog << " - Meany [mm] = " << fMeanY << endl;
180 *fLog << " - atg(y/x) [deg] = " << atg << endl;
181}
182
183// --------------------------------------------------------------------------
184//
185// Print the hillas Parameters to *fLog, depending on the geometry in
186// units of deg
187//
188void MHillas::Print(const MGeomCam &geom) const
189{
190 const Double_t atg = atan2(fMeanX, fMeanY)*kRad2Deg;
191
192 *fLog << all;
193 *fLog << GetDescriptor() << endl;
194 *fLog << " - Length [deg] = " << fLength*geom.GetConvMm2Deg() << endl;
195 *fLog << " - Width [deg] = " << fWidth *geom.GetConvMm2Deg() << endl;
196 *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl;
197 *fLog << " - Size [phe] = " << fSize << endl;
198 *fLog << " - Meanx [deg] = " << fMeanX *geom.GetConvMm2Deg() << endl;
199 *fLog << " - Meany [deg] = " << fMeanY *geom.GetConvMm2Deg() << endl;
200 *fLog << " - atg(y/x) [deg] = " << atg << endl;
201}
202
203// --------------------------------------------------------------------------
204//
205// Paint the ellipse corresponding to the parameters
206//
207void MHillas::Paint(Option_t *opt)
208{
209 static const Float_t sOffsetW=40.0;
210 static const Float_t sOffsetL=266.0;
211
212 if (fLength<0 || fWidth<0)
213 return;
214
215 const Float_t l = fLength+sOffsetL;
216 const Float_t w = fWidth +sOffsetW;
217
218 TLine line;
219 line.SetLineWidth(1);
220 line.SetLineColor(kRed);
221
222 // Length line
223 line.PaintLine(-l*fCosDelta+fMeanX, -l*fSinDelta+fMeanY,
224 +l*fCosDelta+fMeanX, +l*fSinDelta+fMeanY);
225
226
227 // Width line
228 line.PaintLine(+w*fSinDelta+fMeanX, -w*fCosDelta+fMeanY,
229 -w*fSinDelta+fMeanX, +w*fCosDelta+fMeanY);
230
231 // Rough estimate for disp
232 const Double_t p = 362*(1-fWidth/fLength);
233
234 // Disp line
235 line.PaintLine( p*fCosDelta+20*fSinDelta+fMeanX, p*fSinDelta-20*fCosDelta+fMeanY,
236 p*fCosDelta-20*fSinDelta+fMeanX, p*fSinDelta+20*fCosDelta+fMeanY);
237 line.PaintLine(-p*fCosDelta+20*fSinDelta+fMeanX, -p*fSinDelta-20*fCosDelta+fMeanY,
238 -p*fCosDelta-20*fSinDelta+fMeanX, -p*fSinDelta+20*fCosDelta+fMeanY);
239
240 TEllipse e(fMeanX, fMeanY, fLength, fWidth, 0, 360, fDelta*kRad2Deg+180);
241 e.SetLineWidth(2);
242 e.SetLineColor(kGreen);
243 e.Paint();
244}
245
246// --------------------------------------------------------------------------
247//
248// Calculate the image parameters from a Cherenkov photon event
249// assuming Cher.photons/pixel and pixel coordinates are given
250// In case you don't call Calc from within an eventloop make sure, that
251// you call the Reset member function before.
252// Returns:
253// 0 no error
254// 1 number of pixels in event == 0 (special MC events)
255// 2 size==0
256// 3 number of used pixel < 3
257// 4 CorrXY == 0
258//
259Int_t MHillas::Calc(const MGeomCam &geom, const MSignalCam &evt, Int_t island)
260{
261 const UInt_t numpix = evt.GetNumPixels();
262
263 //
264 // sanity check 1 (special MC events)
265 //
266 if (numpix==0)
267 return 1;
268
269 //
270 // calculate mean value of pixel coordinates and fSize
271 // -----------------------------------------------------
272 //
273 // Because this are only simple sums of roughly 600 values
274 // with an accuracy less than three digits there should not
275 // be any need to use double precision (Double_t) for the
276 // calculation of fMeanX, fMeanY and fSize.
277 //
278 fMeanX = 0;
279 fMeanY = 0;
280 fSize = 0;
281
282 UInt_t numused = 0;
283
284 for (UInt_t i=0; i<numpix; i++)
285 {
286 MSignalPix &pix = evt[i];
287 if (!pix.IsPixelUsed())
288 continue;
289
290 if (island>=0 && pix.GetIdxIsland()!=island)
291 continue;
292
293 const MGeom &gpix = geom[i];
294
295 const Float_t nphot = pix.GetNumPhotons();
296
297 fSize += nphot; // [counter]
298 fMeanX += nphot * gpix.GetX(); // [mm]
299 fMeanY += nphot * gpix.GetY(); // [mm]
300
301 numused++;
302 }
303
304 //
305 // sanity check 2
306 //
307 if (fSize==0)
308 return 2;
309
310 fMeanX /= fSize; // [mm]
311 fMeanY /= fSize; // [mm]
312
313 //
314 // sanity check 3
315 //
316 if (numused<3)
317 return 3;
318
319 //
320 // calculate 2nd moments
321 // ---------------------
322 //
323 // To make sure that the more complicated sum is evaluated as
324 // accurate as possible (because it is needed for more
325 // complicated calculations (see below) we calculate it using
326 // double prcision.
327 //
328 Double_t corrxx=0; // [m^2]
329 Double_t corrxy=0; // [m^2]
330 Double_t corryy=0; // [m^2]
331
332 for (UInt_t i=0; i<numpix; i++)
333 {
334 const MSignalPix &pix = evt[i];
335 if (!pix.IsPixelUsed())
336 continue;
337
338 if (island>=0 && pix.GetIdxIsland()!=island)
339 continue;
340
341 const MGeom &gpix = geom[i];
342
343 const Float_t dx = gpix.GetX() - fMeanX; // [mm]
344 const Float_t dy = gpix.GetY() - fMeanY; // [mm]
345
346 const Float_t nphot = pix.GetNumPhotons(); // [#phot]
347
348 corrxx += nphot * dx*dx; // [mm^2]
349 corrxy += nphot * dx*dy; // [mm^2]
350 corryy += nphot * dy*dy; // [mm^2]
351 }
352
353 //
354 // calculate the basic Hillas parameters: orientation and size of axes
355 // -------------------------------------------------------------------
356 //
357 // fDelta is the angle between the major axis of the ellipse and
358 // the x axis. It is independent of the position of the ellipse
359 // in the camera it has values between -pi/2 and pi/2 degrees
360 //
361 // Rem: I tested replacing sqrt() by hypot() but they exactly
362 // consume the same amount of time
363 //
364 const Double_t d0 = corryy - corrxx;
365 const Double_t d1 = corrxy*2;
366 const Double_t d2 = TMath::Sqrt(d0*d0 + d1*d1) + d0; // [0
367
368 const Double_t tand = d2==0 ? 0 : d2 / d1; // Force 0 also if d1==0
369 const Double_t tand2 = tand*tand;
370
371 const Double_t s2 = tand2+1;
372 const Double_t s = TMath::Sqrt(s2);
373
374 // Set default for the case in which the image is symmetric on the y-axis
375 fDelta = TMath::Pi()/2;
376 fCosDelta = 0;
377 fSinDelta = 1;
378
379 Double_t axis1 = corryy;
380 Double_t axis2 = corrxx;
381
382 // This are all cases in which the image is not symmetric on the y-axis
383 if (d1!=0 || d2==0)
384 {
385 fDelta = TMath::ATan(tand);
386
387 fCosDelta = 1.0 /s; // need these in derived classes
388 fSinDelta = tand/s; // like MHillasExt
389
390 axis1 = (tand2*corryy + d2 + corrxx)/s2;
391 axis2 = (tand2*corrxx - d2 + corryy)/s2;
392 }
393
394 //
395 // fLength^2 is the second moment along the major axis of the distribution
396 // fWidth^2 is the second moment along the minor axis of the distribution
397 //
398 // From the algorithm we get: fWidth <= fLength is always true
399 //
400 // very small numbers can get negative by rounding
401 //
402 fLength = axis1<0 ? 0 : TMath::Sqrt(axis1/fSize); // [mm]
403 fWidth = axis2<0 ? 0 : TMath::Sqrt(axis2/fSize); // [mm]
404
405 SetReadyToSave();
406
407 return 0;
408}
409
410// --------------------------------------------------------------------------
411//
412// This function is ment for special usage, please never try to set
413// values via this function
414//
415void MHillas::Set(const TArrayF &arr)
416{
417 if (arr.GetSize() != 6)
418 return;
419
420 fLength = arr.At(0); // [mm] major axis of ellipse
421 fWidth = arr.At(1); // [mm] minor axis of ellipse
422 fDelta = arr.At(2); // [rad] angle of major axis with x-axis
423 fSize = arr.At(3); // [#CerPhot] sum of content of all pixels (number of Cherenkov photons)
424 fMeanX = arr.At(4); // [mm] x-coordinate of center of ellipse
425 fMeanY = arr.At(5); // [mm] y-coordinate of center of ellipse
426}
Note: See TracBrowser for help on using the repository browser.