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

Last change on this file since 4433 was 4415, checked in by stark, 20 years ago
*** empty log message ***
File size: 31.2 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 // Magnitude up to which the stars are loaded from the catalog
282 fAstroCamera->SetLimMag(6);
283 // Radius of FOV around the source position to load the stars
284 fAstroCamera->SetRadiusFOV(3);
285 // Source position
286 fAstroCamera->SetRaDec(ra, dec);
287 // Catalog to load (here: Bright Star Catalog V5)
288 fAstroCamera->ReadBSC("../catalogs/bsc5.dat");
289 // Obersavatory and time to also get local coordinate information
290 fAstroCamera->SetObservatory(magic1);
291 // fRes contains the corrected source position
292 fRes = new TVector3;
293
294 fStars = (MStarLocalCam*)pList->FindCreateObj(AddSerialNumber("MStarLocalCam"));
295 if (!fStars)
296 {
297 *fLog << err << AddSerialNumber("MStarLocalCam") << " cannot be created ... aborting" << endl;
298 return kFALSE;
299 }
300
301 fMinDCForStars = 1.*fMaxNumIntegratedEvents; //[uA]
302
303 // Initialize the TMinuit object
304
305 TMinuit *gMinuit = new TMinuit(fNumVar); //initialize TMinuit with a maximum of params
306 gMinuit->SetFCN(fcn);
307
308 Double_t arglist[2*fNumVar];//[10];
309 Int_t ierflg = 0;
310
311 arglist[0] = 1;
312 gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
313 arglist[0] = fMinuitPrintOutLevel;
314 gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
315
316 // Set MParList object pointer to allow mimuit to access internally
317 gMinuit->SetObjectFit(pList);
318
319 return kTRUE;
320}
321
322Int_t MFindStars::Process()
323{
324 UInt_t numPixels = fGeomCam->GetNumPixels();
325 TArrayC origPixelsUsed;
326 origPixelsUsed.Set(numPixels);
327 if (fNumIntegratedEvents >= fMaxNumIntegratedEvents)
328 {
329 //Fist delete the previus stars in the list
330 fAstroCamera->GetCatList()->Delete();
331 // Time for which to get the star position
332 fAstroCamera->SetTime(*fTimeCurr); // FIX ME: time should be averaged
333 //Clone the catalog due to the validity range of the instance
334 TObject *o = fAstroCamera->Clone();
335 o->SetBit(kCanDelete);
336 o->Draw("mirror");
337
338 fAstroCamera->StarPosInCamera();
339
340 if (fDrive)
341 {
342// Float_t ra = fDrive->GetRa();
343// Float_t dec = fDrive->GetDec();
344
345// fAstro.SetRaDec(ra, dec);
346// fAstro.SetGuiActive();
347
348// fAstro.FillStarList();
349 }
350 else
351 {
352 //Fist delete the previus stars in the list
353 fStars->GetList()->Delete();
354
355 for (UInt_t pix=1; pix<numPixels; pix++)
356 {
357 if (fDisplay.IsUsed(pix))
358 origPixelsUsed[pix]=(Char_t)kTRUE;
359 else
360 origPixelsUsed[pix]=(Char_t)kFALSE;
361 }
362
363 if (DCPedestalCalc())
364 {
365
366 Float_t innermin = fStars->GetInnerPedestalDC()+fDCTailCut*fStars->GetInnerPedestalRMSDC();
367 Float_t outermin = fStars->GetOuterPedestalDC()+fDCTailCut*fStars->GetOuterPedestalRMSDC();
368 fMinDCForStars = innermin>outermin?innermin:outermin;
369 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars = " << fMinDCForStars << endl;
370
371 // Find the star candidats searching the most brights pairs of pixels
372 Float_t maxPixelDC;
373 MGeomPix maxPixel;
374
375 while(FindPixelWithMaxDC(maxPixelDC, maxPixel))
376 {
377
378 MStarLocalPos *starpos = new MStarLocalPos;
379 starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
380 starpos->SetCalcValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
381 starpos->SetFitValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2,0.,1);
382 fStars->GetList()->Add(starpos);
383
384 ShadowStar(starpos);
385
386 }
387
388 fDisplay.SetUsed(origPixelsUsed);
389 }
390
391 }
392
393 //loop to extract position of stars on the camera
394 if (fStars->GetList()->GetSize() == 0)
395 {
396 *fLog << err << GetName() << " No stars candidates in the camera." << endl;
397 return kCONTINUE;
398 }
399 else
400 *fLog << inf << GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
401
402
403 for (UInt_t pix=1; pix<numPixels; pix++)
404 {
405 if (fDisplay.IsUsed(pix))
406 origPixelsUsed[pix]=(Char_t)kTRUE;
407 else
408 origPixelsUsed[pix]=(Char_t)kFALSE;
409 }
410
411 TIter Next(fStars->GetList());
412 MStarLocalPos* star;
413 while ((star=(MStarLocalPos*)Next()))
414 {
415 FindStar(star);
416 ShadowStar(star);
417 }
418
419 //Get the corrected source position
420 TVector3 *res = new TVector3;
421 if (CorrSourcePos(res))
422 fRes = res;
423 else
424 fRes->SetXYZ(-1000, -1000, -1000);
425 //After finding stars reset all vairables
426 fDisplay.Reset();
427 fStars->GetDisplay().Reset(); //FIXME: Put this display just in the container
428 fDisplay.SetUsed(origPixelsUsed);
429 fNumIntegratedEvents=0;
430 }
431
432 for (UInt_t pix=1; pix<numPixels; pix++)
433 {
434 if (fDisplay.IsUsed(pix))
435 origPixelsUsed[pix]=(Char_t)kTRUE;
436 else
437 origPixelsUsed[pix]=(Char_t)kFALSE;
438
439 }
440
441 fDisplay.AddCamContent(*fCurr);
442 fNumIntegratedEvents++;
443 fDisplay.SetUsed(origPixelsUsed);
444
445
446 return kTRUE;
447}
448
449Int_t MFindStars::PostProcess()
450{
451 return kTRUE;
452}
453
454void MFindStars::SetBlindPixels(TArrayS blindpixels)
455{
456 Int_t npix = blindpixels.GetSize();
457
458 for (Int_t idx=0; idx<npix; idx++)
459 {
460 fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;
461 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "MFindStars::SetBlindPixels fDisplay.IsUsed(" <<blindpixels[idx] << ") kFALSE" << endl;
462 }
463
464 fDisplay.SetUsed(fPixelsUsed);
465}
466
467Bool_t MFindStars::DCPedestalCalc()
468{
469 //-------------------------------------------------------------
470 // save pointer to the MINUIT object for optimizing the supercuts
471 // because it will be overwritten
472 // when fitting the alpha distribution in MHFindSignificance
473 TMinuit *savePointer = gMinuit;
474 //-------------------------------------------------------------
475
476 UInt_t numPixels = fGeomCam->GetNumPixels();
477 Float_t ped;
478 Float_t rms;
479
480 TH1F **dchist = new TH1F*[2];
481 for (UInt_t i=0; i<2; i++)
482 dchist[i] = new TH1F("","",26,0.4*fMaxNumIntegratedEvents,3.*fMaxNumIntegratedEvents);
483
484 for (UInt_t pix=1; pix<=lastInnerPixel; pix++)
485 dchist[0]->Fill(fDisplay.GetBinContent(pix+1));
486 for (UInt_t pix=lastInnerPixel+1; pix<numPixels; pix++)
487 dchist[1]->Fill(fDisplay.GetBinContent(pix+1));
488
489 // inner/outer pixels
490 for (UInt_t i=0; i<2; i++)
491 {
492 Float_t nummaxprobdc = dchist[i]->GetBinContent(dchist[i]->GetMaximumBin());
493 Float_t maxprobdc = dchist[i]->GetBinCenter(dchist[i]->GetMaximumBin());
494 UInt_t bin = dchist[i]->GetMaximumBin();
495 do
496 {
497 bin++;
498 }
499 while(dchist[i]->GetBinContent(bin)/nummaxprobdc > 0.5);
500 Float_t halfmaxprobdc = dchist[i]->GetBinCenter(bin);
501
502 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] ";
503 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << dchist[i]->GetBinContent(bin) << "]" << endl;
504
505 Float_t rmsguess = TMath::Abs(maxprobdc-halfmaxprobdc);
506 Float_t min = maxprobdc-3*rmsguess;
507 min = (min<0.?0.:min);
508 Float_t max = maxprobdc+3*rmsguess;
509
510 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
511
512 TF1 func("func","gaus",min,max);
513 func.SetParameters(nummaxprobdc, maxprobdc, rmsguess);
514
515 dchist[i]->Fit("func","QR0");
516
517 UInt_t aproxnumdegrees = 6*(bin-dchist[i]->GetMaximumBin());
518 Float_t chiq = func.GetChisquare();
519 ped = func.GetParameter(1);
520 rms = func.GetParameter(2);
521
522 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
523
524 Int_t numsigmas = 5;
525 Axis_t minbin = ped-numsigmas*rms/dchist[i]->GetBinWidth(1);
526 minbin=minbin<1?1:minbin;
527 Axis_t maxbin = ped+numsigmas*rms/dchist[i]->GetBinWidth(1);
528 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << dchist[i]->Integral((int)minbin,(int)maxbin) << endl;
529
530 //Check results from the fit are consistent
531 if (ped < 0. || rms < 0.)
532 {
533 *fLog << dbg << "dchist[i]->GetEntries()" << dchist[i]->GetEntries();
534// TCanvas *c1 = new TCanvas("c2","c2",500,800);
535// dchist[i]->Draw();
536// c1->Print("dchist.ps");
537// delete c1;
538// exit(1);
539 }
540 else if (TMath::Abs(ped-maxprobdc) > rmsguess || rms > rmsguess)
541 {
542 *fLog << warn << GetName() << " Pedestal DC fit give non consistent results for " << (i==0?"Inner":"Outer") << "pixels." << endl;
543 *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
544 *fLog << warn << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
545 ped = maxprobdc;
546 rms = rmsguess/1.175; // FWHM=2.35*rms
547 }
548
549 if (i == 0)
550 {
551 fStars->SetInnerPedestalDC(ped);
552 fStars->SetInnerPedestalRMSDC(rms);
553 }
554 else
555 {
556 fStars->SetOuterPedestalDC(ped);
557 fStars->SetOuterPedestalRMSDC(rms);
558 }
559
560
561
562 }
563
564
565 for (UInt_t i=0; i<2; i++)
566 delete dchist[i];
567 delete [] dchist;
568
569 //=================================================================
570
571 // reset gMinuit to the MINUIT object for optimizing the supercuts
572 gMinuit = savePointer;
573 //-------------------------------------------
574
575 if (fStars->GetInnerPedestalDC() < 0. || fStars->GetInnerPedestalRMSDC() < 0. || fStars->GetOuterPedestalDC() < 0. || fStars->GetOuterPedestalRMSDC() < 0.)
576 return kFALSE;
577
578 return kTRUE;
579}
580
581Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
582{
583 UInt_t numPixels = fGeomCam->GetNumPixels();
584
585// Find the two close pixels with the maximun dc
586 UInt_t maxPixIdx[2];
587
588 maxDC = 0;
589
590 for (UInt_t pix=1; pix<numPixels; pix++)
591 {
592 if(fDisplay.IsUsed(pix))
593 {
594 Float_t dc[2];
595 dc[0] = fDisplay.GetBinContent(pix+1);
596 if (dc[0] < fMinDCForStars)
597 continue;
598
599 MGeomPix &g = (*fGeomCam)[pix];
600 Int_t numNextNeighbors = g.GetNumNeighbors();
601
602 Float_t dcsum;
603 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
604 {
605 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
606 if(fDisplay.IsUsed(swneighbor))
607 {
608 dc[1] = fDisplay.GetBinContent(swneighbor+1);
609 if (dc[1] < fMinDCForStars)
610 continue;
611
612 dcsum = dc[0] + dc[1];
613
614 if(dcsum > maxDC*2)
615 {
616 if(dc[0]>=dc[1])
617 {
618 maxPixIdx[0] = pix;
619 maxPixIdx[1] = swneighbor;
620 maxDC = dc[0];
621 }
622 else
623 {
624 maxPixIdx[1] = pix;
625 maxPixIdx[0] = swneighbor;
626 maxDC = dc[1];
627 }
628 }
629 }
630 }
631 }
632 }
633
634 if (maxDC == 0)
635 {
636 *fLog << warn << " No found pixels with maximum dc" << endl;
637 return kFALSE;
638 }
639
640 maxPix = (*fGeomCam)[maxPixIdx[0]];
641
642 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;
643
644 return kTRUE;
645}
646
647Bool_t MFindStars::FindStar(MStarLocalPos* star)
648{
649
650 UInt_t numPixels = fGeomCam->GetNumPixels();
651 Float_t innerped = fStars->GetInnerPedestalDC();
652 Float_t outerped = fStars->GetOuterPedestalDC();
653
654 TArrayC origPixelsUsed;
655 origPixelsUsed.Set(numPixels);
656
657 for (UInt_t pix=1; pix<numPixels; pix++)
658 {
659 if (fDisplay.IsUsed(pix))
660 origPixelsUsed[pix]=(Char_t)kTRUE;
661 else
662 origPixelsUsed[pix]=(Char_t)kFALSE;
663 }
664
665 Float_t expX = star->GetXExp();
666 Float_t expY = star->GetYExp();
667
668 Float_t max=0;
669 UInt_t pixmax=0;
670 Float_t meanX=0;
671 Float_t meanY=0;
672 Float_t meanSqX=0;
673 Float_t meanSqY=0;
674 Float_t sumCharge=0;
675 UInt_t usedInnerPx=0;
676 UInt_t usedOuterPx=0;
677
678 const Float_t meanPresition = 3.; //[mm]
679 const UInt_t maxNumIterations = 10;
680 UInt_t numIterations = 0;
681
682 do
683 {
684 // First define a area of interest around the expected position of the star
685 for (UInt_t pix=1; pix<numPixels; pix++)
686 {
687
688 Float_t pixXpos=(*fGeomCam)[pix].GetX();
689 Float_t pixYpos=(*fGeomCam)[pix].GetY();
690 Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
691 (pixYpos-expY)*(pixYpos-expY));
692
693 if ((dist < fRingInterest) && fDisplay.IsUsed(pix))
694 fPixelsUsed[pix]=(Char_t)kTRUE;
695 else
696 fPixelsUsed[pix]=(Char_t)kFALSE;
697 }
698
699 fDisplay.SetUsed(fPixelsUsed);
700
701// determine mean x and mean y
702 usedInnerPx=0;
703 usedOuterPx=0;
704 for(UInt_t pix=0; pix<numPixels; pix++)
705 {
706 if(fDisplay.IsUsed(pix))
707 {
708 pix>lastInnerPixel?usedOuterPx++:usedInnerPx++;
709
710 Float_t charge = fDisplay.GetBinContent(pix+1);
711 Float_t pixXpos = (*fGeomCam)[pix].GetX();
712 Float_t pixYpos = (*fGeomCam)[pix].GetY();
713
714 if (charge>max)
715 {
716 max=charge;
717 pixmax=pix;
718 }
719
720 meanX += charge*pixXpos;
721 meanY += charge*pixYpos;
722 meanSqX += charge*pixXpos*pixXpos;
723 meanSqY += charge*pixYpos*pixYpos;
724 sumCharge += charge;
725 }
726 }
727
728 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " usedInnerPx " << usedInnerPx << " usedOuterPx " << usedOuterPx << endl;
729
730 meanX /= sumCharge;
731 meanY /= sumCharge;
732 meanSqX /= sumCharge;
733 meanSqY /= sumCharge;
734
735 expX = meanX;
736 expY = meanY;
737
738 if (++numIterations > maxNumIterations)
739 {
740 *fLog << warn << GetName() << "Mean calculation not converge after " << maxNumIterations << " iterations" << endl;
741 break;
742 }
743
744 }while(TMath::Abs(meanX-expX) > meanPresition || TMath::Abs(meanY-expY) > meanPresition);
745
746 Float_t rmsX = (meanSqX - meanX*meanX) - fRingInterest*fRingInterest/12;
747 Float_t rmsY = (meanSqY - meanY*meanY) - fRingInterest*fRingInterest/12;
748
749 if ( rmsX > 0)
750 rmsX = TMath::Sqrt(rmsX);
751 else
752 {
753 *fLog << warn << " MFindStars::FindStar negative rmsX² " << rmsX << endl;
754 *fLog << warn << " meanSqX " << meanSqX << " meanX " << meanX << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge << endl;
755 rmsX = 0.;
756 }
757
758 if ( rmsY > 0)
759 rmsY = TMath::Sqrt(rmsY);
760 else
761 {
762 *fLog << warn << " MFindStars::FindStar negative rmsY² " << rmsY << endl;
763 *fLog << warn << " meanSqY " << meanSqY << " meanY " << meanY << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge<< endl;
764 rmsY = 0.;
765 }
766
767 // Substrack pedestal DC
768 sumCharge-= (usedInnerPx*innerped+usedOuterPx*outerped)/(usedInnerPx+usedOuterPx);
769 max-=pixmax>lastInnerPixel?outerped:innerped;
770
771
772 star->SetCalcValues(sumCharge,max,meanX,meanY,rmsX,rmsY);
773
774 if (rmsX <= 0. || rmsY <= 0.)
775 return kFALSE;
776
777
778// fit the star spot using TMinuit
779
780
781 for (UInt_t pix=1; pix<numPixels; pix++)
782 if (fDisplay.IsUsed(pix))
783 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl;
784
785 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars " << fMinDCForStars << " pixmax>lastInnerPixel?outerped:innerped " << (pixmax>lastInnerPixel?outerped:innerped) << " fMaxNumIntegratedEvents " << fMaxNumIntegratedEvents << endl;
786
787 //Initialate variables for fit
788 fVinit[0] = 0;
789 fVinit[1] = max;
790 fLimlo[1] = fMinDCForStars-(pixmax>lastInnerPixel?outerped:innerped);
791 fLimup[1] = 30*fMaxNumIntegratedEvents-(pixmax>lastInnerPixel?outerped:innerped);
792 fVinit[2] = meanX;
793 fVinit[3] = rmsX;
794 fVinit[4] = meanY;
795 fVinit[5] = rmsY;
796 fVinit[6] = 0;
797 //Init steps
798 for(Int_t i=0; i<fNumVar; i++)
799 {
800 if (fVinit[i] != 0)
801 fStep[i] = TMath::Abs(fVinit[i]/sqrt2);
802 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fVinit[" << i << "] " << fVinit[i];
803 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fStep[" << i << "] " << fStep[i];
804 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimlo[" << i << "] " << fLimlo[i];
805 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimup[" << i << "] " << fLimup[i] << endl;
806 }
807 //
808
809 // -------------------------------------------
810 // call MINUIT
811
812 Double_t arglist[10];
813 Int_t ierflg = 0;
814
815 for (Int_t i=0; i<fNumVar; i++)
816 gMinuit->mnparm(i, fVname[i], fVinit[i], fStep[i], fLimlo[i], fLimup[i], ierflg);
817
818 TStopwatch clock;
819 clock.Start();
820
821// Now ready for minimization step
822 arglist[0] = 500;
823 arglist[1] = 1.;
824 gMinuit->mnexcm(fMethod, arglist ,2,ierflg);
825
826 clock.Stop();
827
828 if(fMinuitPrintOutLevel>=0)
829 {
830 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "Time spent for the minimization in MINUIT : " << endl;;
831 clock.Print();
832 }
833
834 Double_t integratedCharge;
835 Double_t bkFit, bkFitError;
836 Double_t maxFit, maxFitError;
837 Double_t meanXFit, meanXFitError;
838 Double_t sigmaMinorAxis, sigmaMinorAxisError;
839 Double_t meanYFit, meanYFitError;
840 Double_t sigmaMajorAxis, sigmaMajorAxisError;
841 Double_t phiFit, phiFitError;
842 Float_t chisquare = GetChisquare();
843 Int_t dregrees = GetDegreesofFreedom()-fNumVar;
844
845 if (!ierflg)
846 {
847 gMinuit->GetParameter(0,bkFit, bkFitError);
848 gMinuit->GetParameter(1,maxFit, maxFitError);
849 gMinuit->GetParameter(2,meanXFit,meanXFitError);
850 gMinuit->GetParameter(3,sigmaMinorAxis,sigmaMinorAxisError);
851 gMinuit->GetParameter(4,meanYFit,meanYFitError);
852 gMinuit->GetParameter(5,sigmaMajorAxis,sigmaMajorAxisError);
853 gMinuit->GetParameter(6,phiFit,phiFitError);
854
855 //FIXME: Do the integral properlly
856 integratedCharge = 0.;
857
858
859 }
860 else
861 {
862 bkFit = 0.;
863 maxFit = 0.;
864 meanXFit = 0.;
865 sigmaMinorAxis = 0.;
866 meanYFit = 0.;
867 sigmaMajorAxis = 0.;
868 integratedCharge = 0.;
869 phiFit = 0.;
870
871 *fLog << err << "TMinuit::Call error " << ierflg << endl;
872 }
873
874 star->SetFitValues(integratedCharge,bkFit,maxFit,meanXFit,meanYFit,sigmaMinorAxis,sigmaMajorAxis,phiFit,chisquare,dregrees);
875
876 // reset the display to the starting values
877 fDisplay.SetUsed(origPixelsUsed);
878
879 if (ierflg)
880 return kCONTINUE;
881 return kTRUE;
882}
883
884Bool_t MFindStars::ShadowStar(MStarLocalPos* star)
885{
886 UInt_t numPixels = fGeomCam->GetNumPixels();
887
888// Define an area around the star which will be set unused.
889 UInt_t shadowPx=0;
890 for (UInt_t pix=1; pix<numPixels; pix++)
891 {
892 Float_t pixXpos = (*fGeomCam)[pix].GetX();
893 Float_t pixYpos = (*fGeomCam)[pix].GetY();
894 Float_t starXpos = star->GetMeanX();
895 Float_t starYpos = star->GetMeanY();
896
897 Float_t starSize = 3*star->GetSigmaMajorAxis();
898
899 Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
900 (pixYpos-starYpos)*(pixYpos-starYpos));
901
902 if (dist > starSize && fDisplay.IsUsed(pix))
903 fPixelsUsed[pix]=(Char_t)kTRUE;
904 else
905 {
906 fPixelsUsed[pix]=(Char_t)kFALSE;
907 shadowPx++;
908 }
909 }
910
911 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " shadowPx " << shadowPx << endl;
912
913 fDisplay.SetUsed(fPixelsUsed);
914
915 return kTRUE;
916}
917
918Bool_t MFindStars::CorrSourcePos(TVector3* srcpos)
919{
920 const Int_t numstar = fStars->GetList()->GetSize();
921 if (numstar < 2)
922 {
923 cout <<" Less than 2 stars in list. Triangulation not possible."<<endl;
924 return kFALSE;
925 }
926 Float_t starmag[numstar];
927 Float_t starposx[numstar];
928 Float_t starposy[numstar];
929 Float_t starposcatx[numstar];
930 Float_t starposcaty[numstar];
931 Int_t i = 0;
932 Int_t matches = 0;
933 Int_t maxdist = 150;
934
935 TIter Next(fStars->GetList());
936 MStarLocalPos* star;
937 while ((star=(MStarLocalPos*)Next()))
938 {
939 starmag[i] = star->GetMagCalc();
940 starposx[i] = star->GetMeanX();
941 starposy[i] = star->GetMeanY();
942 i +=1;
943 }
944 for ( i = 0; i < numstar; i++)
945 {
946 starposcatx[i] = -1000;
947 starposcaty[i] = -1000;
948 for( Int_t dist = 10; dist < maxdist; dist += 10)
949 {
950 TIter Next1(fAstroCamera->GetCatList());
951 TVector3 *starcat;
952 while ((starcat=(TVector3*)Next1()) && dist < maxdist)
953 {
954 if ((starcat->X()-starposx[i])*(starcat->X()-starposx[i])+(starcat->Y()-starposy[i])*(starcat->Y()-starposy[i]) < dist*dist)
955 {
956 starposcatx[i] = starcat->X();
957 starposcaty[i] = starcat->Y();
958 dist = maxdist;
959 matches ++;
960 }
961 }
962 }
963 }
964 // take the 2 'best' determined stars: Assumpion that the bright stars are better.
965 Float_t maxmag1 = 0, maxmag2 = 0;
966 Float_t minchi1 = 0, minchi2 = 0;
967 Int_t magi1 = -1, magi2 = -1;
968 Int_t chii1 = -1, chii2 = -1;
969 for ( i = 0; i < numstar; i++)
970 {
971 if (starposcatx[i] != -1000 && starposcaty[i] != -1000)
972 {
973 if (starmag[i] > maxmag2)
974 {
975 maxmag2 = starmag[i];
976 magi2 = i;
977 }
978 if (maxmag2 > maxmag1)
979 {
980 maxmag2 = maxmag1;
981 maxmag1 = starmag[i];
982 magi2 = magi1;
983 magi1 = i;
984 }
985 }
986 }
987
988 if (matches < 2)
989 {
990 cout << " Could not find 2 stars in the camera matching with the stars in the catalog. Triangulation not possible."<<endl;
991 return kFALSE;
992 }
993
994 // Triangulation to get the correct source position:
995 TVector2 *source = new TVector2;
996 TVector2 *star1 = new TVector2;
997 TVector2 *star2 = new TVector2;
998 TVector3 *dist = new TVector3;
999 source->Set(0.,0.); // Assume correct pointing.
1000 star1->Set(starposcatx[magi1],starposcaty[magi1]);
1001 star2->Set(starposcatx[magi2],starposcaty[magi2]);
1002
1003 DistBetweenStars(star1,star2,source,dist);
1004
1005 Float_t dx = starposx[magi1] - starposx[magi2];
1006 Float_t delta = acos(dx/dist->Y());
1007 Float_t alpha = acos((dist->Y()*dist->Y()+dist->Z()*dist->Z()-dist->X()*dist->X())/2/dist->Y()/dist->Y());
1008 Float_t sigma = alpha - delta;
1009
1010 if (source->X()-starposx[magi1] > 0)
1011 srcpos->SetX(starposx[magi1] + dist->Z()*abs(cos(sigma)));
1012 else
1013 srcpos->SetX(starposx[magi1] - dist->Z()*abs(cos(sigma)));
1014 if (source->Y()-starposy[magi1] > 0)
1015 srcpos->SetY(starposy[magi1] + dist->Z()*abs(sin(sigma)));
1016 else
1017 srcpos->SetY(starposy[magi1] - dist->Z()*abs(sin(sigma)));
1018
1019 return kTRUE;
1020}
1021
1022
1023void 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
1029}
Note: See TracBrowser for help on using the repository browser.