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

Last change on this file since 4043 was 4037, checked in by jlopez, 21 years ago
*** empty log message ***
File size: 13.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 <TTimer.h>
33#include <TString.h>
34#include <TFile.h>
35#include <TTree.h>
36#include <TH1F.h>
37#include <TF1.h>
38
39#include "MObservatory.h"
40#include "MAstroCamera.h"
41#include "MMcConfigRunHeader.h"
42
43#include "MLog.h"
44#include "MLogManip.h"
45
46#include "MHCamera.h"
47
48#include "MGeomCam.h"
49#include "MGeomPix.h"
50#include "MCameraDC.h"
51#include "MTime.h"
52#include "MReportDrive.h"
53#include "MStarLocalCam.h"
54#include "MStarLocalPos.h"
55
56#include "MParList.h"
57
58ClassImp(MFindStars);
59using namespace std;
60
61MFindStars::MFindStars(const char *name, const char *title)
62{
63 fName = name ? name : "MFindStars";
64 fTitle = title ? title : "Tool to find stars from DC Currents";
65
66 fNumIntegratedEvents=0;
67 fMaxNumIntegratedEvents = 10;
68 fRingInterest = 125.; //[mm] ~ 0.4 deg
69 fMinDCForStars = 2.*fMaxNumIntegratedEvents; //[uA]
70
71 fPixelsUsed.Set(577);
72 fPixelsUsed.Reset((Char_t)kTRUE);
73}
74
75Int_t MFindStars::PreProcess(MParList *pList)
76{
77
78 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
79
80 if (!fGeomCam)
81 {
82 *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl;
83 return kFALSE;
84 }
85
86 fDisplay.SetGeometry(*fGeomCam);
87
88 fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));
89
90 if (!fCurr)
91 {
92 *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl;
93 return kFALSE;
94 }
95
96 fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents"));
97
98 if (!fTimeCurr)
99 {
100 *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl;
101 return kFALSE;
102 }
103
104 fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
105
106 if (!fDrive)
107 {
108 *fLog << warn << AddSerialNumber("MReportDrive") << " not found ... ignored." << endl;
109 }
110 else
111 {
112 //Initialitation MAstroCamera
113 // Name of a MC file having MGeomCam and MMcConfigRunHeader
114 TString fname = "../data/Gamma_zbin9_90_7_1480to1489_w0.root";
115
116 // Time for which to get the picture
117 // MTime time;
118 // time.Set(2004, 2, 28, 01, 32, 15);
119
120 // Current observatory
121 MObservatory magic1;
122
123 // Right Ascension [h] and declination [deg] of source
124 // Currently 'perfect' pointing is assumed
125 // const Double_t ra = MAstro::Hms2Rad(5, 34, 31.9);
126 // const Double_t dec = MAstro::Dms2Rad(22, 0, 52.0);
127
128 // --------------------------------------------------------------------------
129 // Create camera display from geometry
130 MMcConfigRunHeader *config=0;
131 MGeomCam *geom=0;
132
133 TFile file(fname);
134 TTree *tree = (TTree*)file.Get("RunHeaders");
135 tree->SetBranchAddress("MMcConfigRunHeader", &config);
136 if (tree->GetBranch("MGeomCam"))
137 tree->SetBranchAddress("MGeomCam", &geom);
138 tree->GetEntry(0);
139
140 fAstro.SetMirrors(*config->GetMirrors());
141 fAstro.SetGeom(*geom);
142
143
144 fAstro.SetLimMag(6);
145 fAstro.SetRadiusFOV(3);
146 fAstro.ReadBSC("../data/bsc5.dat");
147
148 fAstro.SetObservatory(magic1);
149 // Time for which to get the picture
150 // MTime time;
151 // time.Set(2004, 2, 28, 01, 32, 15);
152 // fAstro.SetTime(time);
153 fAstro.SetTime(*fTimeCurr);
154 fAstro.SetGuiActive();
155
156 fAstro.SetStarList(fStars->GetList());
157
158 }
159
160
161 fStars = (MStarLocalCam*)pList->FindCreateObj(AddSerialNumber("MStarLocalCam"));
162 if (!fStars)
163 {
164 *fLog << err << AddSerialNumber("MStarLocalCam") << " cannot be created ... aborting" << endl;
165 return kFALSE;
166 }
167
168 return kTRUE;
169}
170
171Int_t MFindStars::Process()
172{
173 if (fNumIntegratedEvents >= fMaxNumIntegratedEvents)
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 UInt_t numPixels = fGeomCam->GetNumPixels();
188 TArrayC origPixelsUsed;
189 origPixelsUsed.Set(numPixels);
190
191 for (UInt_t pix=1; pix<numPixels; pix++)
192 {
193 if (fDisplay.IsUsed(pix))
194 origPixelsUsed[pix]=(Char_t)kTRUE;
195 else
196 origPixelsUsed[pix]=(Char_t)kFALSE;
197 }
198
199 Float_t ped;
200 Float_t rms;
201 DCPedestalCalc(ped, rms);
202 fMinDCForStars = fMinDCForStars>(ped+5*rms)?fMinDCForStars:(ped+5*rms);
203
204 *fLog << dbg << " DC pedestal = " << ped << " pedestal rms = " << rms << endl;
205 *fLog << dbg << " fMinDCForStars " << fMinDCForStars << endl;
206
207 // Find the star candidats searching the most brights pairs of pixels
208 Float_t maxPixelDC;
209 MGeomPix maxPixel;
210
211 while(FindPixelWithMaxDC(maxPixelDC, maxPixel))
212 {
213 *fLog << dbg << "Star candidate maxDC(" << setw(3) << maxPixelDC << " uA) x position(" << setw(3) << maxPixel.GetX() << " mm) x position(" << setw(3) << maxPixel.GetY() << " mm)" << endl;
214
215 MStarLocalPos *starpos = new MStarLocalPos;
216 starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
217 starpos->SetCalcValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
218 starpos->SetFitValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
219 fStars->GetList()->Add(starpos);
220
221 ShadowStar(starpos);
222
223 }
224
225 fDisplay.SetUsed(origPixelsUsed);
226 }
227
228 //loop to extract position of stars on the camera
229 if (fStars->GetList()->GetSize() == 0)
230 {
231 *fLog << err << GetName() << " No stars candidates in the camera." << endl;
232 return kFALSE;
233 }
234 else
235 *fLog << inf << GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
236
237 TIter Next(fStars->GetList());
238 MStarLocalPos* star;
239 while ((star=(MStarLocalPos*)Next()))
240 {
241 FindStar(star);
242 ShadowStar(star);
243 }
244
245 //After finding stars reset all vairables
246 fDisplay.Reset();
247 fNumIntegratedEvents=0;
248 }
249
250 fDisplay.AddCamContent(*fCurr);
251 fNumIntegratedEvents++;
252
253 return kTRUE;
254}
255
256Int_t MFindStars::PostProcess()
257{
258 if(fStars->GetList()->GetSize() != 0)
259 fStars->Print();
260
261 return kTRUE;
262}
263
264void MFindStars::SetBlindPixels(TArrayS blindpixels)
265{
266 Int_t npix = blindpixels.GetSize();
267
268 for (Int_t idx=0; idx<npix; idx++)
269 fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;
270
271 fDisplay.SetUsed(fPixelsUsed);
272}
273
274Bool_t MFindStars::DCPedestalCalc(Float_t &ped, Float_t &rms)
275{
276 UInt_t numPixels = fGeomCam->GetNumPixels();
277
278 TH1F dchist("dchist","",120,0.,30.*fMaxNumIntegratedEvents);
279 for (UInt_t pix=1; pix<numPixels; pix++)
280 dchist.Fill(fDisplay.GetBinContent(pix+1));
281
282 Float_t nummaxprobdc = dchist.GetBinContent(dchist.GetMaximumBin());
283 Float_t maxprobdc = dchist.GetBinCenter(dchist.GetMaximumBin());
284 UInt_t bin = dchist.GetMaximumBin();
285 do
286 {
287 bin++;
288 }
289 while(dchist.GetBinContent(bin)/nummaxprobdc > 0.5);
290 Float_t halfmaxprobdc = dchist.GetBinCenter(bin);
291
292 *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] ";
293 *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << dchist.GetBinContent(bin) << "]" << endl;
294
295 Float_t rmsguess = TMath::Abs(maxprobdc-halfmaxprobdc);
296 Float_t min = maxprobdc-3*rmsguess;
297 min = (min<0.?0.:min);
298 Float_t max = maxprobdc+3*rmsguess;
299
300 *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
301
302 TF1 func("func","gaus",min,max);
303 func.SetParameters(nummaxprobdc, maxprobdc, rmsguess);
304
305 dchist.Fit("func","QR0");
306
307 UInt_t aproxnumdegrees = 6*(bin-dchist.GetMaximumBin());
308 Float_t chiq = func.GetChisquare();
309 ped = func.GetParameter(1);
310 rms = func.GetParameter(2);
311
312 *fLog << dbg << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
313
314 Int_t numsigmas = 5;
315 Axis_t minbin = ped-numsigmas*rms/dchist.GetBinWidth(1);
316 minbin=minbin<1?1:minbin;
317 Axis_t maxbin = ped+numsigmas*rms/dchist.GetBinWidth(1);
318 *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << dchist.Integral((int)minbin,(int)maxbin) << endl;
319 return kTRUE;
320}
321
322Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
323{
324 UInt_t numPixels = fGeomCam->GetNumPixels();
325
326// Find the two close pixels with the maximun dc
327 UInt_t maxPixIdx[2];
328
329 maxDC = 0;
330
331 for (UInt_t pix=1; pix<numPixels; pix++)
332 {
333 if(fDisplay.IsUsed(pix))
334 {
335 Float_t dc[2];
336 dc[0] = fDisplay.GetBinContent(pix+1);
337 if (dc[0] < fMinDCForStars)
338 continue;
339
340 MGeomPix &g = (*fGeomCam)[pix];
341 Int_t numNextNeighbors = g.GetNumNeighbors();
342
343 Float_t dcsum;
344 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
345 {
346 if(fDisplay.IsUsed(pix))
347 {
348 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
349 dc[1] = fDisplay.GetBinContent(swneighbor+1);
350 if (dc[1] < fMinDCForStars)
351 continue;
352
353 dcsum = dc[0] + dc[1];
354
355 if(dcsum > maxDC*2)
356 {
357 if(dc[0]>=dc[1])
358 {
359 maxPixIdx[0] = pix;
360 maxPixIdx[1] = swneighbor;
361 maxDC = dc[0];
362 }
363 else
364 {
365 maxPixIdx[1] = pix;
366 maxPixIdx[0] = swneighbor;
367 maxDC = dc[1];
368 }
369 }
370 }
371 }
372 }
373 }
374
375 if (maxDC == 0)
376 return kFALSE;
377
378 maxPix = (*fGeomCam)[maxPixIdx[0]];
379 return kTRUE;
380}
381
382Bool_t MFindStars::FindStar(MStarLocalPos* star)
383{
384
385 UInt_t numPixels = fGeomCam->GetNumPixels();
386 TArrayC origPixelsUsed;
387 origPixelsUsed.Set(numPixels);
388
389 for (UInt_t pix=1; pix<numPixels; pix++)
390 {
391 if (fDisplay.IsUsed(pix))
392 origPixelsUsed[pix]=(Char_t)kTRUE;
393 else
394 origPixelsUsed[pix]=(Char_t)kFALSE;
395 }
396
397 Float_t expX = star->GetXExp();
398 Float_t expY = star->GetYExp();
399
400 // First define a area of interest around the expected position of the star
401 for (UInt_t pix=1; pix<numPixels; pix++)
402 {
403
404 Float_t pixXpos=(*fGeomCam)[pix].GetX();
405 Float_t pixYpos=(*fGeomCam)[pix].GetY();
406 Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
407 (pixYpos-expY)*(pixYpos-expY));
408
409 if ((dist < fRingInterest) && fDisplay.IsUsed(pix))
410 fPixelsUsed[pix]=(Char_t)kTRUE;
411 else
412 fPixelsUsed[pix]=(Char_t)kFALSE;
413 }
414
415 fDisplay.SetUsed(fPixelsUsed);
416
417 // determine mean x and mean y of the selected px
418 Float_t meanX=0;
419 Float_t meanY=0;
420 Float_t meanSqX=0;
421 Float_t meanSqY=0;
422 Float_t sumCharge=0;
423 UInt_t usedPx=0;
424 for(UInt_t pix=0; pix<numPixels; pix++)
425 {
426 if(fDisplay.IsUsed(pix))
427 {
428 usedPx++;
429
430 Float_t charge = fDisplay.GetBinContent(pix+1);
431 Float_t pixXpos = (*fGeomCam)[pix].GetX();
432 Float_t pixYpos = (*fGeomCam)[pix].GetY();
433
434 meanX += charge*pixXpos;
435 meanY += charge*pixYpos;
436 meanSqX += charge*pixXpos*pixXpos;
437 meanSqY += charge*pixYpos*pixYpos;
438 sumCharge += charge;
439 }
440 }
441
442 *fLog << dbg << " usedPx " << usedPx << endl;
443
444 meanX /= sumCharge;
445 meanY /= sumCharge;
446 meanSqX /= sumCharge;
447 meanSqY /= sumCharge;
448
449 Float_t rmsX = TMath::Sqrt(meanSqX - meanX*meanX);
450 Float_t rmsY = TMath::Sqrt(meanSqY - meanY*meanY);
451
452 star->SetCalcValues(sumCharge,meanX,meanY,rmsX,rmsY);
453
454 fDisplay.SetUsed(origPixelsUsed);
455
456 return kTRUE;
457}
458
459Bool_t MFindStars::ShadowStar(MStarLocalPos* star)
460{
461 UInt_t numPixels = fGeomCam->GetNumPixels();
462
463// Define an area around the star which will be set unused.
464 UInt_t shadowPx=0;
465 for (UInt_t pix=1; pix<numPixels; pix++)
466 {
467 Float_t pixXpos = (*fGeomCam)[pix].GetX();
468 Float_t pixYpos = (*fGeomCam)[pix].GetY();
469 Float_t starXpos = star->GetMeanXCalc();
470 Float_t starYpos = star->GetMeanYCalc();
471
472 Float_t starSize = 3*star->GetSigmaMajorAxisCalc();
473
474 Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
475 (pixYpos-starYpos)*(pixYpos-starYpos));
476
477 if (dist > starSize && fDisplay.IsUsed(pix))
478 fPixelsUsed[pix]=(Char_t)kTRUE;
479 else
480 {
481 fPixelsUsed[pix]=(Char_t)kFALSE;
482 shadowPx++;
483 }
484 }
485
486 *fLog << dbg << " shadowPx " << shadowPx << endl;
487
488 fDisplay.SetUsed(fPixelsUsed);
489
490 return kTRUE;
491}
492
493
494
Note: See TracBrowser for help on using the repository browser.