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

Last change on this file since 4211 was 4071, checked in by jlopez, 21 years ago
*** empty log message ***
File size: 22.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 "MLog.h"
46#include "MLogManip.h"
47
48#include "MHCamera.h"
49
50#include "MGeomCam.h"
51#include "MGeomPix.h"
52#include "MCameraDC.h"
53#include "MTime.h"
54#include "MReportDrive.h"
55#include "MStarLocalCam.h"
56#include "MStarLocalPos.h"
57
58#include "MParList.h"
59
60ClassImp(MFindStars);
61using namespace std;
62
63const Float_t sqrt2 = sqrt(2.);
64const Float_t sqrt3 = sqrt(3.);
65
66//______________________________________________________________________________
67//
68// The 2D gaussian fucntion used to fit the spot of the star
69//
70static Double_t func(float x,float y,Double_t *par)
71{
72 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]));
73 return value;
74}
75
76//______________________________________________________________________________
77//
78// Function used by Minuit to do the fit
79//
80static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
81{
82
83 MFindStars* find = (MFindStars*)gMinuit->GetObjectFit();
84 MHCamera& display = (MHCamera&)find->GetDisplay();
85 Float_t ped = find->GetPedestalDC();
86 Float_t rms = find->GetPedestalRMSDC();
87 MGeomCam& geom = (MGeomCam&)display.GetGeomCam();
88
89 UInt_t numPixels = geom.GetNumPixels();
90
91//calculate chisquare
92 Double_t chisq = 0;
93 Double_t delta;
94 Double_t x,y,z;
95 Double_t errorz = rms; //[uA]
96
97 UInt_t usedPx=0;
98 for (UInt_t pixid=1; pixid<numPixels; pixid++)
99 {
100 z = display.GetBinContent(pixid+1)-ped;
101
102 if (display.IsUsed(pixid))
103 {
104 x = geom[pixid].GetX();
105 y = geom[pixid].GetY();
106
107 if (errorz > 0.0)
108 {
109 usedPx++;
110 delta = (z-func(x,y,par))/errorz;
111 chisq += delta*delta;
112 }
113 else
114 cerr << " TMinuit::fcn errorz[" << pixid << "] " << errorz << endl;
115 }
116 }
117 f = chisq;
118
119 find->SetChisquare(chisq);
120 find->SetDegreesofFreedom(usedPx);
121}
122
123MFindStars::MFindStars(const char *name, const char *title): fNumVar(5)
124{
125 fName = name ? name : "MFindStars";
126 fTitle = title ? title : "Tool to find stars from DC Currents";
127
128 fNumIntegratedEvents=0;
129 fMaxNumIntegratedEvents = 10;
130 fRingInterest = 125.; //[mm] ~ 0.4 deg
131 fDCTailCut = 4;
132
133 fPixelsUsed.Set(577);
134 fPixelsUsed.Reset((Char_t)kTRUE);
135
136 //Fitting(Minuit) initialitation
137 const Float_t pixelSize = 31.5; //[mm]
138 fMinuitPrintOutLevel = -1;
139
140 fVname = new TString[fNumVar];
141 fVinit.Set(fNumVar);
142 fStep.Set(fNumVar);
143 fLimlo.Set(fNumVar);
144 fLimup.Set(fNumVar);
145 fFix.Set(fNumVar);
146
147 fVname[0] = "max";
148 fVinit[0] = 10.*fMaxNumIntegratedEvents;
149 fStep[0] = fVinit[0]/sqrt2;
150 fLimlo[0] = fMinDCForStars;
151 fLimup[0] = 30.*fMaxNumIntegratedEvents;
152 fFix[0] = 0;
153
154 fVname[1] = "meanx";
155 fVinit[1] = 0.;
156 fStep[1] = fVinit[1]/sqrt2;
157 fLimlo[1] = -600.;
158 fLimup[1] = 600.;
159 fFix[1] = 0;
160
161 fVname[2] = "sigmaminor";
162 fVinit[2] = pixelSize;
163 fStep[2] = fVinit[2]/sqrt2;
164 fLimlo[2] = pixelSize/(2*sqrt3);
165 fLimup[2] = 500.;
166 fFix[2] = 0;
167
168 fVname[3] = "meany";
169 fVinit[3] = 0.;
170 fStep[3] = fVinit[3]/sqrt2;
171 fLimlo[3] = -600.;
172 fLimup[3] = 600.;
173 fFix[3] = 0;
174
175 fVname[4] = "sigmamajor";
176 fVinit[4] = pixelSize;
177 fStep[4] = fVinit[4]/sqrt2;
178 fLimlo[4] = pixelSize/(2*sqrt3);
179 fLimup[4] = 500.;
180 fFix[4] = 0;
181
182 fObjectFit = NULL;
183 // fMethod = "SIMPLEX";
184 fMethod = "MIGRAD";
185 // fMethod = "MINIMIZE";
186 fNulloutput = kFALSE;
187}
188
189Int_t MFindStars::PreProcess(MParList *pList)
190{
191
192 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
193
194 if (!fGeomCam)
195 {
196 *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl;
197 return kFALSE;
198 }
199
200 // Initialize camera display with the MGeomCam information
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 fMinDCForStars = 1.5*fMaxNumIntegratedEvents; //[uA]
285
286 // Initialize the TMinuit object
287
288 TMinuit *gMinuit = new TMinuit(fNumVar); //initialize TMinuit with a maximum of params
289 gMinuit->SetFCN(fcn);
290
291 Double_t arglist[10];
292 Int_t ierflg = 0;
293
294 arglist[0] = 1;
295 gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
296 arglist[0] = fMinuitPrintOutLevel;
297 gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
298
299 // Set object pointer to allow mimuit to access internally
300 gMinuit->SetObjectFit(this);
301
302 return kTRUE;
303}
304
305Int_t MFindStars::Process()
306{
307 UInt_t numPixels = fGeomCam->GetNumPixels();
308 TArrayC origPixelsUsed;
309 origPixelsUsed.Set(numPixels);
310
311 if (fNumIntegratedEvents >= fMaxNumIntegratedEvents)
312 {
313
314 if (fDrive)
315 {
316 Float_t ra = fDrive->GetRa();
317 Float_t dec = fDrive->GetDec();
318
319 fAstro.SetRaDec(ra, dec);
320 fAstro.SetGuiActive();
321
322 fAstro.FillStarList();
323 }
324 else
325 {
326 //Fist delete the previus stars in the list
327 fStars->GetList()->Delete();
328 //
329
330 for (UInt_t pix=1; pix<numPixels; pix++)
331 {
332 if (fDisplay.IsUsed(pix))
333 origPixelsUsed[pix]=(Char_t)kTRUE;
334 else
335 origPixelsUsed[pix]=(Char_t)kFALSE;
336 }
337
338 DCPedestalCalc(fPedestalDC, fPedestalRMSDC);
339 fMinDCForStars = fMinDCForStars>(fPedestalDC+fDCTailCut*fPedestalRMSDC)?fMinDCForStars:(fPedestalDC+fDCTailCut*fPedestalRMSDC);
340
341 *fLog << dbg << " DC pedestal = " << fPedestalDC << " pedestal rms = " << fPedestalRMSDC << endl;
342 *fLog << dbg << " fMinDCForStars " << fMinDCForStars << endl;
343
344 // Find the star candidats searching the most brights pairs of pixels
345 Float_t maxPixelDC;
346 MGeomPix maxPixel;
347
348 while(FindPixelWithMaxDC(maxPixelDC, maxPixel))
349 {
350
351 MStarLocalPos *starpos = new MStarLocalPos;
352 starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
353 starpos->SetCalcValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
354 starpos->SetFitValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2,0.,1);
355 fStars->GetList()->Add(starpos);
356
357 ShadowStar(starpos);
358
359 }
360
361 fDisplay.SetUsed(origPixelsUsed);
362 }
363
364 //loop to extract position of stars on the camera
365 if (fStars->GetList()->GetSize() == 0)
366 {
367 *fLog << err << GetName() << " No stars candidates in the camera." << endl;
368 return kFALSE;
369 }
370 else
371 *fLog << inf << GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
372
373
374 for (UInt_t pix=1; pix<numPixels; pix++)
375 {
376 if (fDisplay.IsUsed(pix))
377 origPixelsUsed[pix]=(Char_t)kTRUE;
378 else
379 origPixelsUsed[pix]=(Char_t)kFALSE;
380 }
381
382 TIter Next(fStars->GetList());
383 MStarLocalPos* star;
384 while ((star=(MStarLocalPos*)Next()))
385 {
386 FindStar(star);
387 ShadowStar(star);
388 }
389
390 //After finding stars reset all vairables
391 fDisplay.Reset();
392 fDisplay.SetUsed(origPixelsUsed);
393 fNumIntegratedEvents=0;
394 }
395
396 for (UInt_t pix=1; pix<numPixels; pix++)
397 {
398 if (fDisplay.IsUsed(pix))
399 origPixelsUsed[pix]=(Char_t)kTRUE;
400 else
401 origPixelsUsed[pix]=(Char_t)kFALSE;
402
403 }
404
405 fDisplay.AddCamContent(*fCurr);
406 fNumIntegratedEvents++;
407 fDisplay.SetUsed(origPixelsUsed);
408
409
410 return kTRUE;
411}
412
413Int_t MFindStars::PostProcess()
414{
415 if(fStars->GetList()->GetSize() != 0)
416 fStars->Print();
417
418 return kTRUE;
419}
420
421void MFindStars::SetBlindPixels(TArrayS blindpixels)
422{
423 Int_t npix = blindpixels.GetSize();
424
425 for (Int_t idx=0; idx<npix; idx++)
426 {
427 fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;
428 *fLog << dbg << "MFindStars::SetBlindPixels fDisplay.IsUsed(" <<blindpixels[idx] << ") kFALSE" << endl;
429 }
430
431 fDisplay.SetUsed(fPixelsUsed);
432}
433
434Bool_t MFindStars::DCPedestalCalc(Float_t &ped, Float_t &rms)
435{
436 //-------------------------------------------------------------
437 // save pointer to the MINUIT object for optimizing the supercuts
438 // because it will be overwritten
439 // when fitting the alpha distribution in MHFindSignificance
440 TMinuit *savePointer = gMinuit;
441 //-------------------------------------------------------------
442
443 UInt_t numPixels = fGeomCam->GetNumPixels();
444
445 TH1F dchist("dchist","",120,0.,30.*fMaxNumIntegratedEvents);
446 for (UInt_t pix=1; pix<numPixels; pix++)
447 dchist.Fill(fDisplay.GetBinContent(pix+1));
448
449 Float_t nummaxprobdc = dchist.GetBinContent(dchist.GetMaximumBin());
450 Float_t maxprobdc = dchist.GetBinCenter(dchist.GetMaximumBin());
451 UInt_t bin = dchist.GetMaximumBin();
452 do
453 {
454 bin++;
455 }
456 while(dchist.GetBinContent(bin)/nummaxprobdc > 0.5);
457 Float_t halfmaxprobdc = dchist.GetBinCenter(bin);
458
459 *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] ";
460 *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << dchist.GetBinContent(bin) << "]" << endl;
461
462 Float_t rmsguess = TMath::Abs(maxprobdc-halfmaxprobdc);
463 Float_t min = maxprobdc-3*rmsguess;
464 min = (min<0.?0.:min);
465 Float_t max = maxprobdc+3*rmsguess;
466
467 *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
468
469 TF1 func("func","gaus",min,max);
470 func.SetParameters(nummaxprobdc, maxprobdc, rmsguess);
471
472 dchist.Fit("func","QR0");
473 // Remove the comments if you want to go through the file
474 // event-by-event:
475 // HandleInput();
476
477 UInt_t aproxnumdegrees = 6*(bin-dchist.GetMaximumBin());
478 Float_t chiq = func.GetChisquare();
479 ped = func.GetParameter(1);
480 rms = func.GetParameter(2);
481
482 *fLog << dbg << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
483
484 Int_t numsigmas = 5;
485 Axis_t minbin = ped-numsigmas*rms/dchist.GetBinWidth(1);
486 minbin=minbin<1?1:minbin;
487 Axis_t maxbin = ped+numsigmas*rms/dchist.GetBinWidth(1);
488 *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << dchist.Integral((int)minbin,(int)maxbin) << endl;
489
490 //Check results from the fit are consistent
491 if (TMath::Abs(ped-maxprobdc) > rmsguess || rms > rmsguess)
492 {
493 *fLog << warn << GetName() << " Pedestal DC fit give non consistent results." << endl;
494 *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
495 *fLog << warn << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
496 ped = maxprobdc;
497 rms = rmsguess/1.175; // FWHM=2.35*rms
498 }
499
500 //=================================================================
501
502 // reset gMinuit to the MINUIT object for optimizing the supercuts
503 gMinuit = savePointer;
504 //-------------------------------------------
505
506 return kTRUE;
507}
508
509Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
510{
511 UInt_t numPixels = fGeomCam->GetNumPixels();
512
513// Find the two close pixels with the maximun dc
514 UInt_t maxPixIdx[2];
515
516 maxDC = 0;
517
518 for (UInt_t pix=1; pix<numPixels; pix++)
519 {
520 if(fDisplay.IsUsed(pix))
521 {
522 Float_t dc[2];
523 dc[0] = fDisplay.GetBinContent(pix+1);
524 if (dc[0] < fMinDCForStars)
525 continue;
526
527 MGeomPix &g = (*fGeomCam)[pix];
528 Int_t numNextNeighbors = g.GetNumNeighbors();
529
530 Float_t dcsum;
531 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
532 {
533 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
534 if(fDisplay.IsUsed(swneighbor))
535 {
536 dc[1] = fDisplay.GetBinContent(swneighbor+1);
537 if (dc[1] < fMinDCForStars)
538 continue;
539
540 dcsum = dc[0] + dc[1];
541
542 if(dcsum > maxDC*2)
543 {
544 if(dc[0]>=dc[1])
545 {
546 maxPixIdx[0] = pix;
547 maxPixIdx[1] = swneighbor;
548 maxDC = dc[0];
549 }
550 else
551 {
552 maxPixIdx[1] = pix;
553 maxPixIdx[0] = swneighbor;
554 maxDC = dc[1];
555 }
556 }
557 }
558 }
559 }
560 }
561
562 if (maxDC == 0)
563 return kFALSE;
564
565 maxPix = (*fGeomCam)[maxPixIdx[0]];
566
567 *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;
568
569 return kTRUE;
570}
571
572Bool_t MFindStars::FindStar(MStarLocalPos* star)
573{
574
575 UInt_t numPixels = fGeomCam->GetNumPixels();
576 TArrayC origPixelsUsed;
577 origPixelsUsed.Set(numPixels);
578
579 for (UInt_t pix=1; pix<numPixels; pix++)
580 {
581 if (fDisplay.IsUsed(pix))
582 origPixelsUsed[pix]=(Char_t)kTRUE;
583 else
584 origPixelsUsed[pix]=(Char_t)kFALSE;
585 }
586
587 Float_t expX = star->GetXExp();
588 Float_t expY = star->GetYExp();
589
590 Float_t max=0;
591 Float_t meanX=0;
592 Float_t meanY=0;
593 Float_t meanSqX=0;
594 Float_t meanSqY=0;
595 Float_t sumCharge=0;
596 UInt_t usedPx=0;
597
598 const Float_t meanPresition = 3.; //[mm]
599 const UInt_t maxNumIterations = 10;
600 UInt_t numIterations = 0;
601
602 do
603 {
604 // First define a area of interest around the expected position of the star
605 for (UInt_t pix=1; pix<numPixels; pix++)
606 {
607
608 Float_t pixXpos=(*fGeomCam)[pix].GetX();
609 Float_t pixYpos=(*fGeomCam)[pix].GetY();
610 Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
611 (pixYpos-expY)*(pixYpos-expY));
612
613 if ((dist < fRingInterest) && fDisplay.IsUsed(pix))
614 fPixelsUsed[pix]=(Char_t)kTRUE;
615 else
616 fPixelsUsed[pix]=(Char_t)kFALSE;
617 }
618
619 fDisplay.SetUsed(fPixelsUsed);
620
621// determine mean x and mean y
622 usedPx=0;
623 for(UInt_t pix=0; pix<numPixels; pix++)
624 {
625 if(fDisplay.IsUsed(pix))
626 {
627 usedPx++;
628
629 Float_t charge = fDisplay.GetBinContent(pix+1);
630 Float_t pixXpos = (*fGeomCam)[pix].GetX();
631 Float_t pixYpos = (*fGeomCam)[pix].GetY();
632
633 if (charge>max) max=charge;
634 meanX += charge*pixXpos;
635 meanY += charge*pixYpos;
636 meanSqX += charge*pixXpos*pixXpos;
637 meanSqY += charge*pixYpos*pixYpos;
638 sumCharge += charge;
639 }
640 }
641
642 *fLog << dbg << " usedPx " << usedPx << endl;
643
644 meanX /= sumCharge;
645 meanY /= sumCharge;
646 meanSqX /= sumCharge;
647 meanSqY /= sumCharge;
648
649 expX = meanX;
650 expY = meanY;
651
652 if (++numIterations > maxNumIterations)
653 {
654 *fLog << warn << GetName() << "Mean calculation not converge after " << maxNumIterations << " iterations" << endl;
655 break;
656 }
657
658 }while(TMath::Abs(meanX-expX) > meanPresition || TMath::Abs(meanY-expY) > meanPresition);
659
660 Float_t rmsX = (meanSqX - meanX*meanX) - fRingInterest*fRingInterest/12;
661 Float_t rmsY = (meanSqY - meanY*meanY) - fRingInterest*fRingInterest/12;
662
663 if ( rmsX > 0)
664 rmsX = TMath::Sqrt(rmsX);
665 else
666 {
667 *fLog << warn << " MFindStars::FindStar negative rmsX² " << rmsX << endl;
668 *fLog << warn << " meanSqX " << meanSqX << " meanX " << meanX << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge << endl;
669 rmsX = 0.;
670 }
671
672 if ( rmsY > 0)
673 rmsY = TMath::Sqrt(rmsY);
674 else
675 {
676 *fLog << warn << " MFindStars::FindStar negative rmsY² " << rmsY << endl;
677 *fLog << warn << " meanSqY " << meanSqY << " meanY " << meanY << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge<< endl;
678 rmsY = 0.;
679 }
680
681 // Substrack pedestal DC
682 sumCharge-=fPedestalDC*usedPx;
683 max-=fPedestalDC;
684
685
686 star->SetCalcValues(sumCharge,max,meanX,meanY,rmsX,rmsY);
687
688 if (rmsX <= 0. || rmsY <= 0.)
689 return kFALSE;
690
691
692// fit the star spot using TMinuit
693
694
695 for (UInt_t pix=1; pix<numPixels; pix++)
696 if (fDisplay.IsUsed(pix))
697 *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl;
698
699 //Initialate variables for fit
700 fVinit[0] = max;
701 fLimlo[0] = fMinDCForStars-fPedestalDC;
702 fLimup[0] = 30*fMaxNumIntegratedEvents-fPedestalDC;
703 fVinit[1] = meanX;
704 fVinit[2] = rmsX;
705 fVinit[3] = meanY;
706 fVinit[4] = rmsY;
707 //Init steps
708 for(Int_t i=0; i<fNumVar; i++)
709 if (fVinit[i] != 0)
710 fStep[i] = TMath::Abs(fVinit[i]/sqrt2);
711 //
712
713 // -------------------------------------------
714 // call MINUIT
715
716 Double_t arglist[10];
717 Int_t ierflg = 0;
718
719 for (Int_t i=0; i<fNumVar; i++)
720 gMinuit->mnparm(i, fVname[i], fVinit[i], fStep[i], fLimlo[i], fLimup[i], ierflg);
721
722 TStopwatch clock;
723 clock.Start();
724
725// Now ready for minimization step
726 arglist[0] = 500;
727 arglist[1] = 1.;
728 gMinuit->mnexcm(fMethod, arglist ,2,ierflg);
729
730 clock.Stop();
731
732 if(fMinuitPrintOutLevel>=0)
733 {
734 *fLog << dbg << "Time spent for the minimization in MINUIT : " << endl;;
735 clock.Print();
736 }
737
738 Double_t integratedCharge;
739 Double_t maxFit, maxFitError;
740 Double_t meanXFit, meanXFitError;
741 Double_t sigmaMinorAxis, sigmaMinorAxisError;
742 Double_t meanYFit, meanYFitError;
743 Double_t sigmaMajorAxis, sigmaMajorAxisError;
744 Float_t chisquare = GetChisquare();
745 Int_t dregrees = GetDegreesofFreedom()-fNumVar;
746
747 if (!ierflg)
748 {
749 gMinuit->GetParameter(0,maxFit, maxFitError);
750 gMinuit->GetParameter(1,meanXFit,meanXFitError);
751 gMinuit->GetParameter(2,sigmaMinorAxis,sigmaMinorAxisError);
752 gMinuit->GetParameter(3,meanYFit,meanYFitError);
753 gMinuit->GetParameter(4,sigmaMajorAxis,sigmaMajorAxisError);
754
755 //FIXME: Do the integral properlly
756 integratedCharge = 0.;
757
758
759 }
760 else
761 {
762 maxFit = 0.;
763 meanXFit = 0.;
764 sigmaMinorAxis = 0.;
765 meanYFit = 0.;
766 sigmaMajorAxis = 0.;
767 integratedCharge = 0.;
768
769 *fLog << err << "TMinuit::Call error " << ierflg << endl;
770 }
771
772 star->SetFitValues(integratedCharge,maxFit,meanXFit,meanYFit,sigmaMinorAxis,sigmaMajorAxis,chisquare,dregrees);
773
774 // reset the display to the starting values
775 fDisplay.SetUsed(origPixelsUsed);
776
777 if (ierflg)
778 return kCONTINUE;
779 return kTRUE;
780}
781
782Bool_t MFindStars::ShadowStar(MStarLocalPos* star)
783{
784 UInt_t numPixels = fGeomCam->GetNumPixels();
785
786// Define an area around the star which will be set unused.
787 UInt_t shadowPx=0;
788 for (UInt_t pix=1; pix<numPixels; pix++)
789 {
790 Float_t pixXpos = (*fGeomCam)[pix].GetX();
791 Float_t pixYpos = (*fGeomCam)[pix].GetY();
792 Float_t starXpos = star->GetMeanX();
793 Float_t starYpos = star->GetMeanY();
794
795 Float_t starSize = 3*star->GetSigmaMajorAxis();
796
797 Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
798 (pixYpos-starYpos)*(pixYpos-starYpos));
799
800 if (dist > starSize && fDisplay.IsUsed(pix))
801 fPixelsUsed[pix]=(Char_t)kTRUE;
802 else
803 {
804 fPixelsUsed[pix]=(Char_t)kFALSE;
805 shadowPx++;
806 }
807 }
808
809 *fLog << dbg << " shadowPx " << shadowPx << endl;
810
811 fDisplay.SetUsed(fPixelsUsed);
812
813 return kTRUE;
814}
815
816
817
818
819
Note: See TracBrowser for help on using the repository browser.