/* ======================================================================== *\
!
! *
! * 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 <TFile.h>
#include <TTree.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;

MFindStars::MFindStars(const char *name, const char *title)
{
  fName  = name  ? name  : "MFindStars";
  fTitle = title ? title : "Tool to find stars from DC Currents";

  fMaxNumIntegratedEvents = 10;
  fRingInterest = 100.; //[mm] ~ 0.3 deg
  fMinDCForStars = 3.*fMaxNumIntegratedEvents; //[uA]
  
  fPixelsUsed.Set(577);
  fPixelsUsed.Reset((Char_t)kTRUE);
}

Int_t MFindStars::PreProcess(MParList *pList)
{

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

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

    fDisplay.SetGeometry(*fGeomCam);
    
    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;
    }

    return kTRUE;
}

Int_t MFindStars::Process()
{
    if (fNumIntegratedEvents < fMaxNumIntegratedEvents)
    {
      fDisplay.AddCamContent(*fCurr);
      fNumIntegratedEvents++;
    }
    else
    {
      if (fDrive)
        {
          Float_t ra  = fDrive->GetRa();
          Float_t dec = fDrive->GetDec();
          
          fAstro.SetRaDec(ra, dec);
          fAstro.SetGuiActive();
          
          fAstro.FillStarList();
        }
      else
        {
          Float_t maxPixelDC;
          MGeomPix maxPixel;

          MHCamera cam = fDisplay;
          cam.PrintInfo();

          // Find the star candidats searching the most brights pairs of pixels
         while(1)
            {
              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;
              
              if (maxPixelDC<fMinDCForStars)
                break;
              
              MStarLocalPos *starpos = new MStarLocalPos;
              starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
              starpos->SetCalcValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
              fStars->GetList()->Add(starpos);

              ShadowStar(starpos);
              
            }
         fDisplay = cam;
         *fLog << dbg <<  GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
        }
      
      //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;
        }
          
      TIter Next(fStars->GetList());
      MStarLocalPos* star;
      while ((star=(MStarLocalPos*)Next()))
        {
          FindStar(star);
          ShadowStar(star);
        }
      
      
      //After finding stars reset all vairables
      fDisplay.Reset();
      fNumIntegratedEvents=0;
    }

  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;

    fDisplay.SetUsed(fPixelsUsed);
}

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);
	    
	    MGeomPix &g = (*fGeomCam)[pix];
	    Int_t numNextNeighbors = g.GetNumNeighbors();
	    
	    Float_t dcsum;
	    for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
	    {
		if(fDisplay.IsUsed(pix))
		{
		    UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
		    dc[1] = fDisplay.GetBinContent(swneighbor+1);
		    
		    
		    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];
                        }
		    }	
		}
	    }
	}
    }

    maxPix = (*fGeomCam)[maxPixIdx[0]];
    return kTRUE;
}

Bool_t MFindStars::FindStar(MStarLocalPos* star)
{    

  MHCamera cam = fDisplay;
  
    Float_t expX = star->GetXExp();
    Float_t expY = star->GetYExp();

    UInt_t numPixels = fGeomCam->GetNumPixels();

// 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 && cam.IsUsed(pix))
	  fPixelsUsed[pix]=(Char_t)kTRUE;
    }

    cam.SetUsed(fPixelsUsed);


    // determine mean x and mean y of the selected px
    Float_t meanX=0;
    Float_t meanY=0;
    Float_t meanSqX=0;
    Float_t meanSqY=0;
    Float_t sumCharge=0;
    UInt_t usedPx=0;	
    for(UInt_t pix=0; pix<numPixels; pix++)
    {
	if(cam.IsUsed(pix))
	{
	    usedPx++;

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

	    meanX+=charge*pixXpos;
	    meanY+=charge*pixYpos;
	    meanSqX+=charge*pixXpos*pixXpos;
	    meanSqY+=charge*pixYpos*pixYpos;
	    sumCharge+=charge;
//	    cam.ResetUsed(i); // this px must not be used again!
	} //for... use	
    }

    star->SetCalcValues(sumCharge,meanX,meanY,meanSqX,meanSqY);
    
    return kTRUE;
}

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

// Define an area around the star which will be set unused.
    for (UInt_t pix=1; pix<numPixels; pix++)
    {

	Float_t pixXpos  = (*fGeomCam)[pix].GetX();
	Float_t pixYpos  = (*fGeomCam)[pix].GetY();
        Float_t starXpos = star->GetMeanXCalc();
        Float_t starYpos = star->GetMeanYCalc();
        
	Float_t starSize = 2*star->GetSigmaMajorAxisCalc();
        
	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;
    }

    fDisplay.SetUsed(fPixelsUsed);

    return kTRUE;
}



