source: trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.cc@ 4800

Last change on this file since 4800 was 4410, checked in by jlopez, 20 years ago
*** empty log message ***
File size: 10.4 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): Roger Firpo 04/2004 <mailto:jlopez@ifae.es>
19! Author(s): Javier López 05/2004 <mailto:jlopez@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MSrcPosFromStars
29//
30// Task to set the position of the source using the positions of the stars
31// in the field of view of the source
32//
33// Output Containers:
34// MSrcPosCam
35//
36//////////////////////////////////////////////////////////////////////////////
37
38
39#include "MSrcPosFromStars.h"
40
41#include <TList.h>
42#include <TH2F.h>
43#include <TF2.h>
44#include <TTimer.h>
45#include <TString.h>
46#include <TCanvas.h>
47
48#include "MSrcPosCam.h"
49#include "MStarLocalCam.h"
50#include "MStarLocalPos.h"
51
52#include "MLog.h"
53#include "MLogManip.h"
54
55#include "MParList.h"
56
57ClassImp(MSrcPosFromStars);
58
59using namespace std;
60
61static const TString gsDefName = "MSrcPosFromStars";
62static const TString gsDefTitle = "task to calculate the position of the source using the position of stars";
63
64// -------------------------------------------------------------------------
65//
66// Default constructor.
67//
68MSrcPosFromStars::MSrcPosFromStars(const char *name, const char *title)
69 : fStars(NULL)
70{
71 fName = name ? name : gsDefName.Data();
72 fTitle = title ? title : gsDefTitle.Data();
73
74 fDistances.Set(0);
75
76}
77
78// -------------------------------------------------------------------------
79//
80Int_t MSrcPosFromStars::PreProcess(MParList *pList)
81{
82
83 if(!MSrcPlace::PreProcess(pList))
84 return kFALSE;
85
86 fStars = (MStarLocalCam*)pList->FindObject(AddSerialNumber("MStarLocalCam"));
87 if (!fStars)
88 {
89 *fLog << err << AddSerialNumber("MStarLocalCam") << " not found ... aborting" << endl;
90 return kFALSE;
91 }
92
93 fNumStars = fDistances.GetSize();
94
95 return kTRUE;
96}
97
98// --------------------------------------------------------------------------
99//
100//
101static Bool_t HandleInput()
102{
103 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
104 while (1)
105 {
106 //
107 // While reading the input process gui events asynchronously
108 //
109 timer.TurnOn();
110 cout << "Type 'q' to exit, <return> to go on: " << endl;
111 TString input;
112 input = getchar();
113 timer.TurnOff();
114
115 if (input=="q\n")
116 return kFALSE;
117
118 if (input=="\n")
119 return kTRUE;
120 };
121
122 return kFALSE;
123}
124
125Int_t MSrcPosFromStars::ComputeNewSrcPosition()
126{
127
128 if (fNumStars == 0)
129 {
130 if (fStars->GetNumStars() > 0)
131 {
132
133 //Look for the star closer to the center of the camera
134 TIter Next(fStars->GetList());
135 MStarLocalPos* star;
136 Float_t mindist = 600; //mm
137 UInt_t starnum = 0;
138 Int_t select = -1;
139 while ((star=(MStarLocalPos*)Next()))
140 {
141 Float_t dist = TMath::Sqrt(star->GetMeanX()*star->GetMeanX() +
142 star->GetMeanY()*star->GetMeanY());
143 if (dist < mindist)
144 {
145 mindist = dist;
146 select = starnum;
147 }
148
149 starnum++;
150 }
151
152 if (select < 0)
153 {
154 *fLog << err << "Not found star closer to center" << endl;
155 return kFALSE;
156 }
157
158 MStarLocalPos& selecStar = (*fStars)[select];
159
160 if (selecStar.GetChiSquareNdof() > 0. && selecStar.GetChiSquareNdof() < 10.)
161 {
162
163 Float_t selecStarPosX = selecStar.GetMeanX();
164 Float_t selecStarPosY = selecStar.GetMeanY();
165
166 GetOutputSrcPosCam()->SetXY(selecStarPosX,selecStarPosY);
167 }
168 }
169 }
170 else if (fStars->GetNumStars() >= fNumStars)
171 {
172
173 Float_t diameterInnerPixel = 30; //[mm]
174 Float_t diameterOuterPixel = 60; //[mm]
175
176 Double_t probability = 1.;
177
178 // Resolution and computing time are proportional to bins^2
179 const Int_t bins = 200;
180 Double_t max_x_mm = 600;
181 Double_t min_x_mm = -max_x_mm;
182 Double_t max_y_mm = max_x_mm;
183 Double_t min_y_mm = -max_x_mm;
184
185 // bins to mmrees
186 const Double_t bin2mm = 2 * max_x_mm / bins;
187
188 // First run, to find the maximum peak and define the area
189 TH2F *hgrid = new TH2F("hgrid", "", bins, min_x_mm, max_x_mm, bins, min_y_mm, max_y_mm);
190
191 for (Int_t ny = 1; ny <= bins; ny++)
192 for (Int_t nx = 1; nx <= bins; nx++)
193 hgrid->SetBinContent(nx, ny, 1.);
194
195 for (UInt_t numstar = 0; numstar < fNumStars; numstar++)
196 {
197 probability = 1;
198
199 MStarLocalPos& star = (*fStars)[numstar];
200
201 if (star.GetChiSquareNdof() > 0. && star.GetChiSquareNdof() < 10.)
202 {
203
204 Float_t starPosX = star.GetMeanX();
205 Float_t starPosY = star.GetMeanY();
206 Float_t starDist = Dist(0.,0.,starPosX,starPosY);
207 Float_t starSigma = (starDist<350.?diameterInnerPixel:diameterOuterPixel);
208
209// *fLog << dbg << "Star[" << numstar << "] pos (" << starPosX << "," << starPosY << ") dist " << starDist << " size " << starSigma << endl;
210
211 if (starSigma != 0)
212 {
213 for (Int_t ny = 1; ny < bins + 1; ny++)
214 {
215 for (Int_t nx = 0; nx < bins + 1; nx++)
216 {
217 Float_t dist = Dist(min_x_mm + nx * bin2mm, starPosX, min_y_mm + ny * bin2mm, starPosY);
218 Float_t prob = Prob(dist, fDistances[numstar], starSigma);
219
220// if ( prob > 0.8)
221// *fLog << dbg << " x " << min_x_mm + nx * bin2mm << " y " << min_y_mm + ny * bin2mm << " dist " << dist << " stardist " << fDistances[numstar] << " prob " << prob << endl;
222
223 probability = hgrid->GetBinContent(nx, ny)*prob;
224 hgrid->SetBinContent(nx, ny, probability);
225
226 }
227 }
228 }
229 else probability *= 1;
230
231 }
232 }
233
234 // Finding the peak
235 Double_t peak[2] = {0,0};
236 Double_t peak_value = 0;
237
238 for (Int_t ny = 0; ny < bins + 1; ny++)
239 {
240 for (Int_t nx = 0; nx < bins + 1; nx++)
241 {
242 if (hgrid->GetBinContent(nx, ny) > peak_value)
243 {
244 peak_value = hgrid->GetBinContent(nx, ny);
245 peak[0] = min_x_mm + nx * bin2mm;
246 peak[1] = min_y_mm + ny * bin2mm;
247 }
248 }
249 }
250
251 *fLog << dbg << "The peak is at (x, y) = (" << peak[0] << ", " << peak[1] << ") mm" << endl;
252
253
254// TCanvas c1;
255// hgrid->Draw("lego");
256// if(!HandleInput())
257// exit(1);
258
259
260 // Here we define a small area surrounding the peak to avoid wasting time and resolution anywhere else
261
262 min_x_mm = peak[0] - diameterInnerPixel;
263 max_x_mm = peak[0] + diameterInnerPixel;
264 min_y_mm = peak[1] - diameterInnerPixel;
265 max_y_mm = peak[1] + diameterInnerPixel;
266
267 const Double_t xbin2mm = (max_x_mm - min_x_mm) / bins;
268 const Double_t ybin2mm = (max_y_mm - min_y_mm) / bins;
269
270 // Other run centered in the peak for more precision
271 TH2F *hagrid = new TH2F("hagrid", "", bins, min_x_mm, max_x_mm, bins, min_y_mm, max_y_mm);
272
273 for (Int_t ny = 1; ny <= bins; ny++)
274 for (Int_t nx = 1; nx <= bins; nx++)
275 hagrid->SetBinContent(nx, ny, 1.);
276
277
278 for (UInt_t numstar = 0; numstar < fNumStars; numstar++)
279 {
280 probability = 1;
281
282 MStarLocalPos& star = (*fStars)[numstar];
283
284 if (star.GetChiSquareNdof() > 0. && star.GetChiSquareNdof() < 10.)
285 {
286
287 Float_t starPosX = star.GetMeanX();
288 Float_t starPosY = star.GetMeanY();
289 Float_t starDist = Dist(0.,0.,starPosX,starPosY);
290 Float_t starSigma = (starDist<350.?diameterInnerPixel:diameterOuterPixel);
291
292 if (starSigma != 0)
293 {
294 for (Int_t ny = 0; ny < bins; ny++)
295 {
296 for (Int_t nx = 0; nx < bins; nx++)
297 {
298 Float_t prob = Prob(Dist(min_x_mm + nx * xbin2mm, starPosX, min_y_mm + ny * ybin2mm, starPosY), fDistances[numstar], starSigma);
299
300 probability = hagrid->GetBinContent(nx, ny)*prob;
301 hagrid->SetBinContent(nx, ny, probability);
302 }
303 }
304 }
305 else probability *= 1;
306
307 }
308 }
309
310 // This time we fit the histogram with a 2D gaussian
311 // Although we don't expect our function to be gaussian...
312 TF2 *gauss2d = new TF2("gauss2d","[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))*exp(-(y-[3])*(y-[3])/(2*[4]*[4]))", min_x_mm, max_x_mm, min_y_mm, max_y_mm);
313
314 gauss2d->SetParName(0,"f0");
315 gauss2d->SetParName(1,"meanx");
316 gauss2d->SetParName(2,"sigmax");
317 gauss2d->SetParName(3,"meany");
318 gauss2d->SetParName(4,"sigmay");
319
320 gauss2d->SetParameter(0,10);
321 gauss2d->SetParameter(1,peak[0]);
322 gauss2d->SetParameter(2,0.002);
323 gauss2d->SetParameter(3,peak[1]);
324 gauss2d->SetParameter(4,0.002);
325
326 GetOutputSrcPosCam()->SetXY(gauss2d->GetParameter(1), gauss2d->GetParameter(3));
327
328 delete hgrid;
329 delete hagrid;
330 }
331
332
333 return kTRUE;
334}
335
Note: See TracBrowser for help on using the repository browser.