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

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