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

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