source: tags/Mars-V0.5/manalysis/MHillas.cc

Last change on this file was 1018, checked in by tbretz, 23 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 8.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): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
19! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de)
20!
21! Copyright: MAGIC Software Development, 2000-2001
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27// //
28// MHillas //
29// //
30// Storage Container for the Hillas parameter //
31// //
32// FIXME: - Here everybody should find an explanation of the parameters //
33// - using boooleans for fIsPixelUsed, fIsPixelCore, ... is rather //
34// slow because you have to loop over all pixels in any loop. //
35// There could be a huge speed improvement using Hash Tables //
36// (linked lists, see THashTable and THashList, too) //
37// //
38/////////////////////////////////////////////////////////////////////////////
39#include "MHillas.h"
40
41#include <math.h>
42#include <fstream.h>
43
44#include <TEllipse.h>
45
46#include "MGeomPix.h"
47#include "MGeomCam.h"
48
49#include "MCerPhotPix.h"
50#include "MCerPhotEvt.h"
51
52#include "MLog.h"
53
54ClassImp(MHillas);
55
56// --------------------------------------------------------------------------
57//
58// Default constructor.
59//
60MHillas::MHillas(const char *name, const char *title) : fEllipse(NULL)
61{
62 fName = name ? name : "MHillas";
63 fTitle = title ? title : "Storage container for Hillas parameter of one event";
64
65 Reset();
66 // FIXME: (intelligent) initialization of values missing
67}
68
69// --------------------------------------------------------------------------
70//
71// Destructor. Deletes the TEllipse if one exists.
72//
73MHillas::~MHillas()
74{
75 Clear();
76}
77
78void MHillas::Reset()
79{
80 fAlpha = 0;
81 fTheta = 0;
82 fWidth = 0;
83 fLength = 0;
84 fSize = 0;
85 fDist = 0;
86
87 Clear();
88}
89
90// --------------------------------------------------------------------------
91//
92// Print the hillas Parameters to *fLog
93//
94void MHillas::Print(Option_t *) const
95{
96 *fLog << "Hillas Parameter:" << endl;
97 *fLog << " - Alpha = " << fabs(fAlpha) << "°" << endl;
98 *fLog << " - Width = " << fWidth << " mm" << endl;
99 *fLog << " - Length = " << fLength << " mm" << endl;
100 *fLog << " - Size = " << fSize << " #CerPhot" << endl;
101 *fLog << " - Dist = " << fDist << " mm" << endl;
102}
103
104/*
105// --------------------------------------------------------------------------
106//
107// call the Paint function of the Ellipse if a TEllipse exists
108//
109void MHillas::Paint(Option_t *)
110{
111 if (!fEllipse)
112 return;
113
114 fEllipse->Paint();
115}
116*/
117
118// --------------------------------------------------------------------------
119//
120// Instead of adding MHillas itself to the Pad
121// (s. AppendPad in TObject) we create an ellipse,
122// which is added to the Pad by it's Draw function
123// You can remove it by deleting the Ellipse Object
124// (s. Clear() )
125//
126void MHillas::Draw(Option_t *opt)
127{
128 Clear();
129
130 fEllipse = new TEllipse(cos(fTheta)*fDist, sin(fTheta)*fDist,
131 fLength, fWidth,
132 0, 360, fTheta*kRad2Deg+fAlpha-180);
133
134 fEllipse->SetLineWidth(2);
135 fEllipse->Draw();
136 //AppendPad(opt);
137
138 /*
139 This is from TH1
140 TString opt = option;
141 opt.ToLower();
142 if (gPad && !opt.Contains("same")) {
143 //the following statement is necessary in case one attempts to draw
144 //a temporary histogram already in the current pad
145 if (TestBit(kCanDelete)) gPad->GetListOfPrimitives()->Remove(this);
146 gPad->Clear();
147 }
148 AppendPad(opt.Data());
149 */
150}
151
152// --------------------------------------------------------------------------
153//
154// If a TEllipse object exists it is deleted
155//
156void MHillas::Clear(Option_t *)
157{
158 if (!fEllipse)
159 return;
160
161 delete fEllipse;
162
163 fEllipse = NULL;
164}
165
166// --------------------------------------------------------------------------
167//
168// Calculate the Hillas parameters from a cerenkov photon event
169// (The calcualtion is some kind of two dimentional statistics)
170//
171// FIXME: MHillas::Calc is rather slow at the moment because it loops
172// unnecessarily over all pixels in all its loops (s.MImgCleanStd)
173// The speed could be improved very much by using Hash-Tables
174// (linked lists, see THashTable and THashList, too)
175//
176Bool_t MHillas::Calc(const MGeomCam &geom, const MCerPhotEvt &evt)
177{
178 const UInt_t nevt = evt.GetNumPixels();
179
180 //
181 // sanity check
182 //
183 if (nevt <= 2)
184 return kFALSE;
185
186 //
187 // calculate mean valu of pixels
188 //
189 float xmean =0;
190 float ymean =0;
191
192 fSize = 0;
193
194 //
195 // FIXME! npix should be retrieved from MCerPhotEvt
196 //
197 UShort_t npix=0;
198 for (UInt_t i=0; i<nevt; i++)
199 {
200 const MCerPhotPix &pix = evt[i];
201
202 if (!pix.IsPixelUsed())
203 continue;
204
205 const MGeomPix &gpix = geom[pix.GetPixId()];
206
207 const float nphot = pix.GetNumPhotons();
208
209 fSize += nphot;
210 xmean += nphot * gpix.GetX(); // [mm]
211 ymean += nphot * gpix.GetY(); // [mm]
212
213 npix++;
214 }
215
216 //
217 // sanity check
218 //
219 if (fSize==0 || npix<=2)
220 return kFALSE;
221
222 xmean /= fSize; // [mm]
223 ymean /= fSize; // [mm]
224
225 //
226 // calculate sdev
227 //
228 float sigmaxx=0;
229 float sigmaxy=0;
230 float sigmayy=0;
231
232 for (UInt_t i=0; i<nevt; i++)
233 {
234 const MCerPhotPix &pix = evt[i];
235
236 if (!pix.IsPixelUsed())
237 continue;
238
239 const MGeomPix &gpix = geom[pix.GetPixId()];
240
241 const float dx = gpix.GetX() - xmean;
242 const float dy = gpix.GetY() - ymean;
243
244 const float nphot = pix.GetNumPhotons();
245
246 sigmaxx += nphot * dx*dx; // [mm^2]
247 sigmaxy += nphot * dx*dy; // [mm^2]
248 sigmayy += nphot * dy*dy; // [mm^2]
249 }
250
251 //
252 // check for orientation
253 //
254 const float theta = atan(sigmaxy/(sigmaxx-sigmayy)*2)/2;
255
256 float c = cos(theta); // [1]
257 float s = sin(theta); // [1]
258
259 //
260 // calculate the length of the two axis
261 //
262 float axis1 = 2.0*c*s*sigmaxy + c*c*sigmaxx + s*s*sigmayy; // [mm^2]
263 float axis2 = -2.0*c*s*sigmaxy + s*s*sigmaxx + c*c*sigmayy; // [mm^2]
264
265 axis1 /= fSize; // [mm^2]
266 axis2 /= fSize; // [mm^2]
267
268 //
269 // check for numerical negatives
270 // (very small number can get negative by chance)
271 //
272 if (axis1 < 0) axis1=0;
273 if (axis2 < 0) axis2=0;
274
275 //
276 // calculate the main Hillas parameters
277 //
278 // fLength, fWidth describes the two axis of the ellipse
279 // fAlpha is the angle between the length-axis and the center
280 // of the camera
281 // fDist is the distance between the center of the camera and the
282 // denter of the ellipse
283 //
284 const int rotation = axis1<axis2;
285
286 fLength = rotation ? sqrt(axis2) : sqrt(axis1); // [mm]
287 fWidth = rotation ? sqrt(axis1) : sqrt(axis2); // [mm]
288
289 const float a = c*xmean + s*ymean;
290 const float b = c*ymean - s*xmean;
291
292 fAlpha = rotation ? atan(a/b) : atan(-b/a); // [rad]
293 fAlpha *= kRad2Deg; // [deg]
294
295 fDist = sqrt(xmean*xmean + ymean*ymean); // [mm]
296
297 fTheta = atan(ymean/xmean); // [rad]
298 if (xmean<0) fTheta += kPI; // [rad]
299
300 SetReadyToSave();
301
302 return kTRUE;
303}
304
305void MHillas::AsciiRead(ifstream &fin)
306{
307 fin >> fAlpha;
308 fin >> fTheta;
309 fin >> fWidth;
310 fin >> fLength;
311 fin >> fSize;
312 fin >> fDist;
313}
314
315void MHillas::AsciiWrite(ofstream &fout) const
316{
317 fout << fAlpha << " ";
318 fout << fTheta << " ";
319 fout << fWidth << " ";
320 fout << fLength << " ";
321 fout << fSize << " ";
322 fout << fDist << endl;
323}
Note: See TracBrowser for help on using the repository browser.