source: trunk/MagicSoft/Mars/mimage/MHillasExt.cc@ 8675

Last change on this file since 8675 was 8619, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 9.1 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-2007
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//
60// WARNING: Before you can use fAsym, fSlope* and fM3* you must
61// multiply by the sign of MHillasSrc::fCosDeltaAlpha
62//
63////////////////////////////////////////////////////////////////////////////
64/*
65 // fAsymna d/(d na) of ( sum(x*q^na)/sum(q^na), sum(y*q^na)/sum(q^na) )
66 // projected onto the major axis
67 // fAsym0 (F-B)/(F+B) along the major axis
68 */
69#include "MHillasExt.h"
70
71#include <TArrayF.h>
72#include <TMarker.h>
73#include <TVector2.h>
74#include <TMinuit.h>
75#include <TVirtualPad.h>
76
77#include "MGeomPix.h"
78#include "MGeomCam.h"
79
80#include "MSignalPix.h"
81#include "MSignalCam.h"
82
83#include "MLog.h"
84#include "MLogManip.h"
85
86#include "MHillas.h"
87
88ClassImp(MHillasExt);
89
90using namespace std;
91
92// -------------------------------------------------------------------------
93//
94// Default constructor.
95//
96MHillasExt::MHillasExt(const char *name, const char *title)
97{
98 fName = name ? name : "MHillasExt";
99 fTitle = title ? title : "Storage container for extended parameter set of one event";
100
101 Reset();
102}
103
104// -------------------------------------------------------------------------
105//
106// Reset all values to 0
107//
108void MHillasExt::Reset()
109{
110 fAsym = 0;
111 fM3Long = 0;
112 fM3Trans = 0;
113 fSlopeLong = 0;
114 fSlopeTrans = 0;
115}
116
117// -------------------------------------------------------------------------
118//
119// calculation of additional parameters based on the camera geometry
120// and the cerenkov photon event
121//
122Int_t MHillasExt::Calc(const MGeomCam &geom, const MSignalCam &evt, const MHillas &hil, Int_t island)
123{
124 //
125 // calculate the additional image parameters
126 // --------------------------------------------
127 //
128 // loop to get third moments along ellipse axes and two max pixels
129 //
130 // For the moments double precision is used to make sure, that
131 // the complex matrix multiplication and sum is evaluated correctly.
132 //
133 Int_t maxpixid = 0;
134 Float_t maxpix = 0;
135
136 // Variables to caluclate time slope
137 // Formula: http://mo.mathematik.uni-stuttgart.de/inhalt/erlaeuterung/erlaeuterung300/
138 UInt_t cnt = 0;
139
140 Double_t sumx = 0;
141 Double_t sumy = 0;
142 Double_t sumt = 0;
143 Double_t sumxy = 0;
144 Double_t sumxt = 0;
145 Double_t sumyt = 0;
146 Double_t sumx2 = 0;
147 Double_t sumy2 = 0;
148 Double_t sumx3 = 0;
149 Double_t sumy3 = 0;
150 Double_t sumx2y = 0;
151 Double_t sumxy2 = 0;
152
153 const UInt_t npix = evt.GetNumPixels();
154 for (UInt_t i=0; i<npix; i++)
155 {
156 const MSignalPix &pix = evt[i];
157 if (!pix.IsPixelUsed())
158 continue;
159
160 if (island>=0 && pix.GetIdxIsland()!=island)
161 continue;
162
163 const MGeomPix &gpix = geom[i];
164
165 const Double_t x = gpix.GetX();
166 const Double_t y = gpix.GetY();
167 const Double_t t = pix.GetArrivalTime();
168
169 // --- time slope ----
170 sumx += x;
171 sumy += y;
172 sumt += t;
173
174 sumx2 += x*x;
175 sumy2 += y*y;
176 sumxy += x*y;
177 sumxt += x*t;
178 sumyt += y*t;
179
180 // --- 3rd moment ---
181 const Double_t dx = x - hil.GetMeanX(); // [mm]
182 const Double_t dy = y - hil.GetMeanY(); // [mm]
183
184 Double_t nphot = pix.GetNumPhotons(); // [1]
185
186 sumx3 += nphot * dx*dx*dx;
187 sumy3 += nphot * dy*dy*dy;
188 sumx2y += nphot * dx*dx*dy;
189 sumx2y += nphot * dx*dy*dy;
190
191 cnt++;
192
193 //
194 // Now we are working on absolute values of nphot, which
195 // must take pixel size into account
196 //
197 nphot *= geom.GetPixRatio(i);
198
199 // --- max pixel ---
200 if (nphot>maxpix)
201 {
202 maxpix = nphot; // [1]
203 maxpixid = i;
204 }
205 }
206
207 const Double_t c = hil.GetCosDelta();
208 const Double_t s = hil.GetSinDelta();
209
210 //
211 // Time slope
212 //
213 const Double_t dxt = c*sumxt + s*sumyt;
214 const Double_t dyt = -s*sumxt + c*sumyt;
215
216 const Double_t dx = c*sumx + s*sumy;
217 const Double_t dy = -s*sumx + c*sumy;
218
219 const Double_t dx2 = c*c*sumx2 + 2*c*s*sumxy + s*s*sumy2;
220 const Double_t dy2 = s*s*sumx2 - 2*c*s*sumxy + c*c*sumy2;
221
222 const Double_t detx = cnt*dx2 - dx*dx;
223 const Double_t dety = cnt*dy2 - dy*dy;
224
225 fSlopeLong = detx==0 ? 0 : (cnt*dxt - sumt*dx)/detx;
226 fSlopeTrans = dety==0 ? 0 : (cnt*dyt - sumt*dy)/dety;
227
228 //
229 // Third moments along axes get normalized
230 //
231 const Double_t m3l = c*c*c*sumx3 + 3*(s*c*c*sumx2y + c*s*s*sumxy2) + s*s*s*sumy3;
232 const Double_t m3t = c*c*c*sumy3 + 3*(s*s*c*sumx2y - s*c*c*sumxy2) - s*s*s*sumx3;
233
234 fM3Long = MMath::Sqrt3(m3l/hil.GetSize()); // [mm]
235 fM3Trans = MMath::Sqrt3(m3t/hil.GetSize()); // [mm]
236
237 //
238 // Asymmetry
239 //
240 const MGeomPix &maxp = geom[maxpixid];
241 fAsym = (hil.GetMeanX()-maxp.GetX())*c + (hil.GetMeanY()-maxp.GetY())*s; // [mm]
242
243 SetReadyToSave();
244
245 return 0;
246}
247
248// --------------------------------------------------------------------------
249//
250// This function is ment for special usage, please never try to set
251// values via this function
252//
253void MHillasExt::Set(const TArrayF &arr)
254{
255 if (arr.GetSize() != 5)
256 return;
257
258 fAsym = arr.At(0);
259 fM3Long = arr.At(1);
260 fM3Trans = arr.At(2);
261 fSlopeLong = arr.At(3);
262 fSlopeTrans = arr.At(4);
263}
264
265// -------------------------------------------------------------------------
266//
267// Print contents of MHillasExt to *fLog.
268//
269void MHillasExt::Print(Option_t *o) const
270{
271 const Bool_t showtrans = TString(o).Contains("trans");
272
273 *fLog << all;
274 *fLog << GetDescriptor() << endl;
275 *fLog << " - Asymmetry [mm] = " << fAsym << endl;
276 *fLog << " - 3.Moment Long [mm] = " << fM3Long << endl;
277 if (showtrans)
278 *fLog << " - 3.Moment Trans [mm] = " << fM3Trans << endl;
279 *fLog << " - Slope Long [mm] = " << fSlopeLong << endl;
280 if (showtrans)
281 *fLog << " - Slope Trans [mm] = " << fSlopeTrans << endl;
282}
283
284// -------------------------------------------------------------------------
285//
286// Print contents of MHillasExt to *fLog, depending on the geometry in
287// units of deg.
288//
289void MHillasExt::Print(const MGeomCam &geom) const
290{
291 *fLog << all;
292 *fLog << GetDescriptor() << endl;
293 *fLog << " - Asymmetry [deg] = " << fAsym *geom.GetConvMm2Deg() << endl;
294 *fLog << " - 3.Moment Long [deg] = " << fM3Long *geom.GetConvMm2Deg() << endl;
295 *fLog << " - 3.Moment Trans [mm] = " << fM3Trans*geom.GetConvMm2Deg() << endl;
296 *fLog << " - Slope Long [mm] = " << fSlopeLong*geom.GetConvMm2Deg() << endl;
297 *fLog << " - Slope Trans [mm] = " << fSlopeTrans*geom.GetConvMm2Deg() << endl;
298}
299
300// -------------------------------------------------------------------------
301//
302// Paint the 3rd moment on top of the shower. Therefor "MHillas" is needed.
303// it is searched via gPad->FindObject. If opt is not IsNull an object with
304// the name opt is searched.
305//
306void MHillasExt::Paint(Option_t *opt)
307{
308 const TString name(opt);
309
310 const MHillas *hil = dynamic_cast<const MHillas*>(gPad->FindObject(name.IsNull() ? "MHillas" : name.Data()));
311 if (!hil)
312 return;
313
314 TVector2 v(fM3Long, fM3Trans);
315 v = v.Rotate(hil->GetDelta()+TMath::Pi());
316 v += hil->GetMean();
317
318 TMarker m;
319 m.SetMarkerColor(kRed);
320 m.SetMarkerStyle(kFullDotLarge);
321 m.PaintMarker(v.X(), v.Y());
322}
Note: See TracBrowser for help on using the repository browser.