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

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