source: trunk/Mars/mimage/MHillasExt.cc@ 20037

Last change on this file since 20037 was 19203, checked in by tbretz, 7 years ago
Still requires numerical protection for a few very unlikely numerical uncertainty cases.
File size: 12.0 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@astro.uni-wuerzburg.de>
19! Author(s): Rudolf Bock 10/2001 <mailto:Rudolf.Bock@cern.ch>
20! Author(s): Wolfgang Wittek 06/2002 <mailto:wittek@mppmu.mpg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2010
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28//
29// MHillasExt
30//
31// Storage Container for extended image parameters
32//
33// extended image parameters
34//
35//
36// Version 1:
37// ----------
38// fConc ratio of sum of two highest pixels over fSize
39// fConc1 ratio of highest pixel over fSize
40// fAsym distance from highest pixel to center, projected onto major axis
41// fM3Long third moment along major axis
42// fM3Trans third moment along minor axis
43//
44// Version 2:
45// ----------
46// fConc removed
47// fConc1 removed
48//
49// Version 3:
50// ----------
51// fMaxDist added. Distance between center and most distant used pixel
52//
53// Version 4:
54// ----------
55// fMaxDist removed.
56// fSlopeLong added
57// fSlopeTrans added
58//
59// Version 5:
60// ----------
61// fTimeSpread added
62// fTimeSpreadWeighted added
63//
64// fSlopeSpread added
65// fSlopeSpreadWeighted added
66//
67//
68// WARNING: Before you can use fAsym, fSlope* and fM3* you must
69// multiply by the sign of MHillasSrc::fCosDeltaAlpha
70//
71////////////////////////////////////////////////////////////////////////////
72/*
73 // fAsymna d/(d na) of ( sum(x*q^na)/sum(q^na), sum(y*q^na)/sum(q^na) )
74 // projected onto the major axis
75 // fAsym0 (F-B)/(F+B) along the major axis
76 */
77#include "MHillasExt.h"
78
79#include <TArrayF.h>
80#include <TMarker.h>
81#include <TVector2.h>
82#include <TMinuit.h>
83#include <TVirtualPad.h>
84
85#include "MGeomPix.h"
86#include "MGeomCam.h"
87
88#include "MSignalPix.h"
89#include "MSignalCam.h"
90
91#include "MLog.h"
92#include "MLogManip.h"
93
94#include "MHillas.h"
95
96ClassImp(MHillasExt);
97
98using namespace std;
99
100// -------------------------------------------------------------------------
101//
102// Default constructor.
103//
104MHillasExt::MHillasExt(const char *name, const char *title)
105{
106 fName = name ? name : "MHillasExt";
107 fTitle = title ? title : "Storage container for extended parameter set of one event";
108
109 Reset();
110}
111
112// -------------------------------------------------------------------------
113//
114// Reset all values to 0
115//
116void MHillasExt::Reset()
117{
118 fAsym = 0;
119 fM3Long = 0;
120 fM3Trans = 0;
121 fSlopeLong = 0;
122 fSlopeTrans = 0;
123
124 fTimeSpread = -1;
125 fTimeSpreadWeighted = -1;
126 fSlopeSpread = -1;
127 fSlopeSpreadWeighted = -1;
128}
129
130// -------------------------------------------------------------------------
131//
132// calculation of additional parameters based on the camera geometry
133// and the cerenkov photon event
134//
135Int_t MHillasExt::Calc(const MGeomCam &geom, const MSignalCam &evt, const MHillas &hil, Int_t island)
136{
137 //
138 // calculate the additional image parameters
139 // --------------------------------------------
140 //
141 // loop to get third moments along ellipse axes and two max pixels
142 //
143 // For the moments double precision is used to make sure, that
144 // the complex matrix multiplication and sum is evaluated correctly.
145 //
146 Int_t maxpixid = 0;
147 Float_t maxpix = 0;
148
149 // Variables to caluclate time slope
150 // Formula: http://mo.mathematik.uni-stuttgart.de/inhalt/erlaeuterung/erlaeuterung300/
151 UInt_t cnt = 0;
152
153 Double_t sumx = 0;
154 Double_t sumy = 0;
155 Double_t sumt = 0;
156 Double_t sumw = 0;
157
158 Double_t sumxt = 0;
159 Double_t sumyt = 0;
160
161 Double_t sumx2 = 0;
162 Double_t sumy2 = 0;
163 Double_t sumxy = 0;
164 Double_t sumt2 = 0;
165
166 Double_t sumxw = 0;
167 Double_t sumyw = 0;
168 Double_t sumtw = 0;
169
170 Double_t sumx2w = 0;
171 Double_t sumy2w = 0;
172 Double_t sumt2w = 0;
173 Double_t sumxyw = 0;
174 Double_t sumxtw = 0;
175 Double_t sumytw = 0;
176
177 Double_t sumdx3w = 0;
178 Double_t sumdy3w = 0;
179
180 Double_t sumdx2dyw = 0;
181 Double_t sumdxdy2w = 0;
182
183 const UInt_t npix = evt.GetNumPixels();
184 for (UInt_t i=0; i<npix; i++)
185 {
186 const MSignalPix &pix = evt[i];
187 if (!pix.IsPixelUsed())
188 continue;
189
190 if (island>=0 && pix.GetIdxIsland()!=island)
191 continue;
192
193 Double_t nphot = pix.GetNumPhotons(); // [#phot]
194 if (nphot<0)
195 continue;
196
197 const MGeom &gpix = geom[i];
198
199 const Double_t x = gpix.GetX();
200 const Double_t y = gpix.GetY();
201 const Double_t t = pix.GetArrivalTime();
202
203 // --- time slope ----
204 sumx += x;
205 sumy += y;
206 sumt += t;
207 sumw += nphot;
208
209 sumxt += x*t;
210 sumyt += y*t;
211
212 sumx2 += x*x;
213 sumy2 += y*y;
214 sumxy += x*y;
215 sumt2 += t*t;
216
217 sumxw += x*nphot;
218 sumyw += y*nphot;
219 sumtw += t*nphot;
220
221 sumx2w += x*x*nphot;
222 sumy2w += y*y*nphot;
223 sumt2w += t*t*nphot;
224
225 sumxyw += x*y*nphot;
226
227 sumxtw += x*t*nphot;
228 sumytw += y*t*nphot;
229
230 // --- 3rd moment ---
231 const Double_t dx = x - hil.GetMeanX(); // [mm]
232 const Double_t dy = y - hil.GetMeanY(); // [mm]
233
234 sumdx3w += dx*dx*dx*nphot;
235 sumdy3w += dy*dy*dy*nphot;
236 sumdx2dyw += dx*dx*dy*nphot;
237 sumdxdy2w += dx*dy*dy*nphot;
238
239 cnt++;
240
241 //
242 // Now we are working on absolute values of nphot, which
243 // must take pixel size into account
244 //
245 nphot *= geom.GetPixRatio(i);
246
247 // --- max pixel ---
248 if (nphot>maxpix)
249 {
250 maxpix = nphot; // [1]
251 maxpixid = i;
252 }
253 }
254
255 const Double_t c = hil.GetCosDelta();
256 const Double_t s = hil.GetSinDelta();
257
258 //
259 // Time slope
260 //
261 const Double_t dxt = c*sumxt + s*sumyt;
262 const Double_t dyt = -s*sumxt + c*sumyt;
263
264 const Double_t dx = c*sumx + s*sumy;
265 const Double_t dy = -s*sumx + c*sumy;
266
267 const Double_t dx2 = c*c*sumx2 + 2*c*s*sumxy + s*s*sumy2;
268 const Double_t dy2 = s*s*sumx2 - 2*c*s*sumxy + c*c*sumy2;
269
270 const Double_t detx = cnt*dx2 - dx*dx;
271 const Double_t dety = cnt*dy2 - dy*dy;
272
273 fSlopeLong = detx==0 ? 0 : (cnt*dxt - sumt*dx)/detx;
274 fSlopeTrans = dety==0 ? 0 : (cnt*dyt - sumt*dy)/dety;
275
276 //
277 // Time spread
278 //
279 const Double_t m = fSlopeLong;
280
281 const Double_t z = dx; //c*sumx + s*sumy;
282 const Double_t zw = c*sumxw + s*sumyw;
283
284 const Double_t zt = dxt; //c*sumxt + s*sumyt;
285 const Double_t ztw = c*sumxtw + s*sumytw;
286
287 const Double_t x2y2 = dx2; //c*c*sumx2 + s*s*sumy2 + 2*c*s*sumxy;
288 const Double_t x2y2w = c*c*sumx2w + s*s*sumy2w + 2*c*s*sumxyw;
289
290 const Double_t sumdt = sumt - m*z;
291 const Double_t sumdtw = sumtw - m*zw;
292
293 const Double_t sumdt2 = sumt2 + m*m*x2y2 - 2*m*zt;
294 const Double_t sumdt2w = sumt2w + m*m*x2y2w - 2*m*ztw;
295
296 const Double_t meant = sumt /cnt;
297 const Double_t meandt = sumdt /cnt;
298 const Double_t meantw = sumtw /sumw;
299 const Double_t meandtw = sumdtw/sumw;
300
301 const Double_t rad_ts = sumt2 /cnt - meant *meant;
302 const Double_t rad_tsw = sumt2w /sumw - meantw*meantw;
303 const Double_t rad_ss = sumdt2 /cnt - meandt *meandt;
304 const Double_t rad_ssw = sumdt2w/sumw - meandtw*meandtw;
305
306 fTimeSpread = rad_ts <0 ? 0 : TMath::Sqrt(rad_ts);
307 fTimeSpreadWeighted = rad_tsw<0 ? 0 : TMath::Sqrt(rad_tsw);
308
309 fSlopeSpread = rad_ss <0 ? 0 : TMath::Sqrt(rad_ss);
310 fSlopeSpreadWeighted = rad_ssw<0 ? 0 : TMath::Sqrt(rad_ssw);
311
312 //
313 // Third moments along axes get normalized
314 //
315 const Double_t m3l = c*c*c*sumdx3w + s*s*s*sumdy3w + 3*(s*c*c*sumdx2dyw + c*s*s*sumdxdy2w);
316 const Double_t m3t = c*c*c*sumdy3w - s*s*s*sumdx3w + 3*(s*s*c*sumdx2dyw - s*c*c*sumdxdy2w);
317
318 fM3Long = MMath::Sqrt3(m3l/sumw); // [mm]
319 fM3Trans = MMath::Sqrt3(m3t/sumw); // [mm]
320
321 //
322 // Asymmetry
323 //
324 const MGeom &maxp = geom[maxpixid];
325 fAsym = (hil.GetMeanX()-maxp.GetX())*c + (hil.GetMeanY()-maxp.GetY())*s; // [mm]
326
327 SetReadyToSave();
328
329 return 0;
330}
331
332// --------------------------------------------------------------------------
333//
334// This function is ment for special usage, please never try to set
335// values via this function
336//
337void MHillasExt::Set(const TArrayF &arr)
338{
339 if (arr.GetSize() != 9)
340 return;
341
342 fAsym = arr.At(0);
343 fM3Long = arr.At(1);
344 fM3Trans = arr.At(2);
345 fSlopeLong = arr.At(3);
346 fSlopeTrans = arr.At(4);
347
348 fTimeSpread = arr.At(5);
349 fTimeSpreadWeighted = arr.At(6);
350 fSlopeSpread = arr.At(7);
351 fSlopeSpreadWeighted = arr.At(8);
352}
353
354// -------------------------------------------------------------------------
355//
356// Print contents of MHillasExt to *fLog.
357//
358void MHillasExt::Print(Option_t *o) const
359{
360 const Bool_t showtrans = TString(o).Contains("trans");
361
362 *fLog << all;
363 *fLog << GetDescriptor() << endl;
364 *fLog << " - Asymmetry [mm] = " << fAsym << endl;
365 *fLog << " - 3.Moment Long [mm] = " << fM3Long << endl;
366 if (showtrans)
367 *fLog << " - 3.Moment Trans [mm] = " << fM3Trans << endl;
368 *fLog << " - Slope Long [mm] = " << fSlopeLong << endl;
369 if (showtrans)
370 *fLog << " - Slope Trans [mm] = " << fSlopeTrans << endl;
371 *fLog << " - Time Spread [ns] = " << fTimeSpread << endl;
372 *fLog << " - Slope Spread [ns] = " << fSlopeSpread << endl;
373 *fLog << " - Time Spread W [ns] = " << fTimeSpreadWeighted << endl;
374 *fLog << " - Slope Spread W [ns] = " << fSlopeSpreadWeighted << endl;
375}
376
377// -------------------------------------------------------------------------
378//
379// Print contents of MHillasExt to *fLog, depending on the geometry in
380// units of deg.
381//
382void MHillasExt::Print(const MGeomCam &geom) const
383{
384 *fLog << all;
385 *fLog << GetDescriptor() << endl;
386 *fLog << " - Asymmetry [deg] = " << fAsym *geom.GetConvMm2Deg() << endl;
387 *fLog << " - 3.Moment Long [deg] = " << fM3Long *geom.GetConvMm2Deg() << endl;
388 *fLog << " - 3.Moment Trans [mm] = " << fM3Trans*geom.GetConvMm2Deg() << endl;
389 *fLog << " - Slope Long [mm] = " << fSlopeLong*geom.GetConvMm2Deg() << endl;
390 *fLog << " - Slope Trans [mm] = " << fSlopeTrans*geom.GetConvMm2Deg() << endl;
391 *fLog << " - Slope Spread [ns] = " << fSlopeSpread << endl;
392 *fLog << " - Time Spread W [ns] = " << fTimeSpreadWeighted << endl;
393 *fLog << " - Slope Spread W [ns] = " << fSlopeSpreadWeighted << endl;
394}
395
396// -------------------------------------------------------------------------
397//
398// Paint the 3rd moment on top of the shower. Therefor "MHillas" is needed.
399// it is searched via gPad->FindObject. If opt is not IsNull an object with
400// the name opt is searched.
401//
402void MHillasExt::Paint(Option_t *opt)
403{
404 const TString name(opt);
405
406 const MHillas *hil = dynamic_cast<const MHillas*>(gPad->FindObject(name.IsNull() ? "MHillas" : name.Data()));
407 if (!hil)
408 return;
409
410 TVector2 v(fM3Long, fM3Trans);
411 v = v.Rotate(hil->GetDelta()+TMath::Pi());
412 v += hil->GetMean();
413
414 TMarker m;
415 m.SetMarkerColor(kRed);
416 m.SetMarkerStyle(kFullDotLarge);
417 m.PaintMarker(v.X(), v.Y());
418}
Note: See TracBrowser for help on using the repository browser.