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

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