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

Last change on this file since 4486 was 4429, checked in by jlopez, 21 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 }
520 else if (TMath::Abs(ped-maxprobdc) > rmsguess || rms > rmsguess)
521 {
522 *fLog << warn << GetName() << " Pedestal DC fit give non consistent results for " << (i==0?"Inner":"Outer") << "pixels." << endl;
523 *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
524 *fLog << warn << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
525 ped = maxprobdc;
526 rms = rmsguess/1.175; // FWHM=2.35*rms
527 }
528
529// TCanvas *c1 = new TCanvas("c2","c2",500,800);
530// dchist[i]->Draw();
531// c1->Print("dchist.ps");
532// getchar();
533// delete c1;
534// exit(1);
535
536 if (i == 0)
537 {
538 fStars->SetInnerPedestalDC(ped);
539 fStars->SetInnerPedestalRMSDC(rms);
540 }
541 else
542 {
543 fStars->SetOuterPedestalDC(ped);
544 fStars->SetOuterPedestalRMSDC(rms);
545 }
546
547
548
549 }
550
551
552 for (UInt_t i=0; i<2; i++)
553 delete dchist[i];
554 delete [] dchist;
555
556 //=================================================================
557
558 // reset gMinuit to the MINUIT object for optimizing the supercuts
559 gMinuit = savePointer;
560 //-------------------------------------------
561
562 if (fStars->GetInnerPedestalDC() < 0. || fStars->GetInnerPedestalRMSDC() < 0. || fStars->GetOuterPedestalDC() < 0. || fStars->GetOuterPedestalRMSDC() < 0.)
563 return kFALSE;
564
565 return kTRUE;
566}
567
568Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
569{
570 UInt_t numPixels = fGeomCam->GetNumPixels();
571
572// Find the two close pixels with the maximun dc
573 UInt_t maxPixIdx[2];
574
575 maxDC = 0;
576
577 for (UInt_t pix=1; pix<numPixels; pix++)
578 {
579 if(fDisplay.IsUsed(pix))
580 {
581 Float_t dc[2];
582 dc[0] = fDisplay.GetBinContent(pix+1);
583 if (dc[0] < fMinDCForStars)
584 continue;
585
586 MGeomPix &g = (*fGeomCam)[pix];
587 Int_t numNextNeighbors = g.GetNumNeighbors();
588
589 Float_t dcsum;
590 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
591 {
592 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
593 if(fDisplay.IsUsed(swneighbor))
594 {
595 dc[1] = fDisplay.GetBinContent(swneighbor+1);
596 if (dc[1] < fMinDCForStars)
597 continue;
598
599 dcsum = dc[0] + dc[1];
600
601 if(dcsum > maxDC*2)
602 {
603 if(dc[0]>=dc[1])
604 {
605 maxPixIdx[0] = pix;
606 maxPixIdx[1] = swneighbor;
607 maxDC = dc[0];
608 }
609 else
610 {
611 maxPixIdx[1] = pix;
612 maxPixIdx[0] = swneighbor;
613 maxDC = dc[1];
614 }
615 }
616 }
617 }
618 }
619 }
620
621 if (maxDC == 0)
622 {
623 *fLog << warn << " No found pixels with maximum dc" << endl;
624 return kFALSE;
625 }
626
627 maxPix = (*fGeomCam)[maxPixIdx[0]];
628
629 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;
630
631 return kTRUE;
632}
633
634Bool_t MFindStars::FindStar(MStarLocalPos* star)
635{
636
637 UInt_t numPixels = fGeomCam->GetNumPixels();
638 Float_t innerped = fStars->GetInnerPedestalDC();
639 Float_t outerped = fStars->GetOuterPedestalDC();
640
641 TArrayC origPixelsUsed;
642 origPixelsUsed.Set(numPixels);
643
644 for (UInt_t pix=1; pix<numPixels; pix++)
645 {
646 if (fDisplay.IsUsed(pix))
647 origPixelsUsed[pix]=(Char_t)kTRUE;
648 else
649 origPixelsUsed[pix]=(Char_t)kFALSE;
650 }
651
652 Float_t expX = star->GetXExp();
653 Float_t expY = star->GetYExp();
654
655 Float_t max=0;
656 UInt_t pixmax=0;
657 Float_t meanX=0;
658 Float_t meanY=0;
659 Float_t meanSqX=0;
660 Float_t meanSqY=0;
661 Float_t sumCharge=0;
662 UInt_t usedInnerPx=0;
663 UInt_t usedOuterPx=0;
664
665 const Float_t meanPresition = 3.; //[mm]
666 const UInt_t maxNumIterations = 10;
667 UInt_t numIterations = 0;
668
669 do
670 {
671 // First define a area of interest around the expected position of the star
672 for (UInt_t pix=1; pix<numPixels; pix++)
673 {
674
675 Float_t pixXpos=(*fGeomCam)[pix].GetX();
676 Float_t pixYpos=(*fGeomCam)[pix].GetY();
677 Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
678 (pixYpos-expY)*(pixYpos-expY));
679
680 if ((dist < fRingInterest) && fDisplay.IsUsed(pix))
681 fPixelsUsed[pix]=(Char_t)kTRUE;
682 else
683 fPixelsUsed[pix]=(Char_t)kFALSE;
684 }
685
686 fDisplay.SetUsed(fPixelsUsed);
687
688// determine mean x and mean y
689 usedInnerPx=0;
690 usedOuterPx=0;
691 for(UInt_t pix=0; pix<numPixels; pix++)
692 {
693 if(fDisplay.IsUsed(pix))
694 {
695 pix>lastInnerPixel?usedOuterPx++:usedInnerPx++;
696
697 Float_t charge = fDisplay.GetBinContent(pix+1);
698 Float_t pixXpos = (*fGeomCam)[pix].GetX();
699 Float_t pixYpos = (*fGeomCam)[pix].GetY();
700
701 if (charge>max)
702 {
703 max=charge;
704 pixmax=pix;
705 }
706
707 meanX += charge*pixXpos;
708 meanY += charge*pixYpos;
709 meanSqX += charge*pixXpos*pixXpos;
710 meanSqY += charge*pixYpos*pixYpos;
711 sumCharge += charge;
712 }
713 }
714
715 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " usedInnerPx " << usedInnerPx << " usedOuterPx " << usedOuterPx << endl;
716
717 meanX /= sumCharge;
718 meanY /= sumCharge;
719 meanSqX /= sumCharge;
720 meanSqY /= sumCharge;
721
722 expX = meanX;
723 expY = meanY;
724
725 if (++numIterations > maxNumIterations)
726 {
727 *fLog << warn << GetName() << "Mean calculation not converge after " << maxNumIterations << " iterations" << endl;
728 break;
729 }
730
731 }while(TMath::Abs(meanX-expX) > meanPresition || TMath::Abs(meanY-expY) > meanPresition);
732
733 Float_t rmsX = (meanSqX - meanX*meanX) - fRingInterest*fRingInterest/12;
734 Float_t rmsY = (meanSqY - meanY*meanY) - fRingInterest*fRingInterest/12;
735
736 if ( rmsX > 0)
737 rmsX = TMath::Sqrt(rmsX);
738 else
739 {
740 *fLog << warn << " MFindStars::FindStar negative rmsX² " << rmsX << endl;
741 *fLog << warn << " meanSqX " << meanSqX << " meanX " << meanX << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge << endl;
742 rmsX = 0.;
743 }
744
745 if ( rmsY > 0)
746 rmsY = TMath::Sqrt(rmsY);
747 else
748 {
749 *fLog << warn << " MFindStars::FindStar negative rmsY² " << rmsY << endl;
750 *fLog << warn << " meanSqY " << meanSqY << " meanY " << meanY << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge<< endl;
751 rmsY = 0.;
752 }
753
754 // Substrack pedestal DC
755 sumCharge-= (usedInnerPx*innerped+usedOuterPx*outerped)/(usedInnerPx+usedOuterPx);
756 max-=pixmax>lastInnerPixel?outerped:innerped;
757
758
759 star->SetCalcValues(sumCharge,max,meanX,meanY,rmsX,rmsY);
760
761 if (rmsX <= 0. || rmsY <= 0.)
762 return kFALSE;
763
764
765// fit the star spot using TMinuit
766
767
768 for (UInt_t pix=1; pix<numPixels; pix++)
769 if (fDisplay.IsUsed(pix))
770 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl;
771
772 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars " << fMinDCForStars << " pixmax>lastInnerPixel?outerped:innerped " << (pixmax>lastInnerPixel?outerped:innerped) << " fMaxNumIntegratedEvents " << fMaxNumIntegratedEvents << endl;
773
774 //Initialate variables for fit
775 fVinit[0] = max;
776 fLimlo[0] = fMinDCForStars-(pixmax>lastInnerPixel?outerped:innerped);
777 fLimup[0] = 30*fMaxNumIntegratedEvents-(pixmax>lastInnerPixel?outerped:innerped);
778 fVinit[1] = meanX;
779 fVinit[2] = rmsX;
780 fVinit[3] = meanY;
781 fVinit[4] = rmsY;
782 //Init steps
783 for(Int_t i=0; i<fNumVar; i++)
784 {
785 if (fVinit[i] != 0)
786 fStep[i] = TMath::Abs(fVinit[i]/sqrt2);
787 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fVinit[" << i << "] " << fVinit[i];
788 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fStep[" << i << "] " << fStep[i];
789 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimlo[" << i << "] " << fLimlo[i];
790 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimup[" << i << "] " << fLimup[i] << endl;
791 }
792 //
793
794 // -------------------------------------------
795 // call MINUIT
796
797 Double_t arglist[10];
798 Int_t ierflg = 0;
799
800 for (Int_t i=0; i<fNumVar; i++)
801 gMinuit->mnparm(i, fVname[i], fVinit[i], fStep[i], fLimlo[i], fLimup[i], ierflg);
802
803 TStopwatch clock;
804 clock.Start();
805
806// Now ready for minimization step
807 arglist[0] = 500;
808 arglist[1] = 1.;
809 gMinuit->mnexcm(fMethod, arglist ,2,ierflg);
810
811 clock.Stop();
812
813 if(fMinuitPrintOutLevel>=0)
814 {
815 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "Time spent for the minimization in MINUIT : " << endl;;
816 clock.Print();
817 }
818
819 Double_t integratedCharge;
820 Double_t maxFit, maxFitError;
821 Double_t meanXFit, meanXFitError;
822 Double_t sigmaMinorAxis, sigmaMinorAxisError;
823 Double_t meanYFit, meanYFitError;
824 Double_t sigmaMajorAxis, sigmaMajorAxisError;
825 Float_t chisquare = GetChisquare();
826 Int_t dregrees = GetDegreesofFreedom()-fNumVar;
827
828 if (!ierflg)
829 {
830 gMinuit->GetParameter(0,maxFit, maxFitError);
831 gMinuit->GetParameter(1,meanXFit,meanXFitError);
832 gMinuit->GetParameter(2,sigmaMinorAxis,sigmaMinorAxisError);
833 gMinuit->GetParameter(3,meanYFit,meanYFitError);
834 gMinuit->GetParameter(4,sigmaMajorAxis,sigmaMajorAxisError);
835
836 //FIXME: Do the integral properlly
837 integratedCharge = 0.;
838
839
840 }
841 else
842 {
843 maxFit = 0.;
844 meanXFit = 0.;
845 sigmaMinorAxis = 0.;
846 meanYFit = 0.;
847 sigmaMajorAxis = 0.;
848 integratedCharge = 0.;
849
850 *fLog << err << "TMinuit::Call error " << ierflg << endl;
851 }
852
853 star->SetFitValues(integratedCharge,maxFit,meanXFit,meanYFit,sigmaMinorAxis,sigmaMajorAxis,chisquare,dregrees);
854
855 // reset the display to the starting values
856 fDisplay.SetUsed(origPixelsUsed);
857
858 if (ierflg)
859 return kCONTINUE;
860 return kTRUE;
861}
862
863Bool_t MFindStars::ShadowStar(MStarLocalPos* star)
864{
865 UInt_t numPixels = fGeomCam->GetNumPixels();
866
867// Define an area around the star which will be set unused.
868 UInt_t shadowPx=0;
869 for (UInt_t pix=1; pix<numPixels; pix++)
870 {
871 Float_t pixXpos = (*fGeomCam)[pix].GetX();
872 Float_t pixYpos = (*fGeomCam)[pix].GetY();
873 Float_t starXpos = star->GetMeanX();
874 Float_t starYpos = star->GetMeanY();
875
876 Float_t starSize = 3*star->GetSigmaMajorAxis();
877
878 Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
879 (pixYpos-starYpos)*(pixYpos-starYpos));
880
881 if (dist > starSize && fDisplay.IsUsed(pix))
882 fPixelsUsed[pix]=(Char_t)kTRUE;
883 else
884 {
885 fPixelsUsed[pix]=(Char_t)kFALSE;
886 shadowPx++;
887 }
888 }
889
890 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " shadowPx " << shadowPx << endl;
891
892 fDisplay.SetUsed(fPixelsUsed);
893
894 return kTRUE;
895}
896
897
898
899
900
Note: See TracBrowser for help on using the repository browser.