/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!   Author(s): Robert Wagner, 2/2004 <mailto:rwagner@mppmu.mpg.de>
!   Author(s): Javier López , 4/2004 <mailto:jlopez@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MFindStars
//
/////////////////////////////////////////////////////////////////////////////
#include "MFindStars.h"

#include <TMinuit.h>
#include <TStopwatch.h>
#include <TTimer.h>
#include <TString.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1F.h>
#include <TF1.h>

#include "MObservatory.h"
#include "MAstroCamera.h"
#include "MMcConfigRunHeader.h"

#include "MLog.h"
#include "MLogManip.h"

#include "MHCamera.h"

#include "MGeomCam.h"
#include "MGeomPix.h"
#include "MCameraDC.h"
#include "MTime.h"
#include "MReportDrive.h"
#include "MStarLocalCam.h"
#include "MStarLocalPos.h"

#include "MParList.h"

ClassImp(MFindStars);
using namespace std;

const Float_t sqrt2 = sqrt(2.);
const Float_t sqrt3 = sqrt(3.);

//______________________________________________________________________________
//
// 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 = rms; //[uA]

    UInt_t usedPx=0;
    for (UInt_t pixid=1; pixid<numPixels; pixid++) 
    {
	z = display.GetBinContent(pixid+1)-ped;

	if (display.IsUsed(pixid))
	{
	    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";
  fTitle = title ? title : "Tool to find stars from DC Currents";

  fNumIntegratedEvents=0;
  fMaxNumIntegratedEvents = 10;
  fRingInterest = 125.; //[mm] ~ 0.4 deg
  fDCTailCut = 4;
  
  fPixelsUsed.Set(577);
  fPixelsUsed.Reset((Char_t)kTRUE);
  
  //Fitting(Minuit) initialitation
  const Float_t pixelSize = 31.5; //[mm]
  fMinuitPrintOutLevel = -1;
  
  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[1]/sqrt2;
  fLimlo[1] = -600.;
  fLimup[1] = 600.;
  fFix[1]   = 0;

  fVname[2] = "sigmaminor";
  fVinit[2] = pixelSize;
  fStep[2]  = fVinit[2]/sqrt2;
  fLimlo[2] = pixelSize/(2*sqrt3);
  fLimup[2] = 500.;
  fFix[2]   = 0;

  fVname[3] = "meany";
  fVinit[3] = 0.;
  fStep[3]  = fVinit[3]/sqrt2;
  fLimlo[3] = -600.;
  fLimup[3] = 600.;
  fFix[3]   = 0;

  fVname[4] = "sigmamajor";
  fVinit[4] = pixelSize;
  fStep[4]  = fVinit[4]/sqrt2;
  fLimlo[4] = pixelSize/(2*sqrt3);
  fLimup[4] = 500.;
  fFix[4]   = 0;

  fObjectFit  = NULL;
  //  fMethod     = "SIMPLEX";
  fMethod     = "MIGRAD";
  //  fMethod     = "MINIMIZE";
  fNulloutput = kFALSE;
}

Int_t MFindStars::PreProcess(MParList *pList)
{

    fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));

    if (!fGeomCam)
    {
      *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl;
      return kFALSE;
    }

    // Initialize camera display with the MGeomCam information
    fDisplay.SetGeometry(*fGeomCam);
    fDisplay.SetUsed(fPixelsUsed);
    
    fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));

    if (!fCurr)
    {
      *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl;
      return kFALSE;
    }

    fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents"));

    if (!fTimeCurr)
    {
      *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl;
      return kFALSE;
    }

    fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));

    if (!fDrive)
      {
        *fLog << warn << AddSerialNumber("MReportDrive") << " not found ... ignored." << endl;
      }
    else
      {
        //Initialitation MAstroCamera
        // Name of a MC file having MGeomCam and MMcConfigRunHeader
        TString fname = "../data/Gamma_zbin9_90_7_1480to1489_w0.root";

        // Time for which to get the picture
        //    MTime time;
        //    time.Set(2004, 2, 28, 01, 32, 15);
        
        // Current observatory
        MObservatory magic1;
        
        // Right Ascension [h] and declination [deg] of source
        // Currently 'perfect' pointing is assumed
        //    const Double_t ra  = MAstro::Hms2Rad(5, 34, 31.9);
        //    const Double_t dec = MAstro::Dms2Rad(22, 0, 52.0);

        // --------------------------------------------------------------------------
        // Create camera display from geometry
        MMcConfigRunHeader *config=0;
        MGeomCam           *geom=0;
        
        TFile file(fname);
        TTree *tree = (TTree*)file.Get("RunHeaders");
        tree->SetBranchAddress("MMcConfigRunHeader", &config);
        if (tree->GetBranch("MGeomCam"))
          tree->SetBranchAddress("MGeomCam", &geom);
        tree->GetEntry(0);
        
        fAstro.SetMirrors(*config->GetMirrors());
        fAstro.SetGeom(*geom);
        
        
        fAstro.SetLimMag(6);
        fAstro.SetRadiusFOV(3);
        fAstro.ReadBSC("../data/bsc5.dat");
        
        fAstro.SetObservatory(magic1);
        // Time for which to get the picture
        //    MTime time;
        //    time.Set(2004, 2, 28, 01, 32, 15);
        //   fAstro.SetTime(time);
        fAstro.SetTime(*fTimeCurr);
        fAstro.SetGuiActive();
        
        fAstro.SetStarList(fStars->GetList());
        
      }
    

    fStars = (MStarLocalCam*)pList->FindCreateObj(AddSerialNumber("MStarLocalCam"));
    if (!fStars)
    {
      *fLog << err << AddSerialNumber("MStarLocalCam") << " cannot be created ... aborting" << endl;
      return kFALSE;
    }

    fMinDCForStars = 1.5*fMaxNumIntegratedEvents; //[uA]

    // Initialize the TMinuit object

    TMinuit *gMinuit = new TMinuit(fNumVar);  //initialize TMinuit with a maximum of params
    gMinuit->SetFCN(fcn);

    Double_t arglist[10];
    Int_t ierflg = 0;

    arglist[0] = 1;
    gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
    arglist[0] = fMinuitPrintOutLevel;
    gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);

    // Set object pointer to allow mimuit to access internally
    gMinuit->SetObjectFit(this);

    return kTRUE;
}

Int_t MFindStars::Process()
{
  UInt_t numPixels = fGeomCam->GetNumPixels();
  TArrayC origPixelsUsed;
  origPixelsUsed.Set(numPixels);

    if (fNumIntegratedEvents >= fMaxNumIntegratedEvents)
    {

      if (fDrive)
        {
          Float_t ra  = fDrive->GetRa();
          Float_t dec = fDrive->GetDec();
          
          fAstro.SetRaDec(ra, dec);
          fAstro.SetGuiActive();
          
          fAstro.FillStarList();
        }
      else
        {
          //Fist delete the previus stars in the list
          fStars->GetList()->Delete();
          //

          for (UInt_t pix=1; pix<numPixels; pix++)
            {
              if (fDisplay.IsUsed(pix))
                origPixelsUsed[pix]=(Char_t)kTRUE;
              else
                  origPixelsUsed[pix]=(Char_t)kFALSE;
            }
          
	  DCPedestalCalc(fPedestalDC, fPedestalRMSDC);
	  fMinDCForStars = fMinDCForStars>(fPedestalDC+fDCTailCut*fPedestalRMSDC)?fMinDCForStars:(fPedestalDC+fDCTailCut*fPedestalRMSDC);

	  *fLog << dbg << " DC pedestal = " << fPedestalDC << " pedestal rms = " << fPedestalRMSDC << endl;
	  *fLog << dbg << " fMinDCForStars " << fMinDCForStars << endl;
	  
          // Find the star candidats searching the most brights pairs of pixels
          Float_t maxPixelDC;
          MGeomPix maxPixel;

	  while(FindPixelWithMaxDC(maxPixelDC, maxPixel))
            {

              MStarLocalPos *starpos = new MStarLocalPos;
              starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
              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);

              ShadowStar(starpos);
              
            }

          fDisplay.SetUsed(origPixelsUsed);
        }
      
      //loop to extract position of stars on the camera
      if (fStars->GetList()->GetSize() == 0)
        {
          *fLog << err << GetName() << " No stars candidates in the camera." << endl;
          return kFALSE;
        }
      else
          *fLog << inf << GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
          

      for (UInt_t pix=1; pix<numPixels; pix++)
        {
          if (fDisplay.IsUsed(pix))
            origPixelsUsed[pix]=(Char_t)kTRUE;
          else
            origPixelsUsed[pix]=(Char_t)kFALSE;
        }

      TIter Next(fStars->GetList());
      MStarLocalPos* star;
      while ((star=(MStarLocalPos*)Next()))
        {
          FindStar(star);
          ShadowStar(star);
        }
      
      //After finding stars reset all vairables
      fDisplay.Reset();
      fDisplay.SetUsed(origPixelsUsed);
      fNumIntegratedEvents=0;
    }

    for (UInt_t pix=1; pix<numPixels; pix++)
      {
        if (fDisplay.IsUsed(pix))
          origPixelsUsed[pix]=(Char_t)kTRUE;
        else
            origPixelsUsed[pix]=(Char_t)kFALSE;
        
      }

    fDisplay.AddCamContent(*fCurr);
    fNumIntegratedEvents++;
    fDisplay.SetUsed(origPixelsUsed);
    

  return kTRUE;
}

Int_t MFindStars::PostProcess()
{
  if(fStars->GetList()->GetSize() != 0)
    fStars->Print();

  return kTRUE;
}

void MFindStars::SetBlindPixels(TArrayS blindpixels)
{
    Int_t npix = blindpixels.GetSize();

    for (Int_t idx=0; idx<npix; idx++)
      {
	fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;
        *fLog << dbg << "MFindStars::SetBlindPixels fDisplay.IsUsed(" <<blindpixels[idx]  << ") kFALSE" << endl;
      }
    
    fDisplay.SetUsed(fPixelsUsed);
}

Bool_t MFindStars::DCPedestalCalc(Float_t &ped, Float_t &rms)
{
    //-------------------------------------------------------------
    // save pointer to the MINUIT object for optimizing the supercuts
    // because it will be overwritten 
    // when fitting the alpha distribution in MHFindSignificance
    TMinuit *savePointer = gMinuit;
    //-------------------------------------------------------------

   UInt_t numPixels = fGeomCam->GetNumPixels();

   TH1F dchist("dchist","",120,0.,30.*fMaxNumIntegratedEvents);
   for (UInt_t pix=1; pix<numPixels; pix++)
       dchist.Fill(fDisplay.GetBinContent(pix+1));

   Float_t nummaxprobdc = dchist.GetBinContent(dchist.GetMaximumBin());
   Float_t maxprobdc = dchist.GetBinCenter(dchist.GetMaximumBin());
   UInt_t bin = dchist.GetMaximumBin();
   do
   {
       bin++;
   }
   while(dchist.GetBinContent(bin)/nummaxprobdc > 0.5);
   Float_t halfmaxprobdc = dchist.GetBinCenter(bin);

   *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] ";
   *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << dchist.GetBinContent(bin) << "]" << endl;

   Float_t rmsguess = TMath::Abs(maxprobdc-halfmaxprobdc);
   Float_t min = maxprobdc-3*rmsguess;
   min = (min<0.?0.:min);
   Float_t max = maxprobdc+3*rmsguess;

   *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;

   TF1 func("func","gaus",min,max);
   func.SetParameters(nummaxprobdc, maxprobdc, rmsguess);
   
   dchist.Fit("func","QR0");
   // Remove the comments if you want to go through the file
   // event-by-event:
   //   HandleInput();

   UInt_t aproxnumdegrees = 6*(bin-dchist.GetMaximumBin());
   Float_t chiq = func.GetChisquare();
   ped = func.GetParameter(1);
   rms = func.GetParameter(2);

   *fLog << dbg << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;

   Int_t numsigmas = 5;
   Axis_t minbin = ped-numsigmas*rms/dchist.GetBinWidth(1);
   minbin=minbin<1?1:minbin;
   Axis_t maxbin = ped+numsigmas*rms/dchist.GetBinWidth(1);
   *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << dchist.Integral((int)minbin,(int)maxbin) << endl;

    //Check results from the fit are consistent
    if (TMath::Abs(ped-maxprobdc) > rmsguess || rms > rmsguess)
      {
        *fLog << warn << GetName() << " Pedestal DC fit give non consistent results." << endl;
        *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
        *fLog << warn << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
        ped = maxprobdc;
        rms = rmsguess/1.175; // FWHM=2.35*rms
      }

    //=================================================================

    // reset gMinuit to the MINUIT object for optimizing the supercuts 
    gMinuit = savePointer;
    //-------------------------------------------

   return kTRUE;
}
    
Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
{
    UInt_t numPixels = fGeomCam->GetNumPixels();

// Find the two close pixels with the maximun dc
    UInt_t maxPixIdx[2];

    maxDC = 0;

    for (UInt_t pix=1; pix<numPixels; pix++)
    {
	if(fDisplay.IsUsed(pix))
	{
	    Float_t dc[2];
	    dc[0] = fDisplay.GetBinContent(pix+1);
	    if (dc[0] < fMinDCForStars)
		continue;

	    MGeomPix &g = (*fGeomCam)[pix];
	    Int_t numNextNeighbors = g.GetNumNeighbors();
	    
	    Float_t dcsum;
	    for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
	    {
              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];
                        }
                      else
                        {
                          maxPixIdx[1] = pix;
                          maxPixIdx[0] = swneighbor;
                          maxDC = dc[1];
                        }
                    }	
                }
            }
        }
    }

    if (maxDC == 0)
	return kFALSE;

    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;
}

Bool_t MFindStars::FindStar(MStarLocalPos* star)
{    

  UInt_t numPixels = fGeomCam->GetNumPixels();
  TArrayC origPixelsUsed;
  origPixelsUsed.Set(numPixels);
  
  for (UInt_t pix=1; pix<numPixels; pix++)
    {
      if (fDisplay.IsUsed(pix))
        origPixelsUsed[pix]=(Char_t)kTRUE;
      else
        origPixelsUsed[pix]=(Char_t)kFALSE;
    }
  
  Float_t expX = star->GetXExp();
  Float_t expY = star->GetYExp();
  
  Float_t max=0;
  Float_t meanX=0;
  Float_t meanY=0;
  Float_t meanSqX=0;
  Float_t meanSqY=0;
  Float_t sumCharge=0;
  UInt_t usedPx=0;	

  const Float_t meanPresition = 3.; //[mm]
  const UInt_t maxNumIterations = 10;
  UInt_t numIterations = 0;

  do
    {
  // First define a area of interest around the expected position of the star
  for (UInt_t pix=1; pix<numPixels; pix++)
    {
      
      Float_t pixXpos=(*fGeomCam)[pix].GetX();
      Float_t pixYpos=(*fGeomCam)[pix].GetY();
      Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
                          (pixYpos-expY)*(pixYpos-expY));
      
      if ((dist < fRingInterest) && fDisplay.IsUsed(pix))
        fPixelsUsed[pix]=(Char_t)kTRUE;
      else
        fPixelsUsed[pix]=(Char_t)kFALSE;
    }
  
    fDisplay.SetUsed(fPixelsUsed);
    
// determine mean x and mean y
    usedPx=0;	
    for(UInt_t pix=0; pix<numPixels; pix++)
    {
	if(fDisplay.IsUsed(pix))
	{
	    usedPx++;

	    Float_t charge  = fDisplay.GetBinContent(pix+1);
	    Float_t pixXpos = (*fGeomCam)[pix].GetX();
	    Float_t pixYpos = (*fGeomCam)[pix].GetY();

            if (charge>max) max=charge;
	    meanX     += charge*pixXpos;
	    meanY     += charge*pixYpos;
	    meanSqX   += charge*pixXpos*pixXpos;
	    meanSqY   += charge*pixYpos*pixYpos;
	    sumCharge += charge;
	}
    }

    *fLog << dbg << " usedPx " << usedPx << endl;

    meanX   /= sumCharge;
    meanY   /= sumCharge;
    meanSqX /= sumCharge;
    meanSqY /= sumCharge;

    expX = meanX;
    expY = meanY;

    if (++numIterations >  maxNumIterations)
      {
        *fLog << warn << GetName() << "Mean calculation not converge after " << maxNumIterations << " iterations" << endl;
        break;
      }
        
    }while(TMath::Abs(meanX-expX) > meanPresition || TMath::Abs(meanY-expY) > meanPresition);
  
    Float_t rmsX = (meanSqX - meanX*meanX) - fRingInterest*fRingInterest/12;
    Float_t rmsY = (meanSqY - meanY*meanY) - fRingInterest*fRingInterest/12;

    if ( rmsX > 0)
      rmsX =  TMath::Sqrt(rmsX);
    else
      {
        *fLog << warn << " MFindStars::FindStar negative rmsX² " << rmsX << endl;
        *fLog << warn << " meanSqX " << meanSqX << " meanX " << meanX << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge << endl;
        rmsX = 0.;
      }

    if ( rmsY > 0)
      rmsY =  TMath::Sqrt(rmsY);
    else
      {
        *fLog << warn << " MFindStars::FindStar negative rmsY² " << rmsY << endl;
        *fLog << warn << " meanSqY " << meanSqY << " meanY " << meanY << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge<< endl;
        rmsY = 0.;
      }
    
    // Substrack pedestal DC
    sumCharge-=fPedestalDC*usedPx;
    max-=fPedestalDC;
    

    star->SetCalcValues(sumCharge,max,meanX,meanY,rmsX,rmsY);

    if (rmsX <= 0. || rmsY <= 0.)
      return kFALSE;
    
    
// fit the star spot using TMinuit

    
    for (UInt_t pix=1; pix<numPixels; pix++)
      if (fDisplay.IsUsed(pix))
        *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl;
  
  //Initialate variables for fit
    fVinit[0] = max;
    fLimlo[0] = fMinDCForStars-fPedestalDC;
    fLimup[0] = 30*fMaxNumIntegratedEvents-fPedestalDC;
    fVinit[1] = meanX;
    fVinit[2] = rmsX;
    fVinit[3] = meanY;
    fVinit[4] = rmsY;
    //Init steps
    for(Int_t i=0; i<fNumVar; i++)
	if (fVinit[i] != 0)
	  fStep[i] = TMath::Abs(fVinit[i]/sqrt2);
    //

    // -------------------------------------------
    // call MINUIT

    Double_t arglist[10];
    Int_t ierflg = 0;

    for (Int_t i=0; i<fNumVar; i++)
      gMinuit->mnparm(i, fVname[i], fVinit[i], fStep[i], fLimlo[i], fLimup[i], ierflg);

    TStopwatch clock;
    clock.Start();

// Now ready for minimization step
    arglist[0] = 500;
    arglist[1] = 1.;
    gMinuit->mnexcm(fMethod, arglist ,2,ierflg);

    clock.Stop();

    if(fMinuitPrintOutLevel>=0)
      {
	*fLog << dbg << "Time spent for the minimization in MINUIT :   " << endl;;
	clock.Print();
      }

    Double_t integratedCharge;
    Double_t maxFit, maxFitError;
    Double_t meanXFit, meanXFitError;
    Double_t sigmaMinorAxis, sigmaMinorAxisError;
    Double_t meanYFit, meanYFitError;
    Double_t sigmaMajorAxis, sigmaMajorAxisError;
    Float_t chisquare = GetChisquare();
    Int_t   dregrees  = GetDegreesofFreedom()-fNumVar;

    if (!ierflg)
      {
        gMinuit->GetParameter(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.;

	*fLog << err << "TMinuit::Call error " << ierflg << endl;
      }
    
    star->SetFitValues(integratedCharge,maxFit,meanXFit,meanYFit,sigmaMinorAxis,sigmaMajorAxis,chisquare,dregrees);
    
    // reset the display to the starting values
    fDisplay.SetUsed(origPixelsUsed);

    if (ierflg)
      return kCONTINUE;
    return kTRUE;
}

Bool_t MFindStars::ShadowStar(MStarLocalPos* star)
{
    UInt_t numPixels = fGeomCam->GetNumPixels();

// Define an area around the star which will be set unused.
    UInt_t shadowPx=0;	
    for (UInt_t pix=1; pix<numPixels; pix++)
    {
	Float_t pixXpos  = (*fGeomCam)[pix].GetX();
	Float_t pixYpos  = (*fGeomCam)[pix].GetY();
        Float_t starXpos = star->GetMeanX();
        Float_t starYpos = star->GetMeanY();
        
	Float_t starSize = 3*star->GetSigmaMajorAxis();
        
	Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
			    (pixYpos-starYpos)*(pixYpos-starYpos));

        if (dist > starSize && fDisplay.IsUsed(pix))
	  fPixelsUsed[pix]=(Char_t)kTRUE;
        else
          {
            fPixelsUsed[pix]=(Char_t)kFALSE;
            shadowPx++;
          }
    }

    *fLog << dbg << " shadowPx " << shadowPx << endl;

    fDisplay.SetUsed(fPixelsUsed);

    return kTRUE;
}





