source: trunk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.cc@ 3900

Last change on this file since 3900 was 3900, checked in by jlopez, 21 years ago
*** empty log message ***
File size: 14.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!
18! Author(s): Javier López , 4/2004 <mailto:jlopez@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24#include "MPSFFitCalc.h"
25
26#include <iostream>
27
28
29#include <TH1D.h>
30#include <TString.h>
31#include <TArrayS.h>
32#include <TArrayI.h>
33#include <TArrayD.h>
34#include <TMinuit.h>
35#include <TStopwatch.h>
36
37#include "MGeomPix.h"
38#include "MGeomCam.h"
39#include "MMinuitInterface.h"
40
41#include "MLog.h"
42#include "MLogManip.h"
43
44#include "MParList.h"
45
46ClassImp(MPSFFitCalc);
47
48using namespace std;
49
50
51const Int_t numVar = 5;
52const Float_t pixelSize = 31.5; //[mm]
53const Float_t sqrt2 = sqrt(2.);
54const Float_t sqrt3 = sqrt(3.);
55const Bool_t usePrintOut = kFALSE;
56const Int_t minuitPrintOut = -1;
57
58UInt_t numPixels;
59Bool_t isPixelUsed[577];
60UInt_t numPixelsUsed;
61Float_t pixelPosX[577];
62Float_t pixelPosY[577];
63Float_t pixelValue[577];
64Float_t ChiSquare;
65
66MPSFFitCalc::MPSFFitCalc(ImgCleanMode_t imgmode, const char *name, const char *title)
67{
68 fName = name ? name : "MPSFFitCalc";
69 fTitle = title ? title : "Task that fits the PSF using the dc readout of the camera.";
70
71 fImgCleanMode = imgmode;
72 fNumRings = 6;
73 fPedLevel = 3.0;
74
75
76// Initialization of the camera dc mask 'fIsPixelused[]'
77 numPixels = 577;
78 for (UInt_t pixid=0; pixid<numPixels; pixid++)
79 isPixelUsed[pixid] = kTRUE;
80
81 fTotalMeanFit.Reset();
82 fNumTotalMeanFits = 0;
83
84 fMaxDC = 30.;
85
86 fPedestalDCHist = new TH1D("ped","",(Int_t)fMaxDC*10,0.,fMaxDC);
87
88 fVname = new TString[numVar];
89 fVinit.Set(numVar);
90 fStep.Set(numVar);
91 fLimlo.Set(numVar);
92 fLimup.Set(numVar);
93 fFix.Set(numVar);
94
95
96 fVname[0] = "max";
97 fVinit[0] = fMaxDC;
98 fStep[0] = fVinit[0]/sqrt2;
99 fLimlo[0] = 1.;
100 fLimup[0] = 40.;
101 fFix[0] = 0;
102
103 fVname[1] = "meanminor";
104 fVinit[1] = 0.;
105 fStep[1] = fVinit[0]/sqrt2;
106 fLimlo[1] = -600.;
107 fLimup[1] = 600.;
108 fFix[1] = 0;
109
110 fVname[2] = "sigmaminor";
111 fVinit[2] = pixelSize;
112 fStep[2] = fVinit[0]/sqrt2;
113 fLimlo[2] = pixelSize/(2*sqrt3);
114 fLimup[2] = 500.;
115 fFix[2] = 0;
116
117 fVname[3] = "meanmajor";
118 fVinit[3] = 0.;
119 fStep[3] = fVinit[0]/sqrt2;
120 fLimlo[3] = -600.;
121 fLimup[3] = 600.;
122 fFix[3] = 0;
123
124 fVname[4] = "sigmamajor";
125 fVinit[4] = pixelSize;
126 fStep[4] = fVinit[0]/sqrt2;
127 fLimlo[4] = pixelSize/(2*sqrt3);
128 fLimup[4] = 500.;
129 fFix[4] = 0;
130
131 fObjectFit = NULL;
132 // fMethod = "SIMPLEX";
133 fMethod = "MIGRAD";
134 fNulloutput = kFALSE;
135}
136
137MPSFFitCalc::~MPSFFitCalc()
138{
139 delete fPedestalDCHist;
140}
141
142//______________________________________________________________________________
143//
144// The 2D gaussian fucntion used to fit the spot of the star
145//
146static Double_t func(float x,float y,Double_t *par)
147{
148 Double_t value=par[0]*exp(-(x-par[1])*(x-par[1])/(2*par[2]*par[2]))*exp(-(y-par[3])*(y-par[3])/(2*par[4]*par[4]));
149 return value;
150}
151
152//______________________________________________________________________________
153//
154// Function used by Minuit to do the fit
155//
156static void fcnPSF(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
157{
158
159//calculate chisquare
160 Double_t chisq = 0;
161 Double_t delta;
162 Double_t x,y,z;
163 Double_t errorz = 0.2; //[uA]
164
165
166 for (UInt_t pixid=1; pixid<numPixels; pixid++)
167 {
168
169 if (isPixelUsed[pixid] && pixelValue[pixid]>0.)
170 {
171 x = pixelPosX[pixid];
172 y = pixelPosY[pixid];
173 z = pixelValue[pixid];
174
175 if (errorz > 0.0)
176 {
177 delta = (z-func(x,y,par))/errorz;
178 chisq += delta*delta;
179 }
180 else
181 cerr << "ERROR fcnPSF: This should never happen errorz[" << pixid << "] " << errorz << endl;
182 }
183 }
184 f = chisq;
185 ChiSquare = chisq;
186
187}
188
189
190void MPSFFitCalc::InitFitVariables()
191{
192
193 for (UInt_t pixid=1; pixid<577; pixid++)
194 pixelValue[pixid]=(Float_t)((*fCamera)[pixid]);
195
196 numPixelsUsed = 0;
197 Float_t totalDC = 0.0;
198
199 fVinit[0] = fMaxDC;
200
201 if(usePrintOut)
202 *fLog << dbg << "Pixels used in the fit \t";
203
204 for (UInt_t idx=0; idx<numPixels; idx++)
205 {
206 if (isPixelUsed[idx])
207 {
208 fVinit[1] += pixelPosX[idx]*pixelValue[idx];
209 fVinit[3] += pixelPosY[idx]*pixelValue[idx];
210 totalDC += pixelValue[idx];
211 numPixelsUsed++;
212
213 if(usePrintOut)
214 *fLog << dbg << ' ' << idx;
215 }
216 }
217 if(usePrintOut)
218 *fLog << dbg << endl;
219
220
221 fVinit[1] /= totalDC;
222 fVinit[3] /= totalDC;
223
224
225 fVinit[2] = pixelSize*sqrt((Float_t)numPixelsUsed)/2;
226 fVinit[4] = pixelSize*sqrt((Float_t)numPixelsUsed)/3;
227
228 //Init steps
229
230 for(Int_t i=0; i<numVar; i++)
231 if (fVinit[i] != 0)
232 fStep[i] = abs(fVinit[i]/sqrt2);
233
234
235 for (UInt_t pixidx=0; pixidx<numPixels; pixidx++)
236 if ( (*fGeomCam)[pixidx].GetSector() == 6)
237 fPedestalDCHist->Fill(pixelValue[pixidx]);
238
239 fPedestalDC = fPedestalDCHist->GetBinCenter(fPedestalDCHist->GetMaximumBin());
240
241 for (UInt_t pixid=1; pixid<577; pixid++)
242 pixelValue[pixid]-=fPedestalDC;
243
244 if(usePrintOut)
245 {
246 *fLog << dbg << "numPixelsUsed \t" << numPixelsUsed << endl;
247 *fLog << dbg << "fPedestalDC \t" << fPedestalDC << endl;
248 *fLog << dbg << "Maximun DC Init. value \t" << fVinit[0] << endl;
249 *fLog << dbg << "Mean Minor Axis Init. value \t" << fVinit[1] << " mm\t";
250 *fLog << dbg << "Sigma Minor Axis Init. value \t" << fVinit[2] << endl;
251 *fLog << dbg << "Mean Major Axis Init. value \t" << fVinit[3] << " mm\t";
252 *fLog << dbg << "Sigma Major Axis Init. value \t" << fVinit[4] << endl;
253 }
254}
255
256void MPSFFitCalc::RingImgClean()
257{
258
259 Bool_t isPixelUsedTmp[577];
260
261 fMaxDC=0;
262 UInt_t maxDCpix[2];
263
264
265// Look for the two neighbor pixels with the maximun signal and set all pixels as unused
266 Float_t dc[2];
267 Float_t dcsum;
268
269// Find the two close pixels with the maximun dc
270 for(UInt_t pixidx=0; pixidx<numPixels; pixidx++)
271 {
272
273 if(isPixelUsed[pixidx])
274 {
275 dc[0] = (Float_t)(*fCamera)[pixidx];
276 isPixelUsedTmp[pixidx] = kFALSE;
277
278 MGeomPix g = (*fGeomCam)[pixidx];
279 Int_t numNextNeighbors = g.GetNumNeighbors();
280
281 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
282 {
283 if(isPixelUsed[pixidx])
284 {
285 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
286 dc[1] = (Float_t)(*fCamera)[swneighbor];
287
288 dcsum = dc[0] + dc[1];
289
290 if(dcsum > fMaxDC*2)
291 {
292 maxDCpix[0] = pixidx;
293 maxDCpix[1] = swneighbor;
294 fMaxDC = dcsum/2;
295 }
296 }
297 }
298 }
299 }
300
301// Those variables are:
302// 1. an array of 2 dimensions
303// 1.1 the first dimension store the ring around the 'maxiun signal pixel'
304// 1.2 the sw numbers of the pixels in this ring
305// 2. an array with the number of pixels in each ring
306 UInt_t isPixelUesdInRing[fNumRings][577];
307 UInt_t numPixelsUsedInRing[fNumRings];
308
309// Store the 'maximun signal pixel in the [0] ring and set it as used
310
311 for (Int_t core=0; core<2; core++)
312 {
313 isPixelUesdInRing[0][0] = maxDCpix[core];
314 numPixelsUsedInRing[0] = 1;
315 isPixelUsedTmp[isPixelUesdInRing[0][0]] = kTRUE;
316
317 UInt_t count;
318
319// Loop over the neighbors of the previus ring and store the sw numbers in the 2D array to be
320// use in the next loop/ring
321 for (UShort_t ring=0; ring<fNumRings-1; ring++)
322 {
323 count = 0; // In this variable we count the pixels we are in each ring
324 for(UInt_t nextPixelUsed=0; nextPixelUsed<numPixelsUsedInRing[ring]; nextPixelUsed++)
325 {
326 MGeomPix g = (*fGeomCam)[isPixelUesdInRing[ring][nextPixelUsed]];
327 Int_t numNextNeighbors = g.GetNumNeighbors();
328 for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
329 {
330 UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
331 if (!isPixelUsedTmp[swneighbor])
332 {
333 isPixelUsedTmp[swneighbor] = kTRUE;
334 isPixelUesdInRing[ring+1][count] = swneighbor;
335 count++;
336 }
337 }
338 numPixelsUsedInRing[ring+1] = count;
339 }
340 }
341
342
343 if(usePrintOut)
344 {
345 for (UInt_t row=0; row<fNumRings; row++)
346 {
347
348 *fLog << dbg << "fIsPixeUsed[" << row << "][" << numPixelsUsedInRing[row] << "] ";
349 for (UInt_t column=0; column<numPixelsUsedInRing[row]; column++)
350 *fLog << dbg << isPixelUesdInRing[row][column] << ' ';
351 *fLog << dbg << endl;
352 }
353 }
354 }
355
356
357 for(UInt_t pixidx=0; pixidx<numPixels; pixidx++)
358 {
359 if(isPixelUsed[pixidx] && isPixelUsedTmp[pixidx])
360 isPixelUsed[pixidx] = kTRUE;
361 else
362 isPixelUsed[pixidx] = kFALSE;
363 }
364
365}
366
367void MPSFFitCalc::RmsImgClean(Float_t pedestal)
368{}
369
370void MPSFFitCalc::MaskImgClean(TArrayS blindpixels)
371{
372
373 Int_t npix = blindpixels.GetSize();
374
375 for (Int_t idx=0; idx<npix; idx++)
376 isPixelUsed[blindpixels[idx]] = kFALSE;
377
378}
379
380// --------------------------------------------------------------------------
381//
382// The PreProcess searches for the following input containers:
383// - MCameraDC
384//
385// The following output containers are also searched and created if
386// they were not found:
387//
388// - MPSFFit
389//
390Int_t MPSFFitCalc::PreProcess(MParList *pList)
391{
392
393 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
394
395 if (!fGeomCam)
396 {
397 *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl;
398 return kFALSE;
399 }
400
401// Initialize positions of the pixels to be used in the minuit minimizations
402
403 for (UInt_t pixid=1; pixid<577; pixid++)
404 {
405 MGeomPix &pix=(*fGeomCam)[pixid];
406 pixelPosX[pixid] = pix.GetX();
407 pixelPosY[pixid] = pix.GetY();
408 }
409
410 fCamera = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));
411
412 if (!fCamera)
413 {
414 *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl;
415 return kFALSE;
416 }
417
418 fFit = (MPSFFit*)pList->FindCreateObj(AddSerialNumber("MPSFFit"));
419 if (!fFit)
420 {
421 *fLog << err << AddSerialNumber("MPSFFit") << " cannot be created ... aborting" << endl;
422 return kFALSE;
423 }
424
425 // Minuit initialization
426
427 TMinuit *gMinuit = new TMinuit(6); //initialize TMinuit with a maximum of 5 params
428 gMinuit->SetFCN(fcnPSF);
429
430 Double_t arglist[10];
431 Int_t ierflg = 0;
432
433 arglist[0] = 1;
434 gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
435 arglist[0] = minuitPrintOut;
436 gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
437
438 return kTRUE;
439
440}
441
442Int_t MPSFFitCalc::Process()
443{
444
445 // ------------------------------------------
446 // image cleaning
447
448 switch(fImgCleanMode)
449 {
450 case kRing:
451 {
452 RingImgClean();
453 break;
454 }
455 case kRms:
456 {
457 RmsImgClean(fPedLevel);
458 break;
459 }
460 case kMask:
461 {
462 MaskImgClean(fBlindPixels);
463 break;
464 }
465 case kCombined:
466 {
467 MaskImgClean(fBlindPixels);
468 RingImgClean();
469 MaskImgClean(fBlindPixels);
470 break;
471 }
472 default:
473 {
474 *fLog << err << "Image Cleaning mode " << fImgCleanMode << " not defined" << endl;
475 return kFALSE;
476 }
477 }
478
479 InitFitVariables();
480
481 // -------------------------------------------
482 // call MINUIT
483
484 Double_t arglist[10];
485 Int_t ierflg = 0;
486
487 for (Int_t i=0; i<numVar; i++)
488 gMinuit->mnparm(i, fVname[i], fVinit[i], fStep[i], fLimlo[i], fLimup[i], ierflg);
489
490 TStopwatch clock;
491 clock.Start();
492
493// Now ready for minimization step
494 arglist[0] = 500;
495 arglist[1] = 1.;
496 gMinuit->mnexcm(fMethod, arglist ,2,ierflg);
497
498 clock.Stop();
499
500 if(usePrintOut)
501 {
502 *fLog << dbg << "Time spent for the minimization in MINUIT : " << endl;;
503 clock.Print();
504 }
505
506 if (ierflg)
507 {
508 *fLog << err << "MMinuitInterface::CallMinuit error " << ierflg << endl;
509 return kCONTINUE;
510 }
511
512 fNumTotalMeanFits++;
513
514 Double_t parm,parmerr;
515
516 gMinuit->GetParameter(0,parm,parmerr);
517 fFit->SetMaximun(parm);
518 fTotalMeanFit.SetMaximun(fTotalMeanFit.GetMaximun()+parm);
519
520 gMinuit->GetParameter(1,parm,parmerr);
521 fFit->SetMeanMinorAxis(parm);
522 fTotalMeanFit.SetMeanMinorAxis(fTotalMeanFit.GetMeanMinorAxis()+parm);
523
524 gMinuit->GetParameter(2,parm,parmerr);
525 fFit->SetSigmaMinorAxis(parm);
526 fTotalMeanFit.SetSigmaMinorAxis(fTotalMeanFit.GetSigmaMinorAxis()+parm);
527
528 gMinuit->GetParameter(3,parm,parmerr);
529 fFit->SetMeanMajorAxis(parm);
530 fTotalMeanFit.SetMeanMajorAxis(fTotalMeanFit.GetMeanMajorAxis()+parm);
531
532 gMinuit->GetParameter(4,parm,parmerr);
533 fFit->SetSigmaMajorAxis(parm);
534 fTotalMeanFit.SetSigmaMajorAxis(fTotalMeanFit.GetSigmaMajorAxis()+parm);
535
536 fFit->SetChisquare(ChiSquare/(numPixelsUsed-numVar));
537 fTotalMeanFit.SetChisquare(fTotalMeanFit.GetChisquare()+ChiSquare/(numPixelsUsed-numVar));
538
539 if(usePrintOut)
540 {
541 fFit->Print();
542 fTotalMeanFit.Print();
543 }
544
545 return kTRUE;
546}
547
548Int_t MPSFFitCalc::PostProcess()
549{
550
551 fTotalMeanFit.SetMaximun(fTotalMeanFit.GetMaximun()/fNumTotalMeanFits);
552 fTotalMeanFit.SetMeanMinorAxis(fTotalMeanFit.GetMeanMinorAxis()/fNumTotalMeanFits);
553 fTotalMeanFit.SetSigmaMinorAxis(fTotalMeanFit.GetSigmaMinorAxis()/fNumTotalMeanFits);
554 fTotalMeanFit.SetMeanMajorAxis(fTotalMeanFit.GetMeanMajorAxis()/fNumTotalMeanFits);
555 fTotalMeanFit.SetSigmaMajorAxis(fTotalMeanFit.GetSigmaMajorAxis()/fNumTotalMeanFits);
556 fTotalMeanFit.SetChisquare(fTotalMeanFit.GetChisquare()/fNumTotalMeanFits);
557
558 fFit->SetMaximun(fTotalMeanFit.GetMaximun());
559 fFit->SetMeanMinorAxis(fTotalMeanFit.GetMeanMinorAxis());
560 fFit->SetSigmaMinorAxis(fTotalMeanFit.GetSigmaMinorAxis());
561 fFit->SetMeanMajorAxis(fTotalMeanFit.GetMeanMajorAxis());
562 fFit->SetSigmaMajorAxis(fTotalMeanFit.GetSigmaMajorAxis());
563 fFit->SetChisquare(fTotalMeanFit.GetChisquare());
564
565 if(usePrintOut)
566 fTotalMeanFit.Print();
567
568 return kTRUE;
569}
Note: See TracBrowser for help on using the repository browser.