source: trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MFindStars.cc@ 5138

Last change on this file since 5138 was 5120, checked in by stamerra, 20 years ago
*** empty log message ***
File size: 37.1 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! Javier López , 4/2004 <mailto:jlopez@ifae.es>
19! Wolfgang Wittek, 8/2004 <mailto:wittek@mppmu.mpg.de>
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#include <TEllipse.h>
43
44
45#include "MObservatory.h"
46#include "MAstroCamera.h"
47#include "MMcConfigRunHeader.h"
48
49#include "MLog.h"
50#include "MLogManip.h"
51
52#include "MHCamera.h"
53
54#include "MGeomCam.h"
55#include "MGeomPix.h"
56#include "MCameraDC.h"
57#include "MTime.h"
58#include "MReportDrive.h"
59#include "MStarCam.h"
60#include "MStarPos.h"
61
62#include "MParList.h"
63#include "MTaskList.h"
64
65ClassImp(MFindStars);
66using namespace std;
67
68const Float_t sqrt2 = sqrt(2.);
69const Float_t sqrt3 = sqrt(3.);
70const UInt_t lastInnerPixel = 396;
71
72
73//______________________________________________________________________________
74//
75// The 2D uncorrelated gaussian function used to fit the spot of the star
76//
77static Double_t func(float x,float y,Double_t *par)
78{
79 Double_t value=par[0]*exp( -(x-par[1])*(x-par[1])/(2*par[2]*par[2]) )
80 *exp( -(y-par[3])*(y-par[3])/(2*par[4]*par[4]) );
81 return value;
82}
83
84//______________________________________________________________________________
85//
86// The 2D correlated gaussian function used to fit the spot of the star
87//
88static Double_t funcCG(float x,float y,Double_t *par)
89{
90 Double_t temp = 1.0-par[5]*par[5];
91 // Double_t value = par[0] / (2.0*TMath::Pi()*par[2]*par[4]*sqrt(temp))
92 Double_t value = par[0]
93 * exp( -0.5/temp * ( (x-par[1])*(x-par[1])/(par[2]*par[2])
94 -2.0*par[5] * (x-par[1])*(y-par[3])/(par[2]*par[4])
95 + (y-par[3])*(y-par[3])/(par[4]*par[4]) ) );
96 return value;
97}
98
99//______________________________________________________________________________
100//
101// Function used by Minuit to do the fit
102//
103static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
104{
105
106 MParList* plist = (MParList*)gMinuit->GetObjectFit();
107 MTaskList* tlist = (MTaskList*)plist->FindObject("MTaskList");
108 MFindStars* find = (MFindStars*)tlist->FindObject("MFindStars");
109 MStarCam* stars =
110 (MStarCam*)plist->FindObject("MStarCam","MStarCam");
111 MGeomCam* geom = (MGeomCam*)plist->FindObject("MGeomCam");
112
113 MHCamera& display = (MHCamera&)find->GetDisplay();
114
115 Float_t innerped = stars->GetInnerPedestalDC();
116 Float_t innerrms = stars->GetInnerPedestalRMSDC();
117 Float_t outerped = stars->GetOuterPedestalDC();
118 Float_t outerrms = stars->GetOuterPedestalRMSDC();
119
120 UInt_t numPixels = geom->GetNumPixels();
121
122
123//calculate chisquare
124 Double_t chisq = 0;
125 Double_t delta;
126 Double_t x,y,z;
127 Double_t errorz=0;
128
129 UInt_t usedPx=0;
130 for (UInt_t pixid=1; pixid<numPixels; pixid++)
131 {
132
133
134 if (display.IsUsed(pixid))
135 {
136 x = (*geom)[pixid].GetX();
137 y = (*geom)[pixid].GetY();
138
139 z = display.GetBinContent(pixid+1)
140 - (pixid>lastInnerPixel?outerped:innerped);
141 errorz=(pixid>lastInnerPixel?outerrms:innerrms);
142
143
144 if (errorz > 0.0)
145 {
146 usedPx++;
147
148 Double_t fu;
149 if (find->GetUseCorrelatedGauss())
150 fu = funcCG(x,y,par);
151 else
152 fu = func(x,y,par);
153
154 delta = (z-fu) / errorz;
155 chisq += delta*delta;
156
157 if (iflag == 3)
158 {
159 gLog << "fcn : usedPx, pixid, content, pedestal,z, fu, errorz, delta = "
160 << usedPx << ", " << pixid << ", "
161 << display.GetBinContent(pixid+1) << ", "
162 << (pixid>lastInnerPixel?outerped:innerped)
163 << ", " << z << ", " << fu << ", " << errorz
164 << ", " << delta << endl;
165 }
166 }
167 else
168 cerr << " TMinuit::fcn errorz[" << pixid << "] " << errorz << endl;
169 }
170 }
171 f = chisq;
172
173 find->SetChisquare(chisq);
174 find->SetDegreesofFreedom(usedPx);
175
176 //gLog << "fcn : chisq, usedPx = " << chisq << ", " << usedPx << endl;
177}
178
179//-------------------------------------------------------------------------
180//
181// Constructor
182//
183
184MFindStars::MFindStars(const char *name, const char *title):
185 fGeomCam(NULL), fCurr(NULL), fTimeCurr(NULL), fDrive(NULL), fStars(NULL),
186 fNumVar(6)
187{
188 fName = name ? name : "MFindStars";
189 fTitle = title ? title : "Tool to find stars from DC Currents";
190
191 // the correlated Gauss function
192 // is fitted by default
193 fNumVar = 6;
194 fUseCorrelatedGauss = kTRUE;
195
196 fNumIntegratedEvents=0;
197 fMaxNumIntegratedEvents = 10;
198 fRingInterest = 125.; //[mm] ~ 0.4 deg
199 fDCTailCut = 4;
200
201 fPixelsUsed.Set(577);
202 fPixelsUsed.Reset((Char_t)kTRUE);
203
204 //Fitting(Minuit) initialitation
205 const Float_t pixelSize = 31.5; //[mm]
206 fMinuitPrintOutLevel = -1;
207 //fMinuitPrintOutLevel = 3;
208
209 fVname = new TString[fNumVar];
210 fVinit.Set(fNumVar);
211 fStep.Set(fNumVar);
212 fLimlo.Set(fNumVar);
213 fLimup.Set(fNumVar);
214 fFix.Set(fNumVar);
215
216 fVname[0] = "max";
217 fVinit[0] = 10.*fMaxNumIntegratedEvents;
218 fStep[0] = fVinit[0]/sqrt2;
219 fLimlo[0] = fMinDCForStars;
220 fLimup[0] = 30.*fMaxNumIntegratedEvents;
221 fFix[0] = 0;
222
223 fVname[1] = "meanx";
224 fVinit[1] = 0.;
225 fStep[1] = fVinit[1]/sqrt2;
226 fLimlo[1] = -600.;
227 fLimup[1] = 600.;
228 fFix[1] = 0;
229
230 fVname[2] = "sigmax";
231 fVinit[2] = pixelSize;
232 fStep[2] = fVinit[2]/sqrt2;
233 fLimlo[2] = pixelSize/(2*sqrt3);
234 fLimup[2] = 500.;
235 fFix[2] = 0;
236
237 fVname[3] = "meany";
238 fVinit[3] = 0.;
239 fStep[3] = fVinit[3]/sqrt2;
240 fLimlo[3] = -600.;
241 fLimup[3] = 600.;
242 fFix[3] = 0;
243
244 fVname[4] = "sigmay";
245 fVinit[4] = pixelSize;
246 fStep[4] = fVinit[4]/sqrt2;
247 fLimlo[4] = pixelSize/(2*sqrt3);
248 fLimup[4] = 500.;
249 fFix[4] = 0;
250
251 if (fUseCorrelatedGauss)
252 {
253 fVname[5] = "xycorr";
254 fVinit[5] = 0.0;
255 fStep[5] = 0.1;
256 fLimlo[5] = -1.0;
257 fLimup[5] = 1.0;
258 fFix[5] = 0;
259 }
260
261 fObjectFit = NULL;
262 // fMethod = "SIMPLEX";
263 fMethod = "MIGRAD";
264 // fMethod = "MINIMIZE";
265 fNulloutput = kFALSE;
266
267 // Set output level
268 // fLog->SetOutputLevel(3); // No dbg messages
269
270 fGeometryFile="";
271 fBSCFile="";
272}
273
274//-------------------------------------------------------------------------
275//
276// Set 'fUseCorrelatedGauss' flag
277//
278// if 'fUseCorrelatedGauss' is TRUE a 2dim Gaussian with correlation
279// will be fitted
280//
281
282void MFindStars::SetUseCorrelatedGauss(Bool_t usecorrgauss)
283{
284 fUseCorrelatedGauss = usecorrgauss;
285
286 if (usecorrgauss)
287 fNumVar = 6;
288 else
289 fNumVar = 5;
290}
291
292//-------------------------------------------------------------------------
293//
294// PreProcess
295//
296
297Int_t MFindStars::PreProcess(MParList *pList)
298{
299 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
300
301 if (!fGeomCam)
302 {
303 *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl;
304 return kFALSE;
305 }
306
307 // Initialize camera display with the MGeomCam information
308 fDisplay.SetGeometry(*fGeomCam,"FindStarsDisplay","FindStarsDisplay");
309 fDisplay.SetUsed(fPixelsUsed);
310
311 fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));
312
313 if (!fCurr)
314 {
315 *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl;
316 return kFALSE;
317 }
318
319 fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents"));
320
321 if (!fTimeCurr)
322 {
323 *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl;
324 return kFALSE;
325 }
326
327 fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
328
329 if (!fDrive)
330 {
331
332 *fLog << warn << AddSerialNumber("MReportDrive") << " not found ... ignored." << endl;
333
334 }
335 else
336 {
337 MObservatory magic1;
338
339 MMcConfigRunHeader *config=0;
340 MGeomCam *geom=0;
341
342 TFile file(fGeometryFile);
343 TTree *tree = dynamic_cast<TTree*>(file.Get("RunHeaders"));
344 tree->SetBranchAddress("MMcConfigRunHeader.", &config);
345 if (tree->GetBranch("MGeomCam.")) tree->SetBranchAddress("MGeomCam.", &geom);
346
347 TClonesArray mirrorArray = *config->GetMirrors();
348 fAstro.SetMirrors(*config->GetMirrors());
349 fAstro.SetGeom(*geom);
350 fAstro.ReadBSC(fBSCFile);
351 fAstro.SetObservatory(magic1);
352
353 }
354
355
356 fStars = (MStarCam*)pList->FindCreateObj(AddSerialNumber("MStarCam"),"MStarCam");
357 if (!fStars)
358 {
359 *fLog << err << AddSerialNumber("MStarCam") << " cannot be created ... aborting" << endl;
360 return kFALSE;
361 }
362
363 fMinDCForStars = 1.*fMaxNumIntegratedEvents; //[uA]
364
365 // Initialize the TMinuit object
366
367 TMinuit *gMinuit = new TMinuit(fNumVar); //initialize TMinuit with a maximum of params
368 gMinuit->SetFCN(fcn);
369
370 Double_t arglist[10];
371 Int_t ierflg = 0;
372
373 arglist[0] = 1;
374 gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
375 arglist[0] = fMinuitPrintOutLevel;
376 gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
377
378 // suppress warnings
379 Int_t errWarn;
380 Double_t tmpwar = 0;
381 gMinuit->mnexcm("SET NOW", &tmpwar, 0, errWarn);
382
383
384 // Set MParList object pointer to allow mimuit to access internally
385 gMinuit->SetObjectFit(pList);
386 return kTRUE;
387}
388
389//------------------------------------------------------------------------
390//
391// Process :
392//
393//
394//
395//
396
397Int_t MFindStars::Process()
398{
399 UInt_t numPixels = fGeomCam->GetNumPixels();
400 TArrayC origPixelsUsed;
401 origPixelsUsed.Set(numPixels);
402
403 if (fNumIntegratedEvents >= fMaxNumIntegratedEvents) {
404
405 //Fist delete the previous stars in the list
406 fStars->GetList()->Delete();
407
408 if (fDrive) {
409
410 //If drive information is provided we take RaDec info
411 //from the drive and let the star list fill by the astrocamera.
412
413 //FIXME: rwagner: Doesn't work as expected
414 //FIXME: rwagner: For the time being set manually.
415 //fAstro.SetRaDec(fDrive->GetRa(), fDrive->GetDec());
416
417 fAstro.SetTime(*fTimeCurr);
418 fAstro.SetGuiActive();
419 fAstro.FillStarList(fStars->GetList());
420
421 cout << "Number of Stars added by astrocamera is "
422 <<fStars->GetList()->GetSize() << endl;
423
424 MStarPos* starpos;
425 TIter Next(fStars->GetList());
426 while ((starpos=(MStarPos*)Next())) {
427 //starpos->SetCalcValues(40,40,starpos->GetXExp(), starpos->GetYExp(),
428 // fRingInterest/2, fRingInterest/2,
429 // 0., 0., 0., 0.);
430 //starpos->SetFitValues (40,40,starpos->GetXExp(), starpos->GetYExp(),
431 // fRingInterest/2, fRingInterest/2, 0.,
432 // 0., 0., 0., 0., -1);
433 }
434
435 for (UInt_t pix=1; pix<numPixels; pix++) {
436 if (fDisplay.IsUsed(pix))
437 origPixelsUsed[pix]=(Char_t)kTRUE;
438 else
439 origPixelsUsed[pix]=(Char_t)kFALSE;
440 }
441
442 DCPedestalCalc();
443
444 }
445 else
446 {
447
448 cout << "No drive information available: Will not use a star catalog to identify stars."<< endl;
449
450 for (UInt_t pix=1; pix<numPixels; pix++) {
451 if (fDisplay.IsUsed(pix))
452 origPixelsUsed[pix]=(Char_t)kTRUE;
453 else
454 origPixelsUsed[pix]=(Char_t)kFALSE;
455 }
456
457 if (DCPedestalCalc()) {
458
459 Float_t innermin = fStars->GetInnerPedestalDC()+fDCTailCut*fStars->GetInnerPedestalRMSDC();
460 Float_t outermin = fStars->GetOuterPedestalDC()+fDCTailCut*fStars->GetOuterPedestalRMSDC();
461 fMinDCForStars = innermin>outermin?innermin:outermin;
462 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars = " << fMinDCForStars << endl;
463
464 // Find the star candidates searching the most brights pairs of pixels
465 Float_t maxPixelDC;
466 MGeomPix maxPixel;
467
468 while(FindPixelWithMaxDC(maxPixelDC, maxPixel)) {
469 MStarPos *starpos = new MStarPos;
470
471 starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
472
473 //starpos->SetCalcValues(maxPixelDC, maxPixelDC,
474 // maxPixel.GetX(), maxPixel.GetY(),
475 // fRingInterest/2, fRingInterest/2,
476 // 0., 0., 0., 0.);
477
478 //starpos->SetFitValues(maxPixelDC, maxPixelDC,
479 // maxPixel.GetX(), maxPixel.GetY(),
480 // fRingInterest/2, fRingInterest/2, 0.,
481 // 0., 0., 0., 0., -1);
482
483 fStars->GetList()->Add(starpos);
484 ShadowStar(starpos);
485 }
486 fDisplay.SetUsed(origPixelsUsed);
487 }
488 }
489
490 //Show the stars found
491 fStars->Print("namepossizechierr");
492
493 //loop to extract position of stars on the camera
494 if (fStars->GetList()->GetSize() == 0) {
495 *fLog << err << GetName() << "No stars candidates in the camera." << endl;
496 return kCONTINUE;
497 }
498 else
499 {
500 *fLog << inf << GetName() << " found " << fStars->GetList()->GetSize()
501 << " star candidates in the camera." << endl;
502 }
503
504 for (UInt_t pix=1; pix<numPixels; pix++) {
505 if (fDisplay.IsUsed(pix))
506 origPixelsUsed[pix]=(Char_t)kTRUE;
507 else
508 origPixelsUsed[pix]=(Char_t)kFALSE;
509 }
510
511 TIter Next(fStars->GetList());
512 MStarPos* star;
513 while ((star=(MStarPos*)Next())) {
514 FindStar(star);
515 ShadowStar(star);
516 }
517
518 //Show the stars found: Here it is interesting to see what FindStars
519 //made out of the positions we gave to it.
520 if (fMinuitPrintOutLevel>=0)
521 fStars->Print("namepossizchierr");
522
523 //After finding stars reset all variables
524 fDisplay.Reset();
525 fStars->GetDisplay().Reset(); //FIXME: Put this display just in the container
526 fDisplay.SetUsed(origPixelsUsed);
527 fNumIntegratedEvents=0;
528 }
529
530 for (UInt_t pix=1; pix<numPixels; pix++) {
531 if (fDisplay.IsUsed(pix))
532 origPixelsUsed[pix]=(Char_t)kTRUE;
533 else
534 origPixelsUsed[pix]=(Char_t)kFALSE;
535
536 }
537
538 fDisplay.AddCamContent(*fCurr);
539 fNumIntegratedEvents++;
540 fDisplay.SetUsed(origPixelsUsed);
541
542 return kTRUE;
543}
544
545Int_t MFindStars::PostProcess()
546{
547 return kTRUE;
548}
549
550void MFindStars::SetBlindPixels(TArrayS blindpixels)
551{
552 Int_t npix = blindpixels.GetSize();
553
554 for (Int_t idx=0; idx<npix; idx++)
555 {
556 fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;
557 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "MFindStars::SetBlindPixels fDisplay.IsUsed(" <<blindpixels[idx] << ") kFALSE" << endl;
558 }
559
560 fDisplay.SetUsed(fPixelsUsed);
561}
562
563//-------------------------------------------------------------------------
564//
565// DCPedestalCalc :
566//
567//
568//
569//
570
571Bool_t MFindStars::DCPedestalCalc()
572{
573 //-------------------------------------------------------------
574 // save pointer to the MINUIT object for optimizing the supercuts
575 // because it will be overwritten
576 // when fitting the alpha distribution in MHFindSignificance
577 TMinuit *savePointer = gMinuit;
578 //-------------------------------------------------------------
579
580 UInt_t numPixels = fGeomCam->GetNumPixels();
581 Float_t ped;
582 Float_t rms;
583
584 //TH1F **dchist = new TH1F*[2];
585 TH1F *dchist[2];
586 TString tit;
587 TString nam;
588 for (UInt_t i=0; i<2; i++)
589 {
590 nam = i;
591 tit = "dchist[";
592 tit += i;
593 tit += "]";
594 dchist[i] = new TH1F(nam, tit,
595 26,0.4*fMaxNumIntegratedEvents,3.*fMaxNumIntegratedEvents);
596 }
597
598 for (UInt_t pix=1; pix<=lastInnerPixel; pix++)
599 dchist[0]->Fill(fDisplay.GetBinContent(pix+1));
600 for (UInt_t pix=lastInnerPixel+1; pix<numPixels; pix++)
601 dchist[1]->Fill(fDisplay.GetBinContent(pix+1));
602
603 // inner/outer pixels
604 for (UInt_t i=0; i<2; i++)
605 {
606 Float_t nummaxprobdc = dchist[i]->GetBinContent(dchist[i]->GetMaximumBin());
607 Float_t maxprobdc = dchist[i]->GetBinCenter(dchist[i]->GetMaximumBin());
608 UInt_t bin = dchist[i]->GetMaximumBin();
609 do
610 {
611 bin++;
612 }
613 while(dchist[i]->GetBinContent(bin)/nummaxprobdc > 0.5);
614 Float_t halfmaxprobdc = dchist[i]->GetBinCenter(bin);
615
616 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] ";
617 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << dchist[i]->GetBinContent(bin) << "]" << endl;
618
619 Float_t rmsguess = TMath::Abs(maxprobdc-halfmaxprobdc);
620 Float_t min = maxprobdc-3*rmsguess;
621 min = (min<0.?0.:min);
622 Float_t max = maxprobdc+3*rmsguess;
623
624 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
625
626 TF1 func("func","gaus",min,max);
627 func.SetParameters(nummaxprobdc, maxprobdc, rmsguess);
628
629 dchist[i]->Fit("func","QR0");
630
631 UInt_t aproxnumdegrees = 6*(bin-dchist[i]->GetMaximumBin());
632 Float_t chiq = func.GetChisquare();
633 ped = func.GetParameter(1);
634 rms = func.GetParameter(2);
635
636 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
637
638 Int_t numsigmas = 5;
639 Axis_t minbin = ped-numsigmas*rms/dchist[i]->GetBinWidth(1);
640 minbin=minbin<1?1:minbin;
641 Axis_t maxbin = ped+numsigmas*rms/dchist[i]->GetBinWidth(1);
642 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << dchist[i]->Integral((int)minbin,(int)maxbin) << endl;
643
644 //Check results from the fit are consistent
645 if (ped < 0. || rms < 0.)
646 {
647 *fLog << dbg << "dchist[i]->GetEntries()" << dchist[i]->GetEntries();
648// TCanvas *c1 = new TCanvas("c2","c2",500,800);
649// dchist[i]->Draw();
650// c1->Print("dchist.ps");
651// delete c1;
652// exit(1);
653 }
654 else if (TMath::Abs(ped-maxprobdc) > rmsguess || rms > rmsguess)
655 {
656 *fLog << warn << GetName() << " Pedestal DC fit give non consistent results for " << (i==0?"Inner":"Outer") << "pixels." << endl;
657 *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
658 *fLog << warn << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
659 ped = maxprobdc;
660 rms = rmsguess/1.175; // FWHM=2.35*rms
661 }
662
663 if (i == 0)
664 {
665 fStars->SetInnerPedestalDC(ped);
666 fStars->SetInnerPedestalRMSDC(rms);
667 }
668 else
669 {
670 fStars->SetOuterPedestalDC(ped);
671 fStars->SetOuterPedestalRMSDC(rms);
672 }
673 }
674
675
676
677 for (UInt_t i=0; i<2; i++)
678 {
679 delete dchist[i];
680 }
681 //delete [] dchist;
682
683 //=================================================================
684
685 // reset gMinuit to the MINUIT object for optimizing the supercuts
686 gMinuit = savePointer;
687 //-------------------------------------------
688
689 if (fStars->GetInnerPedestalDC() < 0. || fStars->GetInnerPedestalRMSDC() < 0. || fStars->GetOuterPedestalDC() < 0. || fStars->GetOuterPedestalRMSDC() < 0.)
690 return kFALSE;
691
692 return kTRUE;
693}
694
695Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
696{
697 UInt_t numPixels = fGeomCam->GetNumPixels();
698
699// Find the two close pixels with the maximun dc
700 UInt_t maxPixIdx[2];
701
702 maxDC = 0;
703
704 for (UInt_t pix=1; pix<numPixels; pix++)
705 {
706 if(fDisplay.IsUsed(pix))
707 {
708 Float_t dc[2];
709 dc[0] = fDisplay.GetBinContent(pix+1);
710 if (dc[0] < fMinDCForStars)
711 continue;
712
713 MGeomPix &g = (*fGeomCam)[pix];
714 Int_t numNextNeighbors = g.GetNumNeighbors();
715
716 Float_t dcsum;
717 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
718 {
719 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
720 if(fDisplay.IsUsed(swneighbor))
721 {
722 dc[1] = fDisplay.GetBinContent(swneighbor+1);
723 if (dc[1] < fMinDCForStars)
724 continue;
725
726 dcsum = dc[0] + dc[1];
727
728 if(dcsum > maxDC*2)
729 {
730 if(dc[0]>=dc[1])
731 {
732 maxPixIdx[0] = pix;
733 maxPixIdx[1] = swneighbor;
734 maxDC = dc[0];
735 }
736 else
737 {
738 maxPixIdx[1] = pix;
739 maxPixIdx[0] = swneighbor;
740 maxDC = dc[1];
741 }
742 }
743 }
744 }
745 }
746 }
747
748 if (maxDC == 0)
749 {
750 *fLog << warn << " No found pixels with maximum dc" << endl;
751 return kFALSE;
752 }
753
754 maxPix = (*fGeomCam)[maxPixIdx[0]];
755
756 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;
757
758 return kTRUE;
759}
760
761
762//----------------------------------------------------------------------------
763//
764// FindStar
765//
766//
767//
768
769Bool_t MFindStars::FindStar(MStarPos* star)
770{
771 UInt_t numPixels = fGeomCam->GetNumPixels();
772 Float_t innerped = fStars->GetInnerPedestalDC();
773 Float_t outerped = fStars->GetOuterPedestalDC();
774
775 TArrayC origPixelsUsed;
776 origPixelsUsed.Set(numPixels);
777
778 for (UInt_t pix=1; pix<numPixels; pix++)
779 {
780 if (fDisplay.IsUsed(pix))
781 origPixelsUsed[pix]=(Char_t)kTRUE;
782 else
783 origPixelsUsed[pix]=(Char_t)kFALSE;
784 }
785
786 Float_t expX = star->GetXExp();
787 Float_t expY = star->GetYExp();
788
789 Float_t max =0;
790 UInt_t pixmax =0;
791 Float_t meanX =0;
792 Float_t meanY =0;
793 Float_t meanSqX =0;
794 Float_t meanSqY =0;
795 Float_t sumCharge =0;
796 Float_t sumSqCharge =0;
797 UInt_t usedInnerPx =0;
798 UInt_t usedOuterPx =0;
799 Float_t factor = 1.0;
800
801 Float_t meanXold = 1.e10;
802 Float_t meanYold = 1.e10;
803
804 const Float_t meanPresition = 3.; //[mm]
805 const UInt_t maxNumIterations = 10;
806 UInt_t numIterations = 0;
807
808 //-------------------- start of iteration loop -----------------------
809 for (UInt_t it=0; it<maxNumIterations; it++)
810 {
811 //rwagner: Need to find px with highest dc and need to assume it is
812 // somewhere near the center of the star
813 //after that we need to REDEFINE our roi! because expected pos could be
814 // quite off that px!
815
816 // 1.) Initial roi around expected position
817
818 for (UInt_t pix=1; pix<numPixels; pix++)
819 {
820
821 Float_t pixXpos=(*fGeomCam)[pix].GetX();
822 Float_t pixYpos=(*fGeomCam)[pix].GetY();
823 Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
824 (pixYpos-expY)*(pixYpos-expY));
825
826 if ((dist < fRingInterest) && fDisplay.IsUsed(pix))
827 fPixelsUsed[pix]=(Char_t)kTRUE;
828 else
829 fPixelsUsed[pix]=(Char_t)kFALSE;
830 }
831 fDisplay.SetUsed(fPixelsUsed);
832
833 // 2.) Find px with highest dc in that region
834
835 max = 0.0;
836 pixmax = 0;
837 for(UInt_t pix=0; pix<numPixels; pix++)
838 if(fDisplay.IsUsed(pix))
839 {
840 Float_t charge = fDisplay.GetBinContent(pix+1);
841 if (charge>max)
842 {
843 max=charge;
844 pixmax=pix;
845 }
846 }
847
848 // 3.) make it new center
849
850 expX = (*fGeomCam)[pixmax].GetX();
851 expY = (*fGeomCam)[pixmax].GetY();
852
853 for (UInt_t pix=1; pix<numPixels; pix++)
854 fPixelsUsed[pix]=(Char_t)kTRUE;
855 fDisplay.SetUsed(fPixelsUsed);
856
857 // First define a area of interest around the expected position of the star
858 for (UInt_t pix=1; pix<numPixels; pix++)
859 {
860 Float_t pixXpos=(*fGeomCam)[pix].GetX();
861 Float_t pixYpos=(*fGeomCam)[pix].GetY();
862 Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
863 (pixYpos-expY)*(pixYpos-expY));
864
865 if ((dist < fRingInterest) && fDisplay.IsUsed(pix))
866 fPixelsUsed[pix]=(Char_t)kTRUE;
867 else
868 fPixelsUsed[pix]=(Char_t)kFALSE;
869 }
870
871 fDisplay.SetUsed(fPixelsUsed);
872
873 // determine mean x and mean y
874 max = 0.0;
875 pixmax = 0;
876 meanX = 0.0;
877 meanY = 0.0;
878 meanSqX = 0.0;
879 meanSqY = 0.0;
880 sumCharge = 0.0;
881 sumSqCharge = 0.0;
882 usedInnerPx = 0;
883 usedOuterPx = 0;
884
885 for(UInt_t pix=0; pix<numPixels; pix++)
886 {
887 if(fDisplay.IsUsed(pix))
888 {
889 pix>lastInnerPixel?usedOuterPx++:usedInnerPx++;
890
891 Float_t charge = fDisplay.GetBinContent(pix+1);
892 Float_t pixXpos = (*fGeomCam)[pix].GetX();
893 Float_t pixYpos = (*fGeomCam)[pix].GetY();
894
895 if (charge>max)
896 {
897 max=charge;
898 pixmax=pix;
899 }
900
901 meanX += charge*pixXpos;
902 meanY += charge*pixYpos;
903 meanSqX += charge*pixXpos*pixXpos;
904 meanSqY += charge*pixYpos*pixYpos;
905 sumCharge += charge;
906 sumSqCharge += charge*charge;
907 }
908 }
909
910 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " usedInnerPx " << usedInnerPx << " usedOuterPx " << usedOuterPx << endl;
911
912 meanX /= sumCharge;
913 meanY /= sumCharge;
914 meanSqX /= sumCharge;
915 meanSqY /= sumCharge;
916 factor = sqrt( sumSqCharge/(sumCharge*sumCharge) );
917
918 // stop iteration when (meanX, meanY) becomes stable
919 if ( TMath::Abs(meanX-meanXold) < meanPresition &&
920 TMath::Abs(meanY-meanYold) < meanPresition )
921 break;
922
923 meanXold = meanX;
924 meanYold = meanY;
925 numIterations++;
926 }
927 // this statement was commented out because the condition
928 // was never fulfilled :
929 //while(TMath::Abs(meanX-expX) > meanPresition ||
930 // TMath::Abs(meanY-expY) > meanPresition);
931 //-------------------- end of iteration loop -----------------------
932
933 if (numIterations >= maxNumIterations)
934 {
935 *fLog << warn << GetName() << "Mean calculation not converged after "
936 << maxNumIterations << " iterations" << endl;
937 }
938
939
940 //Float_t rmsX = (meanSqX - meanX*meanX) - fRingInterest*fRingInterest/12;
941 //Float_t rmsY = (meanSqY - meanY*meanY) - fRingInterest*fRingInterest/12;
942
943 Float_t rmsX = (meanSqX - meanX*meanX);
944 Float_t rmsY = (meanSqY - meanY*meanY);
945
946 if ( rmsX > 0)
947 rmsX = TMath::Sqrt(rmsX);
948 else
949 {
950 *fLog << warn << " MFindStars::FindStar negative rmsX² " << rmsX << endl;
951 *fLog << warn << " meanSqX " << meanSqX << " meanX " << meanX << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge << endl;
952 rmsX = 0.;
953 }
954
955 if ( rmsY > 0)
956 rmsY = TMath::Sqrt(rmsY);
957 else
958 {
959 *fLog << warn << " MFindStars::FindStar negative rmsY² " << rmsY << endl;
960 *fLog << warn << " meanSqY " << meanSqY << " meanY " << meanY << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge<< endl;
961 rmsY = 0.;
962 }
963
964 Float_t deltameanX = factor * rmsX;
965 Float_t deltameanY = factor * rmsY;
966
967
968 // Substrack pedestal DC
969 sumCharge-= (usedInnerPx*innerped+usedOuterPx*outerped)/(usedInnerPx+usedOuterPx);
970 max-=pixmax>lastInnerPixel?outerped:innerped;
971
972
973 star->SetCalcValues( sumCharge, max, meanX, meanY, rmsX, rmsY,
974 0., deltameanX*deltameanX, 0., deltameanY*deltameanY);
975
976 if (rmsX <= 0. || rmsY <= 0.)
977 return kFALSE;
978
979
980// fit the star spot using TMinuit
981
982
983 for (UInt_t pix=1; pix<numPixels; pix++)
984 if (fDisplay.IsUsed(pix))
985 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl;
986
987 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars " << fMinDCForStars << " pixmax>lastInnerPixel?outerped:innerped " << (pixmax>lastInnerPixel?outerped:innerped) << " fMaxNumIntegratedEvents " << fMaxNumIntegratedEvents << endl;
988
989 //Initialize variables for fit
990 fVinit[0] = max;
991 fLimlo[0] = fMinDCForStars-(pixmax>lastInnerPixel?outerped:innerped);
992 fLimup[0] = 30*fMaxNumIntegratedEvents-(pixmax>lastInnerPixel?outerped:innerped);
993 fVinit[1] = meanX;
994 fVinit[2] = rmsX;
995 fVinit[3] = meanY;
996 fVinit[4] = rmsY;
997 //Init steps
998 for(Int_t i=0; i<fNumVar; i++)
999 {
1000 if (fVinit[i] != 0)
1001 fStep[i] = TMath::Abs(fVinit[i]/sqrt2);
1002 else
1003 fStep[i] = 0.1;
1004
1005 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fVinit[" << i << "] " << fVinit[i];
1006 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fStep[" << i << "] " << fStep[i];
1007 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimlo[" << i << "] " << fLimlo[i];
1008 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimup[" << i << "] " << fLimup[i] << endl;
1009 }
1010
1011
1012 // -------------------------------------------
1013 // call MINUIT
1014
1015 Double_t arglist[10];
1016 Int_t ierflg = 0;
1017
1018 for (Int_t i=0; i<fNumVar; i++)
1019 gMinuit->mnparm(i, fVname[i], fVinit[i], fStep[i], fLimlo[i], fLimup[i], ierflg);
1020
1021 TStopwatch clock;
1022 clock.Start();
1023
1024// Now ready for minimization step
1025 arglist[0] = 500;
1026 arglist[1] = 1.;
1027 gMinuit->mnexcm(fMethod, arglist ,2,ierflg);
1028
1029 // call fcn with iflag = 3
1030 //arglist[0] = 3;
1031 //Int_t ierflg3;
1032 //gMinuit->mnexcm("CALL", arglist, 1, ierflg3);
1033
1034 clock.Stop();
1035
1036 if(fMinuitPrintOutLevel>=0)
1037 {
1038 if (fMinuitPrintOutLevel>=0) *fLog << dbg << "Time spent for the minimization in MINUIT : " << endl;;
1039 clock.Print();
1040 }
1041
1042
1043 //---------- for the uncorrelated Gauss fit (start) ---------
1044 if (!fUseCorrelatedGauss)
1045 {
1046 Double_t integratedCharge;
1047 Double_t maxFit, maxFitError;
1048 Double_t meanXFit, meanXFitError;
1049 Double_t sigmaX, sigmaXError;
1050 Double_t meanYFit, meanYFitError;
1051 Double_t sigmaY, sigmaYError;
1052 Float_t chisquare = GetChisquare();
1053 Int_t degrees = GetDegreesofFreedom()-fNumVar;
1054
1055 if (!ierflg)
1056 {
1057 gMinuit->GetParameter(0, maxFit, maxFitError);
1058 gMinuit->GetParameter(1, meanXFit, meanXFitError);
1059 gMinuit->GetParameter(2, sigmaX, sigmaXError);
1060 gMinuit->GetParameter(3, meanYFit, meanYFitError);
1061 gMinuit->GetParameter(4, sigmaY, sigmaYError);
1062
1063 //FIXME: Do the integral properlly
1064 integratedCharge = 0.;
1065 }
1066 else
1067 {
1068 maxFit = 0.;
1069 meanXFit = 0.;
1070 sigmaX = 0.;
1071 meanYFit = 0.;
1072 sigmaY = 0.;
1073 integratedCharge = 0.;
1074
1075 *fLog << err << "TMinuit::Call error " << ierflg << endl;
1076 }
1077
1078 // get error matrix
1079 Double_t matrix[5][5];
1080 gMinuit->mnemat(&matrix[0][0],5);
1081
1082 star->SetFitValues(integratedCharge,maxFit, meanXFit, meanYFit,
1083 sigmaX, sigmaY, 0.0,
1084 matrix[1][1], matrix[1][3], matrix[3][3],
1085 chisquare, degrees);
1086 }
1087 //---------- for the uncorrelated Gauss fit (end) ---------
1088
1089 //---------- for the correlated Gauss fit (start) ---------
1090 if (fUseCorrelatedGauss)
1091 {
1092 TArrayD fValues(fNumVar);
1093 TArrayD fErrors(fNumVar);
1094
1095 const static Int_t fNdim = 6;
1096 Double_t matrix[fNdim][fNdim];
1097
1098 Double_t fmin = 0.0;
1099 Double_t fedm = 0.0;
1100 Double_t errdef = 0.0;
1101 Int_t npari = 0;
1102 Int_t nparx = 0;
1103 Int_t istat = 0;
1104 gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1105 if (fMinuitPrintOutLevel>=0)
1106 {
1107 *fLog << dbg
1108 << "Status of Correlated Gauss fit to the DC currents : " << endl;
1109 }
1110 if (fMinuitPrintOutLevel>=0)
1111 {
1112 *fLog << dbg
1113 << " fmin, fedm, errdef, npari, nparx, istat = "
1114 << fmin << ", " << fedm << ", " << errdef << ", " << npari
1115 << ", " << nparx << ", " << istat << endl;
1116 }
1117
1118 if (!ierflg)
1119 {
1120 for (Int_t k=0; k<fNumVar; k++)
1121 gMinuit->GetParameter( k, fValues[k], fErrors[k] );
1122 }
1123 else
1124 {
1125 *fLog << err << "Correlated Gauss fit failed : ierflg, istat = "
1126 << ierflg << ", " << istat << endl;
1127 }
1128
1129 Float_t charge = 0.0;
1130 Float_t chi2 = fmin;
1131 Int_t ndof = GetDegreesofFreedom()-fNumVar;
1132
1133 //get error matrix
1134 if (istat >= 1)
1135 {
1136 gMinuit->mnemat( &matrix[0][0], fNdim);
1137
1138 // copy covariance matrix into a matrix which includes also the fixed
1139 // parameters
1140 TString name;
1141 Double_t bnd1, bnd2, val, err;
1142 Int_t jvarbl;
1143 Int_t kvarbl;
1144 for (Int_t j=0; j<fNumVar; j++)
1145 {
1146 if (gMinuit)
1147 gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);
1148
1149 for (Int_t k=0; k<fNumVar; k++)
1150 {
1151 if (gMinuit)
1152 gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);
1153
1154 matrix[j][k] = jvarbl==0 || kvarbl==0 ? 0 : matrix[jvarbl-1][kvarbl-1];
1155 }
1156 }
1157 }
1158 else
1159 {
1160 // error matrix was not calculated, construct it
1161 if (fMinuitPrintOutLevel>=0)
1162 {
1163 *fLog << warn
1164 << "error matrix is not defined, construct it" << endl;
1165 }
1166 for (Int_t j=0; j<fNumVar; j++)
1167 {
1168 for (Int_t k=0; k<fNumVar; k++)
1169 matrix[j][k] = 0;
1170
1171 matrix[j][j] = fErrors[j]*fErrors[j];
1172 }
1173 }
1174
1175 if (fMinuitPrintOutLevel>=0)
1176 {
1177 *fLog << "Before calling SetFitValues : " << endl;
1178 *fLog << "fValues = " << fValues[0] << ", " << fValues[1] << ", "
1179 << fValues[2] << ", " << fValues[3] << ", " << fValues[4]
1180 << ", " << fValues[5] << endl;
1181
1182 *fLog << "fErrors = " << fErrors[0] << ", " << fErrors[1] << ", "
1183 << fErrors[2] << ", " << fErrors[3] << ", " << fErrors[4]
1184 << ", " << fErrors[5] << endl;
1185
1186 *fLog << "error matrix = " << matrix[0][0] << " " << matrix[0][1]
1187 << " " << matrix[0][2] << " " << matrix[0][3]
1188 << " " << matrix[0][4] << " " << matrix[0][5]
1189 << endl;
1190 *fLog << " " << matrix[1][0] << " " << matrix[1][1]
1191 << " " << matrix[1][2] << " " << matrix[1][3]
1192 << " " << matrix[1][4] << " " << matrix[1][5]
1193 << endl;
1194 *fLog << " " << matrix[2][0] << " " << matrix[2][1]
1195 << " " << matrix[2][2] << " " << matrix[2][3]
1196 << " " << matrix[2][4] << " " << matrix[2][5]
1197 << endl;
1198 *fLog << " " << matrix[3][0] << " " << matrix[3][1]
1199 << " " << matrix[3][2] << " " << matrix[3][3]
1200 << " " << matrix[3][4] << " " << matrix[3][5]
1201 << endl;
1202 *fLog << " " << matrix[4][0] << " " << matrix[4][1]
1203 << " " << matrix[4][2] << " " << matrix[4][3]
1204 << " " << matrix[4][4] << " " << matrix[4][5]
1205 << endl;
1206 *fLog << " " << matrix[5][0] << " " << matrix[5][1]
1207 << " " << matrix[5][2] << " " << matrix[5][3]
1208 << " " << matrix[5][4] << " " << matrix[5][5]
1209 << endl;
1210 }
1211
1212 star->SetFitValues(charge, fValues[0], fValues[1], fValues[3],
1213 fValues[2], fValues[4], fValues[5],
1214 matrix[1][1], matrix[1][3],matrix[3][3],
1215 chi2, ndof);
1216 }
1217
1218 //---------- for the correlated Gauss fit (end) ---------
1219
1220
1221 // reset the display to the starting values
1222 fDisplay.SetUsed(origPixelsUsed);
1223
1224 if (ierflg)
1225 return kCONTINUE;
1226 return kTRUE;
1227}
1228
1229Bool_t MFindStars::ShadowStar(MStarPos* star)
1230{
1231 UInt_t numPixels = fGeomCam->GetNumPixels();
1232
1233// Define an area around the star which will be set unused.
1234 UInt_t shadowPx=0;
1235 for (UInt_t pix=1; pix<numPixels; pix++)
1236 {
1237 Float_t pixXpos = (*fGeomCam)[pix].GetX();
1238 Float_t pixYpos = (*fGeomCam)[pix].GetY();
1239 Float_t starXpos = star->GetMeanX();
1240 Float_t starYpos = star->GetMeanY();
1241
1242 Float_t starSize = 3*sqrt( star->GetSigmaX()*star->GetSigmaX()
1243 + star->GetSigmaY()*star->GetSigmaY() );
1244
1245 Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
1246 (pixYpos-starYpos)*(pixYpos-starYpos));
1247
1248 if (dist > starSize && fDisplay.IsUsed(pix))
1249 fPixelsUsed[pix]=(Char_t)kTRUE;
1250 else
1251 {
1252 fPixelsUsed[pix]=(Char_t)kFALSE;
1253 shadowPx++;
1254 }
1255 }
1256
1257 if (fMinuitPrintOutLevel>=0) *fLog << dbg << " shadowPx " << shadowPx << endl;
1258
1259 fDisplay.SetUsed(fPixelsUsed);
1260
1261 return kTRUE;
1262}
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
Note: See TracBrowser for help on using the repository browser.