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

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