source: trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc@ 9435

Last change on this file since 9435 was 9376, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 9.5 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): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de>
19! Author(s): Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2008
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MMuonSearchPar
29//
30// Storage Container for muon
31//
32// This class is the container for muon parameters. The calculation
33// of Radius and center of the ring is done here.
34// Muons are searched by fitting the image with a circle.
35//
36// In order to use further information of muons such as the width of arcs,
37// the size of the fraction of the ring and the muons size, use the
38// infomation stored in
39//
40// MMuonCalibPar.
41//
42// The information will be available by using the task of
43//
44// MMuonCalibParCalc.
45//
46// Version 2:
47// ----------
48// + Float_t fTime; // Mean arrival time of core pixels
49// + Float_t fTimeRms; // Rms of arrival time of core pixels
50//
51// Input Containers:
52// MGeomCam
53// MHillas
54// MSignalCam
55//
56/////////////////////////////////////////////////////////////////////////////
57#include "MMuonSearchPar.h"
58
59#include <TMinuit.h>
60#include <TEllipse.h>
61
62#include "MLog.h"
63#include "MLogManip.h"
64
65#include "MHillas.h"
66
67#include "MGeomCam.h"
68#include "MGeomPix.h"
69
70#include "MSignalPix.h"
71#include "MSignalCam.h"
72
73using namespace std;
74
75ClassImp(MMuonSearchPar);
76
77// --------------------------------------------------------------------------
78//
79// Default constructor.
80//
81MMuonSearchPar::MMuonSearchPar(const char *name, const char *title)
82{
83 fName = name ? name : "MMuonSearchPar";
84 fTitle = title ? title : "Parameters to find Muons";
85}
86
87// --------------------------------------------------------------------------
88//
89void MMuonSearchPar::Reset()
90{
91 fRadius = -1;
92 fDeviation = -1;
93 fCenterX = 0;
94 fCenterY = 0;
95 fTime = 0;
96 fTimeRms = -1;
97}
98
99// --------------------------------------------------------------------------
100//
101// return TMath::Hypot(fCenterX, fCenterY);
102//
103Float_t MMuonSearchPar::GetDist() const
104{
105 return TMath::Hypot(fCenterX, fCenterY);
106}
107
108// --------------------------------------------------------------------------
109//
110// This is a wrapper function to have direct access to the data members
111// in the function calculating the minimization value.
112//
113void MMuonSearchPar::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
114{
115 const MMuonSearchPar *optim = (MMuonSearchPar*)gMinuit->GetObjectFit();
116 f = optim->Fcn(par);
117}
118
119// --------------------------------------------------------------------------
120//
121// This function gives you the ring radius fitted best to the camera image
122// and its RMS for the input position.
123//
124Double_t MMuonSearchPar::Fcn(Double_t *par) const
125{
126 const Int_t entries = fSignal.GetSize();
127
128 Double_t meanr=0;
129 Double_t devr =0;
130 Double_t sums =0;
131
132 // It seems that the loop is easy enough for a compiler optimization.
133 // Using pointer arithmetics doesn't improve the speed of the fit.
134 for (Int_t i=0; i<entries; i++ )
135 {
136 const Double_t dx = fX[i]-par[0];
137 const Double_t dy = fY[i]-par[1];
138
139 const Double_t sq = dx*dx + dy*dy;
140
141 sums += fSignal[i];
142 meanr += fSignal[i] * TMath::Sqrt(sq);
143 devr += fSignal[i] * sq;
144 }
145
146 par[2] = meanr/sums;
147 par[3] = devr/sums - par[2]*par[2];
148
149 return par[3];
150}
151
152// --------------------------------------------------------------------------
153//
154// This function finds the center position of the circle which gives minimum
155// RMS of the fit, changing the center position of the circle.
156//
157void MMuonSearchPar::CalcMinimumDeviation(const MGeomCam &geom,
158 const MSignalCam &evt,
159 Double_t &x, Double_t &y,
160 Double_t &sigma, Double_t &radius)
161{
162 // ------- Make a temporaray copy of the signal ---------
163 // ------- and calculate arrival time parameters --------
164 const Int_t n = geom.GetNumPixels();
165
166 fSignal.Set(n);
167 fX.Set(n);
168 fY.Set(n);
169
170 Int_t p=0;
171 Int_t q=0;
172
173 Double_t mean=0;
174 Double_t sq =0;
175
176 for (int i=0; i<n; i++)
177 {
178 const MSignalPix &pix = evt[i];
179 if (!pix.IsPixelUsed())
180 continue;
181
182 fSignal[p] = pix.GetNumPhotons();
183
184 fX[p] = geom[i].GetX();
185 fY[p] = geom[i].GetY();
186 p++;
187
188 //timing
189 if (!pix.IsPixelCore())
190 continue;
191
192 mean += pix.GetArrivalTime();
193 sq += pix.GetArrivalTime()*pix.GetArrivalTime();
194 q++;
195 }
196
197 mean /= q;
198 sq /= q;
199
200 fTime = mean;
201 fTimeRms = TMath::Sqrt(sq-mean*mean);
202
203 fSignal.Set(p);
204
205
206 // ----------------- Setup and call minuit -------------------
207 const Float_t delta = 30.; // 3 mm (1/10 of an inner pixel size) Step to move.
208 //const Double_t r = geom.GetMaxRadius()*2;
209
210 // Save gMinuit
211 TMinuit *minsave = gMinuit;
212
213 // Initialize TMinuit with 4 parameters
214 TMinuit minuit;
215 minuit.SetPrintLevel(-1); // Switch off printing
216 minuit.Command("set nowarn"); // Switch off warning
217 // Define Parameters
218 minuit.DefineParameter(0, "x", x, delta, 0, 0);//-r, r);
219 minuit.DefineParameter(1, "y", y, delta, 0, 0);//-r, r);
220 minuit.DefineParameter(2, "r", 0, 1, 0, 0);
221 minuit.DefineParameter(3, "sigma", 0, 1, 0, 0);
222 // Fix return parameters
223 minuit.FixParameter(2);
224 minuit.FixParameter(3);
225 // Setup Minuit for 'this' and use fit function fcn
226 minuit.SetObjectFit(this);
227 minuit.SetFCN(fcn);
228
229 // Perform Simplex minimization
230 Int_t er;
231 Double_t tmp[2] = { 0, 0 };
232 minuit.mnexcm("simplex", tmp, 2, er);
233
234 // Get resulting parameters
235 Double_t error;
236 minuit.GetParameter(0, x, error);
237 minuit.GetParameter(1, y, error);
238 minuit.GetParameter(2, radius, error);
239 minuit.GetParameter(3, sigma, error);
240
241 sigma = sigma>0 ? TMath::Sqrt(sigma) : 0;
242
243 gMinuit = minsave;
244}
245
246// --------------------------------------------------------------------------
247//
248// Calculation of muon parameters
249//
250void MMuonSearchPar::Calc(const MGeomCam &geom, const MSignalCam &evt,
251 const MHillas &hillas)
252{
253 Double_t x = hillas.GetMeanX();
254 Double_t y = hillas.GetMeanY();
255
256 // -------------------------------------------------
257 // Keiichi suggested trying to precalculate the Muon
258 // center a bit better, but it neither improves the
259 // fit result nor the speed
260 //
261 // const Float_t tmpr = 300.; // assume that the temporal cherenkov angle is 1 deg. (300 mm).
262 //
263 // const Double_t a = TMath::Tan(hillas.GetDelta());
264 //
265 // const Double_t dx = a/TMath::Sqrt(tmpr+a*a)/3.;
266 // const Double_t dy = -tmpr/TMath::Sqrt(1+a*a)/3.;
267 //
268 // Double_t par1[] = { x+dx, y+dy, 0, 0 };
269 // Double_t par2[] = { x-dx, y-dy, 0, 0 };
270 //
271 // const Double_t dev1 = MMuonSearchPar::Fcn(par1);
272 // const Double_t dev2 = MMuonSearchPar::Fcn(par2);
273 //
274 // if (dev1<dev2)
275 // {
276 // x += dx;
277 // y += dy;
278 // }
279 // else
280 // {
281 // x -= dx;
282 // y -= dy;
283 // }
284 // -------------------------------------------------
285
286 Double_t sigma, rad;
287
288 // find the best fit.
289 CalcMinimumDeviation(geom, evt, x, y, sigma, rad);
290
291 fCenterX = x;
292 fCenterY = y;
293 fRadius = rad;
294 fDeviation = sigma;
295
296 //SetReadyToSave();
297}
298
299void MMuonSearchPar::Print(Option_t *) const
300{
301 *fLog << all;
302 *fLog << GetDescriptor() << endl;
303 *fLog << " - Est. Radius [mm] = " << fRadius << endl;
304 *fLog << " - Deviation [mm] = " << fDeviation << endl;
305 *fLog << " - Center Pos. X [mm] = " << fCenterX << endl;
306 *fLog << " - Center Pos. Y [mm] = " << fCenterY << endl;
307}
308
309void MMuonSearchPar::Print(const MGeomCam &geom, Option_t *) const
310{
311 *fLog << all;
312 *fLog << GetDescriptor() << endl;
313 *fLog << " - Est. Radius [deg] = " << fRadius*geom.GetConvMm2Deg() << endl;
314 *fLog << " - Deviation [deg] = " << fDeviation*geom.GetConvMm2Deg() << endl;
315 *fLog << " - Center Pos. X [deg] = " << fCenterX*geom.GetConvMm2Deg() << endl;
316 *fLog << " - Center Pos. Y [deg] = " << fCenterY*geom.GetConvMm2Deg() << endl;
317}
318
319// --------------------------------------------------------------------------
320//
321// Paint the ellipse corresponding to the parameters
322//
323void MMuonSearchPar::Paint(Option_t *opt)
324{
325 //if (fRadius<180 || fRadius>400 || fDeviation>45)
326 // return;
327
328 TEllipse e1(fCenterX, fCenterY, fRadius-fDeviation, fRadius-fDeviation);
329 TEllipse e2(fCenterX, fCenterY, fRadius+fDeviation, fRadius+fDeviation);
330 e1.SetLineWidth(1);
331 e2.SetLineWidth(1);
332 e1.SetLineColor(kYellow);
333 e2.SetLineColor(kYellow);
334 e1.Paint();
335 e2.Paint();
336}
Note: See TracBrowser for help on using the repository browser.