source: trunk/MagicSoft/Mars/mtemp/mmpi/MFindStars.cc@ 8782

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