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

Last change on this file since 13736 was 9860, checked in by tbretz, 15 years ago
Use the definition of the biased weighted sample variance to calculate the weighted time spread. This includes the weights (like in the third moment) linearly, not squared as before.
File size: 11.8 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 const MGeom &gpix = geom[i];
194
195 const Double_t x = gpix.GetX();
196 const Double_t y = gpix.GetY();
197 const Double_t t = pix.GetArrivalTime();
198
199 Double_t nphot = pix.GetNumPhotons(); // [1]
200
201 // --- time slope ----
202 sumx += x;
203 sumy += y;
204 sumt += t;
205 sumw += nphot;
206
207 sumxt += x*t;
208 sumyt += y*t;
209
210 sumx2 += x*x;
211 sumy2 += y*y;
212 sumxy += x*y;
213 sumt2 += t*t;
214
215 sumxw += x*nphot;
216 sumyw += y*nphot;
217 sumtw += t*nphot;
218
219 sumx2w += x*x*nphot;
220 sumy2w += y*y*nphot;
221 sumt2w += t*t*nphot;
222
223 sumxyw += x*y*nphot;
224
225 sumxtw += x*t*nphot;
226 sumytw += y*t*nphot;
227
228 // --- 3rd moment ---
229 const Double_t dx = x - hil.GetMeanX(); // [mm]
230 const Double_t dy = y - hil.GetMeanY(); // [mm]
231
232 sumdx3w += dx*dx*dx*nphot;
233 sumdy3w += dy*dy*dy*nphot;
234 sumdx2dyw += dx*dx*dy*nphot;
235 sumdxdy2w += dx*dy*dy*nphot;
236
237 cnt++;
238
239 //
240 // Now we are working on absolute values of nphot, which
241 // must take pixel size into account
242 //
243 nphot *= geom.GetPixRatio(i);
244
245 // --- max pixel ---
246 if (nphot>maxpix)
247 {
248 maxpix = nphot; // [1]
249 maxpixid = i;
250 }
251 }
252
253 const Double_t c = hil.GetCosDelta();
254 const Double_t s = hil.GetSinDelta();
255
256 //
257 // Time slope
258 //
259 const Double_t dxt = c*sumxt + s*sumyt;
260 const Double_t dyt = -s*sumxt + c*sumyt;
261
262 const Double_t dx = c*sumx + s*sumy;
263 const Double_t dy = -s*sumx + c*sumy;
264
265 const Double_t dx2 = c*c*sumx2 + 2*c*s*sumxy + s*s*sumy2;
266 const Double_t dy2 = s*s*sumx2 - 2*c*s*sumxy + c*c*sumy2;
267
268 const Double_t detx = cnt*dx2 - dx*dx;
269 const Double_t dety = cnt*dy2 - dy*dy;
270
271 fSlopeLong = detx==0 ? 0 : (cnt*dxt - sumt*dx)/detx;
272 fSlopeTrans = dety==0 ? 0 : (cnt*dyt - sumt*dy)/dety;
273
274 //
275 // Time spread
276 //
277 const Double_t m = fSlopeLong;
278
279 const Double_t z = dx; //c*sumx + s*sumy;
280 const Double_t zw = c*sumxw + s*sumyw;
281
282 const Double_t zt = dxt; //c*sumxt + s*sumyt;
283 const Double_t ztw = c*sumxtw + s*sumytw;
284
285 const Double_t x2y2 = dx2; //c*c*sumx2 + s*s*sumy2 + 2*c*s*sumxy;
286 const Double_t x2y2w = c*c*sumx2w + s*s*sumy2w + 2*c*s*sumxyw;
287
288 const Double_t sumdt = sumt - m*z;
289 const Double_t sumdtw = sumtw - m*zw;
290
291 const Double_t sumdt2 = sumt2 + m*m*x2y2 - 2*m*zt;
292 const Double_t sumdt2w = sumt2w + m*m*x2y2w - 2*m*ztw;
293
294 const Double_t meant = sumt /cnt;
295 const Double_t meandt = sumdt /cnt;
296 const Double_t meantw = sumtw /sumw;
297 const Double_t meandtw = sumdtw/sumw;
298
299 fTimeSpread = TMath::Sqrt(sumt2 /cnt - meant *meant);
300 fTimeSpreadWeighted = TMath::Sqrt(sumt2w/sumw - meantw*meantw);
301
302 fSlopeSpread = TMath::Sqrt(sumdt2 /cnt - meandt *meandt);
303 fSlopeSpreadWeighted = TMath::Sqrt(sumdt2w/sumw - meandtw*meandtw);
304
305 //
306 // Third moments along axes get normalized
307 //
308 const Double_t m3l = c*c*c*sumdx3w + s*s*s*sumdy3w + 3*(s*c*c*sumdx2dyw + c*s*s*sumdxdy2w);
309 const Double_t m3t = c*c*c*sumdy3w - s*s*s*sumdx3w + 3*(s*s*c*sumdx2dyw - s*c*c*sumdxdy2w);
310
311 fM3Long = MMath::Sqrt3(m3l/sumw); // [mm]
312 fM3Trans = MMath::Sqrt3(m3t/sumw); // [mm]
313
314 //
315 // Asymmetry
316 //
317 const MGeom &maxp = geom[maxpixid];
318 fAsym = (hil.GetMeanX()-maxp.GetX())*c + (hil.GetMeanY()-maxp.GetY())*s; // [mm]
319
320 SetReadyToSave();
321
322 return 0;
323}
324
325// --------------------------------------------------------------------------
326//
327// This function is ment for special usage, please never try to set
328// values via this function
329//
330void MHillasExt::Set(const TArrayF &arr)
331{
332 if (arr.GetSize() != 9)
333 return;
334
335 fAsym = arr.At(0);
336 fM3Long = arr.At(1);
337 fM3Trans = arr.At(2);
338 fSlopeLong = arr.At(3);
339 fSlopeTrans = arr.At(4);
340
341 fTimeSpread = arr.At(5);
342 fTimeSpreadWeighted = arr.At(6);
343 fSlopeSpread = arr.At(7);
344 fSlopeSpreadWeighted = arr.At(8);
345}
346
347// -------------------------------------------------------------------------
348//
349// Print contents of MHillasExt to *fLog.
350//
351void MHillasExt::Print(Option_t *o) const
352{
353 const Bool_t showtrans = TString(o).Contains("trans");
354
355 *fLog << all;
356 *fLog << GetDescriptor() << endl;
357 *fLog << " - Asymmetry [mm] = " << fAsym << endl;
358 *fLog << " - 3.Moment Long [mm] = " << fM3Long << endl;
359 if (showtrans)
360 *fLog << " - 3.Moment Trans [mm] = " << fM3Trans << endl;
361 *fLog << " - Slope Long [mm] = " << fSlopeLong << endl;
362 if (showtrans)
363 *fLog << " - Slope Trans [mm] = " << fSlopeTrans << endl;
364 *fLog << " - Time Spread [ns] = " << fTimeSpread << endl;
365 *fLog << " - Slope Spread [ns] = " << fSlopeSpread << endl;
366 *fLog << " - Time Spread W [ns] = " << fTimeSpreadWeighted << endl;
367 *fLog << " - Slope Spread W [ns] = " << fSlopeSpreadWeighted << endl;
368}
369
370// -------------------------------------------------------------------------
371//
372// Print contents of MHillasExt to *fLog, depending on the geometry in
373// units of deg.
374//
375void MHillasExt::Print(const MGeomCam &geom) const
376{
377 *fLog << all;
378 *fLog << GetDescriptor() << endl;
379 *fLog << " - Asymmetry [deg] = " << fAsym *geom.GetConvMm2Deg() << endl;
380 *fLog << " - 3.Moment Long [deg] = " << fM3Long *geom.GetConvMm2Deg() << endl;
381 *fLog << " - 3.Moment Trans [mm] = " << fM3Trans*geom.GetConvMm2Deg() << endl;
382 *fLog << " - Slope Long [mm] = " << fSlopeLong*geom.GetConvMm2Deg() << endl;
383 *fLog << " - Slope Trans [mm] = " << fSlopeTrans*geom.GetConvMm2Deg() << endl;
384 *fLog << " - Slope Spread [ns] = " << fSlopeSpread << endl;
385 *fLog << " - Time Spread W [ns] = " << fTimeSpreadWeighted << endl;
386 *fLog << " - Slope Spread W [ns] = " << fSlopeSpreadWeighted << endl;
387}
388
389// -------------------------------------------------------------------------
390//
391// Paint the 3rd moment on top of the shower. Therefor "MHillas" is needed.
392// it is searched via gPad->FindObject. If opt is not IsNull an object with
393// the name opt is searched.
394//
395void MHillasExt::Paint(Option_t *opt)
396{
397 const TString name(opt);
398
399 const MHillas *hil = dynamic_cast<const MHillas*>(gPad->FindObject(name.IsNull() ? "MHillas" : name.Data()));
400 if (!hil)
401 return;
402
403 TVector2 v(fM3Long, fM3Trans);
404 v = v.Rotate(hil->GetDelta()+TMath::Pi());
405 v += hil->GetMean();
406
407 TMarker m;
408 m.SetMarkerColor(kRed);
409 m.SetMarkerStyle(kFullDotLarge);
410 m.PaintMarker(v.X(), v.Y());
411}
Note: See TracBrowser for help on using the repository browser.