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

Last change on this file since 4056 was 4054, checked in by jlopez, 22 years ago
*** empty log message ***
File size: 20.9 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 <TMinuit.h>
33#include <TStopwatch.h>
34#include <TTimer.h>
35#include <TString.h>
36#include <TFile.h>
37#include <TTree.h>
38#include <TH1F.h>
39#include <TF1.h>
40
41#include "MObservatory.h"
42#include "MAstroCamera.h"
43#include "MMcConfigRunHeader.h"
44
45#include "MMinuitInterface.h"
46
47#include "MLog.h"
48#include "MLogManip.h"
49
50#include "MHCamera.h"
51
52#include "MGeomCam.h"
53#include "MGeomPix.h"
54#include "MCameraDC.h"
55#include "MTime.h"
56#include "MReportDrive.h"
57#include "MStarLocalCam.h"
58#include "MStarLocalPos.h"
59
60#include "MParList.h"
61
62ClassImp(MFindStars);
63using namespace std;
64
65const Float_t sqrt2 = sqrt(2.);
66const Float_t sqrt3 = sqrt(3.);
67
68//______________________________________________________________________________
69//
70// The 2D gaussian fucntion used to fit the spot of the star
71//
72static Double_t func(float x,float y,Double_t *par)
73{
74 Double_t value=par[0]*exp(-(x-par[1])*(x-par[1])/(2*par[2]*par[2]))*exp(-(y-par[3])*(y-par[3])/(2*par[4]*par[4]));
75 return value;
76}
77
78//______________________________________________________________________________
79//
80// Function used by Minuit to do the fit
81//
82static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
83{
84
85 MFindStars* find = (MFindStars*)gMinuit->GetObjectFit();
86 MHCamera* display = (MHCamera*)find->GetDisplay();
87 Float_t ped = find->GetPedestalDC();
88 Float_t rms = find->GetPedestalRMSDC();
89 MGeomCam& geom = (MGeomCam&)display->GetGeomCam();
90
91 UInt_t numPixels = geom.GetNumPixels();
92
93//calculate chisquare
94 Double_t chisq = 0;
95 Double_t delta;
96 Double_t x,y,z;
97 Double_t errorz = rms; //[uA]
98
99 UInt_t usedPx=0;
100 for (UInt_t pixid=1; pixid<numPixels; pixid++)
101 {
102 z = display->GetBinContent(pixid+1)-ped;
103
104 if (display->IsUsed(pixid) && z > 0.)
105 {
106 x = geom[pixid].GetX();
107 y = geom[pixid].GetY();
108
109 if (errorz > 0.0)
110 {
111 usedPx++;
112 delta = (z-func(x,y,par))/errorz;
113 chisq += delta*delta;
114 }
115 else
116 cerr << " TMinuit::fcn errorz[" << pixid << "] " << errorz << endl;
117 }
118 }
119 f = chisq;
120
121 find->SetChisquare(chisq);
122 find->SetDegreesofFreedom(usedPx);
123}
124
125MFindStars::MFindStars(const char *name, const char *title): fNumVar(5)
126{
127 fName = name ? name : "MFindStars";
128 fTitle = title ? title : "Tool to find stars from DC Currents";
129
130 fNumIntegratedEvents=0;
131 fMaxNumIntegratedEvents = 10;
132 fRingInterest = 125.; //[mm] ~ 0.4 deg
133 fMinDCForStars = 2.*fMaxNumIntegratedEvents; //[uA]
134
135 fPixelsUsed.Set(577);
136 fPixelsUsed.Reset((Char_t)kTRUE);
137
138 //Fitting(Minuit) initialitation
139 const Float_t pixelSize = 31.5; //[mm]
140
141 fVname = new TString[fNumVar];
142 fVinit.Set(fNumVar);
143 fStep.Set(fNumVar);
144 fLimlo.Set(fNumVar);
145 fLimup.Set(fNumVar);
146 fFix.Set(fNumVar);
147
148 fVname[0] = "max";
149 fVinit[0] = 10.*fMaxNumIntegratedEvents;
150 fStep[0] = fVinit[0]/sqrt2;
151 fLimlo[0] = fMinDCForStars;
152 fLimup[0] = 30.*fMaxNumIntegratedEvents;
153 fFix[0] = 0;
154
155 fVname[1] = "meanx";
156 fVinit[1] = 0.;
157 fStep[1] = fVinit[0]/sqrt2;
158 fLimlo[1] = -600.;
159 fLimup[1] = 600.;
160 fFix[1] = 0;
161
162 fVname[2] = "sigmaminor";
163 fVinit[2] = pixelSize;
164 fStep[2] = fVinit[0]/sqrt2;
165 fLimlo[2] = pixelSize/(2*sqrt3);
166 fLimup[2] = 500.;
167 fFix[2] = 0;
168
169 fVname[3] = "meany";
170 fVinit[3] = 0.;
171 fStep[3] = fVinit[0]/sqrt2;
172 fLimlo[3] = -600.;
173 fLimup[3] = 600.;
174 fFix[3] = 0;
175
176 fVname[4] = "sigmamajor";
177 fVinit[4] = pixelSize;
178 fStep[4] = fVinit[0]/sqrt2;
179 fLimlo[4] = pixelSize/(2*sqrt3);
180 fLimup[4] = 500.;
181 fFix[4] = 0;
182
183 fObjectFit = NULL;
184 // fMethod = "SIMPLEX";
185 // fMethod = "MIGRAD";
186 fMethod = "MINIMIZE";
187 fNulloutput = kFALSE;
188}
189
190Int_t MFindStars::PreProcess(MParList *pList)
191{
192
193 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
194
195 if (!fGeomCam)
196 {
197 *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl;
198 return kFALSE;
199 }
200
201 fDisplay.SetGeometry(*fGeomCam);
202 fDisplay.SetUsed(fPixelsUsed);
203
204 fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));
205
206 if (!fCurr)
207 {
208 *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl;
209 return kFALSE;
210 }
211
212 fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents"));
213
214 if (!fTimeCurr)
215 {
216 *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl;
217 return kFALSE;
218 }
219
220 fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
221
222 if (!fDrive)
223 {
224 *fLog << warn << AddSerialNumber("MReportDrive") << " not found ... ignored." << endl;
225 }
226 else
227 {
228 //Initialitation MAstroCamera
229 // Name of a MC file having MGeomCam and MMcConfigRunHeader
230 TString fname = "../data/Gamma_zbin9_90_7_1480to1489_w0.root";
231
232 // Time for which to get the picture
233 // MTime time;
234 // time.Set(2004, 2, 28, 01, 32, 15);
235
236 // Current observatory
237 MObservatory magic1;
238
239 // Right Ascension [h] and declination [deg] of source
240 // Currently 'perfect' pointing is assumed
241 // const Double_t ra = MAstro::Hms2Rad(5, 34, 31.9);
242 // const Double_t dec = MAstro::Dms2Rad(22, 0, 52.0);
243
244 // --------------------------------------------------------------------------
245 // Create camera display from geometry
246 MMcConfigRunHeader *config=0;
247 MGeomCam *geom=0;
248
249 TFile file(fname);
250 TTree *tree = (TTree*)file.Get("RunHeaders");
251 tree->SetBranchAddress("MMcConfigRunHeader", &config);
252 if (tree->GetBranch("MGeomCam"))
253 tree->SetBranchAddress("MGeomCam", &geom);
254 tree->GetEntry(0);
255
256 fAstro.SetMirrors(*config->GetMirrors());
257 fAstro.SetGeom(*geom);
258
259
260 fAstro.SetLimMag(6);
261 fAstro.SetRadiusFOV(3);
262 fAstro.ReadBSC("../data/bsc5.dat");
263
264 fAstro.SetObservatory(magic1);
265 // Time for which to get the picture
266 // MTime time;
267 // time.Set(2004, 2, 28, 01, 32, 15);
268 // fAstro.SetTime(time);
269 fAstro.SetTime(*fTimeCurr);
270 fAstro.SetGuiActive();
271
272 fAstro.SetStarList(fStars->GetList());
273
274 }
275
276
277 fStars = (MStarLocalCam*)pList->FindCreateObj(AddSerialNumber("MStarLocalCam"));
278 if (!fStars)
279 {
280 *fLog << err << AddSerialNumber("MStarLocalCam") << " cannot be created ... aborting" << endl;
281 return kFALSE;
282 }
283
284 return kTRUE;
285}
286
287Int_t MFindStars::Process()
288{
289 UInt_t numPixels = fGeomCam->GetNumPixels();
290 TArrayC origPixelsUsed;
291 origPixelsUsed.Set(numPixels);
292
293 if (fNumIntegratedEvents >= fMaxNumIntegratedEvents)
294 {
295
296 if (fDrive)
297 {
298 Float_t ra = fDrive->GetRa();
299 Float_t dec = fDrive->GetDec();
300
301 fAstro.SetRaDec(ra, dec);
302 fAstro.SetGuiActive();
303
304 fAstro.FillStarList();
305 }
306 else
307 {
308 //Fist delete the previus stars in the list
309 fStars->GetList()->Delete();
310 //
311
312 for (UInt_t pix=1; pix<numPixels; pix++)
313 {
314 if (fDisplay.IsUsed(pix))
315 origPixelsUsed[pix]=(Char_t)kTRUE;
316 else
317 origPixelsUsed[pix]=(Char_t)kFALSE;
318 }
319
320 DCPedestalCalc(fPedestalDC, fPedestalRMSDC);
321 fMinDCForStars = fMinDCForStars>(fPedestalDC+5*fPedestalRMSDC)?fMinDCForStars:(fPedestalDC+5*fPedestalRMSDC);
322
323 *fLog << dbg << " DC pedestal = " << fPedestalDC << " pedestal rms = " << fPedestalRMSDC << endl;
324 *fLog << dbg << " fMinDCForStars " << fMinDCForStars << endl;
325
326 // Find the star candidats searching the most brights pairs of pixels
327 Float_t maxPixelDC;
328 MGeomPix maxPixel;
329
330 while(FindPixelWithMaxDC(maxPixelDC, maxPixel))
331 {
332
333 MStarLocalPos *starpos = new MStarLocalPos;
334 starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
335 starpos->SetCalcValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
336 starpos->SetFitValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2,0.,1);
337 fStars->GetList()->Add(starpos);
338
339 ShadowStar(starpos);
340
341 }
342
343 fDisplay.SetUsed(origPixelsUsed);
344 }
345
346 //loop to extract position of stars on the camera
347 if (fStars->GetList()->GetSize() == 0)
348 {
349 *fLog << err << GetName() << " No stars candidates in the camera." << endl;
350 return kFALSE;
351 }
352 else
353 *fLog << inf << GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
354
355
356 for (UInt_t pix=1; pix<numPixels; pix++)
357 {
358 if (fDisplay.IsUsed(pix))
359 origPixelsUsed[pix]=(Char_t)kTRUE;
360 else
361 origPixelsUsed[pix]=(Char_t)kFALSE;
362 }
363
364 TIter Next(fStars->GetList());
365 MStarLocalPos* star;
366 while ((star=(MStarLocalPos*)Next()))
367 {
368 FindStar(star);
369 ShadowStar(star);
370 }
371
372 //After finding stars reset all vairables
373 fDisplay.Reset();
374 fDisplay.SetUsed(origPixelsUsed);
375 fNumIntegratedEvents=0;
376 }
377
378 for (UInt_t pix=1; pix<numPixels; pix++)
379 {
380 if (fDisplay.IsUsed(pix))
381 origPixelsUsed[pix]=(Char_t)kTRUE;
382 else
383 origPixelsUsed[pix]=(Char_t)kFALSE;
384
385 }
386
387 fDisplay.AddCamContent(*fCurr);
388 fNumIntegratedEvents++;
389 fDisplay.SetUsed(origPixelsUsed);
390
391
392 return kTRUE;
393}
394
395Int_t MFindStars::PostProcess()
396{
397 if(fStars->GetList()->GetSize() != 0)
398 fStars->Print();
399
400 return kTRUE;
401}
402
403void MFindStars::SetBlindPixels(TArrayS blindpixels)
404{
405 Int_t npix = blindpixels.GetSize();
406
407 for (Int_t idx=0; idx<npix; idx++)
408 {
409 fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;
410 *fLog << dbg << "MFindStars::SetBlindPixels fDisplay.IsUsed(" <<blindpixels[idx] << ") kFALSE" << endl;
411 }
412
413 fDisplay.SetUsed(fPixelsUsed);
414}
415
416Bool_t MFindStars::DCPedestalCalc(Float_t &ped, Float_t &rms)
417{
418 //-------------------------------------------------------------
419 // save pointer to the MINUIT object for optimizing the supercuts
420 // because it will be overwritten
421 // when fitting the alpha distribution in MHFindSignificance
422 TMinuit *savePointer = gMinuit;
423 //-------------------------------------------------------------
424
425 UInt_t numPixels = fGeomCam->GetNumPixels();
426
427 TH1F dchist("dchist","",120,0.,30.*fMaxNumIntegratedEvents);
428 for (UInt_t pix=1; pix<numPixels; pix++)
429 dchist.Fill(fDisplay.GetBinContent(pix+1));
430
431 Float_t nummaxprobdc = dchist.GetBinContent(dchist.GetMaximumBin());
432 Float_t maxprobdc = dchist.GetBinCenter(dchist.GetMaximumBin());
433 UInt_t bin = dchist.GetMaximumBin();
434 do
435 {
436 bin++;
437 }
438 while(dchist.GetBinContent(bin)/nummaxprobdc > 0.5);
439 Float_t halfmaxprobdc = dchist.GetBinCenter(bin);
440
441 *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] ";
442 *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << dchist.GetBinContent(bin) << "]" << endl;
443
444 Float_t rmsguess = TMath::Abs(maxprobdc-halfmaxprobdc);
445 Float_t min = maxprobdc-3*rmsguess;
446 min = (min<0.?0.:min);
447 Float_t max = maxprobdc+3*rmsguess;
448
449 *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
450
451 TF1 func("func","gaus",min,max);
452 func.SetParameters(nummaxprobdc, maxprobdc, rmsguess);
453
454 dchist.Fit("func","QR0");
455
456 UInt_t aproxnumdegrees = 6*(bin-dchist.GetMaximumBin());
457 Float_t chiq = func.GetChisquare();
458 ped = func.GetParameter(1);
459 rms = func.GetParameter(2);
460
461 *fLog << dbg << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
462
463 Int_t numsigmas = 5;
464 Axis_t minbin = ped-numsigmas*rms/dchist.GetBinWidth(1);
465 minbin=minbin<1?1:minbin;
466 Axis_t maxbin = ped+numsigmas*rms/dchist.GetBinWidth(1);
467
468 *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << dchist.Integral((int)minbin,(int)maxbin) << endl;
469
470 //Check results from the fit are consistent
471 if (TMath::Abs(ped-maxprobdc) > rmsguess || rms > rmsguess)
472 {
473 *fLog << warn << GetName() << " Pedestal DC fit give non consistent results." << endl;
474 *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
475 *fLog << warn << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
476 ped = maxprobdc;
477 rms = rmsguess/1.175; // FWHM=2.35*rms
478 }
479
480 //=================================================================
481 // reset gMinuit to the MINUIT object for optimizing the supercuts
482 gMinuit = savePointer;
483 //-------------------------------------------
484
485 return kTRUE;
486}
487
488Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
489{
490 UInt_t numPixels = fGeomCam->GetNumPixels();
491
492// Find the two close pixels with the maximun dc
493 UInt_t maxPixIdx[2];
494
495 maxDC = 0;
496
497 for (UInt_t pix=1; pix<numPixels; pix++)
498 {
499 if(fDisplay.IsUsed(pix))
500 {
501 Float_t dc[2];
502 dc[0] = fDisplay.GetBinContent(pix+1);
503 if (dc[0] < fMinDCForStars)
504 continue;
505
506 MGeomPix &g = (*fGeomCam)[pix];
507 Int_t numNextNeighbors = g.GetNumNeighbors();
508
509 Float_t dcsum;
510 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
511 {
512 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
513 if(fDisplay.IsUsed(swneighbor))
514 {
515 dc[1] = fDisplay.GetBinContent(swneighbor+1);
516 if (dc[1] < fMinDCForStars)
517 continue;
518
519 dcsum = dc[0] + dc[1];
520
521 if(dcsum > maxDC*2)
522 {
523 if(dc[0]>=dc[1])
524 {
525 maxPixIdx[0] = pix;
526 maxPixIdx[1] = swneighbor;
527 maxDC = dc[0];
528 }
529 else
530 {
531 maxPixIdx[1] = pix;
532 maxPixIdx[0] = swneighbor;
533 maxDC = dc[1];
534 }
535 }
536 }
537 }
538 }
539 }
540
541 if (maxDC == 0)
542 return kFALSE;
543
544 maxPix = (*fGeomCam)[maxPixIdx[0]];
545
546 *fLog << dbg << "Star candidate maxDC(" << setw(3) << maxDC << " uA) x position(" << setw(3) << maxPix.GetX() << " mm) x position(" << setw(3) << maxPix.GetY() << " mm) swnumber(" << maxPixIdx[0] << ")" << endl;
547
548 return kTRUE;
549}
550
551Bool_t MFindStars::FindStar(MStarLocalPos* star)
552{
553
554 UInt_t numPixels = fGeomCam->GetNumPixels();
555 TArrayC origPixelsUsed;
556 origPixelsUsed.Set(numPixels);
557
558 for (UInt_t pix=1; pix<numPixels; pix++)
559 {
560 if (fDisplay.IsUsed(pix))
561 origPixelsUsed[pix]=(Char_t)kTRUE;
562 else
563 origPixelsUsed[pix]=(Char_t)kFALSE;
564 }
565
566 Float_t expX = star->GetXExp();
567 Float_t expY = star->GetYExp();
568
569 // First define a area of interest around the expected position of the star
570 for (UInt_t pix=1; pix<numPixels; pix++)
571 {
572
573 Float_t pixXpos=(*fGeomCam)[pix].GetX();
574 Float_t pixYpos=(*fGeomCam)[pix].GetY();
575 Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
576 (pixYpos-expY)*(pixYpos-expY));
577
578 if ((dist < fRingInterest) && fDisplay.IsUsed(pix))
579 fPixelsUsed[pix]=(Char_t)kTRUE;
580 else
581 fPixelsUsed[pix]=(Char_t)kFALSE;
582 }
583
584 fDisplay.SetUsed(fPixelsUsed);
585
586// determine mean x and mean y
587 Float_t max=0;
588 Float_t meanX=0;
589 Float_t meanY=0;
590 Float_t meanSqX=0;
591 Float_t meanSqY=0;
592 Float_t sumCharge=0;
593 UInt_t usedPx=0;
594 for(UInt_t pix=0; pix<numPixels; pix++)
595 {
596 if(fDisplay.IsUsed(pix))
597 {
598 usedPx++;
599
600 Float_t charge = fDisplay.GetBinContent(pix+1);
601 Float_t pixXpos = (*fGeomCam)[pix].GetX();
602 Float_t pixYpos = (*fGeomCam)[pix].GetY();
603
604 if (charge>max) max=charge;
605 meanX += charge*pixXpos;
606 meanY += charge*pixYpos;
607 meanSqX += charge*pixXpos*pixXpos;
608 meanSqY += charge*pixYpos*pixYpos;
609 sumCharge += charge;
610 }
611 }
612
613 *fLog << dbg << " usedPx " << usedPx << endl;
614
615 meanX /= sumCharge;
616 meanY /= sumCharge;
617 meanSqX /= sumCharge;
618 meanSqY /= sumCharge;
619
620 //Substract pedestal from sumCharge
621 sumCharge -= usedPx*fPedestalDC;
622
623 Float_t rmsX = TMath::Sqrt(meanSqX - meanX*meanX);
624 Float_t rmsY = TMath::Sqrt(meanSqY - meanY*meanY);
625
626 star->SetCalcValues(sumCharge,max,meanX,meanY,rmsX,rmsY);
627
628
629// fit the star spot using TMinuit
630
631 for (UInt_t pix=1; pix<numPixels; pix++)
632 if (fDisplay.IsUsed(pix))
633 *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl;
634
635 //Initialate variables for fit
636 fVinit[0] = star->GetMaxCalc()-fPedestalDC;
637 fLimlo[0] = fMinDCForStars-fPedestalDC;
638 fLimup[0] = fLimup[0]-fPedestalDC;
639 fVinit[1] = meanX;
640 fVinit[2] = rmsX;
641 fVinit[3] = meanY;
642 fVinit[4] = rmsY;
643 //Init steps
644 for(Int_t i=0; i<fNumVar; i++)
645 if (fVinit[i] != 0)
646 fStep[i] = TMath::Abs(fVinit[i]/sqrt2);
647 //
648
649 TStopwatch clock;
650 clock.Start();
651
652 *fLog << dbg << " before calling CallMinuit" << endl;
653
654 MMinuitInterface inter;
655 Bool_t rc = inter.CallMinuit(fcn, fVname,
656 fVinit, fStep, fLimlo, fLimup, fFix,
657 this, fMethod, fNulloutput);
658
659 *fLog << dbg << "after calling CallMinuit" << endl;
660 *fLog << dbg << "Time spent for the minimization in MINUIT : " << endl;;
661 clock.Stop();
662 clock.Print();
663
664 Double_t integratedCharge;
665 Double_t maxFit, maxFitError;
666 Double_t meanXFit, meanXFitError;
667 Double_t sigmaMinorAxis, sigmaMinorAxisError;
668 Double_t meanYFit, meanYFitError;
669 Double_t sigmaMajorAxis, sigmaMajorAxisError;
670 Float_t chisquare = GetChisquare();
671 Int_t dregrees = GetDegreesofFreedom()-fNumVar;
672
673 if (rc)
674 {
675 gMinuit->GetParameter(0,maxFit, maxFitError);
676 gMinuit->GetParameter(1,meanXFit,meanXFitError);
677 gMinuit->GetParameter(2,sigmaMinorAxis,sigmaMinorAxisError);
678 gMinuit->GetParameter(3,meanYFit,meanYFitError);
679 gMinuit->GetParameter(4,sigmaMajorAxis,sigmaMajorAxisError);
680
681 //FIXME: Do the integral properlly
682 integratedCharge = 0.;
683
684
685 }
686 else
687 {
688 maxFit = 0.;
689 meanXFit = 0.;
690 sigmaMinorAxis = 0.;
691 meanYFit = 0.;
692 sigmaMajorAxis = 0.;
693 integratedCharge = 0.;
694 }
695
696
697
698 star->SetFitValues(integratedCharge,maxFit,meanXFit,meanYFit,sigmaMinorAxis,sigmaMajorAxis,chisquare,dregrees);
699
700 // reset the display to the starting values
701 fDisplay.SetUsed(origPixelsUsed);
702
703 return rc;
704}
705
706Bool_t MFindStars::ShadowStar(MStarLocalPos* star)
707{
708 UInt_t numPixels = fGeomCam->GetNumPixels();
709
710// Define an area around the star which will be set unused.
711 UInt_t shadowPx=0;
712 for (UInt_t pix=1; pix<numPixels; pix++)
713 {
714 Float_t pixXpos = (*fGeomCam)[pix].GetX();
715 Float_t pixYpos = (*fGeomCam)[pix].GetY();
716 Float_t starXpos = star->GetMeanX();
717 Float_t starYpos = star->GetMeanY();
718
719 Float_t starSize = 3*star->GetSigmaMajorAxis();
720
721 Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
722 (pixYpos-starYpos)*(pixYpos-starYpos));
723
724 if (dist > starSize && fDisplay.IsUsed(pix))
725 fPixelsUsed[pix]=(Char_t)kTRUE;
726 else
727 {
728 fPixelsUsed[pix]=(Char_t)kFALSE;
729 shadowPx++;
730 }
731 }
732
733 *fLog << dbg << " shadowPx " << shadowPx << endl;
734
735 fDisplay.SetUsed(fPixelsUsed);
736
737 return kTRUE;
738}
739
740
741
Note: See TracBrowser for help on using the repository browser.