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

Last change on this file since 4002 was 4002, checked in by jlopez, 21 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! 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 while(1)
189 {
190 FindPixelWithMaxDC(maxPixelDC, maxPixel);
191 if (maxPixelDC<fMinDCForStars)
192 break;
193
194 MStarLocalPos *starpos = new MStarLocalPos;
195 starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
196 fStars->GetList()->Add(starpos);
197 }
198 }
199
200 //loop to extract position of stars on the camera
201 TIter Next(fStars->GetList());
202 MStarLocalPos* star;
203 while ((star=(MStarLocalPos*)Next()))
204 {
205 FindStar(star);
206 ShadowStar(star);
207 }
208
209
210 //After finding stars reset all vairables
211 fDisplay.Reset();
212 fNumIntegratedEvents=0;
213 }
214
215 return kTRUE;
216}
217
218Int_t MFindStars::PostProcess()
219{
220 fStars->Print();
221 return kTRUE;
222}
223
224void MFindStars::SetBlindPixels(TArrayS blindpixels)
225{
226 Int_t npix = blindpixels.GetSize();
227
228 for (Int_t idx=0; idx<npix; idx++)
229 fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;
230
231 fDisplay.SetUsed(fPixelsUsed);
232}
233
234Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
235{
236 UInt_t numPixels = fGeomCam->GetNumPixels();
237
238// Find the two close pixels with the maximun dc
239 UInt_t maxPixIdx[2];
240
241 for (UInt_t pix=1; pix<numPixels; pix++)
242 {
243 maxDC = 0;
244
245 if(fDisplay.IsUsed(pix))
246 {
247 Float_t dc[2];
248 dc[0] = fDisplay.GetBinContent(pix+1);
249
250 MGeomPix &g = (*fGeomCam)[pix];
251 Int_t numNextNeighbors = g.GetNumNeighbors();
252
253 Float_t dcsum;
254 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
255 {
256 if(fDisplay.IsUsed(pix))
257 {
258 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
259 dc[1] = fDisplay.GetBinContent(swneighbor+1);
260
261
262 dcsum = dc[0] + dc[1];
263
264 if(dcsum > maxDC*2)
265 {
266 if(dc[0]>=dc[1])
267 {
268 maxPixIdx[0] = pix;
269 maxPixIdx[1] = swneighbor;
270 }
271 else
272 {
273 maxPixIdx[1] = pix;
274 maxPixIdx[0] = swneighbor;
275 }
276 maxDC = dcsum/2;
277 }
278 }
279 }
280 }
281 }
282
283 maxPix = (*fGeomCam)[maxPixIdx[0]];
284 return kTRUE;
285}
286
287Bool_t MFindStars::FindStar(MStarLocalPos* star)
288{
289
290 MHCamera cam = fDisplay;
291
292 Float_t expX = star->GetXExp();
293 Float_t expY = star->GetYExp();
294
295 UInt_t numPixels = fGeomCam->GetNumPixels();
296
297// First define a area of interest around the expected position of the star
298 for (UInt_t pix=1; pix<numPixels; pix++)
299 {
300
301 Float_t pixXpos=(*fGeomCam)[pix].GetX();
302 Float_t pixYpos=(*fGeomCam)[pix].GetY();
303 Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
304 (pixYpos-expY)*(pixYpos-expY));
305
306 if (dist < fRingInterest && cam.IsUsed(pix))
307 fPixelsUsed[pix]=(Char_t)kTRUE;
308 }
309
310 cam.SetUsed(fPixelsUsed);
311
312
313 // determine mean x and mean y of the selected px
314 Float_t meanX=0;
315 Float_t meanY=0;
316 Float_t meanSqX=0;
317 Float_t meanSqY=0;
318 Float_t sumCharge=0;
319 UInt_t usedPx=0;
320 for(UInt_t pix=0; pix<numPixels; pix++)
321 {
322 if(cam.IsUsed(pix))
323 {
324 usedPx++;
325
326 Float_t charge=cam.GetBinContent(pix+1);
327 Float_t pixXpos=(*fGeomCam)[pix].GetX();
328 Float_t pixYpos=(*fGeomCam)[pix].GetY();
329
330 meanX+=charge*pixXpos;
331 meanY+=charge*pixYpos;
332 meanSqX+=charge*pixXpos*pixXpos;
333 meanSqY+=charge*pixYpos*pixYpos;
334 sumCharge+=charge;
335// cam.ResetUsed(i); // this px must not be used again!
336 } //for... use
337 }
338
339 star->SetCalcValues(sumCharge,meanX,meanY,meanSqX,meanSqY);
340
341 return kTRUE;
342}
343
344Bool_t MFindStars::ShadowStar(MStarLocalPos* star)
345{
346 UInt_t numPixels = fGeomCam->GetNumPixels();
347
348// Define an area around the star which will be set unused.
349 for (UInt_t pix=1; pix<numPixels; pix++)
350 {
351
352 Float_t pixXpos = (*fGeomCam)[pix].GetX();
353 Float_t pixYpos = (*fGeomCam)[pix].GetY();
354 Float_t starXpos = star->GetMeanXCalc();
355 Float_t starYpos = star->GetMeanYCalc();
356
357 Float_t starSize = 2*star->GetSigmaMajorAxisCalc();
358
359 Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
360 (pixYpos-starYpos)*(pixYpos-starYpos));
361
362 if (dist > starSize && fDisplay.IsUsed(pix))
363 fPixelsUsed[pix]=(Char_t)kTRUE;
364 }
365
366 fDisplay.SetUsed(fPixelsUsed);
367
368 return kTRUE;
369}
370
371
372
Note: See TracBrowser for help on using the repository browser.