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

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