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 | !
|
---|
22 | ! Copyright: MAGIC Software Development, 2000-2002
|
---|
23 | !
|
---|
24 | !
|
---|
25 | \* ======================================================================== */
|
---|
26 |
|
---|
27 | /////////////////////////////////////////////////////////////////////////////
|
---|
28 | //
|
---|
29 | // MHillas
|
---|
30 | //
|
---|
31 | // Storage Container for image parameters
|
---|
32 | //
|
---|
33 | // basic image parameters
|
---|
34 | //
|
---|
35 | // Version 1:
|
---|
36 | // ----------
|
---|
37 | // fLength major axis of ellipse
|
---|
38 | // fWidth minor axis
|
---|
39 | // fDelta angle of major axis wrt x-axis
|
---|
40 | // fSize total sum of pixels
|
---|
41 | // fMeanX x of center of ellipse
|
---|
42 | // fMeanY y of center of ellipse
|
---|
43 | //
|
---|
44 | // Version 2:
|
---|
45 | // ----------
|
---|
46 | // fNumCorePixels number of pixels called core
|
---|
47 | // fNumUsedPixels number of pixels which survived the cleaning
|
---|
48 | //
|
---|
49 | /////////////////////////////////////////////////////////////////////////////
|
---|
50 | #include "MHillas.h"
|
---|
51 |
|
---|
52 | #include <fstream.h>
|
---|
53 |
|
---|
54 | #include <TEllipse.h>
|
---|
55 |
|
---|
56 | #include "MGeomPix.h"
|
---|
57 | #include "MGeomCam.h"
|
---|
58 |
|
---|
59 | #include "MCerPhotPix.h"
|
---|
60 | #include "MCerPhotEvt.h"
|
---|
61 |
|
---|
62 | #include "MLog.h"
|
---|
63 | #include "MLogManip.h"
|
---|
64 |
|
---|
65 | ClassImp(MHillas);
|
---|
66 |
|
---|
67 | // --------------------------------------------------------------------------
|
---|
68 | //
|
---|
69 | // Default constructor.
|
---|
70 | //
|
---|
71 | MHillas::MHillas(const char *name, const char *title) : fEllipse(NULL)
|
---|
72 | {
|
---|
73 | fName = name ? name : "MHillas";
|
---|
74 | fTitle = title ? title : "Storage container for image parameters of one event";
|
---|
75 |
|
---|
76 | Reset();
|
---|
77 | // FIXME: (intelligent) initialization of values missing
|
---|
78 |
|
---|
79 | fEllipse = new TEllipse;
|
---|
80 | }
|
---|
81 |
|
---|
82 | // --------------------------------------------------------------------------
|
---|
83 | //
|
---|
84 | // Destructor. Deletes the TEllipse if one exists.
|
---|
85 | //
|
---|
86 | MHillas::~MHillas()
|
---|
87 | {
|
---|
88 | Clear();
|
---|
89 | }
|
---|
90 |
|
---|
91 | // --------------------------------------------------------------------------
|
---|
92 | //
|
---|
93 | void MHillas::Reset()
|
---|
94 | {
|
---|
95 | fLength = 0;
|
---|
96 | fWidth = 0;
|
---|
97 | fDelta = 0;
|
---|
98 |
|
---|
99 | fSize = -1;
|
---|
100 | fMeanX = -1;
|
---|
101 | fMeanY = -1;
|
---|
102 |
|
---|
103 | fNumUsedPixels = -1;
|
---|
104 | fNumCorePixels = -1;
|
---|
105 |
|
---|
106 | Clear();
|
---|
107 | }
|
---|
108 |
|
---|
109 | // --------------------------------------------------------------------------
|
---|
110 | //
|
---|
111 | // Print the hillas Parameters to *fLog
|
---|
112 | //
|
---|
113 | void MHillas::Print(Option_t *) const
|
---|
114 | {
|
---|
115 | Double_t atg = atan2(fMeanY, fMeanX)*kRad2Deg;
|
---|
116 |
|
---|
117 | if (atg<0)
|
---|
118 | atg += 180;
|
---|
119 |
|
---|
120 | *fLog << all;
|
---|
121 | *fLog << "Basic Image Parameters (" << GetName() << ")" << endl;
|
---|
122 | *fLog << " - Length [mm] = " << fLength << endl;
|
---|
123 | *fLog << " - Width [mm] = " << fWidth << endl;
|
---|
124 | *fLog << " - Meanx [mm] = " << fMeanX << endl;
|
---|
125 | *fLog << " - Meany [mm] = " << fMeanY << endl;
|
---|
126 | *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl;
|
---|
127 | *fLog << " - atg(y/x) [deg] = " << atg << endl;
|
---|
128 | *fLog << " - Size [1] = " << fSize << " #CherPhot" << endl;
|
---|
129 | }
|
---|
130 |
|
---|
131 | /*
|
---|
132 | // -----------------------------------------------------------
|
---|
133 | //
|
---|
134 | // call the Paint function of the Ellipse if a TEllipse exists
|
---|
135 | //
|
---|
136 | void MHillas::Paint(Option_t *)
|
---|
137 | {
|
---|
138 | fEllipse->SetLineWidth(2);
|
---|
139 | fEllipse->PaintEllipse(fMeanX, fMeanY, fLength, fWidth,
|
---|
140 | 0, 360, fDelta*kRad2Deg+180);
|
---|
141 | }
|
---|
142 | */
|
---|
143 |
|
---|
144 | // --------------------------------------------------------------------------
|
---|
145 | //
|
---|
146 | // Instead of adding MHillas itself to the Pad
|
---|
147 | // (s. AppendPad in TObject) we create an ellipse,
|
---|
148 | // which is added to the Pad by its Draw function
|
---|
149 | // You can remove it by deleting the Ellipse Object
|
---|
150 | // (s. Clear() )
|
---|
151 | //
|
---|
152 | void MHillas::Draw(Option_t *opt)
|
---|
153 | {
|
---|
154 |
|
---|
155 | Clear();
|
---|
156 |
|
---|
157 | fEllipse = new TEllipse(fMeanX, fMeanY, fLength, fWidth,
|
---|
158 | 0, 360, fDelta*kRad2Deg+180);
|
---|
159 |
|
---|
160 | fEllipse->SetLineWidth(2);
|
---|
161 | fEllipse->Draw();
|
---|
162 |
|
---|
163 | /*
|
---|
164 | fEllipse->SetPhimin();
|
---|
165 | fEllipse->SetPhimax();
|
---|
166 | fEllipse->SetR1(fLength);
|
---|
167 | fEllipse->SetR2(fWidth);
|
---|
168 | fEllipse->SetTheta(fDelta*kRad2Deg+180);
|
---|
169 | fEllipse->SetX1(fMeanX);
|
---|
170 | fEllipse->SetY1(fMeanY);
|
---|
171 |
|
---|
172 | fEllipse->SetLineWidth(2);
|
---|
173 | fEllipse->PaintEllipse(fMeanX, fMeanY, fLength, fWidth,
|
---|
174 | 0, 360, fDelta*kRad2Deg+180);
|
---|
175 |
|
---|
176 | AppendPad(opt);
|
---|
177 |
|
---|
178 | // This is from TH1
|
---|
179 | TString opt = option;
|
---|
180 | opt.ToLower();
|
---|
181 | if (gPad && !opt.Contains("same")) {
|
---|
182 | //the following statement is necessary in case one attempts to draw
|
---|
183 | //a temporary histogram already in the current pad
|
---|
184 | if (TestBit(kCanDelete)) gPad->GetListOfPrimitives()->Remove(this);
|
---|
185 | gPad->Clear();
|
---|
186 | }
|
---|
187 | */
|
---|
188 | }
|
---|
189 |
|
---|
190 | // --------------------------------------------------------------------------
|
---|
191 | //
|
---|
192 | // If a TEllipse object exists it is deleted
|
---|
193 | //
|
---|
194 | void MHillas::Clear(Option_t *)
|
---|
195 | {
|
---|
196 | if (!fEllipse)
|
---|
197 | return;
|
---|
198 |
|
---|
199 | delete fEllipse;
|
---|
200 |
|
---|
201 | fEllipse = NULL;
|
---|
202 | }
|
---|
203 |
|
---|
204 |
|
---|
205 | // --------------------------------------------------------------------------
|
---|
206 | //
|
---|
207 | // Calculate the image parameters from a Cherenkov photon event
|
---|
208 | // assuming Cher.photons/pixel and pixel coordinates are given
|
---|
209 | //
|
---|
210 | Bool_t MHillas::Calc(const MGeomCam &geom, const MCerPhotEvt &evt)
|
---|
211 | {
|
---|
212 | // FIXME: MHillas::Calc is rather slow at the moment because it loops
|
---|
213 | // unnecessarily over all pixels in all its loops (s.MImgCleanStd)
|
---|
214 | // The speed could be improved very much by using Hash-Tables
|
---|
215 | // (linked lists, see THashTable and THashList, too)
|
---|
216 | //
|
---|
217 |
|
---|
218 | const UInt_t npixevt = evt.GetNumPixels();
|
---|
219 |
|
---|
220 | //
|
---|
221 | // sanity check
|
---|
222 | //
|
---|
223 | if (npixevt <= 2)
|
---|
224 | return kFALSE;
|
---|
225 |
|
---|
226 | //
|
---|
227 | // calculate mean value of pixel coordinates and fSize
|
---|
228 | // -----------------------------------------------------
|
---|
229 | //
|
---|
230 | fMeanX = 0;
|
---|
231 | fMeanY = 0;
|
---|
232 | fSize = 0;
|
---|
233 |
|
---|
234 | fNumUsedPixels = 0;
|
---|
235 | fNumCorePixels = 0;
|
---|
236 |
|
---|
237 | //
|
---|
238 | // FIXME! npix should be retrieved from MCerPhotEvt
|
---|
239 | //
|
---|
240 | for (UInt_t i=0; i<npixevt; i++)
|
---|
241 | {
|
---|
242 | const MCerPhotPix &pix = evt[i];
|
---|
243 |
|
---|
244 | if (!pix.IsPixelUsed())
|
---|
245 | continue;
|
---|
246 |
|
---|
247 | const MGeomPix &gpix = geom[pix.GetPixId()];
|
---|
248 |
|
---|
249 | const float nphot = pix.GetNumPhotons();
|
---|
250 |
|
---|
251 | fSize += nphot; // [counter]
|
---|
252 | fMeanX += nphot * gpix.GetX(); // [mm]
|
---|
253 | fMeanY += nphot * gpix.GetY(); // [mm]
|
---|
254 |
|
---|
255 | if (pix.IsPixelCore())
|
---|
256 | fNumCorePixels++;
|
---|
257 |
|
---|
258 | fNumUsedPixels++;
|
---|
259 | }
|
---|
260 |
|
---|
261 | //
|
---|
262 | // sanity check
|
---|
263 | //
|
---|
264 | if (fSize==0 || fNumUsedPixels<=2)
|
---|
265 | return kFALSE;
|
---|
266 |
|
---|
267 | fMeanX /= fSize; // [mm]
|
---|
268 | fMeanY /= fSize; // [mm]
|
---|
269 |
|
---|
270 | //
|
---|
271 | // calculate 2nd moments
|
---|
272 | // -------------------
|
---|
273 | //
|
---|
274 | float corrxx=0; // [m^2]
|
---|
275 | float corrxy=0; // [m^2]
|
---|
276 | float corryy=0; // [m^2]
|
---|
277 |
|
---|
278 | for (UInt_t i=0; i<npixevt; i++)
|
---|
279 | {
|
---|
280 | const MCerPhotPix &pix = evt[i];
|
---|
281 |
|
---|
282 | if (!pix.IsPixelUsed())
|
---|
283 | continue;
|
---|
284 |
|
---|
285 | const MGeomPix &gpix = geom[pix.GetPixId()];
|
---|
286 | const float dx = gpix.GetX() - fMeanX; // [mm]
|
---|
287 | const float dy = gpix.GetY() - fMeanY; // [mm]
|
---|
288 |
|
---|
289 | const float nphot = pix.GetNumPhotons(); // [#phot]
|
---|
290 |
|
---|
291 | corrxx += nphot * dx*dx; // [mm^2]
|
---|
292 | corrxy += nphot * dx*dy; // [mm^2]
|
---|
293 | corryy += nphot * dy*dy; // [mm^2]
|
---|
294 | }
|
---|
295 |
|
---|
296 | //
|
---|
297 | // calculate the basic Hillas parameters: orientation and size of axes
|
---|
298 | // -------------------------------------------------------------------
|
---|
299 | //
|
---|
300 | const float d = (corryy - corrxx)/(corrxy*2);
|
---|
301 |
|
---|
302 | fDelta = atan(d + sqrt(d*d + 1));
|
---|
303 |
|
---|
304 | fCosDelta = cos(fDelta); // need these in derived classes
|
---|
305 | fSinDelta = sin(fDelta); // like MHillasExt
|
---|
306 |
|
---|
307 | float axis1 = ( fCosDelta*fSinDelta*corrxy*2 + fCosDelta*fCosDelta*corrxx + fSinDelta*fSinDelta*corryy) / fSize; // [mm^2]
|
---|
308 | float axis2 = (-fCosDelta*fSinDelta*corrxy*2 + fSinDelta*fSinDelta*corrxx + fCosDelta*fCosDelta*corryy) / fSize; // [mm^2]
|
---|
309 |
|
---|
310 | // very small numbers can get negative by rounding
|
---|
311 | if (axis1 < 0) axis1=0;
|
---|
312 | if (axis2 < 0) axis2=0;
|
---|
313 |
|
---|
314 | fLength = sqrt(axis1); // [mm]
|
---|
315 | fWidth = sqrt(axis2); // [mm]
|
---|
316 |
|
---|
317 | SetReadyToSave();
|
---|
318 |
|
---|
319 | return kTRUE;
|
---|
320 | }
|
---|
321 |
|
---|
322 | // --------------------------------------------------------------------------
|
---|
323 | //
|
---|
324 | void MHillas::AsciiRead(ifstream &fin)
|
---|
325 | {
|
---|
326 | fin >> fLength;
|
---|
327 | fin >> fWidth;
|
---|
328 | fin >> fDelta;
|
---|
329 | fin >> fSize;
|
---|
330 | fin >> fMeanX;
|
---|
331 | fin >> fMeanY;
|
---|
332 | }
|
---|
333 |
|
---|
334 | // --------------------------------------------------------------------------
|
---|
335 | /*
|
---|
336 | void MHillas::AsciiWrite(ofstream &fout) const
|
---|
337 | {
|
---|
338 | fout << fLength << " ";
|
---|
339 | fout << fWidth << " ";
|
---|
340 | fout << fDelta << " ";
|
---|
341 | fout << fSize << " ";
|
---|
342 | fout << fMeanX << " ";
|
---|
343 | fout << fMeanY;
|
---|
344 | }
|
---|
345 | */ |
---|