source: trunk/MagicSoft/Mars/mtemp/MFindStars.cc@ 4019

Last change on this file since 4019 was 4011, checked in by jlopez, 22 years ago
*** empty log message ***
File size: 10.6 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! Author(s): Robert Wagner, 2/2004 <mailto:rwagner@mppmu.mpg.de>
18! Author(s): Javier López , 4/2004 <mailto:jlopez@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MFindStars
28//
29/////////////////////////////////////////////////////////////////////////////
30#include "MFindStars.h"
31
32#include <TFile.h>
33#include <TTree.h>
34
35#include "MObservatory.h"
36#include "MAstroCamera.h"
37#include "MMcConfigRunHeader.h"
38
39#include "MLog.h"
40#include "MLogManip.h"
41
42#include "MHCamera.h"
43
44#include "MGeomCam.h"
45#include "MGeomPix.h"
46#include "MCameraDC.h"
47#include "MTime.h"
48#include "MReportDrive.h"
49#include "MStarLocalCam.h"
50#include "MStarLocalPos.h"
51
52#include "MParList.h"
53
54ClassImp(MFindStars);
55using namespace std;
56
57MFindStars::MFindStars(const char *name, const char *title)
58{
59 fName = name ? name : "MFindStars";
60 fTitle = title ? title : "Tool to find stars from DC Currents";
61
62 fMaxNumIntegratedEvents = 10;
63 fRingInterest = 100.; //[mm] ~ 0.3 deg
64 fMinDCForStars = 3.*fMaxNumIntegratedEvents; //[uA]
65
66 fPixelsUsed.Set(577);
67 fPixelsUsed.Reset((Char_t)kTRUE);
68}
69
70Int_t MFindStars::PreProcess(MParList *pList)
71{
72
73 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
74
75 if (!fGeomCam)
76 {
77 *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl;
78 return kFALSE;
79 }
80
81 fDisplay.SetGeometry(*fGeomCam);
82
83 fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));
84
85 if (!fCurr)
86 {
87 *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl;
88 return kFALSE;
89 }
90
91 fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents"));
92
93 if (!fTimeCurr)
94 {
95 *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl;
96 return kFALSE;
97 }
98
99 fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
100
101 if (!fDrive)
102 {
103 *fLog << warn << AddSerialNumber("MReportDrive") << " not found ... ignored." << endl;
104 }
105 else
106 {
107 //Initialitation MAstroCamera
108 // Name of a MC file having MGeomCam and MMcConfigRunHeader
109 TString fname = "../data/Gamma_zbin9_90_7_1480to1489_w0.root";
110
111 // Time for which to get the picture
112 // MTime time;
113 // time.Set(2004, 2, 28, 01, 32, 15);
114
115 // Current observatory
116 MObservatory magic1;
117
118 // Right Ascension [h] and declination [deg] of source
119 // Currently 'perfect' pointing is assumed
120 // const Double_t ra = MAstro::Hms2Rad(5, 34, 31.9);
121 // const Double_t dec = MAstro::Dms2Rad(22, 0, 52.0);
122
123 // --------------------------------------------------------------------------
124 // Create camera display from geometry
125 MMcConfigRunHeader *config=0;
126 MGeomCam *geom=0;
127
128 TFile file(fname);
129 TTree *tree = (TTree*)file.Get("RunHeaders");
130 tree->SetBranchAddress("MMcConfigRunHeader", &config);
131 if (tree->GetBranch("MGeomCam"))
132 tree->SetBranchAddress("MGeomCam", &geom);
133 tree->GetEntry(0);
134
135 fAstro.SetMirrors(*config->GetMirrors());
136 fAstro.SetGeom(*geom);
137
138
139 fAstro.SetLimMag(6);
140 fAstro.SetRadiusFOV(3);
141 fAstro.ReadBSC("../data/bsc5.dat");
142
143 fAstro.SetObservatory(magic1);
144 // Time for which to get the picture
145 // MTime time;
146 // time.Set(2004, 2, 28, 01, 32, 15);
147 // fAstro.SetTime(time);
148 fAstro.SetTime(*fTimeCurr);
149 fAstro.SetGuiActive();
150
151 fAstro.SetStarList(fStars->GetList());
152
153 }
154
155
156 fStars = (MStarLocalCam*)pList->FindCreateObj(AddSerialNumber("MStarLocalCam"));
157 if (!fStars)
158 {
159 *fLog << err << AddSerialNumber("MStarLocalCam") << " cannot be created ... aborting" << endl;
160 return kFALSE;
161 }
162
163 return kTRUE;
164}
165
166Int_t MFindStars::Process()
167{
168 if (fNumIntegratedEvents < fMaxNumIntegratedEvents)
169 {
170 fDisplay.AddCamContent(*fCurr);
171 fNumIntegratedEvents++;
172 }
173 else
174 {
175 if (fDrive)
176 {
177 Float_t ra = fDrive->GetRa();
178 Float_t dec = fDrive->GetDec();
179
180 fAstro.SetRaDec(ra, dec);
181 fAstro.SetGuiActive();
182
183 fAstro.FillStarList();
184 }
185 else
186 {
187 Float_t maxPixelDC;
188 MGeomPix maxPixel;
189
190 MHCamera cam = fDisplay;
191 cam.PrintInfo();
192
193 // Find the star candidats searching the most brights pairs of pixels
194 while(1)
195 {
196 FindPixelWithMaxDC(maxPixelDC, maxPixel);
197 *fLog << dbg << "Star candidate maxDC(" << setw(3) << maxPixelDC << " uA) x position(" << setw(3) << maxPixel.GetX() << " mm) x position(" << setw(3) << maxPixel.GetY() << " mm)" << endl;
198
199 if (maxPixelDC<fMinDCForStars)
200 break;
201
202 MStarLocalPos *starpos = new MStarLocalPos;
203 starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
204 starpos->SetCalcValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
205 fStars->GetList()->Add(starpos);
206
207 ShadowStar(starpos);
208
209 }
210 fDisplay = cam;
211 *fLog << dbg << GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
212 }
213
214 //loop to extract position of stars on the camera
215 if (fStars->GetList()->GetSize() == 0)
216 {
217 *fLog << err << GetName() << " No stars candidates in the camera." << endl;
218 return kFALSE;
219 }
220
221 TIter Next(fStars->GetList());
222 MStarLocalPos* star;
223 while ((star=(MStarLocalPos*)Next()))
224 {
225 FindStar(star);
226 ShadowStar(star);
227 }
228
229
230 //After finding stars reset all vairables
231 fDisplay.Reset();
232 fNumIntegratedEvents=0;
233 }
234
235 return kTRUE;
236}
237
238Int_t MFindStars::PostProcess()
239{
240 if(fStars->GetList()->GetSize() != 0)
241 fStars->Print();
242
243 return kTRUE;
244}
245
246void MFindStars::SetBlindPixels(TArrayS blindpixels)
247{
248 Int_t npix = blindpixels.GetSize();
249
250 for (Int_t idx=0; idx<npix; idx++)
251 fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;
252
253 fDisplay.SetUsed(fPixelsUsed);
254}
255
256Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
257{
258 UInt_t numPixels = fGeomCam->GetNumPixels();
259
260// Find the two close pixels with the maximun dc
261 UInt_t maxPixIdx[2];
262
263 maxDC = 0;
264
265 for (UInt_t pix=1; pix<numPixels; pix++)
266 {
267 if(fDisplay.IsUsed(pix))
268 {
269 Float_t dc[2];
270 dc[0] = fDisplay.GetBinContent(pix+1);
271
272 MGeomPix &g = (*fGeomCam)[pix];
273 Int_t numNextNeighbors = g.GetNumNeighbors();
274
275 Float_t dcsum;
276 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
277 {
278 if(fDisplay.IsUsed(pix))
279 {
280 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
281 dc[1] = fDisplay.GetBinContent(swneighbor+1);
282
283
284 dcsum = dc[0] + dc[1];
285
286 if(dcsum > maxDC*2)
287 {
288 if(dc[0]>=dc[1])
289 {
290 maxPixIdx[0] = pix;
291 maxPixIdx[1] = swneighbor;
292 maxDC = dc[0];
293 }
294 else
295 {
296 maxPixIdx[1] = pix;
297 maxPixIdx[0] = swneighbor;
298 maxDC = dc[1];
299 }
300 }
301 }
302 }
303 }
304 }
305
306 maxPix = (*fGeomCam)[maxPixIdx[0]];
307 return kTRUE;
308}
309
310Bool_t MFindStars::FindStar(MStarLocalPos* star)
311{
312
313 MHCamera cam = fDisplay;
314
315 Float_t expX = star->GetXExp();
316 Float_t expY = star->GetYExp();
317
318 UInt_t numPixels = fGeomCam->GetNumPixels();
319
320// First define a area of interest around the expected position of the star
321 for (UInt_t pix=1; pix<numPixels; pix++)
322 {
323
324 Float_t pixXpos=(*fGeomCam)[pix].GetX();
325 Float_t pixYpos=(*fGeomCam)[pix].GetY();
326 Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
327 (pixYpos-expY)*(pixYpos-expY));
328
329 if (dist < fRingInterest && cam.IsUsed(pix))
330 fPixelsUsed[pix]=(Char_t)kTRUE;
331 }
332
333 cam.SetUsed(fPixelsUsed);
334
335
336 // determine mean x and mean y of the selected px
337 Float_t meanX=0;
338 Float_t meanY=0;
339 Float_t meanSqX=0;
340 Float_t meanSqY=0;
341 Float_t sumCharge=0;
342 UInt_t usedPx=0;
343 for(UInt_t pix=0; pix<numPixels; pix++)
344 {
345 if(cam.IsUsed(pix))
346 {
347 usedPx++;
348
349 Float_t charge=cam.GetBinContent(pix+1);
350 Float_t pixXpos=(*fGeomCam)[pix].GetX();
351 Float_t pixYpos=(*fGeomCam)[pix].GetY();
352
353 meanX+=charge*pixXpos;
354 meanY+=charge*pixYpos;
355 meanSqX+=charge*pixXpos*pixXpos;
356 meanSqY+=charge*pixYpos*pixYpos;
357 sumCharge+=charge;
358// cam.ResetUsed(i); // this px must not be used again!
359 } //for... use
360 }
361
362 star->SetCalcValues(sumCharge,meanX,meanY,meanSqX,meanSqY);
363
364 return kTRUE;
365}
366
367Bool_t MFindStars::ShadowStar(MStarLocalPos* star)
368{
369 UInt_t numPixels = fGeomCam->GetNumPixels();
370
371// Define an area around the star which will be set unused.
372 for (UInt_t pix=1; pix<numPixels; pix++)
373 {
374
375 Float_t pixXpos = (*fGeomCam)[pix].GetX();
376 Float_t pixYpos = (*fGeomCam)[pix].GetY();
377 Float_t starXpos = star->GetMeanXCalc();
378 Float_t starYpos = star->GetMeanYCalc();
379
380 Float_t starSize = 2*star->GetSigmaMajorAxisCalc();
381
382 Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
383 (pixYpos-starYpos)*(pixYpos-starYpos));
384
385 if (dist > starSize && fDisplay.IsUsed(pix))
386 fPixelsUsed[pix]=(Char_t)kTRUE;
387 else
388 fPixelsUsed[pix]=(Char_t)kFALSE;
389 }
390
391 fDisplay.SetUsed(fPixelsUsed);
392
393 return kTRUE;
394}
395
396
397
Note: See TracBrowser for help on using the repository browser.