source: trunk/MagicSoft/Mars/mtemp/meth/MFindStars.cc@ 5135

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