Index: trunk/MagicSoft/Mars/mtemp/MFindStars.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MFindStars.cc (revision 4050) +++ trunk/MagicSoft/Mars/mtemp/MFindStars.cc (revision 4051) @@ -30,4 +30,6 @@ #include "MFindStars.h" +#include +#include #include #include @@ -40,4 +42,6 @@ #include "MAstroCamera.h" #include "MMcConfigRunHeader.h" + +#include "MMinuitInterface.h" #include "MLog.h" @@ -59,5 +63,89 @@ using namespace std; -MFindStars::MFindStars(const char *name, const char *title) +const Float_t sqrt2 = sqrt(2.); +const Float_t sqrt3 = sqrt(3.); + +Bool_t HandleInput() +{ + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + while (1) + { + // + // While reading the input process gui events asynchronously + // + timer.TurnOn(); + cout << "Type 'q' to exit, to go on: " << endl; + char q; + cin >> q; + TString input = q; + timer.TurnOff(); + + if (input=="q\n") + return kFALSE; + + if (input=="\n") + return kTRUE; + }; +} + + +//______________________________________________________________________________ +// +// The 2D gaussian fucntion used to fit the spot of the star +// +static Double_t func(float x,float y,Double_t *par) +{ + 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])); + return value; +} + +//______________________________________________________________________________ +// +// Function used by Minuit to do the fit +// +static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + + MFindStars* find = (MFindStars*)gMinuit->GetObjectFit(); + MHCamera* display = (MHCamera*)find->GetDisplay(); + Float_t ped = find->GetPedestalDC(); + // Float_t rms = find->GetPedestalRMSDC(); + MGeomCam& geom = (MGeomCam&)display->GetGeomCam(); + + UInt_t numPixels = geom.GetNumPixels(); + +//calculate chisquare + Double_t chisq = 0; + Double_t delta; + Double_t x,y,z; + Double_t errorz = 0.2; //[uA] + + UInt_t usedPx=0; + for (UInt_t pixid=1; pixidGetBinContent(pixid+1)-ped; + + if (display->IsUsed(pixid) && z > 0.) + { + x = geom[pixid].GetX(); + y = geom[pixid].GetY(); + + if (errorz > 0.0) + { + usedPx++; + delta = (z-func(x,y,par))/errorz; + chisq += delta*delta; + } + else + cerr << " TMinuit::fcn errorz[" << pixid << "] " << errorz << endl; + } + } + f = chisq; + + find->SetChisquare(chisq); + find->SetDegreesofFreedom(usedPx); +} + +MFindStars::MFindStars(const char *name, const char *title): fNumVar(5) { fName = name ? name : "MFindStars"; @@ -71,4 +159,55 @@ fPixelsUsed.Set(577); fPixelsUsed.Reset((Char_t)kTRUE); + + //Fitting(Minuit) initialitation + const Float_t pixelSize = 31.5; //[mm] + + fVname = new TString[fNumVar]; + fVinit.Set(fNumVar); + fStep.Set(fNumVar); + fLimlo.Set(fNumVar); + fLimup.Set(fNumVar); + fFix.Set(fNumVar); + + fVname[0] = "max"; + fVinit[0] = 10.*fMaxNumIntegratedEvents; + fStep[0] = fVinit[0]/sqrt2; + fLimlo[0] = fMinDCForStars; + fLimup[0] = 30.*fMaxNumIntegratedEvents; + fFix[0] = 0; + + fVname[1] = "meanx"; + fVinit[1] = 0.; + fStep[1] = fVinit[0]/sqrt2; + fLimlo[1] = -600.; + fLimup[1] = 600.; + fFix[1] = 0; + + fVname[2] = "sigmaminor"; + fVinit[2] = pixelSize; + fStep[2] = fVinit[0]/sqrt2; + fLimlo[2] = pixelSize/(2*sqrt3); + fLimup[2] = 500.; + fFix[2] = 0; + + fVname[3] = "meany"; + fVinit[3] = 0.; + fStep[3] = fVinit[0]/sqrt2; + fLimlo[3] = -600.; + fLimup[3] = 600.; + fFix[3] = 0; + + fVname[4] = "sigmamajor"; + fVinit[4] = pixelSize; + fStep[4] = fVinit[0]/sqrt2; + fLimlo[4] = pixelSize/(2*sqrt3); + fLimup[4] = 500.; + fFix[4] = 0; + + fObjectFit = NULL; + // fMethod = "SIMPLEX"; + // fMethod = "MIGRAD"; + fMethod = "MINIMIZE"; + fNulloutput = kFALSE; } @@ -85,4 +224,5 @@ fDisplay.SetGeometry(*fGeomCam); + fDisplay.SetUsed(fPixelsUsed); fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC")); @@ -171,6 +311,11 @@ Int_t MFindStars::Process() { + UInt_t numPixels = fGeomCam->GetNumPixels(); + TArrayC origPixelsUsed; + origPixelsUsed.Set(numPixels); + if (fNumIntegratedEvents >= fMaxNumIntegratedEvents) { + if (fDrive) { @@ -185,7 +330,7 @@ else { - UInt_t numPixels = fGeomCam->GetNumPixels(); - TArrayC origPixelsUsed; - origPixelsUsed.Set(numPixels); + //Fist delete the previus stars in the list + fStars->GetList()->Delete(); + // for (UInt_t pix=1; pix(ped+5*rms)?fMinDCForStars:(ped+5*rms); - - *fLog << dbg << " DC pedestal = " << ped << " pedestal rms = " << rms << endl; + DCPedestalCalc(fPedestalDC, fPedestalRMSDC); + fMinDCForStars = fMinDCForStars>(fPedestalDC+5*fPedestalRMSDC)?fMinDCForStars:(fPedestalDC+5*fPedestalRMSDC); + + *fLog << dbg << " DC pedestal = " << fPedestalDC << " pedestal rms = " << fPedestalRMSDC << endl; *fLog << dbg << " fMinDCForStars " << fMinDCForStars << endl; @@ -211,10 +354,9 @@ while(FindPixelWithMaxDC(maxPixelDC, maxPixel)) { - *fLog << dbg << "Star candidate maxDC(" << setw(3) << maxPixelDC << " uA) x position(" << setw(3) << maxPixel.GetX() << " mm) x position(" << setw(3) << maxPixel.GetY() << " mm)" << endl; - + MStarLocalPos *starpos = new MStarLocalPos; starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY()); - starpos->SetCalcValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2); - starpos->SetFitValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2); + starpos->SetCalcValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2); + starpos->SetFitValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2,0.,1); fStars->GetList()->Add(starpos); @@ -235,19 +377,40 @@ *fLog << inf << GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl; + + for (UInt_t pix=1; pixGetList()); MStarLocalPos* star; while ((star=(MStarLocalPos*)Next())) - { - FindStar(star); - ShadowStar(star); - } + { + FindStar(star); + ShadowStar(star); + } //After finding stars reset all vairables fDisplay.Reset(); + fDisplay.SetUsed(origPixelsUsed); fNumIntegratedEvents=0; } + + for (UInt_t pix=1; pix maxDC*2) - { + UInt_t swneighbor = g.GetNeighbor(nextNeighbor); + if(fDisplay.IsUsed(swneighbor)) + { + dc[1] = fDisplay.GetBinContent(swneighbor+1); + if (dc[1] < fMinDCForStars) + continue; + + dcsum = dc[0] + dc[1]; + + if(dcsum > maxDC*2) + { if(dc[0]>=dc[1]) { - maxPixIdx[0] = pix; - maxPixIdx[1] = swneighbor; - maxDC = dc[0]; + maxPixIdx[0] = pix; + maxPixIdx[1] = swneighbor; + maxDC = dc[0]; } else @@ -367,8 +536,8 @@ maxDC = dc[1]; } - } - } - } - } + } + } + } + } } @@ -377,4 +546,7 @@ maxPix = (*fGeomCam)[maxPixIdx[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; + return kTRUE; } @@ -415,5 +587,6 @@ fDisplay.SetUsed(fPixelsUsed); - // determine mean x and mean y of the selected px +// determine mean x and mean y + Float_t max=0; Float_t meanX=0; Float_t meanY=0; @@ -432,4 +605,5 @@ Float_t pixYpos = (*fGeomCam)[pix].GetY(); + if (charge>max) max=charge; meanX += charge*pixXpos; meanY += charge*pixYpos; @@ -450,9 +624,82 @@ Float_t rmsY = TMath::Sqrt(meanSqY - meanY*meanY); - star->SetCalcValues(sumCharge,meanX,meanY,rmsX,rmsY); + star->SetCalcValues(sumCharge,max,meanX,meanY,rmsX,rmsY); + + +// fit the star spot using TMinuit + + for (UInt_t pix=1; pixGetMaxCalc()-fPedestalDC; + fLimlo[0] = fMinDCForStars-fPedestalDC; + fLimup[0] = fLimup[0]-fPedestalDC; + fVinit[1] = meanX; + fVinit[2] = rmsX; + fVinit[3] = meanY; + fVinit[4] = rmsY; + //Init steps + for(Int_t i=0; iGetParameter(0,maxFit, maxFitError); + gMinuit->GetParameter(1,meanXFit,meanXFitError); + gMinuit->GetParameter(2,sigmaMinorAxis,sigmaMinorAxisError); + gMinuit->GetParameter(3,meanYFit,meanYFitError); + gMinuit->GetParameter(4,sigmaMajorAxis,sigmaMajorAxisError); + + //FIXME: Do the integral properlly + integratedCharge = 0.; + + + } + else + { + maxFit = 0.; + meanXFit = 0.; + sigmaMinorAxis = 0.; + meanYFit = 0.; + sigmaMajorAxis = 0.; + integratedCharge = 0.; + } + + + + star->SetFitValues(integratedCharge,maxFit,meanXFit,meanYFit,sigmaMinorAxis,sigmaMajorAxis,chisquare,dregrees); + + // reset the display to the starting values fDisplay.SetUsed(origPixelsUsed); - return kTRUE; + return rc; } @@ -479,6 +726,6 @@ else { - fPixelsUsed[pix]=(Char_t)kFALSE; - shadowPx++; + fPixelsUsed[pix]=(Char_t)kFALSE; + shadowPx++; } } Index: trunk/MagicSoft/Mars/mtemp/MFindStars.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MFindStars.h (revision 4050) +++ trunk/MagicSoft/Mars/mtemp/MFindStars.h (revision 4051) @@ -51,4 +51,22 @@ Float_t fMinDCForStars; //[uA] + Float_t fPedestalDC; //[ua] + Float_t fPedestalRMSDC; //[ua] + + //Fitting(Minuit) variables + const Int_t fNumVar; + Float_t fTempChisquare; + Int_t fTempDegreesofFreedom; + + TString *fVname; + TArrayD fVinit; + TArrayD fStep; + TArrayD fLimlo; + TArrayD fLimup; + TArrayI fFix; + TObject *fObjectFit; + TString fMethod; + Bool_t fNulloutput; + Bool_t DCPedestalCalc(Float_t &ped, Float_t &rms); Bool_t FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix); @@ -68,5 +86,15 @@ void SetBlindPixels(TArrayS blindpixels); + void SetChisquare(Float_t chi) {fTempChisquare=chi;} + void SetDegreesofFreedom(Int_t free) {fTempDegreesofFreedom=free;} + MHCamera* GetDisplay() { return &fDisplay; } + Float_t GetPedestalDC() { return fPedestalDC; } + Float_t GetPedestalRMSDC() { return fPedestalRMSDC; } + + Float_t GetChisquare() {return fTempChisquare;} + Int_t GetDegreesofFreedom() {return fTempDegreesofFreedom;} + + ClassDef(MFindStars, 0) // Tool to find stars from DC Currents }; Index: trunk/MagicSoft/Mars/mtemp/MStarLocalCam.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MStarLocalCam.cc (revision 4050) +++ trunk/MagicSoft/Mars/mtemp/MStarLocalCam.cc (revision 4051) @@ -79,7 +79,14 @@ } +void MStarLocalCam::Paint(Option_t *o) +{ + TIter Next(fStars); + MStarLocalPos* star; + while ((star=(MStarLocalPos*)Next())) + star->Paint(o); +} + void MStarLocalCam::Print(Option_t *o) const { - //loop to extract position of stars on the camera TIter Next(fStars); MStarLocalPos* star; @@ -88,5 +95,5 @@ { *fLog << inf << "Star[" << starnum << "] info:" << endl; - star->Print(); + star->Print(o); starnum++; } Index: trunk/MagicSoft/Mars/mtemp/MStarLocalCam.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MStarLocalCam.h (revision 4050) +++ trunk/MagicSoft/Mars/mtemp/MStarLocalCam.h (revision 4051) @@ -30,5 +30,6 @@ TList *GetList() const { return fStars; } - void Print(Option_t *o="") const; + void Paint(Option_t *o=NULL); + void Print(Option_t *o=NULL) const; ClassDef(MStarLocalCam, 1) // Storage Container for star positions in the camera Index: trunk/MagicSoft/Mars/mtemp/MStarLocalPos.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MStarLocalPos.cc (revision 4050) +++ trunk/MagicSoft/Mars/mtemp/MStarLocalPos.cc (revision 4051) @@ -53,4 +53,5 @@ fMagCalc = 0.; + fMaxCalc = 0.; fMeanXCalc = 0.; fMeanYCalc = 0.; @@ -61,9 +62,11 @@ fMagFit = 0.; + fMaxFit = 0.; fMeanXFit = 0.; fMeanYFit = 0.; fSigmaMinorAxisFit = 0.; fSigmaMajorAxisFit = 0.; - fChisquarenDof = 0.; + fChiSquare = 0.; + fNdof = 0; } @@ -76,7 +79,8 @@ } -void MStarLocalPos::SetCalcValues(Float_t mag, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis) +void MStarLocalPos::SetCalcValues(Float_t mag, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis) { fMagCalc = mag; + fMaxCalc = max; fMeanXCalc = x; fMeanYCalc = y; @@ -85,24 +89,79 @@ } -void MStarLocalPos::SetFitValues(Float_t mag, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis) +void MStarLocalPos::SetFitValues(Float_t mag, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, Float_t chiSquare, Int_t ndof) { fMagFit = mag; + fMaxFit = max; fMeanXFit = x; fMeanYFit = y; fSigmaMinorAxisFit = sigmaMinorAxis; fSigmaMajorAxisFit = sigmaMajorAxis; + fChiSquare = chiSquare; + fNdof = ndof; } +// -------------------------------------------------------------------------- +// +// Paint the ellipse corresponding to the parameters +// void MStarLocalPos::Paint(Option_t *opt) -{} +{ + //Print a cross in the expected position + + if (fSigmaMinorAxisCalc>0. && fSigmaMajorAxisCalc>0.) + { + TEllipse ecalc(fMeanXCalc, fMeanYCalc, fSigmaMinorAxisCalc, fSigmaMajorAxisCalc, 0, 360, 0); + ecalc.SetLineWidth(2); + ecalc.SetLineColor(kRed); + ecalc.Paint(); + } + + if (fSigmaMinorAxisFit>0. || fSigmaMajorAxisFit>0.) + { + TEllipse efit(fMeanXFit, fMeanYFit, fSigmaMinorAxisFit, fSigmaMajorAxisFit, 0, 360, 0); + efit.SetLineWidth(2); + efit.SetLineColor(kBlack); + efit.Paint(); + } +} void MStarLocalPos::Print(Option_t *opt) const { - *fLog << inf << "Star position:" << endl; - *fLog << inf << " Expected \t X " << setw(4) << fXExp << " mm \tY " << setw(4) << fYExp << " mm" << endl; - *fLog << inf << " Calcultated \t X " << setw(4) << fMeanXCalc << " mm \tY " << setw(4) << fMeanYCalc << " mm" << endl; - *fLog << inf << " Fitted \t X " << setw(4) << fMeanXFit << " mm \tY " << setw(4) << fMeanYFit << " mm" << endl; - *fLog << inf << "Star size:" << endl; - *fLog << inf << " Calcultated \t X " << setw(4) << fSigmaMinorAxisCalc << " mm \tY " << setw(4) << fSigmaMajorAxisCalc << " mm" << endl; - *fLog << inf << " Fitted \t X " << setw(4) << fSigmaMinorAxisFit << " mm \tY " << setw(4) << fSigmaMajorAxisFit << " mm" << endl; + TString o = opt; + + if (o.Contains("mag", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star maginitude:" << endl; + *fLog << inf << " Expected \t" << setw(4) << fMagExp << endl; + *fLog << inf << " Calcultated \t " << setw(4) << fMagCalc << endl; + *fLog << inf << " Fitted \t " << setw(4) << fMagFit << endl; + } + + if (o.Contains("max", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star Maximum:" << endl; + *fLog << inf << " Calcultated \t " << setw(4) << fMaxCalc << " uA" << endl; + *fLog << inf << " Fitted \t " << setw(4) << fMaxFit << " uA" << endl; + } + + if (o.Contains("pos", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star position:" << endl; + *fLog << inf << " Expected \t X " << setw(4) << fXExp << " mm \tY " << setw(4) << fYExp << " mm" << endl; + *fLog << inf << " Calcultated \t X " << setw(4) << fMeanXCalc << " mm \tY " << setw(4) << fMeanYCalc << " mm" << endl; + *fLog << inf << " Fitted \t X " << setw(4) << fMeanXFit << " mm \tY " << setw(4) << fMeanYFit << " mm" << endl; + } + + if (o.Contains("siz", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star size:" << endl; + *fLog << inf << " Calcultated \t X " << setw(4) << fSigmaMinorAxisCalc << " mm \tY " << setw(4) << fSigmaMajorAxisCalc << " mm" << endl; + *fLog << inf << " Fitted \t X " << setw(4) << fSigmaMinorAxisFit << " mm \tY " << setw(4) << fSigmaMajorAxisFit << " mm" << endl; + } + + if (o.Contains("chi", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star Fit Quality:" << endl; + *fLog << inf << " ChiSquare/Ndof \t " << setw(3) << fChiSquare << "/" << fNdof << endl; + } } Index: trunk/MagicSoft/Mars/mtemp/MStarLocalPos.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MStarLocalPos.h (revision 4050) +++ trunk/MagicSoft/Mars/mtemp/MStarLocalPos.h (revision 4051) @@ -19,4 +19,5 @@ Float_t fMagCalc; + Float_t fMaxCalc; //[uA] Float_t fMeanXCalc; //[mm] Float_t fMeanYCalc; //[mm] @@ -27,9 +28,11 @@ Float_t fMagFit; + Float_t fMaxFit; //[uA] Float_t fMeanXFit; //[mm] Float_t fMeanYFit; //[mm] Float_t fSigmaMinorAxisFit; //[mm] Float_t fSigmaMajorAxisFit; //[mm] - Float_t fChisquarenDof; + Float_t fChiSquare; + Int_t fNdof; public: @@ -43,4 +46,5 @@ Float_t GetMagCalc() {return fMagCalc;} + Float_t GetMaxCalc() {return fMaxCalc;} Float_t GetMeanXCalc() {return fMeanXCalc;} Float_t GetMeanYCalc() {return fMeanYCalc;} @@ -49,20 +53,28 @@ Float_t GetMagFit() {return fMagFit;} + Float_t GetMaxFit() {return fMaxFit;} Float_t GetMeanXFit() {return fMeanXFit;} Float_t GetMeanYFit() {return fMeanYFit;} Float_t GetSigmaMinorAxisFit() {return fSigmaMinorAxisFit;} Float_t GetSigmaMajorAxisFit() {return fSigmaMajorAxisFit;} - Float_t GetChisquarenDof() {return fChisquarenDof;} + Float_t GetChiSquare() {return fChiSquare;} + Float_t GetNdof() {return fNdof;} + Float_t GetChiSquareNdof() {return fChiSquare/fNdof;} + + // Float_t GetMeanX(); + // Float_t GetMeanY(); + // Float_t GetSigmaMinorAxis(); + // Float_t GetSigmaMajorAxis(); void Reset(); void SetExpValues(Float_t mag, Float_t x, Float_t y); - void SetCalcValues(Float_t mag, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis); - void SetFitValues(Float_t mag, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis); + void SetCalcValues(Float_t mag, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis); + void SetFitValues(Float_t mag, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, Float_t chi, Int_t ndof); void Paint(Option_t *opt=NULL); void Print(Option_t *opt=NULL) const; - ClassDef(MStarLocalPos, 1) // Container that holds + ClassDef(MStarLocalPos, 1) // Container that holds the star information in the PMT camera };