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

Last change on this file since 9849 was 9849, checked in by tbretz, 14 years ago
Fixed MHillasExt::fTimeSpread and implemented calculation of fSlopeSpread.
File size: 11.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
125// -------------------------------------------------------------------------
126//
127// calculation of additional parameters based on the camera geometry
128// and the cerenkov photon event
129//
130Int_t MHillasExt::Calc(const MGeomCam &geom, const MSignalCam &evt, const MHillas &hil, Int_t island)
131{
132 //
133 // calculate the additional image parameters
134 // --------------------------------------------
135 //
136 // loop to get third moments along ellipse axes and two max pixels
137 //
138 // For the moments double precision is used to make sure, that
139 // the complex matrix multiplication and sum is evaluated correctly.
140 //
141 Int_t maxpixid = 0;
142 Float_t maxpix = 0;
143
144 // Variables to caluclate time slope
145 // Formula: http://mo.mathematik.uni-stuttgart.de/inhalt/erlaeuterung/erlaeuterung300/
146 UInt_t cnt = 0;
147
148 Double_t sumx = 0;
149 Double_t sumy = 0;
150 Double_t sumt = 0;
151 Double_t sumw = 0;
152 Double_t sumxy = 0;
153 Double_t sumxt = 0;
154 Double_t sumyt = 0;
155 Double_t sumxw = 0;
156 Double_t sumyw = 0;
157 Double_t sumtw = 0;
158 Double_t sumx2 = 0;
159 Double_t sumy2 = 0;
160 Double_t sumw2 = 0;
161 Double_t sumt2 = 0;
162 Double_t sumx3 = 0;
163 Double_t sumy3 = 0;
164 Double_t sumxw2 = 0;
165 Double_t sumyw2 = 0;
166 Double_t sumtx2 = 0;
167 Double_t sumty2 = 0;
168 Double_t sumtw2 = 0;
169 Double_t sumx2y = 0;
170 Double_t sumxy2 = 0;
171 Double_t sumtxw2 = 0;
172 Double_t sumtyw2 = 0;
173 Double_t sumt2w2 = 0;
174
175 const UInt_t npix = evt.GetNumPixels();
176 for (UInt_t i=0; i<npix; i++)
177 {
178 const MSignalPix &pix = evt[i];
179 if (!pix.IsPixelUsed())
180 continue;
181
182 if (island>=0 && pix.GetIdxIsland()!=island)
183 continue;
184
185 const MGeom &gpix = geom[i];
186
187 const Double_t x = gpix.GetX();
188 const Double_t y = gpix.GetY();
189 const Double_t t = pix.GetArrivalTime();
190
191 Double_t nphot = pix.GetNumPhotons(); // [1]
192
193 // --- time slope ----
194 sumx += x;
195 sumy += y;
196
197 sumt += t;
198 sumt2 += t*t;
199
200 sumw += nphot;
201 sumw2 += nphot*nphot;
202
203 sumxw += x*nphot;
204 sumyw += x*nphot;
205 sumtw += t*nphot;
206
207 sumtx2 += x*nphot*nphot;
208 sumty2 += y*nphot*nphot;
209 sumtw2 += t*nphot*nphot;
210
211 sumtxw2 += t*x*nphot*nphot;
212 sumtyw2 += t*y*nphot*nphot;
213 sumt2w2 += t*t*nphot*nphot;
214
215 sumx2 += x*x;
216 sumy2 += y*y;
217 sumxy += x*y;
218 sumxt += x*t;
219 sumyt += y*t;
220
221 // --- 3rd moment ---
222 const Double_t dx = x - hil.GetMeanX(); // [mm]
223 const Double_t dy = y - hil.GetMeanY(); // [mm]
224
225 sumx3 += nphot * dx*dx*dx;
226 sumy3 += nphot * dy*dy*dy;
227 sumx2y += nphot * dx*dx*dy;
228 sumx2y += nphot * dx*dy*dy;
229
230 cnt++;
231
232 //
233 // Now we are working on absolute values of nphot, which
234 // must take pixel size into account
235 //
236 nphot *= geom.GetPixRatio(i);
237
238 // --- max pixel ---
239 if (nphot>maxpix)
240 {
241 maxpix = nphot; // [1]
242 maxpixid = i;
243 }
244 }
245
246 const Double_t c = hil.GetCosDelta();
247 const Double_t s = hil.GetSinDelta();
248
249 //
250 // Time spread
251 //
252 fTimeSpread = TMath::Sqrt(sumt2/cnt - sumt*sumt/cnt/cnt);
253 fTimeSpreadWeighted = TMath::Sqrt(sumt2w2 + (sumw2 - 2*sumtw2)*sumtw/sumw)/sumw;
254
255 //
256 // Time slope
257 //
258 const Double_t dxt = c*sumxt + s*sumyt;
259 const Double_t dyt = -s*sumxt + c*sumyt;
260
261 const Double_t dx = c*sumx + s*sumy;
262 const Double_t dy = -s*sumx + c*sumy;
263
264 const Double_t dx2 = c*c*sumx2 + 2*c*s*sumxy + s*s*sumy2;
265 const Double_t dy2 = s*s*sumx2 - 2*c*s*sumxy + c*c*sumy2;
266
267 const Double_t detx = cnt*dx2 - dx*dx;
268 const Double_t dety = cnt*dy2 - dy*dy;
269
270 fSlopeLong = detx==0 ? 0 : (cnt*dxt - sumt*dx)/detx;
271 fSlopeTrans = dety==0 ? 0 : (cnt*dyt - sumt*dy)/dety;
272
273 //
274 // Calculate time spread around longitudinal time slope
275 //
276 const Double_t sumdt = sumt - fSlopeLong*(c*sumx + s*sumy);
277 const Double_t sumdtw = sumtw - fSlopeLong*(c*sumxw + s*sumyw);
278 const Double_t sumdtw2 = sumtw2 - fSlopeLong*(c*sumxw2 + s*sumyw2);
279 const Double_t sumdt2 = sumt2 - 2*fSlopeLong*(c*sumxt + s*sumyt) + fSlopeLong*fSlopeLong * (c*c*sumx2 + s*s*sumy2 + 2*c*s*sumxy);
280 const Double_t sumdt2w2 = sumt2w2 - 2*fSlopeLong*(c*sumtxw2 + s*sumtyw2) + fSlopeLong*fSlopeLong * (c*c*sumx2 + s*s*sumy2 + 2*c*s*sumxy) * sumw2;
281
282 fSlopeSpread = TMath::Sqrt(sumdt2/cnt - sumdt*sumdt/cnt/cnt);
283 fSlopeSpreadWeighted = TMath::Sqrt(sumdt2w2 + (sumw2 - 2*sumdtw2)*sumdtw/sumw)/sumw;
284
285 //
286 // FIXME: Missing Time spread along slope
287 //
288
289 //
290 // Third moments along axes get normalized
291 //
292 const Double_t m3l = c*c*c*sumx3 + 3*(s*c*c*sumx2y + c*s*s*sumxy2) + s*s*s*sumy3;
293 const Double_t m3t = c*c*c*sumy3 + 3*(s*s*c*sumx2y - s*c*c*sumxy2) - s*s*s*sumx3;
294
295 fM3Long = MMath::Sqrt3(m3l/hil.GetSize()); // [mm]
296 fM3Trans = MMath::Sqrt3(m3t/hil.GetSize()); // [mm]
297
298 //
299 // Asymmetry
300 //
301 const MGeom &maxp = geom[maxpixid];
302 fAsym = (hil.GetMeanX()-maxp.GetX())*c + (hil.GetMeanY()-maxp.GetY())*s; // [mm]
303
304 SetReadyToSave();
305
306 return 0;
307}
308
309// --------------------------------------------------------------------------
310//
311// This function is ment for special usage, please never try to set
312// values via this function
313//
314void MHillasExt::Set(const TArrayF &arr)
315{
316 if (arr.GetSize() != 5)
317 return;
318
319 fAsym = arr.At(0);
320 fM3Long = arr.At(1);
321 fM3Trans = arr.At(2);
322 fSlopeLong = arr.At(3);
323 fSlopeTrans = arr.At(4);
324}
325
326// -------------------------------------------------------------------------
327//
328// Print contents of MHillasExt to *fLog.
329//
330void MHillasExt::Print(Option_t *o) const
331{
332 const Bool_t showtrans = TString(o).Contains("trans");
333
334 *fLog << all;
335 *fLog << GetDescriptor() << endl;
336 *fLog << " - Asymmetry [mm] = " << fAsym << endl;
337 *fLog << " - 3.Moment Long [mm] = " << fM3Long << endl;
338 if (showtrans)
339 *fLog << " - 3.Moment Trans [mm] = " << fM3Trans << endl;
340 *fLog << " - Slope Long [mm] = " << fSlopeLong << endl;
341 if (showtrans)
342 *fLog << " - Slope Trans [mm] = " << fSlopeTrans << endl;
343}
344
345// -------------------------------------------------------------------------
346//
347// Print contents of MHillasExt to *fLog, depending on the geometry in
348// units of deg.
349//
350void MHillasExt::Print(const MGeomCam &geom) const
351{
352 *fLog << all;
353 *fLog << GetDescriptor() << endl;
354 *fLog << " - Asymmetry [deg] = " << fAsym *geom.GetConvMm2Deg() << endl;
355 *fLog << " - 3.Moment Long [deg] = " << fM3Long *geom.GetConvMm2Deg() << endl;
356 *fLog << " - 3.Moment Trans [mm] = " << fM3Trans*geom.GetConvMm2Deg() << endl;
357 *fLog << " - Slope Long [mm] = " << fSlopeLong*geom.GetConvMm2Deg() << endl;
358 *fLog << " - Slope Trans [mm] = " << fSlopeTrans*geom.GetConvMm2Deg() << endl;
359}
360
361// -------------------------------------------------------------------------
362//
363// Paint the 3rd moment on top of the shower. Therefor "MHillas" is needed.
364// it is searched via gPad->FindObject. If opt is not IsNull an object with
365// the name opt is searched.
366//
367void MHillasExt::Paint(Option_t *opt)
368{
369 const TString name(opt);
370
371 const MHillas *hil = dynamic_cast<const MHillas*>(gPad->FindObject(name.IsNull() ? "MHillas" : name.Data()));
372 if (!hil)
373 return;
374
375 TVector2 v(fM3Long, fM3Trans);
376 v = v.Rotate(hil->GetDelta()+TMath::Pi());
377 v += hil->GetMean();
378
379 TMarker m;
380 m.SetMarkerColor(kRed);
381 m.SetMarkerStyle(kFullDotLarge);
382 m.PaintMarker(v.X(), v.Y());
383}
Note: See TracBrowser for help on using the repository browser.